Introduction

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)
attach(taiwan_real_estate)
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 = T,
                                                                        right = F))

Descriptive Statistics

Prepare a heatmap with correlation coefficients on it:

library(corrplot)
## corrplot 0.92 loaded
M<-cor(taiwan_real_estate[,1:6])
corrplot(M, method = 'number')

Draw a scatter plot of n_convenience vs. price_twd_msq:

Draw a scatter plot of house_age_years vs. price_twd_msq:

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:

Summarize to calculate the mean, sd, median etc. house price/area by house age:

Simple model

Run a linear regression of price_twd_msq vs. best, but only 1 predictor:

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():

What do they mean?

Discuss model accuracy:

Model diagnostics:

#par(mfrow = c(2, 2))
#plot(mdl_price_vs_age) # ?your name of the model

The four plots show…

Create the diagnostic plots using ggfortify:

#library(ggfortify)
#autoplot(model_lm1) # ?your name of the model

Outliers and high levarage points:

#plot(model_lm1, 5)  # ?your name of the model

Influential values:

# Cook's distance
#plot(model_lm1, 4) # ?your name of the model

or just plot all of diagnostic plots together:

#autoplot(model_lm1, which = 1:6, label.size = 3) # ?your name of the model

Discussion:

Multiple Regression Model

Test and training set

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.

Variable selection methods

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:

model3 <- lm(taiwan_real_estate$price_twd_msq ~ taiwan_real_estate$house_age_years + taiwan_real_estate$dist_to_mrt_m + taiwan_real_estate$n_convenience + taiwan_real_estate$house_age_cat + taiwan_real_estate$latitude, data = train)

summary(model3)
## 
## Call:
## lm(formula = taiwan_real_estate$price_twd_msq ~ taiwan_real_estate$house_age_years + 
##     taiwan_real_estate$dist_to_mrt_m + taiwan_real_estate$n_convenience + 
##     taiwan_real_estate$house_age_cat + taiwan_real_estate$latitude, 
##     data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.606  -5.002  -1.105   4.400  76.521 
## 
## Coefficients:
##                                           Estimate Std. Error t value Pr(>|t|)
## (Intercept)                             -6.109e+03  1.102e+03  -5.545 5.30e-08
## taiwan_real_estate$house_age_years      -4.288e-01  1.014e-01  -4.227 2.92e-05
## taiwan_real_estate$dist_to_mrt_m        -3.868e-03  4.961e-04  -7.798 5.33e-14
## taiwan_real_estate$n_convenience         1.132e+00  1.884e-01   6.008 4.18e-09
## taiwan_real_estate$house_age_cat[15,30) -5.960e-01  1.557e+00  -0.383   0.7020
## taiwan_real_estate$house_age_cat[30,45]  5.390e+00  2.939e+00   1.834   0.0674
## taiwan_real_estate$latitude              2.464e+02  4.412e+01   5.585 4.28e-08
##                                            
## (Intercept)                             ***
## taiwan_real_estate$house_age_years      ***
## taiwan_real_estate$dist_to_mrt_m        ***
## taiwan_real_estate$n_convenience        ***
## taiwan_real_estate$house_age_cat[15,30)    
## taiwan_real_estate$house_age_cat[30,45] .  
## taiwan_real_estate$latitude             ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.849 on 407 degrees of freedom
## Multiple R-squared:  0.5832, Adjusted R-squared:  0.5771 
## F-statistic: 94.92 on 6 and 407 DF,  p-value: < 2.2e-16

Evaluating multi-collinearity

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.

## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
##                                        GVIF Df GVIF^(1/(2*Df))
## taiwan_real_estate$house_age_years 7.044790  1        2.654202
## taiwan_real_estate$dist_to_mrt_m   2.067441  1        1.437860
## taiwan_real_estate$n_convenience   1.625163  1        1.274819
## taiwan_real_estate$house_age_cat   7.498627  2        1.654800
## taiwan_real_estate$latitude        1.581323  1        1.257507

Finally we test our model on test dataset:

## Warning: 'newdata' had 104 rows but variables found have 414 rows
## Warning in predicted_values - actual_values: longer object length is not a
## multiple of shorter object length

## Warning in predicted_values - actual_values: longer object length is not a
## multiple of shorter object length
## Warning in actual_values - predicted_values: longer object length is not a
## multiple of shorter object length
## [1] "MSE:  250.400057644815"
## [1] "MAE:  12.7279078779504"
## [1] "R-squared:  -5.89487311354564"

Interpret results…

Variable selection using best subset regression

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.

## Warning: package 'leaps' was built under R version 4.3.3

Variable selection using stepwise regression

## Start:  AIC=-26558.51
## house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     n_convenience + latitude + longitude + price_twd_msq + house_age_cat
## Start:  AIC=-26558.51
## house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     n_convenience + latitude + longitude + price_twd_msq + house_age_cat
## Warning: attempting model selection on an essentially perfect fit is nonsense
##                   Df Sum of Sq   RSS      AIC
## - n_convenience    1         0     0 -28119.3
## - dist_to_mrt_m    1         0     0 -27646.4
## - latitude         1         0     0 -27297.5
## - house_age_years  1         0     0 -26708.9
## - longitude        1         0     0 -26679.1
## - house_age_cat    2         0     0 -26561.3
## <none>                             0 -26558.5
## - price_twd_msq    1     31866 31866   1814.2
## 
## Step:  AIC=-28119.32
## house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     latitude + longitude + price_twd_msq + house_age_cat
## Warning: attempting model selection on an essentially perfect fit is nonsense
##                   Df Sum of Sq   RSS      AIC
## - house_age_cat    2         0     0 -28122.5
## <none>                             0 -28119.3
## - longitude        1         0     0 -27159.0
## - dist_to_mrt_m    1         0     0 -26912.9
## - latitude         1         0     0 -26252.8
## - house_age_years  1         0     0 -26057.0
## - price_twd_msq    1     34677 34677   1847.2
## 
## Step:  AIC=-28122.46
## house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     latitude + longitude + price_twd_msq
## Warning: attempting model selection on an essentially perfect fit is nonsense
##                   Df Sum of Sq   RSS      AIC
## <none>                             0 -28122.5
## - longitude        1         0     0 -27162.6
## - dist_to_mrt_m    1         0     0 -26916.6
## - latitude         1         0     0 -26256.6
## - house_age_years  1         0     0 -26059.6
## - price_twd_msq    1     35794 35794   1856.3
## Start:  AIC=-26558.51
## house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     n_convenience + latitude + longitude + price_twd_msq + house_age_cat
## Warning: attempting model selection on an essentially perfect fit is nonsense
##                   Df Sum of Sq   RSS      AIC
## - n_convenience    1         0     0 -28119.3
## - dist_to_mrt_m    1         0     0 -27646.4
## - latitude         1         0     0 -27297.5
## - house_age_years  1         0     0 -26708.9
## - longitude        1         0     0 -26679.1
## - house_age_cat    2         0     0 -26561.3
## <none>                             0 -26558.5
## - price_twd_msq    1     31866 31866   1814.2
## 
## Step:  AIC=-28119.32
## house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     latitude + longitude + price_twd_msq + house_age_cat
## Warning: attempting model selection on an essentially perfect fit is nonsense

## Warning: attempting model selection on an essentially perfect fit is nonsense
##                   Df Sum of Sq   RSS      AIC
## - house_age_cat    2         0     0 -28122.5
## <none>                             0 -28119.3
## + n_convenience    1         0     0 -28117.3
## - longitude        1         0     0 -27159.0
## - dist_to_mrt_m    1         0     0 -26912.9
## - latitude         1         0     0 -26252.8
## - house_age_years  1         0     0 -26057.0
## - price_twd_msq    1     34677 34677   1847.2
## 
## Step:  AIC=-28122.46
## house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     latitude + longitude + price_twd_msq
## Warning: attempting model selection on an essentially perfect fit is nonsense

## Warning: attempting model selection on an essentially perfect fit is nonsense
##                   Df Sum of Sq   RSS      AIC
## <none>                             0 -28122.5
## + n_convenience    1         0     0 -28120.5
## + house_age_cat    2         0     0 -28119.3
## - longitude        1         0     0 -27162.6
## - dist_to_mrt_m    1         0     0 -26916.6
## - latitude         1         0     0 -26256.6
## - house_age_years  1         0     0 -26059.6
## - price_twd_msq    1     35794 35794   1856.3

Comparing competing models

##                df       AIC
## model_forward  10 -25381.63
## model.backward  7 -26945.58
## model.step      7 -26945.58
## Warning in summary.lm(model.step): essentially perfect fit: summary may be
## unreliable
## 
## Call:
## lm(formula = house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     latitude + longitude + price_twd_msq, data = taiwan_real_estate)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -2.604e-15 -2.840e-16 -7.900e-17  1.310e-16  3.422e-14 
## 
## Coefficients:
##                  Estimate Std. Error   t value Pr(>|t|)    
## (Intercept)     0.000e+00  1.222e-12 0.000e+00        1    
## house_age_years 0.000e+00  8.006e-18 0.000e+00        1    
## dist_to_mrt_m   0.000e+00  1.433e-19 0.000e+00        1    
## latitude        0.000e+00  9.132e-15 0.000e+00        1    
## longitude       0.000e+00  9.654e-15 0.000e+00        1    
## price_twd_msq   1.000e+00  9.321e-18 1.073e+17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.763e-15 on 408 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 4.918e+33 on 5 and 408 DF,  p-value: < 2.2e-16

From Best subset regression and stepwise selection (forward, backward, both), we see that …….

## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Start:  AIC=-26558.51
## house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     n_convenience + latitude + longitude + price_twd_msq + house_age_cat
## Warning: attempting model selection on an essentially perfect fit is nonsense
##                   Df Sum of Sq   RSS      AIC
## - n_convenience    1         0     0 -28119.3
## - dist_to_mrt_m    1         0     0 -27646.4
## - latitude         1         0     0 -27297.5
## - house_age_years  1         0     0 -26708.9
## - longitude        1         0     0 -26679.1
## - house_age_cat    2         0     0 -26561.3
## <none>                             0 -26558.5
## - price_twd_msq    1     31866 31866   1814.2
## 
## Step:  AIC=-28119.32
## house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     latitude + longitude + price_twd_msq + house_age_cat
## Warning: attempting model selection on an essentially perfect fit is nonsense

## Warning: attempting model selection on an essentially perfect fit is nonsense
##                   Df Sum of Sq   RSS      AIC
## - house_age_cat    2         0     0 -28122.5
## <none>                             0 -28119.3
## + n_convenience    1         0     0 -28117.3
## - longitude        1         0     0 -27159.0
## - dist_to_mrt_m    1         0     0 -26912.9
## - latitude         1         0     0 -26252.8
## - house_age_years  1         0     0 -26057.0
## - price_twd_msq    1     34677 34677   1847.2
## 
## Step:  AIC=-28122.46
## house.price.of.unit.area ~ house_age_years + dist_to_mrt_m + 
##     latitude + longitude + price_twd_msq
## Warning: attempting model selection on an essentially perfect fit is nonsense

## Warning: attempting model selection on an essentially perfect fit is nonsense
##                   Df Sum of Sq   RSS      AIC
## <none>                             0 -28122.5
## + n_convenience    1         0     0 -28120.5
## + house_age_cat    2         0     0 -28119.3
## - longitude        1         0     0 -27162.6
## - dist_to_mrt_m    1         0     0 -26916.6
## - latitude         1         0     0 -26256.6
## - house_age_years  1         0     0 -26059.6
## - price_twd_msq    1     35794 35794   1856.3

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.

## Warning in summary.lm(model.step): essentially perfect fit: summary may be
## unreliable
## [1] "R-squared: 1"
## Warning in summary.lm(model.step): essentially perfect fit: summary may be
## unreliable
## [1] "Adjusted R-squared: 1"
## [1] "AIC: -26945.5790111926"
## [1] "BIC: -26917.3979493758"

Finally, we can check the Out-of-sample Prediction or test error (MSPE):

if (!require(boot)) {
  install.packages("boot")
}
## Loading required package: boot
## 
## Attaching package: 'boot'
## The following object is masked from 'package:car':
## 
##     logit
library(boot)
  
set.seed(123)  # For reproducibility
#cv_results <- cv.glm(taiwan_real_estate, model.step, K = 10)
#mspe <- cv_results$delta[1]
#print(paste("Mean Square Prediction Error (MSPE):", mspe))

Cross Validation

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.

# Perform cross-validation
#set.seed(123)  # For reproducibility
#cv_results <- cv.glm(taiwan_real_estate, model.step, K = 10)
#print(cv_results)
# Extract the raw cross-validation estimate of prediction error
#cv_error <- cv_results$delta[1]
#print(paste("Cross-validation prediction error:", cv_error))

Based on AIC criteria …

We need to check out-of-sample MSPE for both models. Based on out-of-sample prediction error, model …