Peer Assessment I
First, let us load the data and necessary packages:
load("ames_train.Rdata")
library(MASS)
library(dplyr)
library(ggplot2)
library(tidyverse)
library(BAS)1
Make a labeled histogram (with 30 bins) of the ages of the houses in the data set, and describe the distribution.
ames_train <- ames_train %>% mutate(age=2018-Year.Built)
ggplot(data=ames_train, aes(x=age)) + geom_histogram(aes(y =..density..), bins=30) + geom_density(col="red") +
geom_vline(aes(xintercept = mean(age)),col='red')+
geom_vline(aes(xintercept = median(age)),col='blue')+
geom_text(data = ames_train, aes( x = (mean(ames_train$age)), y = .015, label = 'mean', col='red'), size = 3, parse = T) +
geom_text(data = ames_train,aes( x = (median(ames_train$age)),y = .020, label = 'median'), col='blue', size = 3, parse = T)
Density distribution shows that it is multimodal model and right skewed.
Skewness describes that houses with more age can be a bit low in number
Summary can be found below
ames_train %>% summarise(mean = mean(age),
median = median(age),
sd = sd(age),
min = min(age),
max = max(age),
IQR = IQR(age),
total = n())## # A tibble: 1 x 7
## mean median sd min max IQR total
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 45.8 43 29.6 8 146 46 10002
The mantra in real estate is “Location, Location, Location!” Make a graphical display that relates a home price to its neighborhood in Ames, Iowa. Which summary statistics are most appropriate to use for determining the most expensive, least expensive, and most heterogeneous (having the most variation in housing price) neighborhoods? Report which neighborhoods these are based on the summary statistics of your choice. Report the value of your chosen summary statistics for these neighborhoods.
The distribution of the price rance as in each neighborhood can be visualized using boxplot
neighborhood_data <- ames_train %>% dplyr::select(price, Neighborhood)
neighborhood_summary <- neighborhood_data %>% group_by(Neighborhood) %>% summarise(mean=mean(price), median=median(price), mode=mode(price), iqr=IQR(price), sd=sd(price), var=var(price), total=n())
ggplot(data=ames_train, aes(x = Neighborhood, y = price)) + geom_boxplot() + labs(title = "Housing prices per Neighborhood", x = 'Neighborhood', y = "Price") + theme(axis.text.x = element_text(angle = 90, hjust = 1))
2.0.0.1 Expensive Neighborhood
It can be found by comparing the median of each neighborhood with the max of the all medians of all neighborhood
expense_neighborhood <- neighborhood_summary[which(neighborhood_summary$median == max(neighborhood_summary$median)),]
expense_neighborhood## # A tibble: 1 x 8
## Neighborhood mean median mode iqr sd var total
## <fct> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <int>
## 1 StoneBr 339316. 340692. numeric 151358 123459. 15242150036. 202.0.0.2 Least Expense Neighborhood
It can be found by comparing the median of each neighborhood with the min of the all medians of all neighborhood
cheap_neighborhood <- neighborhood_summary[which(neighborhood_summary$median == min(neighborhood_summary$median)),]
cheap_neighborhood## # A tibble: 1 x 8
## Neighborhood mean median mode iqr sd var total
## <fct> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <int>
## 1 MeadowV 92947. 85750 numeric 20150 18940. 358715156. 162.0.0.3 Heterogeneous Neighborhood
It can be found by comparing the median of each neighborhood with the min of the all medians of all neighborhood
het_neighborhood <- neighborhood_summary[which(neighborhood_summary$sd == max(neighborhood_summary$sd)),]
het_neighborhood## # A tibble: 1 x 8
## Neighborhood mean median mode iqr sd var total
## <fct> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <int>
## 1 StoneBr 339316. 340692. numeric 151358 123459. 15242150036. 20So Expense Neighborhood is StoneBr, least expensive neighborhood is MeadowV and heterogeneous neighborhood is StoneBr
3
Which variable has the largest number of missing values? Explain why it makes sense that there are so many missing values for this variable.
max_na = ""
na_count = 0
for(col in names(ames_train)){
dat = as.vector(summary(ames_train[col]))
for(value in dat){
if(startsWith(value, "NA's")){
tempcount = as.numeric(trimws(unlist(strsplit(value,":"))[2]))
if(tempcount > na_count){
na_count = tempcount
max_na = col
}
}
}
}
max_na## [1] "Pool.QC"Pool.QC (pool quality) has the higest number of NA’s (i.e. house with no pool). As Iowa has extreme winters and changing seasons, it is accepted that many houses won’t have a pool in IOWA.
4
We want to predict the natural log of the home prices. Candidate explanatory variables are lot size in square feet (Lot.Area), slope of property (Land.Slope), original construction date (Year.Built), remodel date (Year.Remod.Add), and the number of bedrooms above grade (Bedroom.AbvGr). Pick a model selection or model averaging method covered in the Specialization, and describe how this method works. Then, use this method to find the best multiple regression model for predicting the natural log of the home prices.
ames_train <- ames_train %>% mutate(log_price=log(price))
model <- bas.lm(log_price ~ Lot.Area + Land.Slope + Year.Built + Year.Remod.Add + Bedroom.AbvGr, prior="BIC", modelprior=uniform(), data=ames_train)
summary(model)## P(B != 0 | Y) model 1 model 2 model 3
## Intercept 1.0000000 1.0000 1.0000000 1.000000e+00
## Lot.Area 1.0000000 1.0000 1.0000000 1.000000e+00
## Land.SlopeMod 0.9039993 1.0000 0.0000000 1.000000e+00
## Land.SlopeSev 0.9039993 1.0000 0.0000000 1.000000e+00
## Year.Built 1.0000000 1.0000 1.0000000 1.000000e+00
## Year.Remod.Add 1.0000000 1.0000 1.0000000 1.000000e+00
## Bedroom.AbvGr 1.0000000 1.0000 1.0000000 0.000000e+00
## BF NA 1.0000 0.1061955 5.347284e-13
## PostProbs NA 0.9040 0.0960000 0.000000e+00
## R2 NA 0.5625 0.5544000 5.338000e-01
## dim NA 7.0000 5.0000000 6.000000e+00
## logmarg NA -2198.1673 -2200.4098178 -2.226424e+03
## model 4 model 5
## Intercept 1.000000e+00 1.000000e+00
## Lot.Area 1.000000e+00 0.000000e+00
## Land.SlopeMod 0.000000e+00 1.000000e+00
## Land.SlopeSev 0.000000e+00 1.000000e+00
## Year.Built 1.000000e+00 1.000000e+00
## Year.Remod.Add 1.000000e+00 1.000000e+00
## Bedroom.AbvGr 0.000000e+00 1.000000e+00
## BF 6.841724e-14 2.396619e-17
## PostProbs 0.000000e+00 0.000000e+00
## R2 5.254000e-01 5.244000e-01
## dim 4.000000e+00 6.000000e+00
## logmarg -2.228480e+03 -2.236437e+03We use Bayesian Model Averaging to simulate all possible models. The posterior probability of each model is based on its BIC score which rewards models for accuracy of prediction but penalizes them for each additional variable.
image(model,rotate = FALSE)
The highest probability model is shown in the image above. It uses all the variables to predict log(price).
5
Which home has the largest squared residual in the previous analysis (Question 4)? Looking at all the variables in the data set, can you explain why this home stands out from the rest (what factors contribute to the high squared residual and why are those factors relevant)?
model_predict <- lm(log_price ~ Lot.Area + Land.Slope + Year.Built + Year.Remod.Add + Bedroom.AbvGr, data=ames_train)
summary(model_predict)##
## Call:
## lm(formula = log_price ~ Lot.Area + Land.Slope + Year.Built +
## Year.Remod.Add + Bedroom.AbvGr, data = ames_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0878 -0.1651 -0.0211 0.1657 0.9945
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.371e+01 8.574e-01 -15.996 < 2e-16 ***
## Lot.Area 1.028e-05 1.106e-06 9.296 < 2e-16 ***
## Land.SlopeMod 1.384e-01 4.991e-02 2.773 0.00565 **
## Land.SlopeSev -4.567e-01 1.514e-01 -3.016 0.00263 **
## Year.Built 6.049e-03 3.788e-04 15.968 < 2e-16 ***
## Year.Remod.Add 6.778e-03 5.468e-04 12.395 < 2e-16 ***
## Bedroom.AbvGr 8.686e-02 1.077e-02 8.063 2.12e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.279 on 993 degrees of freedom
## Multiple R-squared: 0.5625, Adjusted R-squared: 0.5598
## F-statistic: 212.8 on 6 and 993 DF, p-value: < 2.2e-16plot(model_predict,which=1, main="House #428 has highest residual magnitude")
plot(model_predict,which=2, main="Residual Plot")
From the plot, it is obvious that House #428 has highest residual magnitude
ames_train[428,] %>% select(PID, price, log_price, Lot.Area, Land.Slope, Year.Built, Year.Remod.Add, Bedroom.AbvGr)## # A tibble: 1 x 8
## PID price log_price Lot.Area Land.Slope Year.Built Year.Remod.Add
## <int> <int> <dbl> <int> <fct> <int> <int>
## 1 9.02e8 12789 9.46 9656 Gtl 1923 1970
## # ... with 1 more variable: Bedroom.AbvGr <int>Predict the price
predict(model_predict, ames_train[428,])## 1
## 11.54419Diff : 9.46 (original log price) vs 11.54419 (predicted log price)
The difference in original price can be calculated as below
exp(9.46)## [1] 12835.88vs
exp(11.54419)## [1] 103175.8It is the comparison between 103175.8 (predicted) and 12789 (original) price. The difference is huge 90,386.8.
This property has the lowest price in the entire data set. The difference in prediction is due to variables which are not included in the model. Eg:
| Variables | No of occurrences |
|---|---|
Overall.Qual and Overall.Cond are both “poor” |
2 out of 1000 houses |
Sale.Condition is abnormal which is true |
61 out of 1000 houses. |
It is obvious that poor quality, poor condition, and abnormal sales condition will impact price in a very strong and negative way.
6
Use the same model selection method you chose in Question 4 to again find the best multiple regression model to predict the natural log of home prices, but this time replacing Lot.Area with log(Lot.Area). Do you arrive at a model including the same set of predictors?
ames_train <- ames_train %>% mutate(log_larea=log(Lot.Area))
model_larea <- bas.lm(log_price ~ log_larea + Land.Slope + Year.Built + Year.Remod.Add + Bedroom.AbvGr, prior="BIC", modelprior=uniform(), data=ames_train)
summary(model_larea)## P(B != 0 | Y) model 1 model 2 model 3
## Intercept 1.00000000 1.0000 1.000000e+00 1.000000e+00
## log_larea 1.00000000 1.0000 1.000000e+00 1.000000e+00
## Land.SlopeMod 0.02318949 0.0000 1.000000e+00 0.000000e+00
## Land.SlopeSev 0.02318949 0.0000 1.000000e+00 0.000000e+00
## Year.Built 1.00000000 1.0000 1.000000e+00 1.000000e+00
## Year.Remod.Add 1.00000000 1.0000 1.000000e+00 1.000000e+00
## Bedroom.AbvGr 0.99998607 1.0000 1.000000e+00 0.000000e+00
## BF NA 1.0000 2.374021e-02 1.412793e-05
## PostProbs NA 0.9768 2.320000e-02 0.000000e+00
## R2 NA 0.6031 6.056000e-01 5.913000e-01
## dim NA 5.0000 7.000000e+00 4.000000e+00
## logmarg NA -2142.5652 -2.146306e+03 -2.153733e+03
## model 4 model 5
## Intercept 1.000000e+00 1.000000e+00
## log_larea 1.000000e+00 1.000000e+00
## Land.SlopeMod 1.000000e+00 0.000000e+00
## Land.SlopeSev 1.000000e+00 0.000000e+00
## Year.Built 1.000000e+00 1.000000e+00
## Year.Remod.Add 1.000000e+00 0.000000e+00
## Bedroom.AbvGr 0.000000e+00 1.000000e+00
## BF 1.340374e-07 2.253065e-33
## PostProbs 0.000000e+00 0.000000e+00
## R2 5.931000e-01 5.355000e-01
## dim 6.000000e+00 4.000000e+00
## logmarg -2.158390e+03 -2.217738e+03image(model_larea,rotate = FALSE)
We didn’t arrive at the same model as before, but only with the variables log_larea, Year.Built, Year.Remod.Add and Bedroom.AbvGr (excluding Land.Slope). And the R squared increased from 0.5625 to 0.6031 indicating the improvement in the prediction model
7
Do you think it is better to log transform Lot.Area, in terms of assumptions for linear regression? Make graphs of the predicted values of log home price versus the true values of log home price for the regression models selected for Lot.Area and log(Lot.Area). Referencing these two plots, provide a written support that includes a quantitative justification for your answer in the first part of question 7.
model_lprice_larea <- lm(log_price ~ log_larea + Year.Built + Year.Remod.Add + Bedroom.AbvGr, data=ames_train)
summary(model_lprice_larea)##
## Call:
## lm(formula = log_price ~ log_larea + Year.Built + Year.Remod.Add +
## Bedroom.AbvGr, data = ames_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.14609 -0.15825 -0.01477 0.15354 1.01578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.557e+01 8.213e-01 -18.964 < 2e-16 ***
## log_larea 2.471e-01 1.654e-02 14.935 < 2e-16 ***
## Year.Built 5.964e-03 3.604e-04 16.547 < 2e-16 ***
## Year.Remod.Add 6.765e-03 5.197e-04 13.017 < 2e-16 ***
## Bedroom.AbvGr 5.726e-02 1.054e-02 5.434 6.94e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2655 on 995 degrees of freedom
## Multiple R-squared: 0.6031, Adjusted R-squared: 0.6015
## F-statistic: 377.9 on 4 and 995 DF, p-value: < 2.2e-16ames_train$fitted_lprice <- model_predict$fitted.values
ames_train$fitted_lprice_larea <- model_lprice_larea$fitted.valuesames_train %>% ggplot(aes(y=log_price,x=fitted_lprice)) + geom_point() + geom_smooth(method="lm") + ggtitle("Log(price) vs Lot.Area")
ames_train %>% ggplot(aes(y=log_price,x=fitted_lprice_larea)) + geom_point() + geom_smooth(method="lm") + ggtitle("Log(price) vs Log(Lot.Area)")
First plot shows that there are very influential points on the right and that affects the slope of the model. So log of the Lot.Area (i.e.) log(Lot.Area) reduces the influence of such points.
As the area is right skewed, please observe the below 2 plots.
ames_train %>% ggplot(aes(x=Lot.Area,y=log_price)) + geom_smooth(method="lm") + geom_point() + ggtitle("Lot.Area vs log(price)")
ames_train %>% ggplot(aes(x=log_larea,y=log_price)) + geom_smooth(method="lm") + geom_point() + ggtitle("log(Lot.Area) vs log(price)")
It is obvious that Log(Lot.Area) is better than Lot.Area for handling influential points