Boston 8 - Decision Tree (Regression) with Bagging, Raondom Forest and Boosting with 5-fold Cross Validation
Nao Mimoto - Dept. of Statistics : The University of Akron – Class Web Page
- 1. Load Data and Prepare
- 2. Data Separation
- 3. Decision Tree has high variance (CV.train fit only)
- 4. Decision Tree Growing using rpart(). It has color plotting.
- 5. Decision Tree (Grow and Prune using tree())
- 6. Bagging (Train and Test only)
- 7. Radnom Forest (Train and Test)
- 8. Boosting (Train and Test)
- 9. OLS for comparison
- CODE ONLY
1. Load Data and Prepare
library(MASS) # install.packages("MASS")
library(tidyverse) # install.packages("tidyverse")
Boston <- as_tibble(Boston)
Boston## # A tibble: 506 × 14
## crim zn indus chas nox rm age dis rad tax ptratio black
## <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.00632 18 2.31 0 0.538 6.58 65.2 4.09 1 296 15.3 397.
## 2 0.0273 0 7.07 0 0.469 6.42 78.9 4.97 2 242 17.8 397.
## 3 0.0273 0 7.07 0 0.469 7.18 61.1 4.97 2 242 17.8 393.
## 4 0.0324 0 2.18 0 0.458 7.00 45.8 6.06 3 222 18.7 395.
## 5 0.0690 0 2.18 0 0.458 7.15 54.2 6.06 3 222 18.7 397.
## 6 0.0298 0 2.18 0 0.458 6.43 58.7 6.06 3 222 18.7 394.
## 7 0.0883 12.5 7.87 0 0.524 6.01 66.6 5.56 5 311 15.2 396.
## 8 0.145 12.5 7.87 0 0.524 6.17 96.1 5.95 5 311 15.2 397.
## 9 0.211 12.5 7.87 0 0.524 5.63 100 6.08 5 311 15.2 387.
## 10 0.170 12.5 7.87 0 0.524 6.00 85.9 6.59 5 311 15.2 387.
## # ℹ 496 more rows
## # ℹ 2 more variables: lstat <dbl>, medv <dbl># Rename "medv" column as "resp" to streamline analysis.
Boston2 <- Boston %>% rename(resp=medv) %>% relocate(resp)
Boston2## # A tibble: 506 × 14
## resp crim zn indus chas nox rm age dis rad tax ptratio
## <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 24 0.00632 18 2.31 0 0.538 6.58 65.2 4.09 1 296 15.3
## 2 21.6 0.0273 0 7.07 0 0.469 6.42 78.9 4.97 2 242 17.8
## 3 34.7 0.0273 0 7.07 0 0.469 7.18 61.1 4.97 2 242 17.8
## 4 33.4 0.0324 0 2.18 0 0.458 7.00 45.8 6.06 3 222 18.7
## 5 36.2 0.0690 0 2.18 0 0.458 7.15 54.2 6.06 3 222 18.7
## 6 28.7 0.0298 0 2.18 0 0.458 6.43 58.7 6.06 3 222 18.7
## 7 22.9 0.0883 12.5 7.87 0 0.524 6.01 66.6 5.56 5 311 15.2
## 8 27.1 0.145 12.5 7.87 0 0.524 6.17 96.1 5.95 5 311 15.2
## 9 16.5 0.211 12.5 7.87 0 0.524 5.63 100 6.08 5 311 15.2
## 10 18.9 0.170 12.5 7.87 0 0.524 6.00 85.9 6.59 5 311 15.2
## # ℹ 496 more rows
## # ℹ 2 more variables: black <dbl>, lstat <dbl>## # A tibble: 506 × 14
## resp crim zn indus chas nox rm age dis rad tax ptratio
## <dbl> <dbl> <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 24 0.00632 18 2.31 0 0.538 6.58 65.2 4.09 1 296 15.3
## 2 21.6 0.0273 0 7.07 0 0.469 6.42 78.9 4.97 2 242 17.8
## 3 34.7 0.0273 0 7.07 0 0.469 7.18 61.1 4.97 2 242 17.8
## 4 33.4 0.0324 0 2.18 0 0.458 7.00 45.8 6.06 3 222 18.7
## 5 36.2 0.0690 0 2.18 0 0.458 7.15 54.2 6.06 3 222 18.7
## 6 28.7 0.0298 0 2.18 0 0.458 6.43 58.7 6.06 3 222 18.7
## 7 22.9 0.0883 12.5 7.87 0 0.524 6.01 66.6 5.56 5 311 15.2
## 8 27.1 0.145 12.5 7.87 0 0.524 6.17 96.1 5.95 5 311 15.2
## 9 16.5 0.211 12.5 7.87 0 0.524 5.63 100 6.08 5 311 15.2
## 10 18.9 0.170 12.5 7.87 0 0.524 6.00 85.9 6.59 5 311 15.2
## # ℹ 496 more rows
## # ℹ 2 more variables: black <dbl>, lstat <dbl># Importance List by Chisq test
# top list:
# 3 4 6 10
# "zn" "indus" "nox" "rad"
# bottom list
# c(7, 8, 9, 13)
# "rm" "age" "dis" "black"
#1 resp (dbl) <- used to be "medv". Response variable (Y).
#2 crim (dbl)
#3 zn (dbl)
#4 indus (dbl)
#5 chas (factor)
#6 nox (dbl)
#7 rm (dbl)
#8 age (dbl)
#9 dis (dbl)
#10 rad (dbl)
#11 tax (dbl)
#12 ptratio (dbl)
#13 black (dbl)
#14 lstat (dbl)
# check to see if Boston3 has NA
Orig <- Boston3
#- Check for N/A in data. Remove if there's any.
summary(Orig)## resp crim zn indus chas
## Min. : 5.00 Min. : 0.00632 Min. : 0.00 Min. : 0.46 0:471
## 1st Qu.:17.02 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1: 35
## Median :21.20 Median : 0.25651 Median : 0.00 Median : 9.69
## Mean :22.53 Mean : 3.61352 Mean : 11.36 Mean :11.14
## 3rd Qu.:25.00 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10
## Max. :50.00 Max. :88.97620 Max. :100.00 Max. :27.74
## 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
## Min. : 1.73
## 1st Qu.: 6.95
## Median :11.36
## Mean :12.65
## 3rd Qu.:16.95
## Max. :37.97## [1] 0## [1] 506 14## [1] 506 142. Data Separation
# Load the file from my website
source("https://nmimoto.github.io/R/ML-00.txt") # load CreateCV()
# Use this function to create CV dataset
my.seed=3231
CreateCV(Orig, numFolds=5, seed=my.seed)##
## number of folds: 5
## number of obs: 506
## training data: 84 x 5 = 420
## size of test data: 86
##
## Global variables created:
## Train.set Train.resp
## Test.set Test.resp
## CV.train[[k]] CV.train.resp[[k]]
## CV.valid[[k]] CV.valid.resp[[k]]# Output (all as Data Frame):
# Train.set Train.resp
# Test.set Test.resp
# CV.train[[[[k]]]] CV.train.resp[[[[k]]]]
# CV.valid[[[[k]]]] CV.valid.resp[[[[k]]]]
#--- plot Training and Test for visualization
fold = 2
plot(CV.train[[fold]]$lstat, CV.train[[fold]]$resp, xlab="lstat", ylab="Response (medv)")
lines(CV.valid[[fold]]$lstat, CV.valid[[fold]]$resp, type="p", col="red", pch=19)
3. Decision Tree has high variance (CV.train fit only)
library(tree) # install.packages('tree', repos='https://cran.case.edu/')
tree31 <- tree(resp ~., CV.train[[1]])
summary(tree31)##
## Regression tree:
## tree(formula = resp ~ ., data = CV.train[[1]])
## Variables actually used in tree construction:
## [1] "lstat" "rm" "tax" "ptratio" "nox" "crim"
## Number of terminal nodes: 11
## Residual mean deviance: 13.93 = 4527 / 325
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -14.5800 -2.0430 -0.2444 0.0000 1.9540 29.2500
##
## Regression tree:
## tree(formula = resp ~ ., data = CV.train[[2]])
## Variables actually used in tree construction:
## [1] "lstat" "rm" "rad" "nox"
## Number of terminal nodes: 8
## Residual mean deviance: 16.3 = 5348 / 328
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -17.90000 -2.35400 -0.06714 0.00000 1.87400 18.77000
4. Decision Tree Growing using rpart(). It has color plotting.
### (Train set fitting only)
library(rpart) # install.packages('rpart', repos='https://cran.case.edu/')
library(rpart.plot) # install.packages('rpart.plot', repos='https://cran.case.edu/')
tree00<- rpart(resp~., data=Train.set)
summary(tree00)## Call:
## rpart(formula = resp ~ ., data = Train.set)
## n= 420
##
## CP nsplit rel error xerror xstd
## 1 0.43840439 0 1.0000000 1.0045898 0.09130839
## 2 0.15578179 1 0.5615956 0.6598170 0.06458250
## 3 0.09303331 2 0.4058138 0.4595255 0.05273658
## 4 0.04500192 3 0.3127805 0.4244977 0.05183411
## 5 0.02718020 4 0.2677786 0.3809280 0.04840369
## 6 0.02576447 5 0.2405984 0.3624618 0.04778535
## 7 0.01782710 6 0.2148339 0.3586712 0.05168727
## 8 0.01000000 7 0.1970068 0.3432229 0.05180577
##
## Variable importance
## rm lstat indus dis age ptratio tax nox crim rad
## 31 22 10 7 6 6 5 5 5 2
## zn
## 2
##
## Node number 1: 420 observations, complexity param=0.4384044
## mean=22.55214, MSE=86.47402
## left son=2 (341 obs) right son=3 (79 obs)
## Primary splits:
## rm < 6.797 to the left, improve=0.4384044, (0 missing)
## lstat < 9.725 to the right, improve=0.4362905, (0 missing)
## indus < 6.66 to the right, improve=0.2451984, (0 missing)
## ptratio < 19.65 to the right, improve=0.2285231, (0 missing)
## nox < 0.6695 to the right, improve=0.2176077, (0 missing)
## Surrogate splits:
## lstat < 5.155 to the right, agree=0.883, adj=0.380, (0 split)
## indus < 3.985 to the right, agree=0.850, adj=0.203, (0 split)
## ptratio < 14.55 to the right, agree=0.845, adj=0.177, (0 split)
## zn < 81.25 to the left, agree=0.824, adj=0.063, (0 split)
## crim < 0.023425 to the right, agree=0.821, adj=0.051, (0 split)
##
## Node number 2: 341 observations, complexity param=0.1557818
## mean=19.58856, MSE=39.93022
## left son=4 (140 obs) right son=5 (201 obs)
## Primary splits:
## lstat < 14.4 to the right, improve=0.4155235, (0 missing)
## nox < 0.6695 to the right, improve=0.2669802, (0 missing)
## crim < 9.280695 to the right, improve=0.2375049, (0 missing)
## age < 77.5 to the right, improve=0.1836006, (0 missing)
## ptratio < 19.9 to the right, improve=0.1757154, (0 missing)
## Surrogate splits:
## age < 83.6 to the right, agree=0.812, adj=0.543, (0 split)
## nox < 0.5765 to the right, agree=0.774, adj=0.450, (0 split)
## dis < 2.23995 to the left, agree=0.774, adj=0.450, (0 split)
## indus < 16.57 to the right, agree=0.768, adj=0.436, (0 split)
## tax < 434.5 to the right, agree=0.760, adj=0.414, (0 split)
##
## Node number 3: 79 observations, complexity param=0.09303331
## mean=35.3443, MSE=85.82829
## left son=6 (54 obs) right son=7 (25 obs)
## Primary splits:
## rm < 7.445 to the left, improve=0.4983287, (0 missing)
## lstat < 5.185 to the right, improve=0.3055814, (0 missing)
## ptratio < 19.45 to the right, improve=0.2503916, (0 missing)
## crim < 2.742235 to the right, improve=0.2101518, (0 missing)
## rad < 16 to the right, improve=0.2101518, (0 missing)
## Surrogate splits:
## lstat < 3.99 to the right, agree=0.785, adj=0.32, (0 split)
## ptratio < 14.75 to the right, agree=0.722, adj=0.12, (0 split)
## indus < 1.215 to the right, agree=0.709, adj=0.08, (0 split)
## zn < 81.25 to the left, agree=0.696, adj=0.04, (0 split)
## chas splits as LR, agree=0.696, adj=0.04, (0 split)
##
## Node number 4: 140 observations, complexity param=0.0271802
## mean=14.70786, MSE=17.82972
## left son=8 (60 obs) right son=9 (80 obs)
## Primary splits:
## crim < 6.91188 to the right, improve=0.3954712, (0 missing)
## tax < 567.5 to the right, improve=0.3580768, (0 missing)
## nox < 0.607 to the right, improve=0.3564383, (0 missing)
## dis < 1.94 to the left, improve=0.3396412, (0 missing)
## lstat < 18.905 to the right, improve=0.3196136, (0 missing)
## Surrogate splits:
## rad < 16 to the right, agree=0.900, adj=0.767, (0 split)
## tax < 567.5 to the right, agree=0.871, adj=0.700, (0 split)
## nox < 0.657 to the right, agree=0.779, adj=0.483, (0 split)
## dis < 2.202 to the left, agree=0.729, adj=0.367, (0 split)
## ptratio < 20.15 to the right, agree=0.729, adj=0.367, (0 split)
##
## Node number 5: 201 observations, complexity param=0.04500192
## mean=22.98806, MSE=27.17508
## left son=10 (194 obs) right son=11 (7 obs)
## Primary splits:
## dis < 1.5604 to the right, improve=0.2992259, (0 missing)
## lstat < 4.695 to the right, improve=0.1828713, (0 missing)
## rm < 6.543 to the left, improve=0.1335235, (0 missing)
## crim < 4.866945 to the left, improve=0.1063103, (0 missing)
## nox < 0.6275 to the left, improve=0.1030935, (0 missing)
## Surrogate splits:
## crim < 8.053285 to the left, agree=0.98, adj=0.429, (0 split)
## age < 99.45 to the left, agree=0.98, adj=0.429, (0 split)
## lstat < 3.745 to the right, agree=0.97, adj=0.143, (0 split)
##
## Node number 6: 54 observations, complexity param=0.02576447
## mean=30.89444, MSE=43.9009
## left son=12 (8 obs) right son=13 (46 obs)
## Primary splits:
## lstat < 11.315 to the right, improve=0.3947200, (0 missing)
## ptratio < 19.45 to the right, improve=0.2453645, (0 missing)
## crim < 5.189025 to the right, improve=0.2207138, (0 missing)
## rad < 16 to the right, improve=0.2050469, (0 missing)
## tax < 534.5 to the right, improve=0.2050469, (0 missing)
## Surrogate splits:
## crim < 7.393425 to the right, agree=0.944, adj=0.625, (0 split)
## rad < 16 to the right, agree=0.926, adj=0.500, (0 split)
## tax < 534.5 to the right, agree=0.926, adj=0.500, (0 split)
## indus < 15.015 to the right, agree=0.907, adj=0.375, (0 split)
## nox < 0.659 to the right, agree=0.907, adj=0.375, (0 split)
##
## Node number 7: 25 observations
## mean=44.956, MSE=41.23606
##
## Node number 8: 60 observations
## mean=11.64167, MSE=11.9821
##
## Node number 9: 80 observations
## mean=17.0075, MSE=9.875944
##
## Node number 10: 194 observations, complexity param=0.0178271
## mean=22.44639, MSE=12.36465
## left son=20 (163 obs) right son=21 (31 obs)
## Primary splits:
## rm < 6.543 to the left, improve=0.26991810, (0 missing)
## lstat < 7.685 to the right, improve=0.25300660, (0 missing)
## nox < 0.5125 to the right, improve=0.10822600, (0 missing)
## ptratio < 18.75 to the right, improve=0.08753569, (0 missing)
## tax < 298 to the right, improve=0.06818074, (0 missing)
## Surrogate splits:
## lstat < 4.905 to the right, agree=0.871, adj=0.194, (0 split)
## crim < 0.01733 to the right, agree=0.845, adj=0.032, (0 split)
##
## Node number 11: 7 observations
## mean=38, MSE=204.1457
##
## Node number 12: 8 observations
## mean=20.9125, MSE=52.61609
##
## Node number 13: 46 observations
## mean=32.63043, MSE=22.04299
##
## Node number 20: 163 observations
## mean=21.64969, MSE=8.876978
##
## Node number 21: 31 observations
## mean=26.63548, MSE=9.817128
5. Decision Tree (Grow and Prune using tree())
#--- a. Growing the tree
library(tree) # install.packages('tree', repos='https://cran.case.edu/')
tree1 = tree(resp~., Train.set)
summary(tree1)##
## Regression tree:
## tree(formula = resp ~ ., data = Train.set)
## Variables actually used in tree construction:
## [1] "rm" "lstat" "dis" "crim" "ptratio"
## Number of terminal nodes: 9
## Residual mean deviance: 13.58 = 5580 / 411
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -17.6800 -2.2180 0.1701 0.0000 2.2820 17.6800## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 420 36320.0 22.55
## 2) rm < 6.797 341 13620.0 19.59
## 4) lstat < 14.4 201 5462.0 22.99
## 8) dis < 1.38485 5 390.7 45.58 *
## 9) dis > 1.38485 196 2454.0 22.41
## 18) rm < 6.543 165 1493.0 21.62 *
## 19) rm > 6.543 31 304.3 26.64 *
## 5) lstat > 14.4 140 2496.0 14.71
## 10) crim < 6.91188 80 790.1 17.01 *
## 11) crim > 6.91188 60 718.9 11.64 *
## 3) rm > 6.797 79 6780.0 35.34
## 6) rm < 7.445 54 2371.0 30.89
## 12) crim < 7.39342 49 1134.0 32.32 *
## 13) crim > 7.39342 5 166.9 16.96 *
## 7) rm > 7.445 25 1031.0 44.96
## 14) ptratio < 17.9 20 269.2 47.08 *
## 15) ptratio > 17.9 5 312.7 36.48 *
# Check the training fit and test prediction
Train.fitted = predict(tree1, type="vector")
Test.pred = predict(tree1, Test.set, type="vector")
# Training fit and Test pred plot
plot(Train.fitted, as.matrix(Train.resp), xlab="Fitted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
GrowOnly <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
GrowOnly## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 3.644889 0.8463676 3.285225 0.857898# RMSE Rsquare
#resp 4.605124 0.6852545
# RMSE Rsquare
#1 4.946511 0.6353647
#----------------------
#--- b. Pruning the tree
set.seed(my.seed)
cv.for.pruning = cv.tree(tree1, FUN=prune.tree, K=5) #5-fold CV
#use FUN= prune.misclass if you are doing classification tree
names(cv.for.pruning)## [1] "size" "dev" "k" "method"
# size is the number of terminal nodes
# dev is the Av. CV error rate
# k is the pruning parameter (alpha)
cv.for.pruning## $size
## [1] 9 8 7 6 5 4 3 2 1
##
## $dev
## [1] 9508.553 10049.027 10335.583 11955.665 11955.665 16318.459 16318.459
## [8] 23056.966 36386.824
##
## $k
## [1] -Inf 449.0161 656.9466 987.1600 1069.9094 2617.0803 3378.8850
## [8] 5657.8527 15922.4477
##
## $method
## [1] "deviance"
##
## attr(,"class")
## [1] "prune" "tree.sequence"# In this case, it's probably better not to purne. If you want to prune, use below.
pruned2 = prune.tree(tree1, best=5) # prune to 5-node tree
plot(pruned2)
text(pruned2, pretty=0, cex=1)
# For pruned tree, check the training fit and test prediction
Train.fitted = predict(pruned2, type="vector")
Test.pred = predict(pruned2, Test.set, type="vector")
# Training fit and Test pred plot
plot(Train.fitted, as.matrix(Train.resp), xlab="Fitted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
GrowAndPrune <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
GrowAndPrune## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 4.562482 0.7592775 4.604073 0.7229465## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 3.644889 0.8463676 3.285225 0.8578986. Bagging (Train and Test only)
library(ipred) # install.packages('ipred', repos='https://cran.case.edu/')
set.seed(my.seed)
treeBag01 = bagging(resp~., data=Train.set, nbagg=500,
coob = TRUE,
control = rpart.control(minsplit = 2, cp = 0))
treeBag01##
## Bagging regression trees with 500 bootstrap replications
##
## Call: bagging.data.frame(formula = resp ~ ., data = Train.set, nbagg = 500,
## coob = TRUE, control = rpart.control(minsplit = 2, cp = 0))
##
## Out-of-bag estimate of root mean squared error: 3.3053# Check the training fit and test prediction
Train.fitted = predict(treeBag01, type="vector")
Test.pred = predict(treeBag01, Test.set, type="vector")
# Training fit and Test pred plot
plot(Train.fitted, as.matrix(Train.resp), xlab="Fitted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
Bagging <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
Bagging## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 3.305279 0.8746225 2.443435 0.92143487. Radnom Forest (Train and Test)
###--- change mtry to 6
library(randomForest)
set.seed(my.seed)
treeRF03 = randomForest(resp~., data=Train.set, mtry=6, ntree=25, importance=TRUE)
treeRF03##
## Call:
## randomForest(formula = resp ~ ., data = Train.set, mtry = 6, ntree = 25, importance = TRUE)
## Type of random forest: regression
## Number of trees: 25
## No. of variables tried at each split: 6
##
## Mean of squared residuals: 11.22322
## % Var explained: 87.02## %IncMSE IncNodePurity
## crim 3.0099614 1935.3244
## zn -0.7357398 44.4304
## indus 2.2977240 1119.4098
## chas 1.5305541 158.2590
## nox 4.2976837 3266.9023
## rm 8.6405574 11904.5049
## age 3.1929801 744.2673
## dis 4.4251748 1672.1709
## rad 1.6021103 128.2961
## tax 2.3290421 990.3573
## ptratio 3.3289742 1214.5150
## black 0.5856600 652.0155
## lstat 9.2190670 13354.8619
# Check the training fit and test prediction
Train.fitted = predict(treeRF03, type="response")
Test.pred = predict(treeRF03, Test.set, type="response")
# Training fit and Test pred plot
plot(Train.fitted, as.matrix(Train.resp), xlab="Fitted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
library(caret) # install.packages("caret")
RandomForest <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
RandomForest## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 3.350107 0.8720183 2.55398 0.91327898. Boosting (Train and Test)
set.seed(my.seed)
treeBT01 = gbm::gbm(resp~., data=Train.set, distribution="gaussian", n.trees=5000,
interaction.depth=4)
summary(treeBT01)
## var rel.inf
## lstat lstat 35.3445087
## rm rm 32.1945704
## dis dis 9.9111121
## crim crim 4.4479677
## nox nox 3.9730781
## black black 3.2990865
## age age 2.9813599
## ptratio ptratio 2.5558200
## tax tax 2.4240095
## indus indus 1.2208156
## rad rad 0.7748624
## chas chas 0.6783217
## zn zn 0.1944874
# Check the training fit and test prediction
Train.fitted = predict(treeBT01, type="response")
Test.pred = predict(treeBT01, Test.set, type="response")
# Training fit and Test pred plot
plot(Train.fitted, as.matrix(Train.resp), xlab="Fitted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
library(caret) # install.packages("caret")
Boosting <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
Boosting## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 0.01087816 0.9999986 2.629869 0.90824199. OLS for comparison
##
## Call:
## lm(formula = resp ~ ., data = Train.set)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.8051 -2.7992 -0.6311 1.8135 25.1832
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.001807 5.720082 6.644 9.87e-11 ***
## crim -0.110956 0.035144 -3.157 0.00171 **
## zn 0.044613 0.016287 2.739 0.00643 **
## indus -0.015933 0.070523 -0.226 0.82137
## chas1 2.665636 0.977653 2.727 0.00668 **
## nox -17.394630 4.348874 -4.000 7.53e-05 ***
## rm 3.613494 0.460419 7.848 3.78e-14 ***
## age 0.001465 0.015248 0.096 0.92352
## dis -1.546831 0.230786 -6.702 6.87e-11 ***
## rad 0.298638 0.073530 4.061 5.85e-05 ***
## tax -0.010363 0.004251 -2.438 0.01521 *
## ptratio -0.933205 0.146661 -6.363 5.35e-10 ***
## black 0.007865 0.003030 2.596 0.00979 **
## lstat -0.548968 0.056673 -9.687 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.911 on 406 degrees of freedom
## Multiple R-squared: 0.7304, Adjusted R-squared: 0.7218
## F-statistic: 84.62 on 13 and 406 DF, p-value: < 2.2e-16Train.fitted=predict(Reg01)
Test.pred =predict(Reg01, newdata=Test.set)
plot(Train.fitted, as.matrix(Train.resp), xlab="Predicted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
library(caret) # install.packages("caret")
OLS <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
OLS## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 4.828203 0.7304214 3.957827 0.8075347## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 4.828203 0.7304214 3.957827 0.8075347## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 3.644889 0.8463676 3.285225 0.857898## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 4.562482 0.7592775 4.604073 0.7229465## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 3.305279 0.8746225 2.443435 0.9214348## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 3.350107 0.8720183 2.55398 0.9132789## tr.RMSE tr.Rsquare test.RMSE test.Rsquare
## resp 0.01087816 0.9999986 2.629869 0.9082419CODE ONLY
##--------------------------------------
## 1. Load Data and Prepare
library(MASS) # install.packages("MASS")
library(tidyverse) # install.packages("tidyverse")
Boston <- as_tibble(Boston)
Boston
# Rename "medv" column as "resp" to streamline analysis.
Boston2 <- Boston %>% rename(resp=medv) %>% relocate(resp)
Boston2
# turn "chas" column into 0/1 factor
Boston3 <- Boston2 %>% mutate( chas=as.factor(chas) )
Boston3
# Importance List by Chisq test
# top list:
# 3 4 6 10
# "zn" "indus" "nox" "rad"
# bottom list
# c(7, 8, 9, 13)
# "rm" "age" "dis" "black"
#1 resp (dbl) <- used to be "medv". Response variable (Y).
#2 crim (dbl)
#3 zn (dbl)
#4 indus (dbl)
#5 chas (factor)
#6 nox (dbl)
#7 rm (dbl)
#8 age (dbl)
#9 dis (dbl)
#10 rad (dbl)
#11 tax (dbl)
#12 ptratio (dbl)
#13 black (dbl)
#14 lstat (dbl)
# check to see if Boston3 has NA
Orig <- Boston3
#- Check for N/A in data. Remove if there's any.
summary(Orig)
sum(is.na(Orig))
# If there is na in the data, run below
dim(Orig)
Orig <- Orig %>% na.omit()
dim(Orig)
##----------------------------------------
## 2. Data Separation
# Load the file from my website
source("https://nmimoto.github.io/R/ML-00.txt") # load CreateCV()
# Use this function to create CV dataset
my.seed=3231
CreateCV(Orig, numFolds=5, seed=my.seed)
# Output (all as Data Frame):
# Train.set Train.resp
# Test.set Test.resp
# CV.train[[[[k]]]] CV.train.resp[[[[k]]]]
# CV.valid[[[[k]]]] CV.valid.resp[[[[k]]]]
#--- plot Training and Test for visualization
fold = 2
plot(CV.train[[fold]]$lstat, CV.train[[fold]]$resp, xlab="lstat", ylab="Response (medv)")
lines(CV.valid[[fold]]$lstat, CV.valid[[fold]]$resp, type="p", col="red", pch=19)
##----------------------------------------
## 3. Decision Tree has high variance (CV.train fit only)
library(tree) # install.packages('tree', repos='https://cran.case.edu/')
tree31 <- tree(resp ~., CV.train[[1]])
summary(tree31)
plot(tree31)
text(tree31, pretty=1, cex=.7)
#---------
tree32 <- tree(resp ~., CV.train[[2]])
summary(tree32)
plot(tree32)
text(tree32, pretty=0, cex=.7)
##----------------------------------------
## 4. Decision Tree Growing using rpart(). It has color plotting.
### (Train set fitting only)
library(rpart) # install.packages('rpart', repos='https://cran.case.edu/')
library(rpart.plot) # install.packages('rpart.plot', repos='https://cran.case.edu/')
tree00<- rpart(resp~., data=Train.set)
summary(tree00)
rpart.plot(tree00) # now you can plot in color
plot(tree00) # old way of plotting
text(tree00)
##----------------------------------------
## 5. Decision Tree (Grow and Prune using tree())
#--- a. Growing the tree
library(tree) # install.packages('tree', repos='https://cran.case.edu/')
tree1 = tree(resp~., Train.set)
summary(tree1)
tree1
plot(tree1)
text(tree1, pretty=0, cex=1)
# Check the training fit and test prediction
Train.fitted = predict(tree1, type="vector")
Test.pred = predict(tree1, Test.set, type="vector")
# Training fit and Test pred plot
plot(Train.fitted, as.matrix(Train.resp), xlab="Fitted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
GrowOnly <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
GrowOnly
# RMSE Rsquare
#resp 4.605124 0.6852545
# RMSE Rsquare
#1 4.946511 0.6353647
#----------------------
#--- b. Pruning the tree
set.seed(my.seed)
cv.for.pruning = cv.tree(tree1, FUN=prune.tree, K=5) #5-fold CV
#use FUN= prune.misclass if you are doing classification tree
names(cv.for.pruning)
plot(cv.for.pruning$size, cv.for.pruning$dev, type="b")
plot(cv.for.pruning$k, cv.for.pruning$dev, type="b")
# size is the number of terminal nodes
# dev is the Av. CV error rate
# k is the pruning parameter (alpha)
cv.for.pruning
# In this case, it's probably better not to purne. If you want to prune, use below.
pruned2 = prune.tree(tree1, best=5) # prune to 5-node tree
plot(pruned2)
text(pruned2, pretty=0, cex=1)
# For pruned tree, check the training fit and test prediction
Train.fitted = predict(pruned2, type="vector")
Test.pred = predict(pruned2, Test.set, type="vector")
# Training fit and Test pred plot
plot(Train.fitted, as.matrix(Train.resp), xlab="Fitted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
GrowAndPrune <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
GrowAndPrune
GrowOnly
##----------------------------------------
## 6. Bagging (Train and Test only)
library(ipred) # install.packages('ipred', repos='https://cran.case.edu/')
set.seed(my.seed)
treeBag01 = bagging(resp~., data=Train.set, nbagg=500,
coob = TRUE,
control = rpart.control(minsplit = 2, cp = 0))
treeBag01
# Check the training fit and test prediction
Train.fitted = predict(treeBag01, type="vector")
Test.pred = predict(treeBag01, Test.set, type="vector")
# Training fit and Test pred plot
plot(Train.fitted, as.matrix(Train.resp), xlab="Fitted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
Bagging <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
Bagging
##----------------------------------------
## 7. Radnom Forest (Train and Test)
###--- change mtry to 6
library(randomForest)
set.seed(my.seed)
treeRF03 = randomForest(resp~., data=Train.set, mtry=6, ntree=25, importance=TRUE)
treeRF03
importance (treeRF03)
varImpPlot (treeRF03)
# Check the training fit and test prediction
Train.fitted = predict(treeRF03, type="response")
Test.pred = predict(treeRF03, Test.set, type="response")
# Training fit and Test pred plot
plot(Train.fitted, as.matrix(Train.resp), xlab="Fitted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
library(caret) # install.packages("caret")
RandomForest <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
RandomForest
##----------------------------------------
## 8. Boosting (Train and Test)
set.seed(my.seed)
treeBT01 = gbm::gbm(resp~., data=Train.set, distribution="gaussian", n.trees=5000,
interaction.depth=4)
summary(treeBT01)
plot(treeBT01, i="rm")
plot(treeBT01, i="lstat")
# Check the training fit and test prediction
Train.fitted = predict(treeBT01, type="response")
Test.pred = predict(treeBT01, Test.set, type="response")
# Training fit and Test pred plot
plot(Train.fitted, as.matrix(Train.resp), xlab="Fitted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
library(caret) # install.packages("caret")
Boosting <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
Boosting
##----------------------------------------
## 9. OLS for comparison
Reg01 = lm(resp~., data=Train.set)
summary(Reg01)
Train.fitted=predict(Reg01)
Test.pred =predict(Reg01, newdata=Test.set)
plot(Train.fitted, as.matrix(Train.resp), xlab="Predicted", ylab="Actual")
lines(Test.pred, as.matrix(Test.resp), type="p", col="red", pch=19)
abline(0,1, col="red")
library(caret) # install.packages("caret")
OLS <- data.frame(
tr.RMSE = caret::RMSE(Train.fitted, as.matrix(Train.resp)),
tr.Rsquare = caret::R2(Train.fitted, Train.resp),
test.RMSE = caret::RMSE(Test.pred, as.matrix(Test.resp)),
test.Rsquare = caret::R2(Test.pred, Test.resp)
)
OLS
OLS
GrowOnly
GrowAndPrune
Bagging
RandomForest
Boosting