#problem 1
rm(list=ls())
##This is analysis of California housing data
#install.packages("tree")
library(tree)
## Warning: package 'tree' was built under R version 4.5.3
#install.packages("randomForest")
#install.packages("gbm")
#compare methods for course project
#California housing example from 1990
#each row corresponds to a neighborhood
dat.miss<-read.csv("C:/Users/Sabuj Ganguly/OneDrive/Documents/Ph.D 3rd sem/Regression/Housing data.csv")
#Study missing values
#next line sweeps data frame for NA values
sum(is.na(dat.miss)) #207 missing values in set, a tiny proportion of records
## [1] 207
missingcount<-apply(is.na(dat.miss),1, sum)
#A sophisticated strategy to assess these missing values and account for them in the analysis would
#be worthwhile, but I'll just use casewise deletion (i.e., remove rows w missing values) to keep
#the exposition on the tree-based modeling.
dat<-dat.miss[missingcount==0,]
dat$ocean_proximity<-as.factor(dat$ocean_proximity)
#The outcome variable is median_house_value
#Get a good look at the data dictionary to understand scale these variables are measured on
#https://www.kaggle.com/datasets/camnugent/california-housing-prices?resource=download
#Create a training and holdout set - let's use 80% training 20% holdout
N<-dim(dat)[1]
n<-round(.8*N)
set.seed(73841)
train.ind<-sort(sample(1:N, n, replace=FALSE))
test.ind<- seq(1,N)[-train.ind]
house.t<-dat[train.ind,] #training data
house.h<-dat[test.ind,] #holdout data
#Articulate the modeling plan
#I plan to fit a regression tree and use cost-complexity pruning, bagging, random forests, and
#boosting to train models on the training set, then assess performance on the holdout set.
#To understand how spatial information (latitude and longitude) perform, we conduct the analysis
#both including and excluding these terms and using only those terms
#Let's also consider the implications of e.g., log transforming outcomes
#Once I start analyzing holdout performance, I won't develop new training strategies
#otherwise I am implicitly training on holdout data
#A neat idea might be to make a spatial map of the data then show where lat and long splits
#are (using a tree that only includes latitude and longitude)
#Fit a regression tree to the training data
#Fit a classification tree to the training data
house.tree<-tree(median_house_value~.,data=house.t) #six other observations were omitted.
graphics.off()
par(mfrow = c(1,1))
plot(house.tree)
text(house.tree, pretty = 0)
summary(house.tree)
##
## Regression tree:
## tree(formula = median_house_value ~ ., data = house.t)
## Variables actually used in tree construction:
## [1] "median_income" "ocean_proximity" "longitude" "latitude"
## Number of terminal nodes: 9
## Residual mean deviance: 5.487e+09 = 8.964e+13 / 16340
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -308900 -44520 -12590 0 32260 408800
#Demonstrate that deviance, residual mean deviance are SSE, MSE respectively
ls(summary(house.tree))
## [1] "call" "dev" "df" "residuals" "size" "type"
## [7] "used"
#compute predicted values on the training set
yhat.house.tree.t<-predict(house.tree)
y.t<-house.t$median_house_value
SSE.house.tree.t<-sum((y.t-yhat.house.tree.t)^2)
MSE.house.tree.t<-mean((y.t-yhat.house.tree.t)^2)
SSE.house.tree.t
## [1] 8.963569e+13
MSE.house.tree.t
## [1] 5483646588
#SSE and MSE are the deviance and residual mean deviance in tree output
###############################QUIZ#################################
#SSE and MSE seem "intuitively gigantic" - why do you think that is?
####################################################################
#What if I fit a tree using only spatial information from latitude, longitude
latlon.tree<-tree(median_house_value~latitude + longitude,data=house.t)
plot(latlon.tree)
text(latlon.tree)
summary(latlon.tree)
##
## Regression tree:
## tree(formula = median_house_value ~ latitude + longitude, data = house.t)
## Number of terminal nodes: 13
## Residual mean deviance: 7.399e+09 = 1.208e+14 / 16330
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -351700 -53570 -16320 0 38860 390400
#Much larger training MSE
#in-sample predictions for latlong tree
yhat.latlon.tree.t<-predict(latlon.tree)
#Let's set up a larger tree - do cost complexity pruning
house.big.tree<-tree(median_house_value~.,data=house.t,control=tree.control(nobs = 16346, minsize = 5, mindev=0.0005))
summary(house.big.tree)
##
## Regression tree:
## tree(formula = median_house_value ~ ., data = house.t, control = tree.control(nobs = 16346,
## minsize = 5, mindev = 5e-04))
## Number of terminal nodes: 114
## Residual mean deviance: 3.176e+09 = 5.156e+13 / 16230
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -311800 -32050 -7455 0 23490 421300
plot(house.big.tree)
text(house.big.tree) # not so readable
#in-sample predictions for larger tree
yhat.house.big.tree.t<-predict(house.big.tree)
#Let's cost complexity prune the much larger tree
house.big.tree.cv<- cv.tree(house.big.tree)
par(mfrow = c(1, 2))
plot(house.big.tree.cv$size, house.big.tree.cv$dev, type = "b")
plot(house.big.tree.cv$k, house.big.tree.cv$dev, type = "b")
house.big.tree.cv
## $size
## [1] 114 113 112 111 110 109 108 107 106 105 104 100 99 98 97 94 93 92 91
## [20] 89 87 86 84 83 82 81 80 77 76 75 74 73 72 70 69 68 67 66
## [39] 63 62 61 60 59 58 56 54 53 52 51 49 46 45 42 41 40 39 38
## [58] 37 36 33 32 31 29 28 26 23 22 21 20 18 17 16 15 14 13 11
## [77] 10 9 8 7 6 5 4 3 2 1
##
## $dev
## [1] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [6] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [11] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [16] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [21] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [26] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [31] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [36] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [41] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [46] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [51] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [56] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [61] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [66] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [71] 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13 9.184166e+13
## [76] 9.184166e+13 9.184166e+13 9.184166e+13 9.331721e+13 9.925797e+13
## [81] 1.003572e+14 1.014041e+14 1.096327e+14 1.221233e+14 1.472852e+14
## [86] 2.157197e+14
##
## $k
## [1] -Inf 1.079980e+11 1.110458e+11 1.116917e+11 1.127354e+11
## [6] 1.150011e+11 1.205640e+11 1.208964e+11 1.218157e+11 1.219461e+11
## [11] 1.238059e+11 1.247160e+11 1.249612e+11 1.268959e+11 1.280067e+11
## [16] 1.288015e+11 1.308031e+11 1.317577e+11 1.340931e+11 1.342989e+11
## [21] 1.359147e+11 1.378865e+11 1.400819e+11 1.410063e+11 1.421355e+11
## [26] 1.423857e+11 1.426139e+11 1.482063e+11 1.535622e+11 1.596828e+11
## [31] 1.604444e+11 1.723854e+11 1.757871e+11 1.857710e+11 1.869755e+11
## [36] 1.889183e+11 1.939476e+11 1.956700e+11 1.983004e+11 1.985889e+11
## [41] 1.990274e+11 1.993447e+11 2.067590e+11 2.280778e+11 2.432833e+11
## [46] 2.579686e+11 2.725111e+11 2.735504e+11 2.894363e+11 2.969075e+11
## [51] 3.046837e+11 3.154270e+11 3.170222e+11 3.399202e+11 3.689565e+11
## [56] 3.796910e+11 3.851877e+11 3.926888e+11 4.022388e+11 4.155287e+11
## [61] 4.233193e+11 4.243845e+11 4.270863e+11 5.257213e+11 6.505707e+11
## [66] 6.656490e+11 7.064992e+11 7.103840e+11 8.289254e+11 8.876273e+11
## [71] 8.981074e+11 9.473146e+11 1.091823e+12 1.123989e+12 1.373380e+12
## [76] 1.510702e+12 1.705109e+12 1.854143e+12 2.187748e+12 2.807915e+12
## [81] 2.887209e+12 3.211974e+12 8.133910e+12 1.273596e+13 2.752198e+13
## [86] 6.657973e+13
##
## $method
## [1] "deviance"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
min(house.big.tree.cv$dev)
## [1] 9.184166e+13
which(house.big.tree.cv$dev==min(house.big.tree.cv$dev))
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
## [76] 76 77 78
house.big.tree.cv$size[which(house.big.tree.cv$dev==min(house.big.tree.cv$dev))]
## [1] 114 113 112 111 110 109 108 107 106 105 104 100 99 98 97 94 93 92 91
## [20] 89 87 86 84 83 82 81 80 77 76 75 74 73 72 70 69 68 67 66
## [39] 63 62 61 60 59 58 56 54 53 52 51 49 46 45 42 41 40 39 38
## [58] 37 36 33 32 31 29 28 26 23 22 21 20 18 17 16 15 14 13 11
## [77] 10 9
#prune down and see how the size 16 tree works
house.prune <- prune.tree(house.big.tree, best = 16)
par(mfrow=c(1,2))
plot(house.tree) #original tree
plot(house.prune) #pruned tree
#check in-sample error rate of pruned tree
summary(house.prune)
##
## Regression tree:
## snip.tree(tree = house.big.tree, nodes = c(183L, 80L, 41L, 21L,
## 8L, 15L, 27L, 182L, 12L, 9L, 23L, 14L, 81L, 26L, 90L, 44L))
## Variables actually used in tree construction:
## [1] "median_income" "ocean_proximity" "longitude"
## [4] "latitude" "housing_median_age"
## Number of terminal nodes: 16
## Residual mean deviance: 4.866e+09 = 7.947e+13 / 16330
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -343600 -41970 -9986 0 32250 408800
#in-sample predictions for pruned tree
yhat.house.prune.t<-predict(house.prune)
#obtain the rest of the in-sample MSEs
MSE.latlon.tree.t<-mean((y.t-yhat.latlon.tree.t)^2)
MSE.house.big.tree.t<-mean((y.t-yhat.house.big.tree.t)^2)
MSE.house.prune.t<-mean((y.t-yhat.house.prune.t)^2)
par(mfrow=c(2,2))
plot(yhat.house.tree.t,y.t);abline(0,1)
plot(yhat.latlon.tree.t,y.t);abline(0,1)
plot(yhat.house.prune.t,y.t);abline(0,1)
plot(yhat.house.big.tree.t,y.t);abline(0,1)
########################################################################################
#############################Holdout analysis###########################################
########################################################################################
#obtain holdout predictions and MSEs
yhat.house.tree.h<-predict(house.tree,newdata=house.h)
yhat.latlon.tree.h<-predict(latlon.tree,newdata=house.h)
yhat.house.big.tree.h<-predict(house.big.tree,newdata=house.h)
yhat.house.prune.h<-predict(house.prune,newdata=house.h)
par(mfrow=c(2,2))
#holdout response
y.h<-house.h$median_house_value
#holdout MSE
MSE.house.tree.h<-mean((y.h-yhat.house.tree.h)^2)
MSE.latlon.tree.h<-mean((y.h-yhat.latlon.tree.h)^2)
MSE.house.big.tree.h<-mean((y.h-yhat.house.big.tree.h)^2)
MSE.house.prune.h<-mean((y.h-yhat.house.prune.h)^2)
#plot some MSE performance metrics as a function of modeling strategy
lab<-c('Default','Lat/Lon','Large','Pruned')
lab.h.vec<-c('MSE.house.tree.h','MSE.latlon.tree.h','MSE.house.big.tree.h','MSE.house.prune.h')
MSE.h.vec<-c(MSE.house.tree.h,MSE.latlon.tree.h,MSE.house.big.tree.h,MSE.house.prune.h)
MSE.t.vec<-c(MSE.house.tree.t,MSE.latlon.tree.t,MSE.house.big.tree.t,MSE.house.prune.t)
par(las=2,mfrow=c(1,1))
plot(MSE.h.vec,axes=FALSE,xlab="Tree type",ylab="MSE",ylim=c(3e9,8e9))
axis(1,at=1:4,labels=lab)
axis(2,at=c(3e9,4e9,5e9,6e9,7e9,8e9))
points(MSE.t.vec,pch=16)
legend(3,7e9,legend=c("Out-of-sample",'In-sample'),pch=c(1,16) )
#This is the second time we see a tree much larger than the default beat both the default and
#pruned trees
#Let's keep these comparisons in mind as we add bagging, random forests, and boosting into the mix
par(mfrow=c(2,2),las=1)
plot(yhat.house.tree.h,y.h);abline(0,1)
plot(yhat.latlon.tree.h,y.h);abline(0,1)
plot(yhat.house.prune.h,y.h);abline(0,1)
plot(yhat.house.big.tree.h,y.h);abline(0,1)
###################################################################################################
#################################Bagging, Random forests ##########################################
###################################################################################################
#Note the randomForest package can be used for bagging by setting number of predictors used in tree
#to the total number of candidate predictors that are being entertained
library(randomForest)
## Warning: package 'randomForest' was built under R version 4.5.3
## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
#Note ntree can be increased - kept it at 100 for LIVE demo
house.bag <- randomForest(median_house_value~.,data=house.t, mtry = 9, importance = TRUE,ntree=100)
class(house.bag)
## [1] "randomForest.formula" "randomForest"
#bagging training predictions
yhat.bag.t <- predict(house.bag)
MSE.bag.t<-mean((y.t-yhat.bag.t)^2)
#bagging holdout predictions
yhat.bag.h <- predict(house.bag,newdata=house.h)
MSE.bag.h<-mean((y.h-yhat.bag.h)^2)
#Fit a random forest - achieved by reducing the number of candidate variables allowed in each tree
house.rf <- randomForest(median_house_value~.,data=house.t, mtry = 3, importance = TRUE,ntree=100)
#random forest training predictions
yhat.rf.t<-predict(house.rf)
MSE.rf.t<-mean((yhat.rf.t - y.t)^2)
#random forest holdout predictions
yhat.rf.h<-predict(house.rf,newdata=house.h)
MSE.rf.h<-mean((yhat.rf.h - y.h)^2)
#Update the performance plot
lab<-c('Default','Lat/Lon','Large','Pruned','Bagging','Random \n Forest')
MSE.h.vec<-c(MSE.house.tree.h,MSE.latlon.tree.h,MSE.house.big.tree.h,MSE.house.prune.h,
MSE.bag.h,MSE.rf.h)
MSE.t.vec<-c(MSE.house.tree.t,MSE.latlon.tree.t,MSE.house.big.tree.t,MSE.house.prune.t,
MSE.bag.t,MSE.rf.t)
par(mfrow=c(1,1),las=2)
plot(MSE.h.vec,axes=FALSE,xlab="Tree type",ylab="MSE",ylim=c(2e9,8e9))
axis(1,at=1:6,labels=lab)
axis(2,at=c(2e9,3e9,4e9,5e9,6e9,7e9,8e9))
points(MSE.t.vec,pch=16)
legend(3,7e9,legend=c('In-sample',"Out-of-sample"),pch=c(1,16) )
text(5.5,3e9,"Why do in- and out-of-sample results \n seem so similar for bagging, RF?",
col='red')
text(5.5,2e9,"Hint: The *in-sample* rates are OOB.",
col='red')
#optional lines
lines(1:6,MSE.h.vec)
lines(1:6,MSE.t.vec)
#assess the performance of the predictors
#This is a nice video about this: https://www.youtube.com/watch?v=8h3H0j2f24I&t=614s
importance(house.bag)
## %IncMSE IncNodePurity
## longitude 44.38520 2.299859e+13
## latitude 45.00571 2.216335e+13
## housing_median_age 68.22429 1.110990e+13
## total_rooms 29.72186 4.958221e+12
## total_bedrooms 23.34638 4.679110e+12
## population 34.83529 6.802496e+12
## households 20.03476 3.794495e+12
## median_income 164.57931 1.036603e+14
## ocean_proximity 58.51955 3.316998e+13
varImpPlot(house.bag)
#for regression tree, node impurity is residual sum of squares
###################################################################################################
#################################Boosting############### ##########################################
###################################################################################################
library(gbm)
## Warning: package 'gbm' was built under R version 4.5.3
## Loaded gbm 2.3.1
## This version of gbm is no longer under development. Consider transitioning to gbm3, https://github.com/gbm-developers/gbm3
set.seed(23784)
##make ocean proximity a factor for compatibility
house.t$ocean_proximity<-factor(house.t$ocean_proximity)
house.h$ocean_proximity<-factor(house.h$ocean_proximity)
house.boost <- gbm(median_house_value~.,data=house.t,
distribution = "gaussian", n.trees = 5000,
interaction.depth = 4)
par(mar=(c(5, 10, 4, 2) + 0.1),las=1) #make left margin large to read variables in next plot
summary(house.boost)
## var rel.inf
## median_income median_income 42.929846
## ocean_proximity ocean_proximity 12.538580
## longitude longitude 11.001765
## latitude latitude 8.504775
## population population 6.478452
## housing_median_age housing_median_age 5.274108
## total_rooms total_rooms 4.800963
## total_bedrooms total_bedrooms 4.753598
## households households 3.717913
plot(house.boost, i = "median_income")
plot(house.boost, i = "ocean_proximity")
plot(house.boost, i = "population")
plot(house.boost, i = "housing_median_age")
yhat.boost.t <- predict(house.boost, n.trees = 5000)
yhat.boost.h <- predict(house.boost,newdata = house.h, n.trees = 5000)
MSE.boost.t<-mean((yhat.boost.t - y.t)^2)
MSE.boost.h<-mean((yhat.boost.h - y.h)^2)
#Update the performance plot
lab<-c('Default','Lat/Lon','Large','Pruned','Bagging','Random \n Forest','Boosting')
MSE.h.vec<-c(MSE.house.tree.h,MSE.latlon.tree.h,MSE.house.big.tree.h,MSE.house.prune.h,
MSE.bag.h,MSE.rf.h,MSE.boost.h)
MSE.t.vec<-c(MSE.house.tree.t,MSE.latlon.tree.t,MSE.house.big.tree.t,MSE.house.prune.t,
MSE.bag.t,MSE.rf.t,MSE.boost.t)
par(mfrow=c(1,1),las=2)
plot(MSE.h.vec,axes=FALSE,xlab="Tree type",ylab="MSE",ylim=c(1e9,8e9))
axis(1,at=1:7,labels=lab)
axis(2,at=c(1e9,2e9,3e9,4e9,5e9,6e9,7e9,8e9))
points(MSE.t.vec,pch=16)
legend(5,7e9,legend=c('In-sample',"Out-of-sample"),pch=c(16,1) )
#Problem 3
rm(list=ls())
##This is analysis of California housing data
#install.packages("tree")
library(tree)
#install.packages("randomForest")
#install.packages("gbm")
#There is a full couse on spatial statistics - anyone taking it? This might be a great data set to
#compare methods for course project
#California housing example from 1990
#each row corresponds to a neighborhood
dat.miss<-read.csv("C:/Users/Sabuj Ganguly/OneDrive/Documents/Ph.D 3rd sem/Regression/Housing data.csv")
#Study missing values
#next line sweeps data frame for NA values
sum(is.na(dat.miss)) #207 missing values in set, a tiny proportion of records
## [1] 207
missingcount<-apply(is.na(dat.miss),1, sum)
#A sophisticated strategy to assess these missing values and account for them in the analysis would
#be worthwhile, but I'll just use casewise deletion (i.e., remove rows w missing values) to keep
#the exposition on the tree-based modeling.
dat<-dat.miss[missingcount==0,]
dat$ocean_proximity<-as.factor(dat$ocean_proximity)
#LOG TRANSFORM THE OUTCOME
dat$median_house_value<-log(dat$median_house_value)
#The outcome variable is median_house_value
#Get a good look at the data dictionary to understand scale these variables are measured on
#https://www.kaggle.com/datasets/camnugent/california-housing-prices?resource=download
#Create a training and holdout set - let's use 80% training 20% holdout
N<-dim(dat)[1]
n<-round(.8*N)
set.seed(73841)
train.ind<-sort(sample(1:N, n, replace=FALSE))
test.ind<- seq(1,N)[-train.ind]
house.t<-dat[train.ind,] #training data
house.h<-dat[test.ind,] #holdout data
#Articulate the modeling plan
#I plan to fit a regression tree and use cost-complexity pruning, bagging, random forests, and
#boosting to train models on the training set, then assess performance on the holdout set.
#To understand how spatial information (latitude and longitude) perform, we conduct the analysis
#both including and excluding these terms and using only those terms
#Let's also consider the implications of e.g., log transforming outcomes
#Once I start analyzing holdout performance, I won't develop new training strategies
#otherwise I am implicitly training on holdout data
#A neat idea might be to make a spatial map of the data then show where lat and long splits
#are (using a tree that only includes latitude and longitude)
#Fit a regression tree to the training data
#Fit a classification tree to the training data
house.tree<-tree(median_house_value~.,data=house.t)
plot(house.tree)
text(house.tree, pretty = 0)
summary(house.tree)
##
## Regression tree:
## tree(formula = median_house_value ~ ., data = house.t)
## Variables actually used in tree construction:
## [1] "ocean_proximity" "median_income"
## Number of terminal nodes: 8
## Residual mean deviance: 0.1246 = 2036 / 16340
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.16000 -0.24030 -0.02189 0.00000 0.22510 1.91300
#Demonstrate that deviance, residual mean deviance are SSE, MSE respectively
ls(summary(house.tree))
## [1] "call" "dev" "df" "residuals" "size" "type"
## [7] "used"
#compute predicted values on the training set
yhat.house.tree.t<-predict(house.tree)
y.t<-house.t$median_house_value
SSE.house.tree.t<-sum((y.t-yhat.house.tree.t)^2)
MSE.house.tree.t<-mean((y.t-yhat.house.tree.t)^2)
SSE.house.tree.t
## [1] 2035.796
MSE.house.tree.t
## [1] 0.124544
#SSE and MSE are the deviance and residual mean deviance in tree output
###############################QUIZ#################################
#SSE and MSE seem "intuitively gigantic" - why do you think that is?
####################################################################
#What if I fit a tree using only spatial information from latitude, longitude
latlon.tree<-tree(median_house_value~latitude + longitude,data=house.t)
plot(latlon.tree)
text(latlon.tree)
summary(latlon.tree)
##
## Regression tree:
## tree(formula = median_house_value ~ latitude + longitude, data = house.t)
## Number of terminal nodes: 13
## Residual mean deviance: 0.1622 = 2649 / 16330
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.537000 -0.247100 -0.003685 0.000000 0.258600 1.839000
#Much larger training MSE
#in-sample predictions for latlong tree
yhat.latlon.tree.t<-predict(latlon.tree)
#Let's set up a larger tree - do cost complexity pruning
house.big.tree<-tree(median_house_value~.,data=house.t,
control=tree.control(nobs = 16346, minsize = 5, mindev=0.0005))
summary(house.big.tree)
##
## Regression tree:
## tree(formula = median_house_value ~ ., data = house.t, control = tree.control(nobs = 16346,
## minsize = 5, mindev = 5e-04))
## Variables actually used in tree construction:
## [1] "ocean_proximity" "median_income" "latitude"
## [4] "longitude" "housing_median_age" "total_rooms"
## [7] "population" "households"
## Number of terminal nodes: 93
## Residual mean deviance: 0.07497 = 1219 / 16250
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.32900 -0.16230 -0.01025 0.00000 0.15000 1.72300
plot(house.big.tree)
text(house.big.tree) # not so readable
#in-sample predictions for larger tree
yhat.house.big.tree.t<-predict(house.big.tree)
#Let's cost complexity prune the much larger tree
house.big.tree.cv<- cv.tree(house.big.tree)
par(mar=c(4,4,2,1))
par(mfrow = c( 2,1))
plot(house.big.tree.cv$size, house.big.tree.cv$dev, type = "b")
plot(house.big.tree.cv$k, house.big.tree.cv$dev, type = "b")
house.big.tree.cv
## $size
## [1] 93 92 91 90 89 88 85 84 83 82 81 80 79 78 77 76 74 73 72 71 70 69 68 64 62
## [26] 61 60 58 57 56 54 53 52 51 50 46 45 44 43 42 41 40 39 36 35 34 33 32 31 30
## [51] 28 27 26 24 22 20 19 18 17 16 15 14 13 12 10 9 8 7 6 5 4 3 2 1
##
## $dev
## [1] 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276
## [9] 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276
## [17] 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276
## [25] 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276
## [33] 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276
## [41] 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276
## [49] 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276
## [57] 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276 2077.276
## [65] 2077.276 2077.276 2077.276 2122.080 2215.129 2316.793 2539.187 2925.677
## [73] 3572.156 5281.063
##
## $k
## [1] -Inf 2.829206 2.888077 2.964482 2.990269 3.014699
## [7] 3.183621 3.204789 3.256758 3.296962 3.418494 3.456331
## [13] 3.476925 3.749916 3.811836 3.849820 3.925207 3.959622
## [19] 3.977917 4.153583 4.427128 4.503697 4.641103 4.742687
## [25] 4.792508 4.972116 5.001589 5.118530 5.248363 5.483189
## [31] 5.590796 5.749868 5.829775 6.314001 6.455388 6.605174
## [37] 6.638390 6.890477 7.017237 7.136849 7.445855 7.661306
## [43] 8.240106 8.384776 8.681937 8.718536 9.182878 9.702217
## [49] 10.054692 10.604465 10.839345 12.239572 13.186718 15.362169
## [55] 15.399442 17.247955 17.500553 17.911542 18.034388 18.789539
## [61] 19.595042 27.990134 28.089134 29.762243 31.179959 43.845505
## [67] 46.426504 61.417668 97.026267 130.895731 187.179778 396.766272
## [73] 662.460304 1708.993470
##
## $method
## [1] "deviance"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
min(house.big.tree.cv$dev)
## [1] 2077.276
which(house.big.tree.cv$dev==min(house.big.tree.cv$dev))
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
## [51] 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
house.big.tree.cv$size[which(house.big.tree.cv$dev==min(house.big.tree.cv$dev))]
## [1] 93 92 91 90 89 88 85 84 83 82 81 80 79 78 77 76 74 73 72 71 70 69 68 64 62
## [26] 61 60 58 57 56 54 53 52 51 50 46 45 44 43 42 41 40 39 36 35 34 33 32 31 30
## [51] 28 27 26 24 22 20 19 18 17 16 15 14 13 12 10 9 8
#prune down and see how the size 16 tree works
house.prune <- prune.tree(house.big.tree, best = 16)
par(mfrow=c(1,2))
plot(house.tree) #original tree
plot(house.prune) #pruned tree
#check in-sample error rate of pruned tree
summary(house.prune)
##
## Regression tree:
## snip.tree(tree = house.big.tree, nodes = c(23L, 211L, 96L, 104L,
## 22L, 10L, 53L, 210L, 97L, 27L, 8L, 9L, 49L, 15L, 14L))
## Variables actually used in tree construction:
## [1] "ocean_proximity" "median_income" "latitude" "longitude"
## Number of terminal nodes: 16
## Residual mean deviance: 0.1089 = 1778 / 16330
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.25700 -0.21030 -0.01775 0.00000 0.21390 1.91300
#in-sample predictions for pruned tree
yhat.house.prune.t<-predict(house.prune)
#obtain the rest of the in-sample MSEs
MSE.latlon.tree.t<-mean((y.t-yhat.latlon.tree.t)^2)
MSE.house.big.tree.t<-mean((y.t-yhat.house.big.tree.t)^2)
MSE.house.prune.t<-mean((y.t-yhat.house.prune.t)^2)
par(mfrow=c(2,2))
plot(yhat.house.tree.t,y.t);abline(0,1)
plot(yhat.latlon.tree.t,y.t);abline(0,1)
plot(yhat.house.prune.t,y.t);abline(0,1)
plot(yhat.house.big.tree.t,y.t);abline(0,1)
########################################################################################
#############################Holdout analysis###########################################
########################################################################################
#obtain holdout predictions and MSEs
yhat.house.tree.h<-predict(house.tree,newdata=house.h)
yhat.latlon.tree.h<-predict(latlon.tree,newdata=house.h)
yhat.house.big.tree.h<-predict(house.big.tree,newdata=house.h)
yhat.house.prune.h<-predict(house.prune,newdata=house.h)
par(mfrow=c(2,2))
#holdout response
y.h<-house.h$median_house_value
#holdout MSE
MSE.house.tree.h<-mean((y.h-yhat.house.tree.h)^2)
MSE.latlon.tree.h<-mean((y.h-yhat.latlon.tree.h)^2)
MSE.house.big.tree.h<-mean((y.h-yhat.house.big.tree.h)^2)
MSE.house.prune.h<-mean((y.h-yhat.house.prune.h)^2)
#plot some MSE performance metrics as a function of modeling strategy
lab<-c('Default','Lat/Lon','Large','Pruned')
lab.h.vec<-c('MSE.house.tree.h','MSE.latlon.tree.h','MSE.house.big.tree.h','MSE.house.prune.h')
MSE.h.vec<-c(MSE.house.tree.h,MSE.latlon.tree.h,MSE.house.big.tree.h,MSE.house.prune.h)
MSE.t.vec<-c(MSE.house.tree.t,MSE.latlon.tree.t,MSE.house.big.tree.t,MSE.house.prune.t)
par(las=2,mfrow=c(1,1))
#changed only because MSE is now on log scale
plot(MSE.h.vec,axes=FALSE,xlab="Tree type",ylab="MSE")
axis(1,at=1:4,labels=lab)
axis(2)
points(MSE.t.vec,pch=16)
legend("topright",legend=c("Out-of-sample",'In-sample'),pch=c(1,16))
#This is the second time we see a tree much larger than the default beat both the default and
#pruned trees
#Let's keep these comparisons in mind as we add bagging, random forests, and boosting into the mix
par(mfrow=c(2,2),las=1)
plot(yhat.house.tree.h,y.h);abline(0,1)
plot(yhat.latlon.tree.h,y.h);abline(0,1)
plot(yhat.house.prune.h,y.h);abline(0,1)
plot(yhat.house.big.tree.h,y.h);abline(0,1)
###################################################################################################
#################################Bagging, Random forests ##########################################
###################################################################################################
#Note the randomForest package can be used for bagging by setting number of predictors used in tree
#to the total number of candidate predictors that are being entertained
library(randomForest)
#Note ntree can be increased - kept it at 100 for LIVE demo
house.bag <- randomForest(median_house_value~.,data=house.t, mtry = 9, importance = TRUE,ntree=100)
class(house.bag)
## [1] "randomForest.formula" "randomForest"
#bagging training predictions
yhat.bag.t <- predict(house.bag)
MSE.bag.t<-mean((y.t-yhat.bag.t)^2)
#bagging holdout predictions
yhat.bag.h <- predict(house.bag,newdata=house.h)
MSE.bag.h<-mean((y.h-yhat.bag.h)^2)
#Fit a random forest - achieved by reducing the number of candidate variables allowed in each tree
house.rf <- randomForest(median_house_value~.,data=house.t, mtry = 3, importance = TRUE,ntree=100)
#random forest training predictions
yhat.rf.t<-predict(house.rf)
MSE.rf.t<-mean((yhat.rf.t - y.t)^2)
#random forest holdout predictions
yhat.rf.h<-predict(house.rf,newdata=house.h)
MSE.rf.h<-mean((yhat.rf.h - y.h)^2)
#Update the performance plot
lab<-c('Default','Lat/Lon','Large','Pruned','Bagging','Random \n Forest')
MSE.h.vec<-c(MSE.house.tree.h,MSE.latlon.tree.h,MSE.house.big.tree.h,MSE.house.prune.h,
MSE.bag.h,MSE.rf.h)
MSE.t.vec<-c(MSE.house.tree.t,MSE.latlon.tree.t,MSE.house.big.tree.t,MSE.house.prune.t,
MSE.bag.t,MSE.rf.t)
par(mfrow=c(1,1),las=2)
#changed only because MSE is now on log scale
plot(MSE.h.vec,axes=FALSE,xlab="Tree type",ylab="MSE")
axis(1,at=1:6,labels=lab)
axis(2)
points(MSE.t.vec,pch=16)
legend("topright",legend=c('In-sample',"Out-of-sample"),pch=c(1,16))
#optional lines
lines(1:6,MSE.h.vec)
lines(1:6,MSE.t.vec)
#assess the performance of the predictors
importance(house.bag)
## %IncMSE IncNodePurity
## longitude 58.89653 483.97190
## latitude 68.30844 559.19126
## housing_median_age 49.54565 183.01287
## total_rooms 29.74484 111.13411
## total_bedrooms 23.30352 97.51177
## population 36.04199 151.89020
## households 20.95493 92.43736
## median_income 185.21959 1806.66099
## ocean_proximity 73.54483 1737.39242
varImpPlot(house.bag)
#for regression tree, node impurity is residual sum of squares
###################################################################################################
#################################Boosting############### ##########################################
###################################################################################################
library(gbm)
set.seed(23784)
##make ocean proximity a factor for compatibility
house.t$ocean_proximity<-factor(house.t$ocean_proximity)
house.h$ocean_proximity<-factor(house.h$ocean_proximity)
house.boost <- gbm(median_house_value~.,data=house.t,
distribution = "gaussian", n.trees = 5000,
interaction.depth = 4)
par(mar=(c(5, 10, 4, 2) + 0.1),las=1) #make left margin large to read variables in next plot
summary(house.boost)
## var rel.inf
## median_income median_income 37.853704
## ocean_proximity ocean_proximity 20.452767
## longitude longitude 9.860508
## latitude latitude 8.340269
## population population 6.388501
## total_bedrooms total_bedrooms 4.854051
## total_rooms total_rooms 4.787926
## households households 4.066814
## housing_median_age housing_median_age 3.395461
plot(house.boost, i = "median_income")
plot(house.boost, i = "ocean_proximity")
plot(house.boost, i = "population")
plot(house.boost, i = "housing_median_age")
yhat.boost.t <- predict(house.boost, n.trees = 5000)
yhat.boost.h <- predict(house.boost,newdata = house.h, n.trees = 5000)
MSE.boost.t<-mean((yhat.boost.t - y.t)^2)
MSE.boost.h<-mean((yhat.boost.h - y.h)^2)
#Update the performance plot
lab<-c('Default','Lat/Lon','Large','Pruned','Bagging','Random \n Forest','Boosting')
MSE.h.vec<-c(MSE.house.tree.h,MSE.latlon.tree.h,MSE.house.big.tree.h,MSE.house.prune.h,
MSE.bag.h,MSE.rf.h,MSE.boost.h)
MSE.t.vec<-c(MSE.house.tree.t,MSE.latlon.tree.t,MSE.house.big.tree.t,MSE.house.prune.t,
MSE.bag.t,MSE.rf.t,MSE.boost.t)
par(mfrow=c(1,1),las=2)
#changed only because MSE is now on log scale
plot(MSE.h.vec,axes=FALSE,xlab="Tree type",ylab="MSE")
axis(1,at=1:7,labels=lab)
axis(2)
points(MSE.t.vec,pch=16)
legend("topright",legend=c('In-sample',"Out-of-sample"),pch=c(16,1))
MSE.house.tree.h
## [1] 0.1235289
MSE.latlon.tree.h
## [1] 0.1662697
MSE.house.big.tree.h
## [1] 0.07936123
MSE.house.prune.h
## [1] 0.1075726
#problem 4
rm(list=ls())
##This is analysis of California housing data using XGBoost
#install.packages("xgboost")
library(xgboost)
## Warning: package 'xgboost' was built under R version 4.5.3
#install.packages("pdp")
library(pdp)
## Warning: package 'pdp' was built under R version 4.5.3
#California housing example from 1990
#each row corresponds to a neighborhood
dat.miss<-read.csv("C:/Users/Sabuj Ganguly/OneDrive/Documents/Ph.D 3rd sem/Regression/Housing data.csv")
#Study missing values
sum(is.na(dat.miss))
## [1] 207
missingcount<-apply(is.na(dat.miss),1,sum)
#same casewise deletion as in class demo
dat<-dat.miss[missingcount==0,]
dat$ocean_proximity<-as.factor(dat$ocean_proximity)
#Create a training and holdout set - same split as original demo
N<-dim(dat)[1]
n<-round(.8*N)
set.seed(73841)
train.ind<-sort(sample(1:N,n,replace=FALSE))
test.ind<-seq(1,N)[-train.ind]
house.t<-dat[train.ind,]
house.h<-dat[test.ind,]
###################################################################################################
#################################XGBoost############################################################
###################################################################################################
set.seed(23784)
##make ocean proximity a factor for compatibility
house.t$ocean_proximity<-factor(house.t$ocean_proximity)
house.h$ocean_proximity<-factor(
house.h$ocean_proximity,
levels=levels(house.t$ocean_proximity)
)
#response variables
y.t<-house.t$median_house_value
y.h<-house.h$median_house_value
#predictor data
x.t<-house.t[,names(house.t)!="median_house_value"]
x.h<-house.h[,names(house.h)!="median_house_value"]
#Fit boosted tree model
house.xgb<-xgboost(
x=x.t,
y=y.t,
objective="reg:squarederror",
nrounds=5000,
max_depth=4,
learning_rate=0.1,
subsample=0.5,
seed=23784,
verbosity=0
)
################################################################################
#Predictions
################################################################################
#training predictions
yhat.xgb.t<-predict(house.xgb,x.t)
#holdout predictions
yhat.xgb.h<-predict(house.xgb,x.h)
################################################################################
#MSE
################################################################################
MSE.xgb.t<-mean((y.t-yhat.xgb.t)^2)
MSE.xgb.h<-mean((y.h-yhat.xgb.h)^2)
MSE.xgb.t
## [1] 204609379
MSE.xgb.h
## [1] 2090843488
################################################################################
#Prediction plots
################################################################################
par(mar=c(4,4,2,1))
par(mfrow=c(1,2))
plot(yhat.xgb.t,y.t,
xlab="Predicted house value",
ylab="Observed house value",
main="Training data")
abline(0,1)
plot(yhat.xgb.h,y.h,
xlab="Predicted house value",
ylab="Observed house value",
main="Holdout data")
abline(0,1)
################################################################################
#Variable importance
################################################################################
importance.xgb<-xgb.importance(model=house.xgb)
importance.xgb
## Feature Gain Cover Frequency
## <char> <num> <num> <num>
## 1: median_income 0.40370146 0.15081632 0.15291341
## 2: ocean_proximity 0.12597806 0.03232186 0.03259277
## 3: longitude 0.11313927 0.12632148 0.14184570
## 4: latitude 0.09589604 0.12606207 0.13317871
## 5: housing_median_age 0.06409742 0.09367756 0.09890408
## 6: population 0.06203043 0.13196627 0.12471517
## 7: total_rooms 0.04804229 0.12345883 0.11840820
## 8: total_bedrooms 0.04770368 0.11022210 0.09996202
## 9: households 0.03941135 0.10515351 0.09747993
par(mfrow=c(1,1))
xgb.plot.importance(
importance.xgb,
measure="Gain"
)
################################################################################
#Partial dependence plots
################################################################################
par(mfrow=c(1,1))
partial(
house.xgb,
pred.var="median_income",
train=x.t,
plot=TRUE
)
partial(
house.xgb,
pred.var="ocean_proximity",
train=x.t,
plot=TRUE
)
partial(
house.xgb,
pred.var="longitude",
train=x.t,
plot=TRUE
)
partial(
house.xgb,
pred.var="latitude",
train=x.t,
plot=TRUE
)
MSE.xgb.t
## [1] 204609379
MSE.xgb.h
## [1] 2090843488
#Problem 5
#Create a training, validation, and holdout set
#60% training, 20% validation, 20% holdout
N<-dim(dat)[1]
n<-round(.8*N)
set.seed(73841)
#First preserve the same 20% holdout as the original demo
train.ind<-sort(sample(1:N, n, replace=FALSE))
test.ind<-seq(1,N)[-train.ind]
dat.trainval<-dat[train.ind,]
house.h<-dat[test.ind,] #FINAL HOLDOUT - do not use for tuning
#Now divide the original 80% into 75% training and 25% validation
#This gives approximately 60% train, 20% validation, 20% holdout overall
N.tv<-dim(dat.trainval)[1]
n.t<-round(.75*N.tv)
set.seed(73841)
train2.ind<-sort(sample(1:N.tv, n.t, replace=FALSE))
valid.ind<-seq(1,N.tv)[-train2.ind]
house.t<-dat.trainval[train2.ind,] #training data
house.v<-dat.trainval[valid.ind,] #validation data
dim(house.t)
## [1] 12260 10
dim(house.v)
## [1] 4086 10
dim(house.h)
## [1] 4087 10
###################################################################################################
#################################Boosting with validation##########################################
###################################################################################################
library(gbm)
set.seed(23784)
##make ocean proximity a factor for compatibility
house.t$ocean_proximity<-factor(house.t$ocean_proximity)
house.v$ocean_proximity<-factor(house.v$ocean_proximity,
levels=levels(house.t$ocean_proximity))
house.h$ocean_proximity<-factor(house.h$ocean_proximity,
levels=levels(house.t$ocean_proximity))
################################################################################
#Values of lambda (shrinkage) and d (interaction depth) to examine
################################################################################
lambda.values<-c(.001,.005,.01,.05,.1)
depth.values<-c(1,2,4,6)
#Create object to store results
boost.results<-expand.grid(lambda=lambda.values,
depth=depth.values)
boost.results$MSE.validation<-NA
################################################################################
#Fit models using TRAINING data and evaluate using VALIDATION data
################################################################################
for(i in 1:nrow(boost.results)){
house.boost.temp<-gbm(median_house_value~.,
data=house.t,
distribution="gaussian",
n.trees=5000,
interaction.depth=boost.results$depth[i],
shrinkage=boost.results$lambda[i],
verbose=FALSE)
yhat.boost.v<-predict(house.boost.temp,
newdata=house.v,
n.trees=5000)
boost.results$MSE.validation[i]<-
mean((house.v$median_house_value-yhat.boost.v)^2)
}
################################################################################
#Look at validation results
################################################################################
boost.results
## lambda depth MSE.validation
## 1 0.001 1 5646950014
## 2 0.005 1 4477610100
## 3 0.010 1 4066813035
## 4 0.050 1 3635148875
## 5 0.100 1 3583255333
## 6 0.001 2 4891999986
## 7 0.005 2 3455059406
## 8 0.010 2 3137608367
## 9 0.050 2 2666642199
## 10 0.100 2 2580214452
## 11 0.001 4 4112133065
## 12 0.005 4 2916478811
## 13 0.010 4 2661045992
## 14 0.050 4 2329019363
## 15 0.100 4 2296502878
## 16 0.001 6 3694825136
## 17 0.005 6 2704765546
## 18 0.010 6 2507253101
## 19 0.050 6 2273969008
## 20 0.100 6 2307942637
#Sort from smallest validation MSE to largest
boost.results<-boost.results[order(boost.results$MSE.validation),]
boost.results
## lambda depth MSE.validation
## 19 0.050 6 2273969008
## 15 0.100 4 2296502878
## 20 0.100 6 2307942637
## 14 0.050 4 2329019363
## 18 0.010 6 2507253101
## 10 0.100 2 2580214452
## 13 0.010 4 2661045992
## 9 0.050 2 2666642199
## 17 0.005 6 2704765546
## 12 0.005 4 2916478811
## 8 0.010 2 3137608367
## 7 0.005 2 3455059406
## 5 0.100 1 3583255333
## 4 0.050 1 3635148875
## 16 0.001 6 3694825136
## 3 0.010 1 4066813035
## 11 0.001 4 4112133065
## 2 0.005 1 4477610100
## 6 0.001 2 4891999986
## 1 0.001 1 5646950014
################################################################################
#Find the best tuning parameters
################################################################################
best.lambda<-boost.results$lambda[1]
best.depth<-boost.results$depth[1]
best.lambda
## [1] 0.05
best.depth
## [1] 6
boost.results$MSE.validation[1]
## [1] 2273969008
################################################################################
#Refit selected model using training + validation data
################################################################################
house.tv<-rbind(house.t,house.v)
house.tv$ocean_proximity<-factor(house.tv$ocean_proximity)
house.h$ocean_proximity<-factor(house.h$ocean_proximity,
levels=levels(house.tv$ocean_proximity))
set.seed(23784)
house.boost.best<-gbm(median_house_value~.,
data=house.tv,
distribution="gaussian",
n.trees=5000,
interaction.depth=best.depth,
shrinkage=best.lambda,
verbose=FALSE)
################################################################################
#Final holdout performance
################################################################################
yhat.boost.h<-predict(house.boost.best,
newdata=house.h,
n.trees=5000)
y.h<-house.h$median_house_value
MSE.boost.h<-mean((y.h-yhat.boost.h)^2)
MSE.boost.h
## [1] 1985024137
sqrt(MSE.boost.h)
## [1] 44553.61
par(mar=(c(5,10,4,2)+0.1),las=1)
summary(house.boost.best)
## var rel.inf
## median_income median_income 44.285063
## ocean_proximity ocean_proximity 13.609492
## longitude longitude 10.623145
## latitude latitude 8.700663
## population population 5.600831
## housing_median_age housing_median_age 5.336120
## total_bedrooms total_bedrooms 4.374389
## total_rooms total_rooms 4.284322
## households households 3.185974
plot(house.boost.best, i="median_income",
n.trees=5000)
plot(house.boost.best, i="ocean_proximity",
n.trees=5000)
plot(house.boost.best, i="population",
n.trees=5000)
plot(house.boost.best, i="housing_median_age",
n.trees=5000)
plot(boost.results$MSE.validation,
pch=16,
xlab="Tuning parameter combination",
ylab="Validation MSE")
axis(1,
at=1:nrow(boost.results),
labels=paste("lambda=",boost.results$lambda,
"\nd=",boost.results$depth),
las=2,
cex.axis=.7)
#problem 6
rm(list=ls())
################################################################################
######################## California Housing Regression ########################
################################################################################
#California housing example from 1990
#each row corresponds to a neighborhood
dat.miss<-read.csv("C:/Users/Sabuj Ganguly/OneDrive/Documents/Ph.D 3rd sem/Regression/Housing data.csv")
#Study missing values
sum(is.na(dat.miss))
## [1] 207
missingcount<-apply(is.na(dat.miss),1,sum)
#same casewise deletion as in original demo
dat<-dat.miss[missingcount==0,]
dat$ocean_proximity<-as.factor(dat$ocean_proximity)
#LOG TRANSFORM THE OUTCOME
dat$median_house_value<-log(dat$median_house_value)
################################################################################
################ Training, validation, and holdout sets ########################
################################################################################
#First preserve 20% as the final holdout set
N<-dim(dat)[1]
n<-round(.8*N)
set.seed(73841)
train.ind<-sort(sample(1:N,n,replace=FALSE))
test.ind<-seq(1,N)[-train.ind]
dat.trainval<-dat[train.ind,]
house.h<-dat[test.ind,] #FINAL HOLDOUT
#Now divide the remaining 80% into training and validation
#75% of 80% = 60% of full data
N.tv<-dim(dat.trainval)[1]
n.t<-round(.75*N.tv)
set.seed(73841)
train2.ind<-sort(sample(1:N.tv,n.t,replace=FALSE))
valid.ind<-seq(1,N.tv)[-train2.ind]
house.t<-dat.trainval[train2.ind,] #about 60%
house.v<-dat.trainval[valid.ind,] #about 20%
dim(house.t)
## [1] 12260 10
dim(house.v)
## [1] 4086 10
dim(house.h)
## [1] 4087 10
################################################################################
######################## Candidate regression variables ########################
################################################################################
candidate.vars<-c("longitude",
"latitude",
"housing_median_age",
"total_rooms",
"total_bedrooms",
"population",
"households",
"median_income",
"ocean_proximity")
################################################################################
######################## Best subset regression ################################
################################################################################
#We will fit every possible nonempty subset of the 9 predictors
#2^9 - 1 = 511 possible models
results<-data.frame(Formula=character(),
Number.variables=numeric(),
Validation.MSE=numeric())
count<-1
for(k in 1:length(candidate.vars)){
combinations<-combn(candidate.vars,k,simplify=FALSE)
for(vars in combinations){
form.text<-paste("median_house_value ~",
paste(vars,collapse=" + "))
fit<-lm(as.formula(form.text),data=house.t)
pred.v<-predict(fit,newdata=house.v)
mse.v<-mean((house.v$median_house_value-pred.v)^2)
results[count,1]<-form.text
results[count,2]<-k
results[count,3]<-mse.v
count<-count+1
}
}
################################################################################
######################## Rank models by validation MSE #########################
################################################################################
results<-results[order(results$Validation.MSE),]
#Top five models
top5<-results[1:5,]
top5
## Formula
## 504 median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + median_income + ocean_proximity
## 478 median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + median_income + ocean_proximity
## 511 median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + households + median_income + ocean_proximity
## 483 median_house_value ~ longitude + latitude + total_rooms + total_bedrooms + population + median_income + ocean_proximity
## 507 median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + households + median_income + ocean_proximity
## Number.variables Validation.MSE
## 504 8 0.1157100
## 478 7 0.1159782
## 511 9 0.1159963
## 483 7 0.1162425
## 507 8 0.1162661
################################################################################
######################## Rashomon effect #######################################
################################################################################
best.validation.mse<-results$Validation.MSE[1]
#Percentage worse than the best model
results$Percent.above.best<-
100*(results$Validation.MSE/best.validation.mse-1)
#Top five including difference from best model
top5<-results[1:5,]
top5
## Formula
## 504 median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + median_income + ocean_proximity
## 478 median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + median_income + ocean_proximity
## 511 median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + households + median_income + ocean_proximity
## 483 median_house_value ~ longitude + latitude + total_rooms + total_bedrooms + population + median_income + ocean_proximity
## 507 median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + households + median_income + ocean_proximity
## Number.variables Validation.MSE Percent.above.best
## 504 8 0.1157100 0.0000000
## 478 7 0.1159782 0.2318017
## 511 9 0.1159963 0.2474092
## 483 7 0.1162425 0.4602232
## 507 8 0.1162661 0.4805898
#How many models are within 1% of the best?
sum(results$Percent.above.best<=1)
## [1] 7
#How many models are within 5% of the best?
sum(results$Percent.above.best<=5)
## [1] 27
#Show models within 1% of best
rashomon.models<-results[results$Percent.above.best<=1,]
rashomon.models
## Formula
## 504 median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + median_income + ocean_proximity
## 478 median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + median_income + ocean_proximity
## 511 median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + households + median_income + ocean_proximity
## 483 median_house_value ~ longitude + latitude + total_rooms + total_bedrooms + population + median_income + ocean_proximity
## 507 median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + households + median_income + ocean_proximity
## 414 median_house_value ~ longitude + latitude + total_bedrooms + population + median_income + ocean_proximity
## 508 median_house_value ~ longitude + latitude + total_rooms + total_bedrooms + population + households + median_income + ocean_proximity
## Number.variables Validation.MSE Percent.above.best
## 504 8 0.1157100 0.0000000
## 478 7 0.1159782 0.2318017
## 511 9 0.1159963 0.2474092
## 483 7 0.1162425 0.4602232
## 507 8 0.1162661 0.4805898
## 414 6 0.1165511 0.7268422
## 508 8 0.1166186 0.7852234
################################################################################
######################## Full regression model #################################
################################################################################
full.formula<-median_house_value ~ longitude + latitude +
housing_median_age + total_rooms + total_bedrooms +
population + households + median_income + ocean_proximity
full.train<-lm(full.formula,data=house.t)
yhat.full.v<-predict(full.train,newdata=house.v)
MSE.full.validation<-
mean((house.v$median_house_value-yhat.full.v)^2)
MSE.full.validation
## [1] 0.1159963
################################################################################
################ Refit top five using training + validation ####################
################################################################################
#Only AFTER model selection do we combine training and validation
house.tv<-rbind(house.t,house.v)
top5$InSample.MSE<-NA
top5$Holdout.MSE<-NA
top5.models<-vector("list",5)
for(i in 1:5){
model.i<-lm(as.formula(top5$Formula[i]),
data=house.tv)
top5.models[[i]]<-model.i
#in-sample predictions
pred.tv<-predict(model.i,newdata=house.tv)
#final holdout predictions
pred.h<-predict(model.i,newdata=house.h)
top5$InSample.MSE[i]<-
mean((house.tv$median_house_value-pred.tv)^2)
top5$Holdout.MSE[i]<-
mean((house.h$median_house_value-pred.h)^2)
}
################################################################################
######################## Refit full regression #################################
################################################################################
full.model<-lm(full.formula,data=house.tv)
yhat.full.t<-predict(full.model,newdata=house.tv)
yhat.full.h<-predict(full.model,newdata=house.h)
MSE.full.t<-
mean((house.tv$median_house_value-yhat.full.t)^2)
MSE.full.h<-
mean((house.h$median_house_value-yhat.full.h)^2)
################################################################################
######################## Regression comparison table ###########################
################################################################################
top5
## Formula
## 504 median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + median_income + ocean_proximity
## 478 median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + median_income + ocean_proximity
## 511 median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + households + median_income + ocean_proximity
## 483 median_house_value ~ longitude + latitude + total_rooms + total_bedrooms + population + median_income + ocean_proximity
## 507 median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + households + median_income + ocean_proximity
## Number.variables Validation.MSE Percent.above.best InSample.MSE Holdout.MSE
## 504 8 0.1157100 0.0000000 0.1091424 0.1068554
## 478 7 0.1159782 0.2318017 0.1092853 0.1066836
## 511 9 0.1159963 0.2474092 0.1088543 0.1066886
## 483 7 0.1162425 0.4602232 0.1099178 0.1076541
## 507 8 0.1162661 0.4805898 0.1089782 0.1065235
regression.results<-data.frame(
Model=c("Regression 1",
"Regression 2",
"Regression 3",
"Regression 4",
"Regression 5",
"Full regression"),
InSample.MSE=c(top5$InSample.MSE,
MSE.full.t),
Holdout.MSE=c(top5$Holdout.MSE,
MSE.full.h)
)
regression.results
## Model InSample.MSE Holdout.MSE
## 1 Regression 1 0.1091424 0.1068554
## 2 Regression 2 0.1092853 0.1066836
## 3 Regression 3 0.1088543 0.1066886
## 4 Regression 4 0.1099178 0.1076541
## 5 Regression 5 0.1089782 0.1065235
## 6 Full regression 0.1088543 0.1066886
################################################################################
######################## Best regression model #################################
################################################################################
#Summary of the model selected using validation MSE
summary(top5.models[[1]])
##
## Call:
## lm(formula = as.formula(top5$Formula[i]), data = house.tv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.41577 -0.19833 -0.00917 0.18997 2.92499
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.521e+00 4.702e-01 -5.363 8.31e-08 ***
## longitude -1.634e-01 5.457e-03 -29.943 < 2e-16 ***
## latitude -1.576e-01 5.403e-03 -29.165 < 2e-16 ***
## housing_median_age 2.546e-03 2.364e-04 10.772 < 2e-16 ***
## total_rooms -1.964e-05 4.247e-06 -4.624 3.79e-06 ***
## total_bedrooms 5.853e-04 2.145e-05 27.285 < 2e-16 ***
## population -1.494e-04 5.033e-06 -29.689 < 2e-16 ***
## median_income 1.717e-01 1.837e-03 93.458 < 2e-16 ***
## ocean_proximityINLAND -3.077e-01 9.408e-03 -32.705 < 2e-16 ***
## ocean_proximityISLAND 4.917e-01 1.909e-01 2.575 0.01002 *
## ocean_proximityNEAR BAY -3.234e-02 1.033e-02 -3.130 0.00175 **
## ocean_proximityNEAR OCEAN -3.370e-02 8.442e-03 -3.992 6.57e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3305 on 16334 degrees of freedom
## Multiple R-squared: 0.6621, Adjusted R-squared: 0.6619
## F-statistic: 2910 on 11 and 16334 DF, p-value: < 2.2e-16
################################################################################
######################## Compare regression models #############################
################################################################################
par(mfrow=c(1,1))
plot(regression.results$Holdout.MSE,
type="b",
xaxt="n",
xlab="Regression model",
ylab="Holdout MSE")
axis(1,
at=1:6,
labels=regression.results$Model,
las=2)
################################################################################
######################## TREE-BASED COMPARISON #################################
################################################################################
library(tree)
library(randomForest)
library(gbm)
################################################################################
#Default regression tree
################################################################################
house.tree<-tree(median_house_value~.,
data=house.tv)
yhat.house.tree.t<-predict(house.tree,newdata=house.tv)
yhat.house.tree.h<-predict(house.tree,newdata=house.h)
MSE.house.tree.t<-
mean((house.tv$median_house_value-yhat.house.tree.t)^2)
MSE.house.tree.h<-
mean((house.h$median_house_value-yhat.house.tree.h)^2)
################################################################################
#Large regression tree
################################################################################
house.big.tree<-tree(median_house_value~.,
data=house.tv,
control=tree.control(
nobs=nrow(house.tv),
minsize=5,
mindev=.0005))
yhat.house.big.tree.t<-predict(house.big.tree,newdata=house.tv)
yhat.house.big.tree.h<-predict(house.big.tree,newdata=house.h)
MSE.house.big.tree.t<-
mean((house.tv$median_house_value-yhat.house.big.tree.t)^2)
MSE.house.big.tree.h<-
mean((house.h$median_house_value-yhat.house.big.tree.h)^2)
################################################################################
#Pruned tree
################################################################################
house.prune<-prune.tree(house.big.tree,best=16)
yhat.house.prune.t<-predict(house.prune,newdata=house.tv)
yhat.house.prune.h<-predict(house.prune,newdata=house.h)
MSE.house.prune.t<-
mean((house.tv$median_house_value-yhat.house.prune.t)^2)
MSE.house.prune.h<-
mean((house.h$median_house_value-yhat.house.prune.h)^2)
################################################################################
#Bagging
################################################################################
set.seed(23784)
house.bag<-randomForest(median_house_value~.,
data=house.tv,
mtry=9,
importance=TRUE,
ntree=100)
#Use newdata here because we want true in-sample fitted performance
yhat.bag.t<-predict(house.bag,newdata=house.tv)
yhat.bag.h<-predict(house.bag,newdata=house.h)
MSE.bag.t<-
mean((house.tv$median_house_value-yhat.bag.t)^2)
MSE.bag.h<-
mean((house.h$median_house_value-yhat.bag.h)^2)
################################################################################
#Random forest
################################################################################
set.seed(23784)
house.rf<-randomForest(median_house_value~.,
data=house.tv,
mtry=3,
importance=TRUE,
ntree=100)
yhat.rf.t<-predict(house.rf,newdata=house.tv)
yhat.rf.h<-predict(house.rf,newdata=house.h)
MSE.rf.t<-
mean((house.tv$median_house_value-yhat.rf.t)^2)
MSE.rf.h<-
mean((house.h$median_house_value-yhat.rf.h)^2)
################################################################################
#Boosting
################################################################################
set.seed(23784)
house.boost<-gbm(median_house_value~.,
data=house.tv,
distribution="gaussian",
n.trees=5000,
interaction.depth=4,
verbose=FALSE)
yhat.boost.t<-predict(house.boost,
newdata=house.tv,
n.trees=5000)
yhat.boost.h<-predict(house.boost,
newdata=house.h,
n.trees=5000)
MSE.boost.t<-
mean((house.tv$median_house_value-yhat.boost.t)^2)
MSE.boost.h<-
mean((house.h$median_house_value-yhat.boost.h)^2)
################################################################################
######################## Final comparison table ################################
################################################################################
tree.results<-data.frame(
Model=c("Default tree",
"Large tree",
"Pruned tree",
"Bagging",
"Random forest",
"Boosting"),
InSample.MSE=c(MSE.house.tree.t,
MSE.house.big.tree.t,
MSE.house.prune.t,
MSE.bag.t,
MSE.rf.t,
MSE.boost.t),
Holdout.MSE=c(MSE.house.tree.h,
MSE.house.big.tree.h,
MSE.house.prune.h,
MSE.bag.h,
MSE.rf.h,
MSE.boost.h)
)
tree.results
## Model InSample.MSE Holdout.MSE
## 1 Default tree 0.124544001 0.12352892
## 2 Large tree 0.074545214 0.07936123
## 3 Pruned tree 0.108756134 0.10757258
## 4 Bagging 0.009743894 0.05030261
## 5 Random forest 0.010975792 0.05107753
## 6 Boosting 0.018905958 0.04736190
#Combine regression and tree-based results
all.results<-rbind(regression.results,
tree.results)
all.results
## Model InSample.MSE Holdout.MSE
## 1 Regression 1 0.109142428 0.10685539
## 2 Regression 2 0.109285290 0.10668357
## 3 Regression 3 0.108854333 0.10668863
## 4 Regression 4 0.109917817 0.10765408
## 5 Regression 5 0.108978180 0.10652349
## 6 Full regression 0.108854333 0.10668863
## 7 Default tree 0.124544001 0.12352892
## 8 Large tree 0.074545214 0.07936123
## 9 Pruned tree 0.108756134 0.10757258
## 10 Bagging 0.009743894 0.05030261
## 11 Random forest 0.010975792 0.05107753
## 12 Boosting 0.018905958 0.04736190
################################################################################
######################## Final performance plot ################################
################################################################################
par(mfrow=c(1,1),las=2)
plot(all.results$Holdout.MSE,
axes=FALSE,
xlab="Model",
ylab="MSE")
axis(1,
at=1:nrow(all.results),
labels=all.results$Model)
axis(2)
points(all.results$InSample.MSE,
pch=16)
legend("topright",
legend=c("Holdout","In-sample"),
pch=c(1,16))
top5[,c("Formula",
"Number.variables",
"Validation.MSE",
"Percent.above.best")]
## Formula
## 504 median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + median_income + ocean_proximity
## 478 median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + median_income + ocean_proximity
## 511 median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + households + median_income + ocean_proximity
## 483 median_house_value ~ longitude + latitude + total_rooms + total_bedrooms + population + median_income + ocean_proximity
## 507 median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + households + median_income + ocean_proximity
## Number.variables Validation.MSE Percent.above.best
## 504 8 0.1157100 0.0000000
## 478 7 0.1159782 0.2318017
## 511 9 0.1159963 0.2474092
## 483 7 0.1162425 0.4602232
## 507 8 0.1162661 0.4805898
top5$Formula[1] # Model 1
## [1] "median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + median_income + ocean_proximity"
top5$Formula[2] # Model 2
## [1] "median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + median_income + ocean_proximity"
top5$Formula[3] # Model 3
## [1] "median_house_value ~ longitude + latitude + housing_median_age + total_rooms + total_bedrooms + population + households + median_income + ocean_proximity"
top5$Formula[4] # Model 4
## [1] "median_house_value ~ longitude + latitude + total_rooms + total_bedrooms + population + median_income + ocean_proximity"
top5$Formula[5] # Model 5
## [1] "median_house_value ~ longitude + latitude + housing_median_age + total_bedrooms + population + households + median_income + ocean_proximity"