Our aim is to predict house values. Before we begin to do any analysis, we should always check whether the dataset has missing value or not, we do so by typing:
taiwan_real_estate <- read.csv("real_estates.csv", row.names=1)
any(is.na(taiwan_real_estate))
## [1] FALSE
Let’s take a look at structure of the data set:
glimpse(taiwan_real_estate)
## Rows: 414
## Columns: 6
## $ house.age <dbl> 32.0, 19.5, 13.3, 13.3, 5.0, 7.1, …
## $ distance.to.the.nearest.MRT.station <dbl> 84.87882, 306.59470, 561.98450, 56…
## $ number.of.convenience.stores <int> 10, 9, 5, 5, 5, 3, 7, 6, 1, 3, 1, …
## $ latitude <dbl> 24.98298, 24.98034, 24.98746, 24.9…
## $ longitude <dbl> 121.5402, 121.5395, 121.5439, 121.…
## $ house.price.of.unit.area <dbl> 37.9, 42.2, 47.3, 54.8, 43.1, 32.1…
Let’s simplify variables’ names:
taiwan_real_estate <- taiwan_real_estate %>%
rename(
house_age_years = house.age,
price_twd_msq = house.price.of.unit.area,
n_convenience = number.of.convenience.stores,
dist_to_mrt_m = distance.to.the.nearest.MRT.station
)
We can also perform binning for “house_age_years”:
#perform binning with specific number of bins
taiwan_real_estate <- taiwan_real_estate %>%
mutate(house_age_cat = cut(house_age_years, breaks=c(0,15,30,45), include.lowest = TRUE, right = FALSE))
Prepare a heatmap with correlation coefficients on it:
M <- cor(taiwan_real_estate[, c("house_age_years", "price_twd_msq", "n_convenience", "dist_to_mrt_m")])
corrplot(M, method = 'number')
Draw a scatter plot of n_convenience vs. price_twd_msq:
ggplot(taiwan_real_estate, aes(x=n_convenience, y=price_twd_msq)) +
geom_point() +
geom_smooth(method="lm", se=FALSE) +
labs(title="Scatter plot of Convenience Stores vs House Price", x="Number of Convenience Stores", y="Price per sqm (TWD)")
## `geom_smooth()` using formula = 'y ~ x'
Draw a scatter plot of house_age_years vs. price_twd_msq:
ggplot(taiwan_real_estate, aes(x=house_age_years, y=price_twd_msq)) +
geom_point() +
geom_smooth(method="lm", se=FALSE) +
labs(title="Scatter plot of House Age vs House Price", x="House Age (years)", y="Price per sqm (TWD)")
## `geom_smooth()` using formula = 'y ~ x'
Draw a scatter plot of distance to nearest MRT station vs. price_twd_msq:
Plot a histogram of price_twd_msq with 10 bins, facet the plot so each house age group gets its own panel:
ggplot(taiwan_real_estate, aes(x=price_twd_msq)) +
geom_histogram(bins=10, fill="blue", color="black", alpha=0.7) +
facet_wrap(~house_age_cat) +
labs(title="Histogram of House Prices", x="Price per sqm (TWD)", y="Count")
Summarize to calculate the mean, sd, median etc. house price/area by house age:
taiwan_real_estate %>%
group_by(house_age_cat) %>%
summarize(
mean_price = mean(price_twd_msq, na.rm = TRUE),
sd_price = sd(price_twd_msq, na.rm = TRUE),
median_price = median(price_twd_msq, na.rm = TRUE)
)
## # A tibble: 3 × 4
## house_age_cat mean_price sd_price median_price
## <fct> <dbl> <dbl> <dbl>
## 1 [0,15) 41.8 14.2 42.6
## 2 [15,30) 32.6 11.4 32.9
## 3 [30,45] 37.7 12.8 38.3
Run a linear regression of price_twd_msq vs. best, but only 1 predictor:
model1 <- lm(price_twd_msq ~ dist_to_mrt_m, data=taiwan_real_estate)
summary(model1)
##
## Call:
## lm(formula = price_twd_msq ~ dist_to_mrt_m, data = taiwan_real_estate)
##
## Residuals:
## Min 1Q Median 3Q Max
## -35.396 -6.007 -1.195 4.831 73.483
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.8514271 0.6526105 70.26 <2e-16 ***
## dist_to_mrt_m -0.0072621 0.0003925 -18.50 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.07 on 412 degrees of freedom
## Multiple R-squared: 0.4538, Adjusted R-squared: 0.4524
## F-statistic: 342.2 on 1 and 412 DF, p-value: < 2.2e-16
We start by displaying the statistical summary of the model using the R function summary():
You can access lots of different aspects of the regression object. To see what’s inside, use names():
names(model1)
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
What do they mean?
Discuss model accuracy:
Model diagnostics:
par(mfrow = c(2, 2))
plot(model1)
The four plots show…
Create the diagnostic plots using ggfortify:
autoplot(model1)
Outliers and high levarage points:
plot(model1, which=5)
Influential values:
# Cook's distance
plot(model1, which=4)
or just plot all of diagnostic plots together:
autoplot(model1, which = 1:6, label.size = 3)
Discussion: The model diagnostics indicate whether our assumptions for linear regression are met. Any deviation from these assumptions could imply that the model may not be the best fit for the data.
We begin by splitting the dataset into two parts, training set and testing set. In this example we will randomly take 75% row in this dataset and put it into the training set, and other 25% row in the testing set:
smp_size <- floor(0.75 * nrow(taiwan_real_estate))
set.seed(12)
train_ind <- sample(seq_len(nrow(taiwan_real_estate)), size=smp_size)
train <- taiwan_real_estate[train_ind, ]
test <- taiwan_real_estate[-train_ind, ]
1st comment: floor() is used to return the largest integer value which is not greater than an individual number, or expression.
2nd comment: set.seed() is used to set the seed of R’s random number generator, this function is used so results from this example can be recreated easily.
Now we have our training set and testing set.
Generally, selecting variables for linear regression is a debatable topic.
There are many methods for variable selecting, namely, forward stepwise selection, backward stepwise selection, etc, some are valid, some are heavily criticized.
I recommend this document: https://www.stat.cmu.edu/~cshalizi/mreg/15/lectures/26/lecture-26.pdf and Gung’s comment: https://stats.stackexchange.com/questions/20836/algorithms-for-automatic-model-selection/20856#20856 if you want to learn more about variable selection process.
If our goal is prediction, it is safer to include all predictors in our model, removing variables without knowing the science behind it usually does more harm than good!!!
We begin to create our multiple linear regression model:
model2 <- lm(price_twd_msq ~ ., data = train)
summary(model2)
##
## Call:
## lm(formula = price_twd_msq ~ ., data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.009 -4.953 -1.296 4.461 75.042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.861e+03 7.542e+03 -0.379 0.704717
## house_age_years -4.209e-01 1.223e-01 -3.442 0.000659 ***
## dist_to_mrt_m -4.558e-03 8.681e-04 -5.251 2.86e-07 ***
## n_convenience 9.826e-01 2.287e-01 4.297 2.34e-05 ***
## latitude 2.505e+02 5.375e+01 4.660 4.74e-06 ***
## longitude -2.755e+01 6.032e+01 -0.457 0.648221
## house_age_cat[15,30) -1.089e+00 1.916e+00 -0.568 0.570124
## house_age_cat[30,45] 5.789e+00 3.577e+00 1.618 0.106655
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.331 on 302 degrees of freedom
## Multiple R-squared: 0.5676, Adjusted R-squared: 0.5575
## F-statistic: 56.62 on 7 and 302 DF, p-value: < 2.2e-16
Discuss the results…
Looking at model summary, we see that variables …. are insignificant, so let’s estimate the model without those variables:
# Assuming 'house_age_years' and 'dist_to_mrt_m' are found insignificant
model3 <- lm(price_twd_msq ~ n_convenience, data = train)
summary(model3)
##
## Call:
## lm(formula = price_twd_msq ~ n_convenience, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -35.573 -7.659 -2.303 6.754 87.252
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.6626 1.1670 23.70 <2e-16 ***
## n_convenience 2.5850 0.2245 11.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.75 on 308 degrees of freedom
## Multiple R-squared: 0.3009, Adjusted R-squared: 0.2986
## F-statistic: 132.6 on 1 and 308 DF, p-value: < 2.2e-16
There are many standards researchers apply for deciding whether a VIF is too large. In some domains, a VIF over 2 is worthy of suspicion. Others set the bar higher, at 5 or 10. Others still will say you shouldn’t pay attention to these at all. Ultimately, the main thing to consider is that small effects are more likely to be “drowned out” by higher VIFs, but this may just be a natural, unavoidable fact with your model.
## Загрузка требуемого пакета: carData
##
## Присоединяю пакет: 'car'
## Следующий объект скрыт от 'package:boot':
##
## logit
## Следующий объект скрыт от 'package:dplyr':
##
## recode
## Следующий объект скрыт от 'package:purrr':
##
## some
## GVIF Df GVIF^(1/(2*Df))
## house_age_years 6.920312 1 2.630649
## dist_to_mrt_m 4.244914 1 2.060319
## n_convenience 1.644746 1 1.282476
## latitude 1.522464 1 1.233882
## longitude 3.014812 1 1.736321
## house_age_cat 7.480049 2 1.653774
Finally we test our model on test dataset:
## [1] 87.09975
Interpret results…
The MSE (Mean Squared Error) on the test set gives us an idea of how well our model is performing in terms of prediction accuracy.
Best subset and stepwise (forward, backward, both) techniques of variable selection can be used to come up with the best linear regression model for the dependent variable medv.
## Subset selection object
## Call: regsubsets.formula(price_twd_msq ~ ., data = train, nbest = 1,
## nvmax = NULL, method = "exhaustive")
## 7 Variables (and intercept)
## Forced in Forced out
## house_age_years FALSE FALSE
## dist_to_mrt_m FALSE FALSE
## n_convenience FALSE FALSE
## latitude FALSE FALSE
## longitude FALSE FALSE
## house_age_cat[15,30) FALSE FALSE
## house_age_cat[30,45] FALSE FALSE
## 1 subsets of each size up to 7
## Selection Algorithm: exhaustive
## house_age_years dist_to_mrt_m n_convenience latitude longitude
## 1 ( 1 ) " " "*" " " " " " "
## 2 ( 1 ) " " "*" "*" " " " "
## 3 ( 1 ) "*" "*" " " "*" " "
## 4 ( 1 ) "*" "*" "*" "*" " "
## 5 ( 1 ) "*" "*" "*" "*" " "
## 6 ( 1 ) "*" "*" "*" "*" " "
## 7 ( 1 ) "*" "*" "*" "*" "*"
## house_age_cat[15,30) house_age_cat[30,45]
## 1 ( 1 ) " " " "
## 2 ( 1 ) " " " "
## 3 ( 1 ) " " " "
## 4 ( 1 ) " " " "
## 5 ( 1 ) " " "*"
## 6 ( 1 ) "*" "*"
## 7 ( 1 ) "*" "*"
## Start: AIC=1638.47
## price_twd_msq ~ 1
##
## Df Sum of Sq RSS AIC
## + dist_to_mrt_m 1 27754.0 33057 1451.5
## + n_convenience 1 18298.4 42512 1529.5
## + longitude 1 17543.7 43267 1535.0
## + latitude 1 16685.2 44125 1541.0
## + house_age_cat 2 5118.7 55692 1615.2
## + house_age_years 1 2037.6 58773 1629.9
## <none> 60811 1638.5
##
## Step: AIC=1451.52
## price_twd_msq ~ dist_to_mrt_m
##
## Df Sum of Sq RSS AIC
## + n_convenience 1 1821.6 31235 1435.9
## + latitude 1 1793.8 31263 1436.2
## + house_age_years 1 1773.7 31283 1436.4
## + house_age_cat 2 1666.5 31390 1439.5
## <none> 33057 1451.5
## + longitude 1 27.3 33029 1453.3
##
## Step: AIC=1435.94
## price_twd_msq ~ dist_to_mrt_m + n_convenience
##
## Df Sum of Sq RSS AIC
## + house_age_years 1 2123.9 29111 1416.1
## + house_age_cat 2 2113.2 29122 1418.2
## + latitude 1 1388.9 29846 1423.8
## <none> 31235 1435.9
## + longitude 1 6.3 31229 1437.9
##
## Step: AIC=1416.11
## price_twd_msq ~ dist_to_mrt_m + n_convenience + house_age_years
##
## Df Sum of Sq RSS AIC
## + latitude 1 1822.97 27288 1398.1
## + house_age_cat 2 871.12 28240 1410.7
## <none> 29111 1416.1
## + longitude 1 32.98 29078 1417.8
##
## Step: AIC=1398.07
## price_twd_msq ~ dist_to_mrt_m + n_convenience + house_age_years +
## latitude
##
## Df Sum of Sq RSS AIC
## + house_age_cat 2 972.81 26315 1390.8
## <none> 27288 1398.1
## + longitude 1 8.02 27280 1400.0
##
## Step: AIC=1390.81
## price_twd_msq ~ dist_to_mrt_m + n_convenience + house_age_years +
## latitude + house_age_cat
##
## Df Sum of Sq RSS AIC
## <none> 26315 1390.8
## + longitude 1 18.161 26297 1392.6
##
## Call:
## lm(formula = price_twd_msq ~ dist_to_mrt_m + n_convenience +
## house_age_years + latitude + house_age_cat, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.913 -5.098 -1.227 4.551 75.387
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.250e+03 1.337e+03 -4.673 4.46e-06 ***
## dist_to_mrt_m -4.274e-03 6.040e-04 -7.076 1.03e-11 ***
## n_convenience 9.871e-01 2.282e-01 4.326 2.06e-05 ***
## house_age_years -4.217e-01 1.221e-01 -3.453 0.000633 ***
## latitude 2.521e+02 5.356e+01 4.708 3.82e-06 ***
## house_age_cat[15,30) -1.024e+00 1.909e+00 -0.536 0.592043
## house_age_cat[30,45] 5.851e+00 3.570e+00 1.639 0.102303
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.319 on 303 degrees of freedom
## Multiple R-squared: 0.5673, Adjusted R-squared: 0.5587
## F-statistic: 66.2 on 6 and 303 DF, p-value: < 2.2e-16
## Start: AIC=1392.6
## price_twd_msq ~ house_age_years + dist_to_mrt_m + n_convenience +
## latitude + longitude + house_age_cat
##
## Df Sum of Sq RSS AIC
## - longitude 1 18.16 26315 1390.8
## <none> 26297 1392.6
## - house_age_cat 2 982.95 27280 1400.0
## - house_age_years 1 1031.51 27329 1402.5
## - n_convenience 1 1607.54 27905 1409.0
## - latitude 1 1891.31 28189 1412.1
## - dist_to_mrt_m 1 2400.71 28698 1417.7
##
## Step: AIC=1390.81
## price_twd_msq ~ house_age_years + dist_to_mrt_m + n_convenience +
## latitude + house_age_cat
##
## Df Sum of Sq RSS AIC
## <none> 26315 1390.8
## - house_age_cat 2 972.8 27288 1398.1
## - house_age_years 1 1035.7 27351 1400.8
## - n_convenience 1 1625.4 27941 1407.4
## - latitude 1 1924.7 28240 1410.7
## - dist_to_mrt_m 1 4348.8 30664 1436.2
##
## Call:
## lm(formula = price_twd_msq ~ house_age_years + dist_to_mrt_m +
## n_convenience + latitude + house_age_cat, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.913 -5.098 -1.227 4.551 75.387
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.250e+03 1.337e+03 -4.673 4.46e-06 ***
## house_age_years -4.217e-01 1.221e-01 -3.453 0.000633 ***
## dist_to_mrt_m -4.274e-03 6.040e-04 -7.076 1.03e-11 ***
## n_convenience 9.871e-01 2.282e-01 4.326 2.06e-05 ***
## latitude 2.521e+02 5.356e+01 4.708 3.82e-06 ***
## house_age_cat[15,30) -1.024e+00 1.909e+00 -0.536 0.592043
## house_age_cat[30,45] 5.851e+00 3.570e+00 1.639 0.102303
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.319 on 303 degrees of freedom
## Multiple R-squared: 0.5673, Adjusted R-squared: 0.5587
## F-statistic: 66.2 on 6 and 303 DF, p-value: < 2.2e-16
## Start: AIC=1638.47
## price_twd_msq ~ 1
##
## Df Sum of Sq RSS AIC
## + dist_to_mrt_m 1 27754.0 33057 1451.5
## + n_convenience 1 18298.4 42512 1529.5
## + longitude 1 17543.7 43267 1535.0
## + latitude 1 16685.2 44125 1541.0
## + house_age_cat 2 5118.7 55692 1615.2
## + house_age_years 1 2037.6 58773 1629.9
## <none> 60811 1638.5
##
## Step: AIC=1451.52
## price_twd_msq ~ dist_to_mrt_m
##
## Df Sum of Sq RSS AIC
## + n_convenience 1 1821.6 31235 1435.9
## + latitude 1 1793.8 31263 1436.2
## + house_age_years 1 1773.7 31283 1436.4
## + house_age_cat 2 1666.5 31390 1439.5
## <none> 33057 1451.5
## + longitude 1 27.3 33029 1453.3
## - dist_to_mrt_m 1 27754.0 60811 1638.5
##
## Step: AIC=1435.94
## price_twd_msq ~ dist_to_mrt_m + n_convenience
##
## Df Sum of Sq RSS AIC
## + house_age_years 1 2123.9 29111 1416.1
## + house_age_cat 2 2113.2 29122 1418.2
## + latitude 1 1388.9 29846 1423.8
## <none> 31235 1435.9
## + longitude 1 6.3 31229 1437.9
## - n_convenience 1 1821.6 33057 1451.5
## - dist_to_mrt_m 1 11277.2 42512 1529.5
##
## Step: AIC=1416.11
## price_twd_msq ~ dist_to_mrt_m + n_convenience + house_age_years
##
## Df Sum of Sq RSS AIC
## + latitude 1 1823.0 27288 1398.1
## + house_age_cat 2 871.1 28240 1410.7
## <none> 29111 1416.1
## + longitude 1 33.0 29078 1417.8
## - house_age_years 1 2123.9 31235 1435.9
## - n_convenience 1 2171.8 31283 1436.4
## - dist_to_mrt_m 1 10564.4 39676 1510.1
##
## Step: AIC=1398.07
## price_twd_msq ~ dist_to_mrt_m + n_convenience + house_age_years +
## latitude
##
## Df Sum of Sq RSS AIC
## + house_age_cat 2 972.8 26315 1390.8
## <none> 27288 1398.1
## + longitude 1 8.0 27280 1400.0
## - n_convenience 1 1706.1 28994 1414.9
## - latitude 1 1823.0 29111 1416.1
## - house_age_years 1 2557.9 29846 1423.8
## - dist_to_mrt_m 1 5494.4 32783 1452.9
##
## Step: AIC=1390.81
## price_twd_msq ~ dist_to_mrt_m + n_convenience + house_age_years +
## latitude + house_age_cat
##
## Df Sum of Sq RSS AIC
## <none> 26315 1390.8
## + longitude 1 18.2 26297 1392.6
## - house_age_cat 2 972.8 27288 1398.1
## - house_age_years 1 1035.7 27351 1400.8
## - n_convenience 1 1625.4 27941 1407.4
## - latitude 1 1924.7 28240 1410.7
## - dist_to_mrt_m 1 4348.8 30664 1436.2
##
## Call:
## lm(formula = price_twd_msq ~ dist_to_mrt_m + n_convenience +
## house_age_years + latitude + house_age_cat, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.913 -5.098 -1.227 4.551 75.387
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.250e+03 1.337e+03 -4.673 4.46e-06 ***
## dist_to_mrt_m -4.274e-03 6.040e-04 -7.076 1.03e-11 ***
## n_convenience 9.871e-01 2.282e-01 4.326 2.06e-05 ***
## house_age_years -4.217e-01 1.221e-01 -3.453 0.000633 ***
## latitude 2.521e+02 5.356e+01 4.708 3.82e-06 ***
## house_age_cat[15,30) -1.024e+00 1.909e+00 -0.536 0.592043
## house_age_cat[30,45] 5.851e+00 3.570e+00 1.639 0.102303
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.319 on 303 degrees of freedom
## Multiple R-squared: 0.5673, Adjusted R-squared: 0.5587
## F-statistic: 66.2 on 6 and 303 DF, p-value: < 2.2e-16
## df AIC
## model_forward 8 2272.556
## model_backward 8 2272.556
## model_step 8 2272.556
##
## Call:
## lm(formula = price_twd_msq ~ dist_to_mrt_m + n_convenience +
## house_age_years + latitude + house_age_cat, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.913 -5.098 -1.227 4.551 75.387
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.250e+03 1.337e+03 -4.673 4.46e-06 ***
## dist_to_mrt_m -4.274e-03 6.040e-04 -7.076 1.03e-11 ***
## n_convenience 9.871e-01 2.282e-01 4.326 2.06e-05 ***
## house_age_years -4.217e-01 1.221e-01 -3.453 0.000633 ***
## latitude 2.521e+02 5.356e+01 4.708 3.82e-06 ***
## house_age_cat[15,30) -1.024e+00 1.909e+00 -0.536 0.592043
## house_age_cat[30,45] 5.851e+00 3.570e+00 1.639 0.102303
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.319 on 303 degrees of freedom
## Multiple R-squared: 0.5673, Adjusted R-squared: 0.5587
## F-statistic: 66.2 on 6 and 303 DF, p-value: < 2.2e-16
From Best subset regression and stepwise selection (forward, backward, both), we see that …….
Plots show that…
Models are compared based on adjusted r square, AIC, BIC criteria for in-sample performance and mean square prediction error (MSPE) for out-of-sample performance.
## [1] 84.88832
## [1] 0.5672572
## [1] 0.558688
## [1] 2272.556
## [1] 2302.449
Finally, we can check the Out-of-sample Prediction or test error (MSPE):
## [1] 55.28947
Please check how function ?cv.glm works.
We will just extract from this object created by cv.glm command - the raw cross-validation estimate of prediction error.
## [1] NaN NaN
Based on AIC criteria …
We need to check out-of-sample MSPE for both models. Based on out-of-sample prediction error, model …