Boston Housing Case Study
###Goal: To Search for the best prediction model using evaluation criteria RMS(P)E for testing Data
##Exploring Boston Housing Data:
###Packages
library(MASS)
library(glmnet)## Warning: package 'glmnet' was built under R version 4.3.2## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 4.3.2## Loaded glmnet 4.1-8library(plotmo)## Warning: package 'plotmo' was built under R version 4.3.2## Loading required package: Formula## Loading required package: plotrix## Warning: package 'plotrix' was built under R version 4.3.2## Loading required package: TeachingDemos## Warning: package 'TeachingDemos' was built under R version 4.3.2library(rpart)
library(rpart.plot)## Warning: package 'rpart.plot' was built under R version 4.3.2library(boot)
data(Boston)
attach(Boston)
dim(Boston)## [1] 506 14str(Boston)## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : int 0 0 0 0 0 0 0 0 0 0 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...summary(Boston)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00##Linear Model
###Random Split Data 80%/20% using #M Number
set.seed(15848936)
sample_index <- sample(nrow(Boston),nrow(Boston)*0.80)
Boston_train <- Boston[sample_index,]
Boston_test <- Boston[-sample_index,]
## Building Linear Model Using all the parameter
set.seed(15848936)
sample_index <- sample(nrow(Boston),nrow(Boston)*0.80)
Boston_train <- Boston[sample_index,]
Boston_test <- Boston[-sample_index,]
##linear model
#Training data
Boston_lm <- lm(medv ~ ., data = Boston_train)
summary(Boston_lm)##
## Call:
## lm(formula = medv ~ ., data = Boston_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.758 -2.752 -0.555 1.866 25.618
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.004841 5.716033 6.824 3.40e-11 ***
## crim -0.139143 0.039070 -3.561 0.000415 ***
## zn 0.049269 0.015171 3.248 0.001264 **
## indus 0.009785 0.070082 0.140 0.889034
## chas 2.153769 0.975097 2.209 0.027771 *
## nox -17.574187 4.269927 -4.116 4.71e-05 ***
## rm 3.544327 0.470513 7.533 3.49e-13 ***
## age -0.001806 0.015088 -0.120 0.904778
## dis -1.515124 0.229429 -6.604 1.31e-10 ***
## rad 0.321394 0.077346 4.155 4.00e-05 ***
## tax -0.012623 0.004348 -2.903 0.003908 **
## ptratio -0.907784 0.148569 -6.110 2.41e-09 ***
## black 0.006198 0.003440 1.802 0.072361 .
## lstat -0.534372 0.058395 -9.151 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.893 on 390 degrees of freedom
## Multiple R-squared: 0.7173, Adjusted R-squared: 0.7079
## F-statistic: 76.11 on 13 and 390 DF, p-value: < 2.2e-16p_lm <- predict(Boston_lm)
In_sample_ASE<-mean((Boston_train$medv - p_lm)^2)
MSE<- (sum((Boston_train$medv - p_lm)^2))/(dim(Boston_train)[1]-13-1)
#Out of sample performance
p_lm_test <- predict(Boston_lm, newdata = Boston_test)
MSPE_out_of_sample <- mean((Boston_test$medv - p_lm_test)^2)
## **For linear Model**:
##MSE
MSE## [1] 23.94085##In sample ASE
In_sample_ASE## [1] 23.11121## Out of Sample MSPE
MSPE_out_of_sample## [1] 17.81081##RMSPE
sqrt(MSPE_out_of_sample)## [1] 4.220285##linear model,with only 11 significant predictors, derived from above model summary
#Training data
Boston_lm <- lm(medv ~ crim+zn+chas+nox+rm+dis+rad+tax+ptratio+black+lstat, data = Boston_train)
summary(Boston_lm)##
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad +
## tax + ptratio + black + lstat, data = Boston_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.7245 -2.7669 -0.5492 1.8744 25.5603
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.020730 5.678406 6.872 2.50e-11 ***
## crim -0.139294 0.038937 -3.577 0.00039 ***
## zn 0.049297 0.014982 3.290 0.00109 **
## chas 2.163859 0.968105 2.235 0.02597 *
## nox -17.546979 3.942928 -4.450 1.12e-05 ***
## rm 3.525591 0.456572 7.722 9.66e-14 ***
## dis -1.514299 0.211342 -7.165 3.88e-12 ***
## rad 0.319545 0.075056 4.257 2.59e-05 ***
## tax -0.012425 0.004033 -3.081 0.00221 **
## ptratio -0.905022 0.146015 -6.198 1.45e-09 ***
## black 0.006136 0.003411 1.799 0.07281 .
## lstat -0.536057 0.054301 -9.872 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.881 on 392 degrees of freedom
## Multiple R-squared: 0.7173, Adjusted R-squared: 0.7093
## F-statistic: 90.4 on 11 and 392 DF, p-value: < 2.2e-16p_lm <- predict(Boston_lm)
In_sample_ASE<-mean((Boston_train$medv - p_lm)^2)
MSE<- (sum((Boston_train$medv - p_lm)^2)/(dim(Boston_train)[1]-11-1))
#Out of sample performance
p_lm_test <- predict(Boston_lm, newdata = Boston_test)
MSPE_out_of_sample <- mean((Boston_test$medv - p_lm_test)^2)
##linear Model:
In_sample_ASE## [1] 23.11322MSE## [1] 23.82076MSPE_out_of_sample## [1] 17.82005##**cross validation**: To comparison in terms of MSPE is more fair, Less subject to the Randomness (boot package used for cross validation)
set.seed(15848936)
sample_index <- sample(nrow(Boston),nrow(Boston)*0.80)
Boston_train <- Boston[sample_index,]
Boston_test <- Boston[-sample_index,]
## split data into folds
boston_glm <- glm(medv~., data = Boston)
cv5 <- cv.glm(boston_glm, data = Boston, K=5)$delta[2]
##split the data into n folds
cv<- cv.glm(boston_glm, data = Boston, K=dim(Boston)[1])$delta[2]
cv5## [1] 23.89094cv## [1] 23.72388##**cross validation**:
## split data into folds
boston_glm <- glm(medv~., data = Boston)
cv5 <- cv.glm(boston_glm, data = Boston, K=5)$delta[2]
##split the data into n folds
cv<- cv.glm(boston_glm, data = Boston, K=dim(Boston)[1])$delta[2]
cv5## [1] 24.23142cv## [1] 23.72388##if we run the model again cv5 changes, but cv remains same.##Regression Tree
## rpart with default cp=0.01
Boston_rpart <- rpart(formula = medv ~ ., data = Boston_train)
prp(Boston_rpart,digits = 4, extra = 1)
## rpart with small cp for bushy tree and then prune using plotcp
Boston_largetree <- rpart(formula = medv ~ ., data = Boston_train, cp = 0.001)
prp(Boston_largetree,digits = 4, extra = 1)
plotcp(Boston_largetree)
# choose the leftmost value for which the value lies below the horizontal line
Boston_rpart <- prune(Boston_largetree, cp=0.0078)
prp(Boston_rpart,digits = 4, extra = 1)
##ASE and MSPE after pruing the model
Boston_train_pred_tree = predict(Boston_rpart)
ASE.rpart<-mean((Boston_train_pred_tree - Boston_train$medv)^2)
Boston_test_pred_tree = predict(Boston_rpart,Boston_test)
mspe.rpart<-mean((Boston_test_pred_tree - Boston_test$medv)^2)
ASE.rpart## [1] 15.25867mspe.rpart## [1] 16.88838##k-NN with optimal k
##knn for classification
library(FNN)## Warning: package 'FNN' was built under R version 4.3.2set.seed(15848936)
sample_index <- sample(nrow(Boston),nrow(Boston)*0.80)
Boston_train <- Boston[sample_index,]
Boston_test <- Boston[-sample_index,]
train.norm <- Boston_train
test.norm <- Boston_test
## normalize numerical predictors to 0-1 scale
## for testing and training dataset
## range [0,1]-standardization ##
cols <- colnames(train.norm[, -14]) #scaling only on p=13 predictors X
for (j in cols) {
train.norm[[j]] <- (train.norm[[j]] - min(Boston_train[[j]])) / (max(Boston_train[[j]]) - min(Boston_train[[j]]))
test.norm[[j]] <- (test.norm[[j]] - min(Boston_train[[j]])) / (max(Boston_train[[j]]) - min(Boston_train[[j]]))
}
## Function scale() may be still used for standardization without using a loop. Please make sure that you are using the scaling parameters from TRAING data only ##
#compute the statistics of scaled train and test data
summary(train.norm)## crim zn indus chas
## Min. :0.0000000 Min. :0.000 Min. :0.0000 Min. :0.00000
## 1st Qu.:0.0008162 1st Qu.:0.000 1st Qu.:0.1695 1st Qu.:0.00000
## Median :0.0026273 Median :0.000 Median :0.3383 Median :0.00000
## Mean :0.0367664 Mean :0.125 Mean :0.3878 Mean :0.07178
## 3rd Qu.:0.0376873 3rd Qu.:0.200 3rd Qu.:0.6466 3rd Qu.:0.00000
## Max. :1.0000000 Max. :1.000 Max. :1.0000 Max. :1.00000
## nox rm age dis
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.00000
## 1st Qu.:0.1296 1st Qu.:0.4508 1st Qu.:0.3950 1st Qu.:0.08714
## Median :0.3086 Median :0.5109 Median :0.7590 Median :0.19184
## Mean :0.3483 Mean :0.5265 Mean :0.6631 Mean :0.24629
## 3rd Qu.:0.4918 3rd Qu.:0.5893 3rd Qu.:0.9361 3rd Qu.:0.37517
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.00000
## rad tax ptratio black
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.1304 1st Qu.:0.1718 1st Qu.:0.4468 1st Qu.:0.9489
## Median :0.1739 Median :0.2729 Median :0.6489 Median :0.9875
## Mean :0.3614 Mean :0.4152 Mean :0.6087 Mean :0.9139
## 3rd Qu.:0.4783 3rd Qu.:0.9141 3rd Qu.:0.8085 3rd Qu.:0.9984
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
## lstat medv
## Min. :0.0000 Min. : 5.00
## 1st Qu.:0.1338 1st Qu.:17.77
## Median :0.2523 Median :21.70
## Mean :0.2931 Mean :23.02
## 3rd Qu.:0.4076 3rd Qu.:26.25
## Max. :1.0000 Max. :50.00summary(test.norm)## crim zn indus chas
## Min. :0.0001018 Min. :0.00000 Min. :0.03372 Min. :0.00000
## 1st Qu.:0.0009933 1st Qu.:0.00000 1st Qu.:0.19795 1st Qu.:0.00000
## Median :0.0039148 Median :0.00000 Median :0.34604 Median :0.00000
## Mean :0.0555069 Mean :0.06863 Mean :0.40558 Mean :0.05882
## 3rd Qu.:0.0561921 3rd Qu.:0.00000 3rd Qu.:0.64663 3rd Qu.:0.00000
## Max. :0.5746830 Max. :0.90000 Max. :1.00000 Max. :1.00000
## nox rm age dis
## Min. :0.01852 Min. :0.1836 Min. :0.0381 Min. :0.004056
## 1st Qu.:0.15329 1st Qu.:0.4356 1st Qu.:0.5247 1st Qu.:0.090478
## Median :0.31481 Median :0.4886 Median :0.8023 Median :0.177295
## Mean :0.35278 Mean :0.5034 Mean :0.7289 Mean :0.226899
## 3rd Qu.:0.55144 3rd Qu.:0.5537 3rd Qu.:0.9475 3rd Qu.:0.320697
## Max. :0.79218 Max. :0.9224 Max. :1.0000 Max. :0.735962
## rad tax ptratio black
## Min. :0.0000 Min. :0.001908 Min. :0.04255 Min. :-0.005578
## 1st Qu.:0.1304 1st Qu.:0.186546 1st Qu.:0.61702 1st Qu.: 0.887374
## Median :0.1739 Median :0.272901 Median :0.75000 Median : 0.982872
## Mean :0.4126 Mean :0.450082 Mean :0.67939 Mean : 0.835058
## 3rd Qu.:1.0000 3rd Qu.:0.914122 3rd Qu.:0.80851 3rd Qu.: 0.997826
## Max. :1.0000 Max. :1.000000 Max. :0.91489 Max. : 1.000000
## lstat medv
## Min. :0.03339 Min. : 5.00
## 1st Qu.:0.17646 1st Qu.:13.95
## Median :0.29484 Median :19.75
## Mean :0.33450 Mean :20.62
## 3rd Qu.:0.45054 3rd Qu.:23.90
## Max. :0.97268 Max. :50.00summary(train.norm$rm)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.4508 0.5109 0.5265 0.5893 1.0000summary(train.norm$lstat)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.1338 0.2523 0.2931 0.4076 1.0000summary(test.norm$rm)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1836 0.4356 0.4886 0.5034 0.5537 0.9224summary(test.norm$lstat)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.03339 0.17646 0.29484 0.33450 0.45054 0.97268#in-sample
Boston.knn.train <- knn.reg(train = train.norm[, 1:13],
test = train.norm[, 1:13],
y = train.norm$medv,
k = 5)
Boston.knn.train <- data.frame(cbind(pred = Boston.knn.train$pred, actual=train.norm$medv))
#calculate in-sample ASE
ASE<-sum((Boston_train$medv-Boston.knn.train$pred)^2)/length(Boston_train$medv)
ASE## [1] 15.05862Boston.knn.reg <- knn.reg(train = train.norm[, 1:13],
test = test.norm[, 1:13],
y = train.norm$medv,
k = 5)
# compute the actual and predicted values and view the first 20 records
Boston.results <- data.frame(cbind(pred = Boston.knn.reg$pred, actual = Boston_test$medv))
# Out-of-Sample Testing
# calculate MSPE (root mean (average) sum squared prediction errors)
MSPE <- sum((Boston_test$medv-Boston.results$pred)^2)/length(Boston_test$medv)
MSPE## [1] 17.93792#Select the best k
##################
set.seed(15848936)
sample_index2 <- sample(nrow(Boston_train),nrow(Boston_train)*0.80)
train2.norm <- train.norm[sample_index2,]
valid.norm <- train.norm[-sample_index2,]
# initialize a data frame with two columns: k and accuracy
RMSE.df <- data.frame(k = seq(1, 30, 1), RMSE.k = rep(0, 30))
# compute knn for different k on validation set
for (i in 1:30) {
knn.reg.pred <- knn.reg(train = train2.norm[, c(1:13)], test = valid.norm[, c(1:13)],
y = train2.norm$medv, k = i)
RMSE.df[i, 2] <- sqrt(sum((valid.norm$medv-knn.reg.pred$pred)^2)/length(valid.norm$medv))
}
RMSE.dfk <- which(RMSE.df[,2] == min(RMSE.df[,2]))
##Optimal K
k## [1] 5#in-sample
Boston.knn.train.best.k <- knn.reg(train = train.norm[, 1:13],
test = train.norm[, 1:13],
y = train.norm$medv,
k = 5)
Boston.knn.best.k <- data.frame(cbind(pred = Boston.knn.train.best.k$pred, actual=train.norm$medv))
#calculate in-sample ASE
ASE<-sum((Boston_train$medv-Boston.knn.best.k$pred)^2)/length(Boston_train$medv)
Boston.knn.reg <- knn.reg(train = train.norm[, 1:13],
test = test.norm[, 1:13],
y = train.norm$medv,
k = 5)
Boston.results <- data.frame(cbind(pred = Boston.knn.reg$pred, actual = Boston_test$medv))
# Out-of-Sample Testing
# calculate MSPE (root mean (average) sum squared prediction errors)
MSPE <- sum((Boston_test$medv-Boston.results$pred)^2)/length(Boston_test$medv)
##Model knn results For optimal K =5
##ASE
ASE## [1] 15.05862##MSPE
MSPE## [1] 17.93792##Random Forests
library(randomForest)## Warning: package 'randomForest' was built under R version 4.3.3## randomForest 4.7-1.1## Type rfNews() to see new features/changes/bug fixes.set.seed(15848936)
sample_index <- sample(nrow(Boston),nrow(Boston)*0.80)
Boston_train <- Boston[sample_index,]
Boston_test <- Boston[-sample_index,]
Boston_rf <- randomForest(medv~.,
data = Boston_train, importance = TRUE)
#Make prediction on training sample
Boston_rf_pred_train<- predict(Boston_rf, Boston_train)
mean((Boston_train$medv-Boston_rf_pred_train)^2) ## [1] 2.271846# in-sample ASE (average/mean squared residuals) = 2.2718 by default
#check variable importance
Boston_rf$importance #rm and lstat are top 2## %IncMSE IncNodePurity
## crim 8.0562285 2006.5178
## zn 0.6230269 338.1159
## indus 6.6178477 2097.3390
## chas 0.4507812 144.8626
## nox 8.6715440 2018.8596
## rm 33.5916404 9066.8001
## age 4.1357630 1032.5357
## dis 7.1259865 2190.6192
## rad 1.4931670 261.7385
## tax 3.8282024 1096.7462
## ptratio 5.8370685 1770.4934
## black 1.2844821 589.3379
## lstat 57.0878785 9621.7533# plot out of bag (OOB) errors for each ntree value from 1 to 500
plot(Boston_rf$mse, type='l', col=2, lwd=2, xlab = "ntree", ylab = "OOB Error")
Boston_rf_pred<- predict(Boston_rf, Boston_test)
mean((Boston_test$medv-Boston_rf_pred)^2) ## [1] 7.969715#MSPE is 7.9697 for test data##Boosting
library(gbm)## Warning: package 'gbm' was built under R version 4.3.3## Loaded gbm 2.1.9## This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3Boston_boost<- gbm(formula = medv~.,
data = Boston_train,
distribution = "gaussian",
n.trees = 100,
shrinkage = 0.1,
interaction.depth = 1) #boosting with defaults
summary(Boston_boost)
#In-sample prediction using Boosting tree
Boston_boost_pred_train<- predict(Boston_boost, Boston_train, n.trees = 100) #use default setting
mean((Boston_train$medv-Boston_boost_pred_train)^2) #In-sample ASE = 12.42303## [1] 12.42303#OUT-of-sample prediction using Boosting tree
Boston_boost_pred_test<- predict(Boston_boost, Boston_test, n.trees = 100)
mean((Boston_test$medv-Boston_boost_pred_test)^2) #out-of-sample MSPE = 9.412584## [1] 9.412584##Generalized Additive Model
library(randomForest)
library(xgboost)## Warning: package 'xgboost' was built under R version 4.3.3library(nlme)
library(mgcv)## This is mgcv 1.8-42. For overview type 'help("mgcv-package")'.set.seed(15848936)
sample_index <- sample(nrow(Boston),nrow(Boston)*0.80)
Boston_train <- Boston[sample_index,]
Boston_test <- Boston[-sample_index,]
Boston_gam <- gam(medv ~ s(crim)+s(zn)+s(indus)+chas+s(nox)
+s(rm)+s(age)+s(dis)+rad+s(tax)+s(ptratio)
+s(black)+s(lstat),data=Boston_train)
summary(Boston_gam)##
## Family: gaussian
## Link function: identity
##
## Formula:
## medv ~ s(crim) + s(zn) + s(indus) + chas + s(nox) + s(rm) + s(age) +
## s(dis) + rad + s(tax) + s(ptratio) + s(black) + s(lstat)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.7363 1.3470 16.137 <2e-16 ***
## chas 0.3867 0.7005 0.552 0.581
## rad 0.1344 0.1436 0.936 0.350
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(crim) 7.704 8.522 7.845 < 2e-16 ***
## s(zn) 1.000 1.000 2.685 0.10218
## s(indus) 4.028 4.911 4.029 0.00307 **
## s(nox) 8.809 8.978 15.517 < 2e-16 ***
## s(rm) 5.756 6.939 22.075 < 2e-16 ***
## s(age) 1.000 1.000 0.636 0.42563
## s(dis) 8.915 8.995 9.176 < 2e-16 ***
## s(tax) 6.435 7.474 7.166 < 2e-16 ***
## s(ptratio) 1.224 1.407 24.764 1.25e-07 ***
## s(black) 2.339 2.889 2.440 0.08158 .
## s(lstat) 6.823 7.937 20.781 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.875 Deviance explained = 89.2%
## GCV = 11.919 Scale est. = 10.236 n = 404plot(Boston_gam, shade=TRUE, seWithMean=TRUE, scale=0, pages=1)
#refit gam on Boston training (nonlinear to linear),
##removing zn, age, ptratio variables from non linearity because the edf is 1 approx
Boston_gam_re <- gam(medv ~ s(crim)+zn+s(indus)+chas+s(nox)
+s(rm)+age+s(dis)+rad+s(tax)+ptratio
+s(black)+s(lstat),data=Boston_train)
summary(Boston_gam_re)##
## Family: gaussian
## Link function: identity
##
## Formula:
## medv ~ s(crim) + zn + s(indus) + chas + s(nox) + s(rm) + age +
## s(dis) + rad + s(tax) + ptratio + s(black) + s(lstat)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.123739 3.088084 11.698 < 2e-16 ***
## zn 0.025936 0.015578 1.665 0.0968 .
## chas 0.386755 0.700650 0.552 0.5813
## age -0.009996 0.012513 -0.799 0.4249
## rad 0.129829 0.143131 0.907 0.3650
## ptratio -0.763947 0.130976 -5.833 1.25e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(crim) 7.670 8.505 7.855 < 2e-16 ***
## s(indus) 4.008 4.891 4.063 0.00314 **
## s(nox) 8.826 8.982 15.673 < 2e-16 ***
## s(rm) 5.811 6.995 22.287 < 2e-16 ***
## s(dis) 8.897 8.994 9.084 < 2e-16 ***
## s(tax) 6.430 7.475 7.174 < 2e-16 ***
## s(black) 2.355 2.910 2.430 0.08089 .
## s(lstat) 6.825 7.939 20.735 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.875 Deviance explained = 89.2%
## GCV = 11.921 Scale est. = 10.244 n = 404plot(Boston_gam_re, shade=TRUE, seWithMean=TRUE, scale=0, pages=1)
#refit gam on Boston training (remove nonsignificant variables)
Boston_gam_final <- gam(medv ~ s(crim)+s(indus)+s(nox)
+s(rm)+s(dis)+s(tax)+ptratio
+s(lstat),data=Boston_train)
summary(Boston_gam_final)##
## Family: gaussian
## Link function: identity
##
## Formula:
## medv ~ s(crim) + s(indus) + s(nox) + s(rm) + s(dis) + s(tax) +
## ptratio + s(lstat)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.5506 2.3458 16.434 < 2e-16 ***
## ptratio -0.8479 0.1277 -6.638 1.2e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(crim) 7.484 8.381 8.021 <2e-16 ***
## s(indus) 3.964 4.821 4.443 0.0013 **
## s(nox) 8.861 8.988 15.794 <2e-16 ***
## s(rm) 5.685 6.868 23.238 <2e-16 ***
## s(dis) 8.898 8.995 10.136 <2e-16 ***
## s(tax) 6.743 7.750 7.061 <2e-16 ***
## s(lstat) 6.738 7.866 26.838 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.872 Deviance explained = 88.8%
## GCV = 11.951 Scale est. = 10.461 n = 404plot(Boston_gam_final, shade=TRUE, seWithMean=TRUE, scale=0, pages=1)
#In-sample fit performance
#Average Sum of Squared Error
Boston_gam.mse.ori <- Boston_gam$dev/nrow(Boston_train) #original model
Boston_gam.mse.adj <- Boston_gam_re$dev/nrow(Boston_train) #adjusted model
Boston_gam.mse.train <- Boston_gam_final$dev/nrow(Boston_train) #final model
Boston_gam.mse.ori;Boston_gam.mse.adj;Boston_gam.mse.train## [1] 8.791109## [1] 8.803315## [1] 9.156565#using the predict() function
Pred_gam_ori <- predict(Boston_gam,Boston_train)#final model
mean((Pred_gam_ori - Boston_train$medv)^2)## [1] 8.791109Pred_gam_adj <- predict(Boston_gam_re,Boston_train)#final model
mean((Pred_gam_adj - Boston_train$medv)^2)## [1] 8.803315Pred_gam_train <- predict(Boston_gam_final,Boston_train)#final model
mean((Pred_gam_train - Boston_train$medv)^2)## [1] 9.156565#out of sample performance
Pred_gam_ori_test <- predict(Boston_gam,Boston_test)#final model
mean((Pred_gam_ori_test - Boston_test$medv)^2)## [1] 9.830322Pred_gam_adj_test <- predict(Boston_gam_re, Boston_test)#final model
mean((Pred_gam_adj_test - Boston_test$medv)^2)## [1] 9.863091Pred_gam_test <- predict(Boston_gam_final,Boston_test)
mean((Pred_gam_test - Boston_test$medv)^2)## [1] 9.49374#AIC
AIC(Boston_gam);AIC(Boston_gam_re);AIC(Boston_gam_final)## [1] 2140.759## [1] 2140.899## [1] 2143.894##Neural Networks
library(neuralnet)## Warning: package 'neuralnet' was built under R version 4.3.3set.seed(15848936)
sample_index <- sample(nrow(Boston),nrow(Boston)*0.80)
Boston_train <- Boston[sample_index,]
Boston_test <- Boston[-sample_index,]
p= dim(Boston)[2]-1
#### Neural networks on (unscaled) original Training data and plot ####
f <- as.formula("medv ~ .")
nn_unscaled <- neuralnet(f,data=Boston_train, hidden=c(5,3), linear.output=T)
summary(nn_unscaled)## Length Class Mode
## call 5 -none- call
## response 404 -none- numeric
## covariate 5252 -none- numeric
## model.list 2 -none- list
## err.fct 1 -none- function
## act.fct 1 -none- function
## linear.output 1 -none- logical
## data 14 data.frame list
## exclude 0 -none- NULL
## net.result 1 -none- list
## weights 1 -none- list
## generalized.weights 1 -none- list
## startweights 1 -none- list
## result.matrix 95 -none- numericplot(nn_unscaled)
# in-sample fit (prediction) # ASE = 81.74783
yhat.train <- predict(object = nn_unscaled, newdata = Boston_train)
ASE_nn <- sum((Boston_train$medv - yhat.train)^2)/nrow(Boston_train)
ASE_nn## [1] 81.74783#### MSPE of the above neural network model for Testing data####
## MSPE = 89.1568
pr_nn <- compute(nn_unscaled, Boston_test[,1:p])
MSPE_nn <- sum((Boston_test$medv - pr_nn$net.result)^2)/nrow(Boston_test) ## MSPE of testing set ##
MSPE_nn## [1] 96.15683## initialize scaling training, testing, and new data frames to originals ##
train.norm <- Boston_train
test.norm <- Boston_test
summary(train.norm)## crim zn indus chas
## Min. : 0.00632 Min. : 0.0 Min. : 0.460 Min. :0.00000
## 1st Qu.: 0.07894 1st Qu.: 0.0 1st Qu.: 5.085 1st Qu.:0.00000
## Median : 0.24008 Median : 0.0 Median : 9.690 Median :0.00000
## Mean : 3.27742 Mean : 12.5 Mean :11.039 Mean :0.07178
## 3rd Qu.: 3.35936 3rd Qu.: 20.0 3rd Qu.:18.100 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.0 Max. :27.740 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4480 1st Qu.:5.914 1st Qu.: 41.25 1st Qu.: 2.088
## Median :0.5350 Median :6.228 Median : 76.60 Median : 3.239
## Mean :0.5543 Mean :6.309 Mean : 67.29 Mean : 3.838
## 3rd Qu.:0.6240 3rd Qu.:6.637 3rd Qu.: 93.80 3rd Qu.: 5.255
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 2.52
## 1st Qu.: 4.000 1st Qu.:277.0 1st Qu.:16.80 1st Qu.:376.74
## Median : 5.000 Median :330.0 Median :18.70 Median :391.95
## Mean : 9.312 Mean :404.5 Mean :18.32 Mean :362.94
## 3rd Qu.:12.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.29
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.730 Min. : 5.00
## 1st Qu.: 6.577 1st Qu.:17.77
## Median :10.875 Median :21.70
## Mean :12.350 Mean :23.02
## 3rd Qu.:16.500 3rd Qu.:26.25
## Max. :37.970 Max. :50.00summary(test.norm)## crim zn indus chas
## Min. : 0.01538 Min. : 0.000 Min. : 1.38 Min. :0.00000
## 1st Qu.: 0.09470 1st Qu.: 0.000 1st Qu.: 5.86 1st Qu.:0.00000
## Median : 0.35462 Median : 0.000 Median : 9.90 Median :0.00000
## Mean : 4.94476 Mean : 6.863 Mean :11.52 Mean :0.05882
## 3rd Qu.: 5.00572 3rd Qu.: 0.000 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :51.13580 Max. :90.000 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3940 Min. :4.519 Min. : 6.60 Min. :1.174
## 1st Qu.:0.4595 1st Qu.:5.835 1st Qu.: 53.85 1st Qu.:2.125
## Median :0.5380 Median :6.111 Median : 80.80 Median :3.079
## Mean :0.5565 Mean :6.188 Mean : 73.67 Mean :3.625
## 3rd Qu.:0.6530 3rd Qu.:6.450 3rd Qu.: 94.90 3rd Qu.:4.656
## Max. :0.7700 Max. :8.375 Max. :100.00 Max. :9.223
## rad tax ptratio black
## Min. : 1.00 Min. :188.0 Min. :13.00 Min. : 0.32
## 1st Qu.: 4.00 1st Qu.:284.8 1st Qu.:18.40 1st Qu.:352.48
## Median : 5.00 Median :330.0 Median :19.65 Median :390.14
## Mean :10.49 Mean :422.8 Mean :18.99 Mean :331.85
## 3rd Qu.:24.00 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.04
## Max. :24.00 Max. :711.0 Max. :21.20 Max. :396.90
## lstat medv
## Min. : 2.940 Min. : 5.00
## 1st Qu.: 8.125 1st Qu.:13.95
## Median :12.415 Median :19.75
## Mean :13.852 Mean :20.62
## 3rd Qu.:18.058 3rd Qu.:23.90
## Max. :36.980 Max. :50.00## normalize all numerical variables (X&Y) to 0-1 scale, range [0,1]-standardization ##
cols <- colnames(train.norm[, ]) #scaling both X and Y
for (j in cols) {
train.norm[[j]] <- (train.norm[[j]] - min(Boston_train[[j]])) / (max(Boston_train[[j]]) - min(Boston_train[[j]]))
test.norm[[j]] <- (test.norm[[j]] - min(Boston_train[[j]])) / (max(Boston_train[[j]]) - min(Boston_train[[j]]))
}
#### Neural networks on (scaled) Training data and plot ####
set.seed(15848936)
nn_scaled <- neuralnet(f,data=train.norm, hidden=c(5,3), linear.output=T)
plot(nn_scaled)
# in-sample fit (prediction) # ASE = 4.4869
pr_nn_scaled_train <- compute(nn_scaled, train.norm[,1:p])
## recover the predicted value back to the original response scale ##
pr_nn_org <- pr_nn_scaled_train$net.result*(max(Boston_train$medv)-min(Boston_train$medv))+min(Boston_train$medv)
train_r <- (train.norm$medv)*(max(Boston_train$medv)-min(Boston_train$medv))+min(Boston_train$medv)
ase.norm<-mean((train_r - pr_nn_org)^2) ##in-sample average sum squared error
ase.norm## [1] 4.486885#### MSPE of the above neural network model for Testing data####
pr_nn_scaled_test <- compute(nn_scaled, test.norm[,1:p])
pr_nn_org_t <- pr_nn_scaled_test$net.result*(max(Boston_train$medv)-min(Boston_train$medv))+min(Boston_train$medv)
test_r <- (test.norm$medv)*(max(Boston_train$medv)-min(Boston_train$medv))+min(Boston_train$medv)
mean((test_r - pr_nn_org_t)^2)## [1] 11.09675## MSPE of testing set = 11.09675