DATA 605 - Final Project
Problem 2
Adam Douglas
5/12/2019
Descriptive and Inferential Statistics.
# Read in the training data
houses <- read_csv("train.csv")# Change character columns to factors
fac <- houses %>% select_if(is.character) %>% colnames()
houses[fac] <- lapply(houses[fac], factor)
# Replace NAs with "No"
houses[fac] <- lapply(houses[fac], fct_explicit_na,'No')
# Replace numeric NAs with 0
houses <- houses %>% mutate_if(is.numeric, ~replace(., is.na(.), 0))Since the SalePrice column will be the target variable, we’ll start there and look at how it is distributed.
# Plot SalePrice
houses %>% ggplot(aes(y=SalePrice)) +
geom_boxplot(outlier.color="firebrick2", outlier.alpha = 0.3) +
scale_y_continuous(labels=comma) +
scale_x_discrete() +
stat_boxplot(geom ='errorbar',width=.3) +
labs(title="Distribution of Home Sale Price",
subtitle="All Homes", y="Price($)",
x="All Homes") + theme_classic()
Looking at the sale price of houses, the values are pretty evenly distributed about $163,000 with some significant outliers well into the $600,000 range.
Let’s get some summary statistics of the other variables in the data set:
summary(houses)## Id MSSubClass MSZoning LotFrontage
## Min. : 1.0 Min. : 20.0 C (all): 10 Min. : 0.00
## 1st Qu.: 365.8 1st Qu.: 20.0 FV : 65 1st Qu.: 42.00
## Median : 730.5 Median : 50.0 RH : 16 Median : 63.00
## Mean : 730.5 Mean : 56.9 RL :1151 Mean : 57.62
## 3rd Qu.:1095.2 3rd Qu.: 70.0 RM : 218 3rd Qu.: 79.00
## Max. :1460.0 Max. :190.0 Max. :313.00
##
## LotArea Street Alley LotShape LandContour
## Min. : 1300 Grvl: 6 Grvl: 50 IR1:484 Bnk: 63
## 1st Qu.: 7554 Pave:1454 Pave: 41 IR2: 41 HLS: 50
## Median : 9478 No :1369 IR3: 10 Low: 36
## Mean : 10517 Reg:925 Lvl:1311
## 3rd Qu.: 11602
## Max. :215245
##
## Utilities LotConfig LandSlope Neighborhood Condition1
## AllPub:1459 Corner : 263 Gtl:1382 NAmes :225 Norm :1260
## NoSeWa: 1 CulDSac: 94 Mod: 65 CollgCr:150 Feedr : 81
## FR2 : 47 Sev: 13 OldTown:113 Artery : 48
## FR3 : 4 Edwards:100 RRAn : 26
## Inside :1052 Somerst: 86 PosN : 19
## Gilbert: 79 RRAe : 11
## (Other):707 (Other): 15
## Condition2 BldgType HouseStyle OverallQual
## Norm :1445 1Fam :1220 1Story :726 Min. : 1.000
## Feedr : 6 2fmCon: 31 2Story :445 1st Qu.: 5.000
## Artery : 2 Duplex: 52 1.5Fin :154 Median : 6.000
## PosN : 2 Twnhs : 43 SLvl : 65 Mean : 6.099
## RRNn : 2 TwnhsE: 114 SFoyer : 37 3rd Qu.: 7.000
## PosA : 1 1.5Unf : 14 Max. :10.000
## (Other): 2 (Other): 19
## OverallCond YearBuilt YearRemodAdd RoofStyle
## Min. :1.000 Min. :1872 Min. :1950 Flat : 13
## 1st Qu.:5.000 1st Qu.:1954 1st Qu.:1967 Gable :1141
## Median :5.000 Median :1973 Median :1994 Gambrel: 11
## Mean :5.575 Mean :1971 Mean :1985 Hip : 286
## 3rd Qu.:6.000 3rd Qu.:2000 3rd Qu.:2004 Mansard: 7
## Max. :9.000 Max. :2010 Max. :2010 Shed : 2
##
## RoofMatl Exterior1st Exterior2nd MasVnrType MasVnrArea
## CompShg:1434 VinylSd:515 VinylSd:504 BrkCmn : 15 Min. : 0.0
## Tar&Grv: 11 HdBoard:222 MetalSd:214 BrkFace:445 1st Qu.: 0.0
## WdShngl: 6 MetalSd:220 HdBoard:207 None :864 Median : 0.0
## WdShake: 5 Wd Sdng:206 Wd Sdng:197 Stone :128 Mean : 103.1
## ClyTile: 1 Plywood:108 Plywood:142 No : 8 3rd Qu.: 164.2
## Membran: 1 CemntBd: 61 CmentBd: 60 Max. :1600.0
## (Other): 2 (Other):128 (Other):136
## ExterQual ExterCond Foundation BsmtQual BsmtCond BsmtExposure
## Ex: 52 Ex: 3 BrkTil:146 Ex:121 Fa: 45 Av:221
## Fa: 14 Fa: 28 CBlock:634 Fa: 35 Gd: 65 Gd:134
## Gd:488 Gd: 146 PConc :647 Gd:618 Po: 2 Mn:114
## TA:906 Po: 1 Slab : 24 TA:649 TA:1311 No:991
## TA:1282 Stone : 6 No: 37 No: 37
## Wood : 3
##
## BsmtFinType1 BsmtFinSF1 BsmtFinType2 BsmtFinSF2
## ALQ:220 Min. : 0.0 ALQ: 19 Min. : 0.00
## BLQ:148 1st Qu.: 0.0 BLQ: 33 1st Qu.: 0.00
## GLQ:418 Median : 383.5 GLQ: 14 Median : 0.00
## LwQ: 74 Mean : 443.6 LwQ: 46 Mean : 46.55
## Rec:133 3rd Qu.: 712.2 Rec: 54 3rd Qu.: 0.00
## Unf:430 Max. :5644.0 Unf:1256 Max. :1474.00
## No : 37 No : 38
## BsmtUnfSF TotalBsmtSF Heating HeatingQC CentralAir
## Min. : 0.0 Min. : 0.0 Floor: 1 Ex:741 N: 95
## 1st Qu.: 223.0 1st Qu.: 795.8 GasA :1428 Fa: 49 Y:1365
## Median : 477.5 Median : 991.5 GasW : 18 Gd:241
## Mean : 567.2 Mean :1057.4 Grav : 7 Po: 1
## 3rd Qu.: 808.0 3rd Qu.:1298.2 OthW : 2 TA:428
## Max. :2336.0 Max. :6110.0 Wall : 4
##
## Electrical 1stFlrSF 2ndFlrSF LowQualFinSF
## FuseA: 94 Min. : 334 Min. : 0 Min. : 0.000
## FuseF: 27 1st Qu.: 882 1st Qu.: 0 1st Qu.: 0.000
## FuseP: 3 Median :1087 Median : 0 Median : 0.000
## Mix : 1 Mean :1163 Mean : 347 Mean : 5.845
## SBrkr:1334 3rd Qu.:1391 3rd Qu.: 728 3rd Qu.: 0.000
## No : 1 Max. :4692 Max. :2065 Max. :572.000
##
## GrLivArea BsmtFullBath BsmtHalfBath FullBath
## Min. : 334 Min. :0.0000 Min. :0.00000 Min. :0.000
## 1st Qu.:1130 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:1.000
## Median :1464 Median :0.0000 Median :0.00000 Median :2.000
## Mean :1515 Mean :0.4253 Mean :0.05753 Mean :1.565
## 3rd Qu.:1777 3rd Qu.:1.0000 3rd Qu.:0.00000 3rd Qu.:2.000
## Max. :5642 Max. :3.0000 Max. :2.00000 Max. :3.000
##
## HalfBath BedroomAbvGr KitchenAbvGr KitchenQual
## Min. :0.0000 Min. :0.000 Min. :0.000 Ex:100
## 1st Qu.:0.0000 1st Qu.:2.000 1st Qu.:1.000 Fa: 39
## Median :0.0000 Median :3.000 Median :1.000 Gd:586
## Mean :0.3829 Mean :2.866 Mean :1.047 TA:735
## 3rd Qu.:1.0000 3rd Qu.:3.000 3rd Qu.:1.000
## Max. :2.0000 Max. :8.000 Max. :3.000
##
## TotRmsAbvGrd Functional Fireplaces FireplaceQu GarageType
## Min. : 2.000 Maj1: 14 Min. :0.000 Ex: 24 2Types : 6
## 1st Qu.: 5.000 Maj2: 5 1st Qu.:0.000 Fa: 33 Attchd :870
## Median : 6.000 Min1: 31 Median :1.000 Gd:380 Basment: 19
## Mean : 6.518 Min2: 34 Mean :0.613 Po: 20 BuiltIn: 88
## 3rd Qu.: 7.000 Mod : 15 3rd Qu.:1.000 TA:313 CarPort: 9
## Max. :14.000 Sev : 1 Max. :3.000 No:690 Detchd :387
## Typ :1360 No : 81
## GarageYrBlt GarageFinish GarageCars GarageArea GarageQual
## Min. : 0 Fin:352 Min. :0.000 Min. : 0.0 Ex: 3
## 1st Qu.:1958 RFn:422 1st Qu.:1.000 1st Qu.: 334.5 Fa: 48
## Median :1977 Unf:605 Median :2.000 Median : 480.0 Gd: 14
## Mean :1869 No : 81 Mean :1.767 Mean : 473.0 Po: 3
## 3rd Qu.:2001 3rd Qu.:2.000 3rd Qu.: 576.0 TA:1311
## Max. :2010 Max. :4.000 Max. :1418.0 No: 81
##
## GarageCond PavedDrive WoodDeckSF OpenPorchSF EnclosedPorch
## Ex: 2 N: 90 Min. : 0.00 Min. : 0.00 Min. : 0.00
## Fa: 35 P: 30 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.00
## Gd: 9 Y:1340 Median : 0.00 Median : 25.00 Median : 0.00
## Po: 7 Mean : 94.24 Mean : 46.66 Mean : 21.95
## TA:1326 3rd Qu.:168.00 3rd Qu.: 68.00 3rd Qu.: 0.00
## No: 81 Max. :857.00 Max. :547.00 Max. :552.00
##
## 3SsnPorch ScreenPorch PoolArea PoolQC
## Min. : 0.00 Min. : 0.00 Min. : 0.000 Ex: 2
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.000 Fa: 2
## Median : 0.00 Median : 0.00 Median : 0.000 Gd: 3
## Mean : 3.41 Mean : 15.06 Mean : 2.759 No:1453
## 3rd Qu.: 0.00 3rd Qu.: 0.00 3rd Qu.: 0.000
## Max. :508.00 Max. :480.00 Max. :738.000
##
## Fence MiscFeature MiscVal MoSold
## GdPrv: 59 Gar2: 2 Min. : 0.00 Min. : 1.000
## GdWo : 54 Othr: 2 1st Qu.: 0.00 1st Qu.: 5.000
## MnPrv: 157 Shed: 49 Median : 0.00 Median : 6.000
## MnWw : 11 TenC: 1 Mean : 43.49 Mean : 6.322
## No :1179 No :1406 3rd Qu.: 0.00 3rd Qu.: 8.000
## Max. :15500.00 Max. :12.000
##
## YrSold SaleType SaleCondition SalePrice
## Min. :2006 WD :1267 Abnorml: 101 Min. : 34900
## 1st Qu.:2007 New : 122 AdjLand: 4 1st Qu.:129975
## Median :2008 COD : 43 Alloca : 12 Median :163000
## Mean :2008 ConLD : 9 Family : 20 Mean :180921
## 3rd Qu.:2009 ConLI : 5 Normal :1198 3rd Qu.:214000
## Max. :2010 ConLw : 5 Partial: 125 Max. :755000
## (Other): 9Typically, square footage is a big factor in a home’s price. We can look at total square footage by creating a new variable and plotting that versus sale prices.
houses$TotalSqFt <- houses$`1stFlrSF` + houses$`2ndFlrSF`houses %>% ggplot(aes(x=TotalSqFt, y=SalePrice)) +
geom_point(alpha=0.2) + scale_y_continuous(labels=comma) +
scale_x_continuous(labels=comma) +
labs(title="Home Sale Price by Above Ground Square Footage",
subtitle="All Homes", y="Price($)",
x="Sq. Ft.") + theme_classic() +
theme(panel.grid.major.y = element_line(colour = "grey90",linetype='dashed'))
Visually, there appears to be a positive correlation between square footage of the main house and the sale price.
What about the finished basement areas?
houses %>% mutate(BsmtFinSF = BsmtFinSF1+BsmtFinSF2) %>%
filter(BsmtFinSF > 0) %>%
ggplot(aes(x=BsmtFinSF, y=SalePrice)) +
geom_point(alpha=0.2) + scale_y_continuous(labels=comma) +
scale_x_continuous(labels=comma) +
labs(title="Home Sale Price by Finished Basement Square Footage",
subtitle="w/ Finished Basement Space", y="Price($)",
x="Sq. Ft.") + theme_classic() +
theme(panel.grid.major.y = element_line(colour = "grey90",linetype='dashed'))
Looking at only houses with finished basement space, there may be a relationship as well.
Finally, looking at garage area, we see a similar relationship:
houses %>% ggplot(aes(x=GarageArea, y=SalePrice)) +
geom_point(alpha=0.2) + scale_y_continuous(labels=comma) +
scale_x_continuous(labels=comma) +
labs(title="Home Sale Price by Garage Square Footage",
subtitle="All Houses", y="Price($)",
x="Sq. Ft.") + theme_classic() +
theme(panel.grid.major.y = element_line(colour = "grey90",linetype='dashed'))
# Subset of variables
housecor <- houses %>% mutate(BsmtFinSF = BsmtFinSF1+BsmtFinSF2) %>% dplyr::select(SalePrice, TotalSqFt, BsmtFinSF, GarageArea)
# Compute correlations
corr <- cor(housecor)
ggcorrplot(corr,lab=TRUE, ggtheme = ggplot2::theme_classic) +
labs(title="Correlation",subtitle="All Houses")
We have some positive correlation between sale price of the houses and total square footage, total finished basement square footage, and garage area. Let’s look at how solid the correlations are:
# Custom Function to nicely display confidence intervals
# and p-values of Correlation Coefficients
corr_hypothesis <- function(x,y,lvl=0.95){
z <- cor.test(x,y,conf.level=lvl)
ci_low <- round(z$conf.int[[1]],3)
ci_high <- round(z$conf.int[[2]],3)
pval <- prettyNum(z$p.value,digits=3)
x.name <- deparse(substitute(x))
if(str_detect(x.name,'\\$')){
x.name <- str_split(x.name,'\\$')[[1]][2]
}
y.name <- deparse(substitute(y))
if(str_detect(y.name,'\\$')){
y.name <- str_split(y.name,'\\$')[[1]][2]
}
paste("Confidence Interval (",x.name," and ",y.name,")",
": ",ci_low," - ", ci_high, " (p-value: ",pval, ")",
sep="")
}# Compute each correlation's CI and p-value
corr_hypothesis(housecor$SalePrice, housecor$TotalSqFt, lvl=0.99)## [1] "Confidence Interval (SalePrice and TotalSqFt): 0.682 - 0.748 (p-value: 1.18e-230)"corr_hypothesis(housecor$SalePrice, housecor$BsmtFinSF, lvl=0.99)## [1] "Confidence Interval (SalePrice and BsmtFinSF): 0.307 - 0.423 (p-value: 1.36e-47)"corr_hypothesis(housecor$SalePrice, housecor$GarageArea, lvl=0.99)## [1] "Confidence Interval (SalePrice and GarageArea): 0.58 - 0.663 (p-value: 5.27e-158)"corr_hypothesis(housecor$BsmtFinSF, housecor$TotalSqFt, lvl=0.99)## [1] "Confidence Interval (BsmtFinSF and TotalSqFt): 0.138 - 0.267 (p-value: 4.7e-15)"corr_hypothesis(housecor$BsmtFinSF, housecor$GarageArea, lvl=0.99)## [1] "Confidence Interval (BsmtFinSF and GarageArea): 0.215 - 0.339 (p-value: 2.22e-27)"corr_hypothesis(housecor$TotalSqFt, housecor$GarageArea, lvl=0.99)## [1] "Confidence Interval (TotalSqFt and GarageArea): 0.426 - 0.529 (p-value: 1.11e-84)"We can see from the hypothesis tests that there is strong evidence for correlation amongst all of these variables. The correlation amongst the independent variables may be of concern when using them in a multiple regression model.
To help prevent familywise error, we used a 99% confidence interval in our test (with only 6 hypothesis tests at \(\alpha = 0.01\) the chances of familywise Type I error are small).
Linear Algebra and Correlation
# Compute precision matrix
precision <- solve(corr)
precision## SalePrice TotalSqFt BsmtFinSF GarageArea
## SalePrice 2.7920278 -1.4727056 -0.4722635 -0.9034946
## TotalSqFt -1.4727056 2.0838074 0.1503380 -0.1222487
## BsmtFinSF -0.4722635 0.1503380 1.1713092 -0.1036110
## GarageArea -0.9034946 -0.1222487 -0.1036110 1.6506838cov(housecor)## SalePrice TotalSqFt BsmtFinSF GarageArea
## SalePrice 6311111264 29680776.55 13855553.07 10589102.52
## TotalSqFt 29680777 271611.42 50389.65 53394.76
## BsmtFinSF 13855553 50389.65 226674.36 28330.71
## GarageArea 10589103 53394.76 28330.71 45712.51corr %*% precision## SalePrice TotalSqFt BsmtFinSF GarageArea
## SalePrice 1.000000e+00 -1.387779e-17 -1.387779e-17 0.000000e+00
## TotalSqFt -3.885781e-16 1.000000e+00 6.938894e-18 -1.110223e-16
## BsmtFinSF -1.665335e-16 6.938894e-17 1.000000e+00 -5.551115e-17
## GarageArea -3.330669e-16 6.938894e-17 -2.775558e-17 1.000000e+00precision %*% corr## SalePrice TotalSqFt BsmtFinSF GarageArea
## SalePrice 1.000000e+00 -4.996004e-16 -2.220446e-16 -5.551115e-16
## TotalSqFt -2.775558e-17 1.000000e+00 6.245005e-17 4.163336e-17
## BsmtFinSF -2.775558e-17 0.000000e+00 1.000000e+00 -4.163336e-17
## GarageArea 2.220446e-16 1.110223e-16 5.551115e-17 1.000000e+00expand(lu(precision))## $L
## 4 x 4 Matrix of class "dtrMatrix" (unitriangular)
## [,1] [,2] [,3] [,4]
## [1,] 1.00000000 . . .
## [2,] -0.52746808 1.00000000 . .
## [3,] -0.16914712 -0.07556672 1.00000000 .
## [4,] -0.32359799 -0.45815781 -0.27831633 1.00000000
##
## $U
## 4 x 4 Matrix of class "dtrMatrix"
## [,1] [,2] [,3] [,4]
## [1,] 2.79202782 -1.47270555 -0.47226347 -0.90349460
## [2,] . 1.30700221 -0.09876586 -0.59881327
## [3,] . . 1.08396382 -0.30168483
## [4,] . . . 1.00000000
##
## $P
## 4 x 4 sparse Matrix of class "pMatrix"
##
## [1,] | . . .
## [2,] . | . .
## [3,] . . | .
## [4,] . . . |Calculus-Based Probability & Statistics
For this part, we will use the LotFrontage variable, as it is right-skewed and above 0 for all values (though there are NA values we need to handle)
# Remove NA values
lotfront <- houses[which(!is.na(houses$LotFrontage)),]$LotFrontage
# Fit an exponential distribution
lotfront.exp <- fitdistr(lotfront,"exponential")
lotfront.exp## rate
## 0.0173540949
## (0.0004541774)Now we compare the two visually using histograms:
# Generate 1000 samples from the exponential distribution
y <- rexp(1000,lotfront.exp$estimate)# Plot observed and simulated values
tibble(lotfront) %>% ggplot(aes(x=lotfront)) +
geom_histogram(aes(y=..density..,fill="Observed"),bins = 50,alpha=0.4) +
geom_histogram(data=tibble(y),aes(x=y,y=..density..,fill="Exponential"),
bins=50,alpha=0.4) +
theme_classic() +
theme(legend.title = element_blank()) +
scale_fill_brewer(palette="Set1") +
labs(title = "Lot Frontage", subtitle = "Observed vs. Modeled (Exponential)",
x="Feet", y="Proportion")
There is some similarity, but the optimal exponential distribution doesn’t fit terribly well with what we’re observing.
Comparing the the \(5^{\text{th}}\) and \(95^{\text{th}}\) percentiles of the exponential distribution with the observed values:
# Exponential
round(qexp(c(0.05,0.95),lotfront.exp$estimate),2)## [1] 2.96 172.62# Observed
quantile(lotfront,probs=c(0.05,0.95))## 5% 95%
## 0 104# Mean of distribution
paste("Sample mean (model): ",round(mean(y),2),sep='')## [1] "Sample mean (model): 55.9"# Mean of observed data
paste("Sample mean (observed): ",round(mean(lotfront),2),sep='')## [1] "Sample mean (observed): 57.62"# 95% confidence interval for sample mean (assuming normality)
z <- qnorm(0.95)
a <- z * sd(lotfront)/sqrt(length(lotfront))
paste("CI for Population Mean: ",round(mean(lotfront - a),2)," - ",
round(mean(lotfront + a),2),sep='')## [1] "CI for Population Mean: 56.13 - 59.12"The model’s sample mean is just outside the 95% confidence interval for the mean. This, along with the disparity in the histograms leads me to believe that this variable is not well modeled by an exponential distribution.
Modeling
# Ensure we have only complete cases (no NAs)
train <- houses %>% dplyr::select(-Id) %>%
filter(complete.cases(.))# Set the seed for random numbers
set.seed(605)
cv <- trainControl(method="cv", number = 10, preProc = c("scale"))
# Fit an ElasticNet Model
houses.glm <- train(SalePrice ~ ., data = train,
method = "glmnet", trControl = cv,
tuneLength=10)# Read in the testing data
test <- read_csv("test.csv")
# Process testing data the same way
# we did the training data
test$TotalSqFt <- test$`1stFlrSF` + test$`2ndFlrSF`
# Change character columns to factors
test[fac] <- lapply(test[fac], factor)
# Replace NAs with "No"
test[fac] <- lapply(test[fac], fct_explicit_na,'No')
# Combine levels in test set and training set
for(c in fac){
levels(test[[c]]) <- fct_c(test[[c]],houses[[c]])}
# Replace numeric NAs with 0
test <- test %>% mutate_if(is.numeric, ~replace(., is.na(.), 0))
# Predict Values
pred <- predict(houses.glm,test)# Export submission
submission <- tibble(ID = test$Id, SalePrice=pred)
write_csv(submission, "submission.csv")This scored a 0.20596 in the official Kaggle scoring (user name: lysanthus)