##This is analysis of email vs spam data set
library(tree)
## Warning: package 'tree' was built under R version 4.5.3
library(kernlab)
data(spam)
class(spam)
## [1] "data.frame"
str(spam)
## 'data.frame': 4601 obs. of 58 variables:
## $ make : num 0 0.21 0.06 0 0 0 0 0 0.15 0.06 ...
## $ address : num 0.64 0.28 0 0 0 0 0 0 0 0.12 ...
## $ all : num 0.64 0.5 0.71 0 0 0 0 0 0.46 0.77 ...
## $ num3d : num 0 0 0 0 0 0 0 0 0 0 ...
## $ our : num 0.32 0.14 1.23 0.63 0.63 1.85 1.92 1.88 0.61 0.19 ...
## $ over : num 0 0.28 0.19 0 0 0 0 0 0 0.32 ...
## $ remove : num 0 0.21 0.19 0.31 0.31 0 0 0 0.3 0.38 ...
## $ internet : num 0 0.07 0.12 0.63 0.63 1.85 0 1.88 0 0 ...
## $ order : num 0 0 0.64 0.31 0.31 0 0 0 0.92 0.06 ...
## $ mail : num 0 0.94 0.25 0.63 0.63 0 0.64 0 0.76 0 ...
## $ receive : num 0 0.21 0.38 0.31 0.31 0 0.96 0 0.76 0 ...
## $ will : num 0.64 0.79 0.45 0.31 0.31 0 1.28 0 0.92 0.64 ...
## $ people : num 0 0.65 0.12 0.31 0.31 0 0 0 0 0.25 ...
## $ report : num 0 0.21 0 0 0 0 0 0 0 0 ...
## $ addresses : num 0 0.14 1.75 0 0 0 0 0 0 0.12 ...
## $ free : num 0.32 0.14 0.06 0.31 0.31 0 0.96 0 0 0 ...
## $ business : num 0 0.07 0.06 0 0 0 0 0 0 0 ...
## $ email : num 1.29 0.28 1.03 0 0 0 0.32 0 0.15 0.12 ...
## $ you : num 1.93 3.47 1.36 3.18 3.18 0 3.85 0 1.23 1.67 ...
## $ credit : num 0 0 0.32 0 0 0 0 0 3.53 0.06 ...
## $ your : num 0.96 1.59 0.51 0.31 0.31 0 0.64 0 2 0.71 ...
## $ font : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num000 : num 0 0.43 1.16 0 0 0 0 0 0 0.19 ...
## $ money : num 0 0.43 0.06 0 0 0 0 0 0.15 0 ...
## $ hp : num 0 0 0 0 0 0 0 0 0 0 ...
## $ hpl : num 0 0 0 0 0 0 0 0 0 0 ...
## $ george : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num650 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ lab : num 0 0 0 0 0 0 0 0 0 0 ...
## $ labs : num 0 0 0 0 0 0 0 0 0 0 ...
## $ telnet : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num857 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ data : num 0 0 0 0 0 0 0 0 0.15 0 ...
## $ num415 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num85 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ technology : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num1999 : num 0 0.07 0 0 0 0 0 0 0 0 ...
## $ parts : num 0 0 0 0 0 0 0 0 0 0 ...
## $ pm : num 0 0 0 0 0 0 0 0 0 0 ...
## $ direct : num 0 0 0.06 0 0 0 0 0 0 0 ...
## $ cs : num 0 0 0 0 0 0 0 0 0 0 ...
## $ meeting : num 0 0 0 0 0 0 0 0 0 0 ...
## $ original : num 0 0 0.12 0 0 0 0 0 0.3 0 ...
## $ project : num 0 0 0 0 0 0 0 0 0 0.06 ...
## $ re : num 0 0 0.06 0 0 0 0 0 0 0 ...
## $ edu : num 0 0 0.06 0 0 0 0 0 0 0 ...
## $ table : num 0 0 0 0 0 0 0 0 0 0 ...
## $ conference : num 0 0 0 0 0 0 0 0 0 0 ...
## $ charSemicolon : num 0 0 0.01 0 0 0 0 0 0 0.04 ...
## $ charRoundbracket : num 0 0.132 0.143 0.137 0.135 0.223 0.054 0.206 0.271 0.03 ...
## $ charSquarebracket: num 0 0 0 0 0 0 0 0 0 0 ...
## $ charExclamation : num 0.778 0.372 0.276 0.137 0.135 0 0.164 0 0.181 0.244 ...
## $ charDollar : num 0 0.18 0.184 0 0 0 0.054 0 0.203 0.081 ...
## $ charHash : num 0 0.048 0.01 0 0 0 0 0 0.022 0 ...
## $ capitalAve : num 3.76 5.11 9.82 3.54 3.54 ...
## $ capitalLong : num 61 101 485 40 40 15 4 11 445 43 ...
## $ capitalTotal : num 278 1028 2259 191 191 ...
## $ type : Factor w/ 2 levels "nonspam","spam": 2 2 2 2 2 2 2 2 2 2 ...
#The variable "type" is the outcome variable
# I'll produce my own training and test sets of the size described in the book
set.seed(3731)
train.ind<-sort(sample(1:4601, 3065, replace=FALSE))
test.ind<- seq(1,4601)[-train.ind]
spam.train<-spam[train.ind,]
spam.test<-spam[test.ind,]
#Fit a classification tree to the training data
spam.tree<-tree(type~.,data=spam.train)
plot(spam.tree)
text(spam.tree)
summary(spam.tree)
##
## Classification tree:
## tree(formula = type ~ ., data = spam.train)
## Variables actually used in tree construction:
## [1] "charDollar" "remove" "charExclamation" "george"
## [5] "hp" "capitalLong" "edu" "our"
## [9] "capitalTotal"
## Number of terminal nodes: 15
## Residual mean deviance: 0.4404 = 1343 / 3050
## Misclassification error rate: 0.07732 = 237 / 3065
conf.mat.train<-table(spam.train$type,predict(spam.tree,type='class'))
conf.mat.train
##
## nonspam spam
## nonspam 1775 78
## spam 159 1053
#How well does this tree perform on the holdout set
pred.tree<-predict(spam.tree, newdata=spam.test,type="class")
conf.mat<-table(spam.test$type,pred.tree)
conf.mat
## pred.tree
## nonspam spam
## nonspam 888 47
## spam 93 508
#misclassification rate
1-sum(diag(conf.mat))/sum(conf.mat)
## [1] 0.09114583
1536-sum(diag(conf.mat)) #Number of incorrect calls
## [1] 140
#Cost complexity pruning the default tree shows misclassification best for full tree
spam.cv <- cv.tree(spam.tree, FUN = prune.misclass)
par(mfrow = c(1, 2))
plot(spam.cv$size, spam.cv$dev, type = "b")
plot(spam.cv$k, spam.cv$dev, type = "b")
spam.cv
## $size
## [1] 15 14 9 8 7 6 5 3 2 1
##
## $dev
## [1] 299 299 309 331 333 337 368 435 644 1212
##
## $k
## [1] -Inf 0.0 7.6 12.0 13.0 14.0 34.0 48.0 177.0 591.0
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
min(spam.cv$dev)
## [1] 299
which(spam.cv$dev==min(spam.cv$dev))
## [1] 1 2
spam.cv$size[2] #size of best tree that is smallest in that class
## [1] 14
#prune down and see how the size 14 tree works
#DOUBLE CHECK best above
spam.prune <- prune.misclass(spam.tree, best = 14)
par(mfrow=c(1,2))
plot(spam.tree) #original tree
plot(spam.prune) #pruned tree
#check in-sample error rate of pruned tree
summary(spam.prune)
##
## Classification tree:
## snip.tree(tree = spam.tree, nodes = 12L)
## Variables actually used in tree construction:
## [1] "charDollar" "remove" "charExclamation" "george"
## [5] "hp" "capitalLong" "edu" "our"
## [9] "capitalTotal"
## Number of terminal nodes: 14
## Residual mean deviance: 0.4563 = 1392 / 3051
## Misclassification error rate: 0.07732 = 237 / 3065
#check out-of-sample error rate of pruned tree
pred.prune<-predict(spam.prune, newdata=spam.test,type="class")
conf.mat.prune<-table(spam.test$type,pred.prune)
conf.mat.prune
## pred.prune
## nonspam spam
## nonspam 888 47
## spam 93 508
#misclassification rate
1-sum(diag(conf.mat.prune))/sum(conf.mat.prune)
## [1] 0.09114583
1536-(sum(diag(conf.mat.prune)))
## [1] 140
#Build a larger tree for cost complexity pruning
spam.tree.ctrl<-tree(type~.,data=spam.train,control=tree.control(nobs = 3065, minsize = 2, mindev=0))
par(mfrow=c(1,1))
plot(spam.tree.ctrl)
text(spam.tree.ctrl)
summary(spam.tree.ctrl)
##
## Classification tree:
## tree(formula = type ~ ., data = spam.train, control = tree.control(nobs = 3065,
## minsize = 2, mindev = 0))
## Variables actually used in tree construction:
## [1] "charDollar" "remove" "charExclamation" "george"
## [5] "hp" "capitalLong" "free" "internet"
## [9] "capitalAve" "you" "edu" "re"
## [13] "capitalTotal" "our" "mail" "all"
## [17] "technology" "will" "your" "charSemicolon"
## [21] "meeting" "address" "num3d" "num650"
## [25] "data" "charRoundbracket" "make" "hpl"
## [29] "over" "email" "charHash" "report"
## [33] "num1999" "business" "receive" "people"
## Number of terminal nodes: 176
## Residual mean deviance: 0.0009597 = 2.773 / 2889
## Misclassification error rate: 0.0003263 = 1 / 3065
#Out of sample misclassification of the large tree
#How well does this tree perform on the holdout set
pred.tree.ctrl<-predict(spam.tree.ctrl, newdata=spam.test,type="class")
conf.mat.ctrl<-table(spam.test$type,pred.tree.ctrl)
conf.mat.ctrl
## pred.tree.ctrl
## nonspam spam
## nonspam 870 65
## spam 71 530
#misclassification rate
1-sum(diag(conf.mat.ctrl))/sum(conf.mat.ctrl)
## [1] 0.08854167
1536-sum(diag(conf.mat.ctrl))
## [1] 136
#Well growing the biggest possible tree overfit relative to the training error rate
#but was superior to the smaller tree in the out-of-sample misclassification rate.
spam.ctrl.cv<- cv.tree(spam.tree.ctrl, FUN = prune.misclass)
par(mfrow = c(1, 2))
plot(spam.ctrl.cv$size, spam.ctrl.cv$dev, type = "b")
plot(spam.ctrl.cv$k, spam.ctrl.cv$dev, type = "b")
spam.ctrl.cv
## $size
## [1] 176 175 170 136 133 125 116 111 104 76 66 64 59 51 39 36 34 30 28
## [20] 24 19 18 16 11 10 9 8 7 6 5 3 2 1
##
## $dev
## [1] 293 293 293 293 293 293 293 293 293 293 293 293 293 293 293
## [16] 293 293 293 293 293 293 296 297 300 296 321 324 328 328 368
## [31] 443 638 1212
##
## $k
## [1] -Inf 0.0000000 0.2000000 0.5000000 0.6666667 0.7500000
## [7] 0.7777778 0.8000000 0.8571429 1.0000000 1.2000000 1.5000000
## [13] 1.6000000 1.7500000 2.0000000 2.3333333 2.5000000 2.7500000
## [19] 3.0000000 3.7500000 4.0000000 5.0000000 6.0000000 7.6000000
## [25] 8.0000000 12.0000000 13.0000000 14.0000000 15.0000000 34.0000000
## [31] 48.0000000 177.0000000 591.0000000
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
min(spam.ctrl.cv$dev)
## [1] 293
which(spam.ctrl.cv$dev==min(spam.ctrl.cv$dev))
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
#be sure to change for current analysis
spam.ctrl.cv$size[1] #Recommends a size __ tree
## [1] 176
#prune down and see how the suize 14 tree works
spam.prune.ctrl <- prune.misclass(spam.tree.ctrl, best = 176)
par(mfrow=c(1,2))
plot(spam.tree) #original tree
plot(spam.prune.ctrl) #pruned tree
#check in-sample error rate of pruned tree
summary(spam.prune.ctrl)
##
## Classification tree:
## tree(formula = type ~ ., data = spam.train, control = tree.control(nobs = 3065,
## minsize = 2, mindev = 0))
## Variables actually used in tree construction:
## [1] "charDollar" "remove" "charExclamation" "george"
## [5] "hp" "capitalLong" "free" "internet"
## [9] "capitalAve" "you" "edu" "re"
## [13] "capitalTotal" "our" "mail" "all"
## [17] "technology" "will" "your" "charSemicolon"
## [21] "meeting" "address" "num3d" "num650"
## [25] "data" "charRoundbracket" "make" "hpl"
## [29] "over" "email" "charHash" "report"
## [33] "num1999" "business" "receive" "people"
## Number of terminal nodes: 176
## Residual mean deviance: 0.0009597 = 2.773 / 2889
## Misclassification error rate: 0.0003263 = 1 / 3065
#check out-of-sample error rate of pruned tree
pred.prune.ctrl<-predict(spam.prune.ctrl, newdata=spam.test,type="class")
conf.mat.prune.ctrl<-table(spam.test$type,pred.prune.ctrl)
conf.mat.prune.ctrl
## pred.prune.ctrl
## nonspam spam
## nonspam 871 64
## spam 71 530
#misclassification rate
1-sum(diag(conf.mat.prune.ctrl))/sum(conf.mat.prune.ctrl)
## [1] 0.08789062
1536-sum(diag(conf.mat.prune.ctrl))
## [1] 135
###Revisit the Gini index version
#spam.tree.Gini<-tree(type~.,data=spam.train,split="gini")
spam.tree.Gini<-tree(type~.,
data=spam.train,
split="gini",
control=tree.control(nobs=3065,
mincut=100,
minsize=200,
mindev=.01))
summary(spam.tree.Gini)
##
## Classification tree:
## tree(formula = type ~ ., data = spam.train, control = tree.control(nobs = 3065,
## mincut = 100, minsize = 200, mindev = 0.01), split = "gini")
## Variables actually used in tree construction:
## [1] "conference" "project" "num415"
## [4] "original" "pm" "technology"
## [7] "num85" "data" "charSquarebracket"
## [10] "report" "charSemicolon" "num1999"
## [13] "hp" "credit" "receive"
## [16] "people" "make" "address"
## [19] "all" "will" "re"
## [22] "charExclamation" "capitalTotal"
## Number of terminal nodes: 24
## Residual mean deviance: 1.051 = 3197 / 3041
## Misclassification error rate: 0.262 = 803 / 3065
summary(spam.tree)
##
## Classification tree:
## tree(formula = type ~ ., data = spam.train)
## Variables actually used in tree construction:
## [1] "charDollar" "remove" "charExclamation" "george"
## [5] "hp" "capitalLong" "edu" "our"
## [9] "capitalTotal"
## Number of terminal nodes: 15
## Residual mean deviance: 0.4404 = 1343 / 3050
## Misclassification error rate: 0.07732 = 237 / 3065
par(mfrow=c(1,1))
plot(spam.tree)
text(spam.tree)
plot(spam.tree.Gini)
text(spam.tree.Gini)
#Compare the various in- and out-of-sample misclassification rates - expectations met?
################################################################################
######################## Bagging ################################################
################################################################################
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.
#There are 57 candidate predictors
#For bagging, allow all predictors to be considered at every split
set.seed(23784)
spam.bag<-randomForest(type~.,
data=spam.train,
mtry=57,
importance=TRUE,
ntree=500)
#training predictions
pred.bag.train<-predict(spam.bag,
newdata=spam.train,
type="class")
conf.mat.bag.train<-table(spam.train$type,
pred.bag.train)
conf.mat.bag.train
## pred.bag.train
## nonspam spam
## nonspam 1852 1
## spam 0 1212
#training misclassification rate
MCR.bag.train<-
1-sum(diag(conf.mat.bag.train))/sum(conf.mat.bag.train)
MCR.bag.train
## [1] 0.0003262643
#holdout predictions
pred.bag<-predict(spam.bag,
newdata=spam.test,
type="class")
conf.mat.bag<-table(spam.test$type,
pred.bag)
conf.mat.bag
## pred.bag
## nonspam spam
## nonspam 899 36
## spam 50 551
#holdout misclassification rate
MCR.bag<-
1-sum(diag(conf.mat.bag))/sum(conf.mat.bag)
MCR.bag
## [1] 0.05598958
#Number incorrectly classified
1536-sum(diag(conf.mat.bag))
## [1] 86
#variable importance
importance(spam.bag)
## nonspam spam MeanDecreaseAccuracy MeanDecreaseGini
## make 3.4375107 7.1340748 7.482192 4.5451238
## address 8.3856032 7.0315143 11.220811 3.4421350
## all 1.7399876 6.5128243 6.859883 5.8280717
## num3d 11.2614372 -0.6113121 9.583824 2.8370080
## our 25.6120408 24.1997680 33.596296 23.5780096
## over 16.1088210 5.2778569 16.357111 4.9756630
## remove 71.8486869 43.8974403 79.650696 180.6049774
## internet 15.1688083 6.2211822 15.190525 8.0687118
## order 11.4536096 1.9803425 12.250713 3.5688217
## mail 5.7479403 5.5055160 8.060951 7.3231617
## receive 17.2489085 2.0130431 17.464684 4.4653853
## will 2.8378200 11.4937870 11.110846 9.2491864
## people -0.9692594 3.4124512 2.245570 3.9128499
## report 10.3119508 8.1193903 12.179344 3.8215417
## addresses 5.3310266 3.0249328 6.265170 1.1292987
## free 36.5417447 30.3690826 44.990199 55.4029369
## business 19.7444486 9.2738823 20.978700 11.3550531
## email 13.0181662 9.0610999 14.245406 9.4070060
## you 14.4237939 15.8286180 22.036126 28.0117604
## credit 12.8634464 -0.4241135 12.880596 2.5884365
## your 15.2606515 15.2803819 20.378356 27.4350838
## font 15.6207333 12.5478591 17.354343 3.4502325
## num000 21.1644025 3.6890746 21.741198 9.3520811
## money 18.3806807 13.6155787 20.588810 16.2026072
## hp 31.4991029 67.7601565 64.548040 65.5995410
## hpl 0.5631187 28.2139028 28.070692 7.2740788
## george 16.9778364 35.7693596 38.520775 28.9517581
## num650 11.2438092 10.7478919 14.596898 7.3190271
## lab -8.9360976 19.5199376 19.372771 1.9241327
## labs 4.2785020 8.2201427 8.687946 2.1191606
## telnet 2.2491732 8.8969068 9.346173 0.9649704
## num857 0.8190727 7.9052724 8.042649 0.5664608
## data -2.5536741 8.9450087 5.547927 4.4758298
## num415 2.6899538 3.6647482 4.089579 0.2738140
## num85 4.6964260 6.2857471 7.123029 1.2681270
## technology 18.0540808 12.3067954 20.195553 7.4682756
## num1999 10.8959235 26.8137012 27.820515 9.9769696
## parts -0.6516205 8.3552941 6.561458 1.2131948
## pm 5.6746546 13.5804438 14.458598 5.5073597
## direct 4.4844411 2.2614370 4.740264 0.6561704
## cs 0.8271131 2.8452786 3.008445 0.5193167
## meeting 5.2054936 33.3063918 33.510794 9.0726166
## original -1.1300318 3.4338742 2.764312 0.4990267
## project 3.5687084 3.7956868 5.058334 2.0685303
## re 10.1173362 15.4814868 18.254379 10.6315605
## edu 26.0358121 80.1141936 79.113022 33.7404122
## table 1.3425112 0.8639931 1.492275 0.1664714
## conference 7.9498170 22.4471542 22.735973 4.0165862
## charSemicolon 8.5157965 8.1002914 11.720263 4.9244527
## charRoundbracket 8.5072320 17.9463037 19.491600 14.0817382
## charSquarebracket 2.3574388 4.5555113 3.972918 1.9406875
## charExclamation 51.0938315 41.6571227 68.471513 304.2808033
## charDollar 43.0037494 35.3829713 53.606872 326.4922243
## charHash 6.1926457 3.1911337 7.443912 2.7632185
## capitalAve 36.7954884 21.9632386 39.753905 68.3362847
## capitalLong 29.2341462 27.5832542 41.734219 46.1432726
## capitalTotal 24.2365427 14.2889067 29.397600 58.2561668
varImpPlot(spam.bag)
################################################################################
######################## Random Forest #########################################
################################################################################
#For classification random forest,
#sqrt(p) is a common number of variables considered at each split
#sqrt(57) is approximately 7.55
set.seed(23784)
spam.rf<-randomForest(type~.,
data=spam.train,
mtry=7,
importance=TRUE,
ntree=500)
#training predictions
pred.rf.train<-predict(spam.rf,
newdata=spam.train,
type="class")
conf.mat.rf.train<-table(spam.train$type,
pred.rf.train)
conf.mat.rf.train
## pred.rf.train
## nonspam spam
## nonspam 1852 1
## spam 11 1201
MCR.rf.train<-
1-sum(diag(conf.mat.rf.train))/sum(conf.mat.rf.train)
MCR.rf.train
## [1] 0.003915171
#holdout predictions
pred.rf<-predict(spam.rf,
newdata=spam.test,
type="class")
conf.mat.rf<-table(spam.test$type,
pred.rf)
conf.mat.rf
## pred.rf
## nonspam spam
## nonspam 910 25
## spam 47 554
MCR.rf<-
1-sum(diag(conf.mat.rf))/sum(conf.mat.rf)
MCR.rf
## [1] 0.046875
#Number incorrectly classified
1536-sum(diag(conf.mat.rf))
## [1] 72
#variable importance
importance(spam.rf)
## nonspam spam MeanDecreaseAccuracy MeanDecreaseGini
## make 4.172762312 7.730498 8.606268 6.4290947
## address 8.775436501 6.989352 10.185876 7.3984319
## all 3.429542992 13.497652 12.122856 12.6517634
## num3d 5.874022376 2.118182 5.905361 1.6734415
## our 17.370010409 21.739281 23.454021 44.9488924
## over 10.425509143 8.677616 12.824353 9.2365711
## remove 34.396460241 28.102982 36.837358 105.2911379
## internet 14.057782133 8.018537 14.480463 15.1322670
## order 8.370562647 6.129955 9.430022 6.6958514
## mail 6.399696975 10.586146 10.791954 11.9998776
## receive 14.498088715 6.968835 15.193493 13.5891028
## will 6.503649346 17.924424 18.668541 14.9806514
## people 3.669521165 7.454563 7.987547 5.8165125
## report 7.887721532 5.645060 9.183213 3.4231855
## addresses 7.237411910 5.752441 8.111511 2.6741580
## free 30.319691024 25.874944 34.191058 96.2002136
## business 15.982405274 11.737043 18.060900 16.9623703
## email 8.688529493 8.775227 11.534552 11.2719564
## you 13.565899965 17.534811 21.589365 36.7619148
## credit 10.476075838 4.319404 10.525087 6.4001370
## your 19.150353774 24.863154 28.951164 85.5409998
## font 10.629564609 8.547643 12.021160 3.7157930
## num000 19.144023742 11.275828 20.222333 35.2129670
## money 16.524854210 13.866646 18.428611 53.4728971
## hp 22.343878699 31.983010 33.900568 55.7070011
## hpl 13.414458885 21.709436 23.222424 26.9939141
## george 17.478682867 26.338498 28.727691 30.0739091
## num650 8.825311841 13.216057 15.344915 8.3044543
## lab -0.760073886 9.238654 9.257735 2.3728877
## labs 5.143684389 9.888046 10.723883 5.8399688
## telnet 3.916863476 8.778323 9.315714 2.7687222
## num857 2.861237911 6.959428 7.326002 1.4336598
## data 3.421868438 8.645324 8.946168 4.2508983
## num415 2.348389027 5.572533 6.056352 1.0793383
## num85 7.396792964 10.900817 12.619649 4.6745368
## technology 12.022905752 10.862526 14.842010 6.0165158
## num1999 17.245118055 23.364987 26.103243 23.4691994
## parts 0.003538345 4.681149 3.588038 0.7465097
## pm 5.695144790 12.052192 12.714622 5.4017764
## direct 5.197500170 2.415312 5.919370 1.7760713
## cs 2.991477790 6.497331 7.019599 1.3724744
## meeting 10.371722194 17.210066 18.679016 7.2812600
## original 1.125585677 6.803357 6.904095 1.6717831
## project 3.953626403 9.158794 9.627615 3.2448685
## re 11.262027858 21.137658 23.386036 14.3692627
## edu 22.139511918 25.827662 29.758170 25.2252029
## table 0.253394991 2.515123 2.396912 0.3314469
## conference 4.741402409 8.485712 9.223826 2.5321774
## charSemicolon 7.934700481 7.790897 11.072341 6.8531019
## charRoundbracket 8.657169191 17.703867 18.257575 17.4500285
## charSquarebracket 8.420872837 8.398282 10.661389 4.5936004
## charExclamation 34.697670170 38.968735 45.600515 182.8879880
## charDollar 32.826689506 28.901154 37.931025 146.1263182
## charHash 7.255782135 6.012007 8.693029 5.0509040
## capitalAve 29.241755112 28.234180 38.750399 89.2816847
## capitalLong 25.383483041 24.178059 33.084177 76.4753412
## capitalTotal 23.536835109 19.603281 30.082176 63.1651356
varImpPlot(spam.rf)
################################################################################
######################## 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
#gbm with Bernoulli loss needs a 0/1 outcome
spam.train.boost<-spam.train
spam.test.boost<-spam.test
spam.train.boost$spam01<-
ifelse(spam.train.boost$type=="spam",1,0)
spam.test.boost$spam01<-
ifelse(spam.test.boost$type=="spam",1,0)
#remove original factor outcome so it is not used as a predictor
spam.train.boost$type<-NULL
spam.test.boost$type<-NULL
set.seed(23784)
spam.boost<-gbm(spam01~.,
data=spam.train.boost,
distribution="bernoulli",
n.trees=5000,
interaction.depth=4,
verbose=FALSE)
summary(spam.boost)
## var rel.inf
## charDollar charDollar 1.929183e+01
## charExclamation charExclamation 1.822695e+01
## remove remove 1.018396e+01
## free free 6.207823e+00
## hp hp 5.271819e+00
## capitalAve capitalAve 4.766057e+00
## your your 4.633751e+00
## capitalTotal capitalTotal 4.235524e+00
## capitalLong capitalLong 3.708443e+00
## edu edu 3.013665e+00
## people people 2.795790e+00
## george george 2.761526e+00
## you you 1.911746e+00
## our our 1.770247e+00
## money money 1.582761e+00
## num000 num000 8.428473e-01
## hpl hpl 8.199242e-01
## will will 7.621925e-01
## num1999 num1999 7.398950e-01
## charRoundbracket charRoundbracket 5.606778e-01
## re re 5.510159e-01
## internet internet 5.459915e-01
## num650 num650 5.170344e-01
## meeting meeting 4.873904e-01
## charSemicolon charSemicolon 4.727749e-01
## receive receive 4.080249e-01
## all all 3.513144e-01
## technology technology 3.092416e-01
## business business 2.902337e-01
## email email 2.597703e-01
## mail mail 2.322399e-01
## make make 2.231470e-01
## report report 2.129757e-01
## over over 1.342410e-01
## project project 1.010395e-01
## order order 9.878013e-02
## num3d num3d 9.808688e-02
## font font 9.155955e-02
## num85 num85 8.565559e-02
## charHash charHash 8.562662e-02
## data data 7.062828e-02
## address address 5.746824e-02
## pm pm 4.582455e-02
## parts parts 4.249470e-02
## charSquarebracket charSquarebracket 4.082814e-02
## conference conference 3.665344e-02
## labs labs 2.119116e-02
## credit credit 1.827047e-02
## lab lab 9.647987e-03
## addresses addresses 6.973161e-03
## direct direct 6.414129e-03
## cs cs 2.986557e-05
## original original 8.272759e-06
## telnet telnet 0.000000e+00
## num857 num857 0.000000e+00
## num415 num415 0.000000e+00
## table table 0.000000e+00
#training predicted probabilities
prob.boost.train<-predict(spam.boost,
newdata=spam.train.boost,
n.trees=5000,
type="response")
#convert probabilities to classes using 0.5 cutoff
pred.boost.train<-ifelse(prob.boost.train>=0.5,
"spam",
"nonspam")
pred.boost.train<-factor(pred.boost.train,
levels=levels(spam.train$type))
conf.mat.boost.train<-table(spam.train$type,
pred.boost.train)
conf.mat.boost.train
## pred.boost.train
## nonspam spam
## nonspam 1852 1
## spam 0 1212
MCR.boost.train<-
1-sum(diag(conf.mat.boost.train))/sum(conf.mat.boost.train)
MCR.boost.train
## [1] 0.0003262643
#holdout predicted probabilities
prob.boost<-predict(spam.boost,
newdata=spam.test.boost,
n.trees=5000,
type="response")
pred.boost<-ifelse(prob.boost>=0.5,
"spam",
"nonspam")
pred.boost<-factor(pred.boost,
levels=levels(spam.test$type))
conf.mat.boost<-table(spam.test$type,
pred.boost)
conf.mat.boost
## pred.boost
## nonspam spam
## nonspam 904 31
## spam 36 565
MCR.boost<-
1-sum(diag(conf.mat.boost))/sum(conf.mat.boost)
MCR.boost
## [1] 0.04361979
#Number incorrectly classified
1536-sum(diag(conf.mat.boost))
## [1] 67
################################################################################
######################## Gini Tree Performance #################################
################################################################################
pred.Gini<-predict(spam.tree.Gini,
newdata=spam.test,
type="class")
conf.mat.Gini<-table(spam.test$type,
pred.Gini)
conf.mat.Gini
## pred.Gini
## nonspam spam
## nonspam 804 131
## spam 248 353
MCR.Gini<-
1-sum(diag(conf.mat.Gini))/sum(conf.mat.Gini)
MCR.Gini
## [1] 0.2467448
################################################################################
######################## Final Comparison ######################################
################################################################################
#Default tree
MCR.tree<-
1-sum(diag(conf.mat))/sum(conf.mat)
#Pruned default tree
MCR.prune<-
1-sum(diag(conf.mat.prune))/sum(conf.mat.prune)
#Large tree
MCR.ctrl<-
1-sum(diag(conf.mat.ctrl))/sum(conf.mat.ctrl)
#Pruned large tree
MCR.prune.ctrl<-
1-sum(diag(conf.mat.prune.ctrl))/sum(conf.mat.prune.ctrl)
results<-data.frame(
Model=c("Default tree",
"Pruned tree",
"Large tree",
"Pruned large tree",
"Gini tree",
"Bagging",
"Random forest",
"Boosting"),
Holdout.Misclassification=c(
MCR.tree,
MCR.prune,
MCR.ctrl,
MCR.prune.ctrl,
MCR.Gini,
MCR.bag,
MCR.rf,
MCR.boost
)
)
results
## Model Holdout.Misclassification
## 1 Default tree 0.09114583
## 2 Pruned tree 0.09114583
## 3 Large tree 0.08854167
## 4 Pruned large tree 0.08789062
## 5 Gini tree 0.24674479
## 6 Bagging 0.05598958
## 7 Random forest 0.04687500
## 8 Boosting 0.04361979
results[order(results$Holdout.Misclassification),]
## Model Holdout.Misclassification
## 8 Boosting 0.04361979
## 7 Random forest 0.04687500
## 6 Bagging 0.05598958
## 4 Pruned large tree 0.08789062
## 3 Large tree 0.08854167
## 1 Default tree 0.09114583
## 2 Pruned tree 0.09114583
## 5 Gini tree 0.24674479
##This is analysis of email vs spam data set
library(tree)
library(kernlab)
data(spam)
class(spam)
## [1] "data.frame"
str(spam)
## 'data.frame': 4601 obs. of 58 variables:
## $ make : num 0 0.21 0.06 0 0 0 0 0 0.15 0.06 ...
## $ address : num 0.64 0.28 0 0 0 0 0 0 0 0.12 ...
## $ all : num 0.64 0.5 0.71 0 0 0 0 0 0.46 0.77 ...
## $ num3d : num 0 0 0 0 0 0 0 0 0 0 ...
## $ our : num 0.32 0.14 1.23 0.63 0.63 1.85 1.92 1.88 0.61 0.19 ...
## $ over : num 0 0.28 0.19 0 0 0 0 0 0 0.32 ...
## $ remove : num 0 0.21 0.19 0.31 0.31 0 0 0 0.3 0.38 ...
## $ internet : num 0 0.07 0.12 0.63 0.63 1.85 0 1.88 0 0 ...
## $ order : num 0 0 0.64 0.31 0.31 0 0 0 0.92 0.06 ...
## $ mail : num 0 0.94 0.25 0.63 0.63 0 0.64 0 0.76 0 ...
## $ receive : num 0 0.21 0.38 0.31 0.31 0 0.96 0 0.76 0 ...
## $ will : num 0.64 0.79 0.45 0.31 0.31 0 1.28 0 0.92 0.64 ...
## $ people : num 0 0.65 0.12 0.31 0.31 0 0 0 0 0.25 ...
## $ report : num 0 0.21 0 0 0 0 0 0 0 0 ...
## $ addresses : num 0 0.14 1.75 0 0 0 0 0 0 0.12 ...
## $ free : num 0.32 0.14 0.06 0.31 0.31 0 0.96 0 0 0 ...
## $ business : num 0 0.07 0.06 0 0 0 0 0 0 0 ...
## $ email : num 1.29 0.28 1.03 0 0 0 0.32 0 0.15 0.12 ...
## $ you : num 1.93 3.47 1.36 3.18 3.18 0 3.85 0 1.23 1.67 ...
## $ credit : num 0 0 0.32 0 0 0 0 0 3.53 0.06 ...
## $ your : num 0.96 1.59 0.51 0.31 0.31 0 0.64 0 2 0.71 ...
## $ font : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num000 : num 0 0.43 1.16 0 0 0 0 0 0 0.19 ...
## $ money : num 0 0.43 0.06 0 0 0 0 0 0.15 0 ...
## $ hp : num 0 0 0 0 0 0 0 0 0 0 ...
## $ hpl : num 0 0 0 0 0 0 0 0 0 0 ...
## $ george : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num650 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ lab : num 0 0 0 0 0 0 0 0 0 0 ...
## $ labs : num 0 0 0 0 0 0 0 0 0 0 ...
## $ telnet : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num857 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ data : num 0 0 0 0 0 0 0 0 0.15 0 ...
## $ num415 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num85 : num 0 0 0 0 0 0 0 0 0 0 ...
## $ technology : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num1999 : num 0 0.07 0 0 0 0 0 0 0 0 ...
## $ parts : num 0 0 0 0 0 0 0 0 0 0 ...
## $ pm : num 0 0 0 0 0 0 0 0 0 0 ...
## $ direct : num 0 0 0.06 0 0 0 0 0 0 0 ...
## $ cs : num 0 0 0 0 0 0 0 0 0 0 ...
## $ meeting : num 0 0 0 0 0 0 0 0 0 0 ...
## $ original : num 0 0 0.12 0 0 0 0 0 0.3 0 ...
## $ project : num 0 0 0 0 0 0 0 0 0 0.06 ...
## $ re : num 0 0 0.06 0 0 0 0 0 0 0 ...
## $ edu : num 0 0 0.06 0 0 0 0 0 0 0 ...
## $ table : num 0 0 0 0 0 0 0 0 0 0 ...
## $ conference : num 0 0 0 0 0 0 0 0 0 0 ...
## $ charSemicolon : num 0 0 0.01 0 0 0 0 0 0 0.04 ...
## $ charRoundbracket : num 0 0.132 0.143 0.137 0.135 0.223 0.054 0.206 0.271 0.03 ...
## $ charSquarebracket: num 0 0 0 0 0 0 0 0 0 0 ...
## $ charExclamation : num 0.778 0.372 0.276 0.137 0.135 0 0.164 0 0.181 0.244 ...
## $ charDollar : num 0 0.18 0.184 0 0 0 0.054 0 0.203 0.081 ...
## $ charHash : num 0 0.048 0.01 0 0 0 0 0 0.022 0 ...
## $ capitalAve : num 3.76 5.11 9.82 3.54 3.54 ...
## $ capitalLong : num 61 101 485 40 40 15 4 11 445 43 ...
## $ capitalTotal : num 278 1028 2259 191 191 ...
## $ type : Factor w/ 2 levels "nonspam","spam": 2 2 2 2 2 2 2 2 2 2 ...
#The variable "type" is the outcome variable
# I'll produce my own training and test sets of the size described in the book
set.seed(3731)
train.ind<-sort(sample(1:4601, 3065, replace=FALSE))
test.ind<-seq(1,4601)[-train.ind]
spam.train<-spam[train.ind,]
spam.test<-spam[test.ind,]
################################################################################
######################## Default Classification Tree ###########################
################################################################################
#Fit a classification tree to the training data
spam.tree<-tree(type~.,data=spam.train)
plot(spam.tree)
text(spam.tree)
summary(spam.tree)
##
## Classification tree:
## tree(formula = type ~ ., data = spam.train)
## Variables actually used in tree construction:
## [1] "charDollar" "remove" "charExclamation" "george"
## [5] "hp" "capitalLong" "edu" "our"
## [9] "capitalTotal"
## Number of terminal nodes: 15
## Residual mean deviance: 0.4404 = 1343 / 3050
## Misclassification error rate: 0.07732 = 237 / 3065
conf.mat.train<-table(spam.train$type,
predict(spam.tree,type='class'))
conf.mat.train
##
## nonspam spam
## nonspam 1775 78
## spam 159 1053
#How well does this tree perform on the holdout set
pred.tree<-predict(spam.tree,
newdata=spam.test,
type="class")
conf.mat<-table(spam.test$type,pred.tree)
conf.mat
## pred.tree
## nonspam spam
## nonspam 888 47
## spam 93 508
#misclassification rate
1-sum(diag(conf.mat))/sum(conf.mat)
## [1] 0.09114583
1536-sum(diag(conf.mat)) #Number of incorrect calls
## [1] 140
################################################################################
######################## Prune Default Tree ####################################
################################################################################
#Cost complexity pruning the default tree shows misclassification best for full tree
spam.cv<-cv.tree(spam.tree,
FUN=prune.misclass)
par(mfrow=c(1,2))
plot(spam.cv$size,
spam.cv$dev,
type="b")
plot(spam.cv$k,
spam.cv$dev,
type="b")
spam.cv
## $size
## [1] 15 14 9 8 7 6 5 3 2 1
##
## $dev
## [1] 299 299 309 331 333 337 368 435 644 1212
##
## $k
## [1] -Inf 0.0 7.6 12.0 13.0 14.0 34.0 48.0 177.0 591.0
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
min(spam.cv$dev)
## [1] 299
which(spam.cv$dev==min(spam.cv$dev))
## [1] 1 2
spam.cv$size[2] #size of best tree that is smallest in that class
## [1] 14
#prune down and see how the size 14 tree works
#DOUBLE CHECK best above
spam.prune<-prune.misclass(spam.tree,
best=14)
par(mfrow=c(1,2))
plot(spam.tree) #original tree
plot(spam.prune) #pruned tree
#check in-sample error rate of pruned tree
summary(spam.prune)
##
## Classification tree:
## snip.tree(tree = spam.tree, nodes = 12L)
## Variables actually used in tree construction:
## [1] "charDollar" "remove" "charExclamation" "george"
## [5] "hp" "capitalLong" "edu" "our"
## [9] "capitalTotal"
## Number of terminal nodes: 14
## Residual mean deviance: 0.4563 = 1392 / 3051
## Misclassification error rate: 0.07732 = 237 / 3065
#check out-of-sample error rate of pruned tree
pred.prune<-predict(spam.prune,
newdata=spam.test,
type="class")
conf.mat.prune<-table(spam.test$type,
pred.prune)
conf.mat.prune
## pred.prune
## nonspam spam
## nonspam 888 47
## spam 93 508
#misclassification rate
1-sum(diag(conf.mat.prune))/sum(conf.mat.prune)
## [1] 0.09114583
1536-sum(diag(conf.mat.prune))
## [1] 140
################################################################################
######################## Large Classification Tree #############################
################################################################################
#Build a larger tree for cost complexity pruning
spam.tree.ctrl<-tree(type~.,
data=spam.train,
control=tree.control(nobs=3065,
minsize=2,
mindev=0))
par(mfrow=c(1,1))
plot(spam.tree.ctrl)
text(spam.tree.ctrl)
summary(spam.tree.ctrl)
##
## Classification tree:
## tree(formula = type ~ ., data = spam.train, control = tree.control(nobs = 3065,
## minsize = 2, mindev = 0))
## Variables actually used in tree construction:
## [1] "charDollar" "remove" "charExclamation" "george"
## [5] "hp" "capitalLong" "free" "internet"
## [9] "capitalAve" "you" "edu" "re"
## [13] "capitalTotal" "our" "mail" "all"
## [17] "technology" "will" "your" "charSemicolon"
## [21] "meeting" "address" "num3d" "num650"
## [25] "data" "charRoundbracket" "make" "hpl"
## [29] "over" "email" "charHash" "report"
## [33] "num1999" "business" "receive" "people"
## Number of terminal nodes: 176
## Residual mean deviance: 0.0009597 = 2.773 / 2889
## Misclassification error rate: 0.0003263 = 1 / 3065
#Out of sample misclassification of the large tree
#How well does this tree perform on the holdout set
pred.tree.ctrl<-predict(spam.tree.ctrl,
newdata=spam.test,
type="class")
conf.mat.ctrl<-table(spam.test$type,
pred.tree.ctrl)
conf.mat.ctrl
## pred.tree.ctrl
## nonspam spam
## nonspam 870 65
## spam 71 530
#misclassification rate
1-sum(diag(conf.mat.ctrl))/sum(conf.mat.ctrl)
## [1] 0.08854167
1536-sum(diag(conf.mat.ctrl))
## [1] 136
#Well growing the biggest possible tree overfit relative to the training error rate
#but was superior to the smaller tree in the out-of-sample misclassification rate.
################################################################################
######################## Prune Large Tree ######################################
################################################################################
spam.ctrl.cv<-cv.tree(spam.tree.ctrl,
FUN=prune.misclass)
par(mfrow=c(1,2))
plot(spam.ctrl.cv$size,
spam.ctrl.cv$dev,
type="b")
plot(spam.ctrl.cv$k,
spam.ctrl.cv$dev,
type="b")
spam.ctrl.cv
## $size
## [1] 176 175 170 136 133 125 116 111 104 76 66 64 59 51 39 36 34 30 28
## [20] 24 19 18 16 11 10 9 8 7 6 5 3 2 1
##
## $dev
## [1] 293 293 293 293 293 293 293 293 293 293 293 293 293 293 293
## [16] 293 293 293 293 293 293 296 297 300 296 321 324 328 328 368
## [31] 443 638 1212
##
## $k
## [1] -Inf 0.0000000 0.2000000 0.5000000 0.6666667 0.7500000
## [7] 0.7777778 0.8000000 0.8571429 1.0000000 1.2000000 1.5000000
## [13] 1.6000000 1.7500000 2.0000000 2.3333333 2.5000000 2.7500000
## [19] 3.0000000 3.7500000 4.0000000 5.0000000 6.0000000 7.6000000
## [25] 8.0000000 12.0000000 13.0000000 14.0000000 15.0000000 34.0000000
## [31] 48.0000000 177.0000000 591.0000000
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
min(spam.ctrl.cv$dev)
## [1] 293
which(spam.ctrl.cv$dev==min(spam.ctrl.cv$dev))
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
#Find the smallest tree among those tied for minimum CV error
best.ind<-which(spam.ctrl.cv$dev==
min(spam.ctrl.cv$dev))
best.size<-min(spam.ctrl.cv$size[best.ind])
best.size
## [1] 19
#prune large tree using selected size
spam.prune.ctrl<-prune.misclass(spam.tree.ctrl,
best=best.size)
par(mfrow=c(1,2))
plot(spam.tree.ctrl) #original large tree
plot(spam.prune.ctrl) #pruned tree
#check in-sample error rate of pruned tree
summary(spam.prune.ctrl)
##
## Classification tree:
## snip.tree(tree = spam.tree.ctrl, nodes = c(33L, 73L, 262L, 10L,
## 76L, 263L, 37L, 260L, 12L, 64L, 72L, 261L))
## Variables actually used in tree construction:
## [1] "charDollar" "remove" "charExclamation" "george"
## [5] "hp" "capitalLong" "edu" "our"
## [9] "num650" "capitalTotal" "free" "business"
## [13] "email"
## Number of terminal nodes: 19
## Residual mean deviance: 0.4056 = 1235 / 3046
## Misclassification error rate: 0.06427 = 197 / 3065
#check out-of-sample error rate of pruned tree
pred.prune.ctrl<-predict(spam.prune.ctrl,
newdata=spam.test,
type="class")
conf.mat.prune.ctrl<-table(spam.test$type,
pred.prune.ctrl)
conf.mat.prune.ctrl
## pred.prune.ctrl
## nonspam spam
## nonspam 882 53
## spam 64 537
#misclassification rate
1-sum(diag(conf.mat.prune.ctrl))/sum(conf.mat.prune.ctrl)
## [1] 0.07617188
1536-sum(diag(conf.mat.prune.ctrl))
## [1] 117
################################################################################
######################## Gini Index Version ####################################
################################################################################
###Revisit the Gini index version
#This unrestricted version reaches maximum depth:
#spam.tree.Gini<-tree(type~.,
# data=spam.train,
# split="gini")
#Use controls to avoid maximum depth error
spam.tree.Gini<-tree(type~.,
data=spam.train,
split="gini",
control=tree.control(nobs=3065,
mincut=100,
minsize=200,
mindev=.01))
summary(spam.tree.Gini)
##
## Classification tree:
## tree(formula = type ~ ., data = spam.train, control = tree.control(nobs = 3065,
## mincut = 100, minsize = 200, mindev = 0.01), split = "gini")
## Variables actually used in tree construction:
## [1] "conference" "project" "num415"
## [4] "original" "pm" "technology"
## [7] "num85" "data" "charSquarebracket"
## [10] "report" "charSemicolon" "num1999"
## [13] "hp" "credit" "receive"
## [16] "people" "make" "address"
## [19] "all" "will" "re"
## [22] "charExclamation" "capitalTotal"
## Number of terminal nodes: 24
## Residual mean deviance: 1.051 = 3197 / 3041
## Misclassification error rate: 0.262 = 803 / 3065
summary(spam.tree)
##
## Classification tree:
## tree(formula = type ~ ., data = spam.train)
## Variables actually used in tree construction:
## [1] "charDollar" "remove" "charExclamation" "george"
## [5] "hp" "capitalLong" "edu" "our"
## [9] "capitalTotal"
## Number of terminal nodes: 15
## Residual mean deviance: 0.4404 = 1343 / 3050
## Misclassification error rate: 0.07732 = 237 / 3065
par(mfrow=c(1,1))
plot(spam.tree)
text(spam.tree)
plot(spam.tree.Gini)
text(spam.tree.Gini)
#Compare the various in- and out-of-sample misclassification rates - expectations met?
################################################################################
######################## Gini Tree Performance #################################
################################################################################
pred.Gini<-predict(spam.tree.Gini,
newdata=spam.test,
type="class")
conf.mat.Gini<-table(spam.test$type,
pred.Gini)
conf.mat.Gini
## pred.Gini
## nonspam spam
## nonspam 804 131
## spam 248 353
MCR.Gini<-
1-sum(diag(conf.mat.Gini))/sum(conf.mat.Gini)
MCR.Gini
## [1] 0.2467448
1536-sum(diag(conf.mat.Gini))
## [1] 379
################################################################################
######################## Bagging ################################################
################################################################################
library(randomForest)
#There are 57 candidate predictors
#For bagging, allow all predictors to be considered at every split
set.seed(23784)
spam.bag<-randomForest(type~.,
data=spam.train,
mtry=57,
importance=TRUE,
ntree=500)
#training predictions
pred.bag.train<-predict(spam.bag,
newdata=spam.train,
type="class")
conf.mat.bag.train<-table(spam.train$type,
pred.bag.train)
conf.mat.bag.train
## pred.bag.train
## nonspam spam
## nonspam 1852 1
## spam 0 1212
#training misclassification rate
MCR.bag.train<-
1-sum(diag(conf.mat.bag.train))/sum(conf.mat.bag.train)
MCR.bag.train
## [1] 0.0003262643
#holdout predictions
pred.bag<-predict(spam.bag,
newdata=spam.test,
type="class")
conf.mat.bag<-table(spam.test$type,
pred.bag)
conf.mat.bag
## pred.bag
## nonspam spam
## nonspam 899 36
## spam 50 551
#holdout misclassification rate
MCR.bag<-
1-sum(diag(conf.mat.bag))/sum(conf.mat.bag)
MCR.bag
## [1] 0.05598958
#Number incorrectly classified
1536-sum(diag(conf.mat.bag))
## [1] 86
#variable importance
importance(spam.bag)
## nonspam spam MeanDecreaseAccuracy MeanDecreaseGini
## make 3.4375107 7.1340748 7.482192 4.5451238
## address 8.3856032 7.0315143 11.220811 3.4421350
## all 1.7399876 6.5128243 6.859883 5.8280717
## num3d 11.2614372 -0.6113121 9.583824 2.8370080
## our 25.6120408 24.1997680 33.596296 23.5780096
## over 16.1088210 5.2778569 16.357111 4.9756630
## remove 71.8486869 43.8974403 79.650696 180.6049774
## internet 15.1688083 6.2211822 15.190525 8.0687118
## order 11.4536096 1.9803425 12.250713 3.5688217
## mail 5.7479403 5.5055160 8.060951 7.3231617
## receive 17.2489085 2.0130431 17.464684 4.4653853
## will 2.8378200 11.4937870 11.110846 9.2491864
## people -0.9692594 3.4124512 2.245570 3.9128499
## report 10.3119508 8.1193903 12.179344 3.8215417
## addresses 5.3310266 3.0249328 6.265170 1.1292987
## free 36.5417447 30.3690826 44.990199 55.4029369
## business 19.7444486 9.2738823 20.978700 11.3550531
## email 13.0181662 9.0610999 14.245406 9.4070060
## you 14.4237939 15.8286180 22.036126 28.0117604
## credit 12.8634464 -0.4241135 12.880596 2.5884365
## your 15.2606515 15.2803819 20.378356 27.4350838
## font 15.6207333 12.5478591 17.354343 3.4502325
## num000 21.1644025 3.6890746 21.741198 9.3520811
## money 18.3806807 13.6155787 20.588810 16.2026072
## hp 31.4991029 67.7601565 64.548040 65.5995410
## hpl 0.5631187 28.2139028 28.070692 7.2740788
## george 16.9778364 35.7693596 38.520775 28.9517581
## num650 11.2438092 10.7478919 14.596898 7.3190271
## lab -8.9360976 19.5199376 19.372771 1.9241327
## labs 4.2785020 8.2201427 8.687946 2.1191606
## telnet 2.2491732 8.8969068 9.346173 0.9649704
## num857 0.8190727 7.9052724 8.042649 0.5664608
## data -2.5536741 8.9450087 5.547927 4.4758298
## num415 2.6899538 3.6647482 4.089579 0.2738140
## num85 4.6964260 6.2857471 7.123029 1.2681270
## technology 18.0540808 12.3067954 20.195553 7.4682756
## num1999 10.8959235 26.8137012 27.820515 9.9769696
## parts -0.6516205 8.3552941 6.561458 1.2131948
## pm 5.6746546 13.5804438 14.458598 5.5073597
## direct 4.4844411 2.2614370 4.740264 0.6561704
## cs 0.8271131 2.8452786 3.008445 0.5193167
## meeting 5.2054936 33.3063918 33.510794 9.0726166
## original -1.1300318 3.4338742 2.764312 0.4990267
## project 3.5687084 3.7956868 5.058334 2.0685303
## re 10.1173362 15.4814868 18.254379 10.6315605
## edu 26.0358121 80.1141936 79.113022 33.7404122
## table 1.3425112 0.8639931 1.492275 0.1664714
## conference 7.9498170 22.4471542 22.735973 4.0165862
## charSemicolon 8.5157965 8.1002914 11.720263 4.9244527
## charRoundbracket 8.5072320 17.9463037 19.491600 14.0817382
## charSquarebracket 2.3574388 4.5555113 3.972918 1.9406875
## charExclamation 51.0938315 41.6571227 68.471513 304.2808033
## charDollar 43.0037494 35.3829713 53.606872 326.4922243
## charHash 6.1926457 3.1911337 7.443912 2.7632185
## capitalAve 36.7954884 21.9632386 39.753905 68.3362847
## capitalLong 29.2341462 27.5832542 41.734219 46.1432726
## capitalTotal 24.2365427 14.2889067 29.397600 58.2561668
varImpPlot(spam.bag)
################################################################################
######################## Random Forest #########################################
################################################################################
#For classification random forest,
#sqrt(p) is a common number of variables considered at each split
#sqrt(57) is approximately 7.55
set.seed(23784)
spam.rf<-randomForest(type~.,
data=spam.train,
mtry=7,
importance=TRUE,
ntree=500)
#training predictions
pred.rf.train<-predict(spam.rf,
newdata=spam.train,
type="class")
conf.mat.rf.train<-table(spam.train$type,
pred.rf.train)
conf.mat.rf.train
## pred.rf.train
## nonspam spam
## nonspam 1852 1
## spam 11 1201
#training misclassification rate
MCR.rf.train<-
1-sum(diag(conf.mat.rf.train))/sum(conf.mat.rf.train)
MCR.rf.train
## [1] 0.003915171
#holdout predictions
pred.rf<-predict(spam.rf,
newdata=spam.test,
type="class")
conf.mat.rf<-table(spam.test$type,
pred.rf)
conf.mat.rf
## pred.rf
## nonspam spam
## nonspam 910 25
## spam 47 554
#holdout misclassification rate
MCR.rf<-
1-sum(diag(conf.mat.rf))/sum(conf.mat.rf)
MCR.rf
## [1] 0.046875
#Number incorrectly classified
1536-sum(diag(conf.mat.rf))
## [1] 72
#variable importance
importance(spam.rf)
## nonspam spam MeanDecreaseAccuracy MeanDecreaseGini
## make 4.172762312 7.730498 8.606268 6.4290947
## address 8.775436501 6.989352 10.185876 7.3984319
## all 3.429542992 13.497652 12.122856 12.6517634
## num3d 5.874022376 2.118182 5.905361 1.6734415
## our 17.370010409 21.739281 23.454021 44.9488924
## over 10.425509143 8.677616 12.824353 9.2365711
## remove 34.396460241 28.102982 36.837358 105.2911379
## internet 14.057782133 8.018537 14.480463 15.1322670
## order 8.370562647 6.129955 9.430022 6.6958514
## mail 6.399696975 10.586146 10.791954 11.9998776
## receive 14.498088715 6.968835 15.193493 13.5891028
## will 6.503649346 17.924424 18.668541 14.9806514
## people 3.669521165 7.454563 7.987547 5.8165125
## report 7.887721532 5.645060 9.183213 3.4231855
## addresses 7.237411910 5.752441 8.111511 2.6741580
## free 30.319691024 25.874944 34.191058 96.2002136
## business 15.982405274 11.737043 18.060900 16.9623703
## email 8.688529493 8.775227 11.534552 11.2719564
## you 13.565899965 17.534811 21.589365 36.7619148
## credit 10.476075838 4.319404 10.525087 6.4001370
## your 19.150353774 24.863154 28.951164 85.5409998
## font 10.629564609 8.547643 12.021160 3.7157930
## num000 19.144023742 11.275828 20.222333 35.2129670
## money 16.524854210 13.866646 18.428611 53.4728971
## hp 22.343878699 31.983010 33.900568 55.7070011
## hpl 13.414458885 21.709436 23.222424 26.9939141
## george 17.478682867 26.338498 28.727691 30.0739091
## num650 8.825311841 13.216057 15.344915 8.3044543
## lab -0.760073886 9.238654 9.257735 2.3728877
## labs 5.143684389 9.888046 10.723883 5.8399688
## telnet 3.916863476 8.778323 9.315714 2.7687222
## num857 2.861237911 6.959428 7.326002 1.4336598
## data 3.421868438 8.645324 8.946168 4.2508983
## num415 2.348389027 5.572533 6.056352 1.0793383
## num85 7.396792964 10.900817 12.619649 4.6745368
## technology 12.022905752 10.862526 14.842010 6.0165158
## num1999 17.245118055 23.364987 26.103243 23.4691994
## parts 0.003538345 4.681149 3.588038 0.7465097
## pm 5.695144790 12.052192 12.714622 5.4017764
## direct 5.197500170 2.415312 5.919370 1.7760713
## cs 2.991477790 6.497331 7.019599 1.3724744
## meeting 10.371722194 17.210066 18.679016 7.2812600
## original 1.125585677 6.803357 6.904095 1.6717831
## project 3.953626403 9.158794 9.627615 3.2448685
## re 11.262027858 21.137658 23.386036 14.3692627
## edu 22.139511918 25.827662 29.758170 25.2252029
## table 0.253394991 2.515123 2.396912 0.3314469
## conference 4.741402409 8.485712 9.223826 2.5321774
## charSemicolon 7.934700481 7.790897 11.072341 6.8531019
## charRoundbracket 8.657169191 17.703867 18.257575 17.4500285
## charSquarebracket 8.420872837 8.398282 10.661389 4.5936004
## charExclamation 34.697670170 38.968735 45.600515 182.8879880
## charDollar 32.826689506 28.901154 37.931025 146.1263182
## charHash 7.255782135 6.012007 8.693029 5.0509040
## capitalAve 29.241755112 28.234180 38.750399 89.2816847
## capitalLong 25.383483041 24.178059 33.084177 76.4753412
## capitalTotal 23.536835109 19.603281 30.082176 63.1651356
varImpPlot(spam.rf)
################################################################################
######################## Boosting ##############################################
################################################################################
library(gbm)
#gbm with Bernoulli loss needs a 0/1 outcome
spam.train.boost<-spam.train
spam.test.boost<-spam.test
spam.train.boost$spam01<-
ifelse(spam.train.boost$type=="spam",1,0)
spam.test.boost$spam01<-
ifelse(spam.test.boost$type=="spam",1,0)
#remove original factor outcome so it is not used as a predictor
spam.train.boost$type<-NULL
spam.test.boost$type<-NULL
set.seed(23784)
spam.boost<-gbm(spam01~.,
data=spam.train.boost,
distribution="bernoulli",
n.trees=5000,
interaction.depth=4,
verbose=FALSE)
summary(spam.boost)
## var rel.inf
## charDollar charDollar 1.929183e+01
## charExclamation charExclamation 1.822695e+01
## remove remove 1.018396e+01
## free free 6.207823e+00
## hp hp 5.271819e+00
## capitalAve capitalAve 4.766057e+00
## your your 4.633751e+00
## capitalTotal capitalTotal 4.235524e+00
## capitalLong capitalLong 3.708443e+00
## edu edu 3.013665e+00
## people people 2.795790e+00
## george george 2.761526e+00
## you you 1.911746e+00
## our our 1.770247e+00
## money money 1.582761e+00
## num000 num000 8.428473e-01
## hpl hpl 8.199242e-01
## will will 7.621925e-01
## num1999 num1999 7.398950e-01
## charRoundbracket charRoundbracket 5.606778e-01
## re re 5.510159e-01
## internet internet 5.459915e-01
## num650 num650 5.170344e-01
## meeting meeting 4.873904e-01
## charSemicolon charSemicolon 4.727749e-01
## receive receive 4.080249e-01
## all all 3.513144e-01
## technology technology 3.092416e-01
## business business 2.902337e-01
## email email 2.597703e-01
## mail mail 2.322399e-01
## make make 2.231470e-01
## report report 2.129757e-01
## over over 1.342410e-01
## project project 1.010395e-01
## order order 9.878013e-02
## num3d num3d 9.808688e-02
## font font 9.155955e-02
## num85 num85 8.565559e-02
## charHash charHash 8.562662e-02
## data data 7.062828e-02
## address address 5.746824e-02
## pm pm 4.582455e-02
## parts parts 4.249470e-02
## charSquarebracket charSquarebracket 4.082814e-02
## conference conference 3.665344e-02
## labs labs 2.119116e-02
## credit credit 1.827047e-02
## lab lab 9.647987e-03
## addresses addresses 6.973161e-03
## direct direct 6.414129e-03
## cs cs 2.986557e-05
## original original 8.272759e-06
## telnet telnet 0.000000e+00
## num857 num857 0.000000e+00
## num415 num415 0.000000e+00
## table table 0.000000e+00
#training predicted probabilities
prob.boost.train<-predict(spam.boost,
newdata=spam.train.boost,
n.trees=5000,
type="response")
#convert probabilities to classes using 0.5 cutoff
pred.boost.train<-ifelse(prob.boost.train>=0.5,
"spam",
"nonspam")
pred.boost.train<-factor(pred.boost.train,
levels=levels(spam.train$type))
conf.mat.boost.train<-table(spam.train$type,
pred.boost.train)
conf.mat.boost.train
## pred.boost.train
## nonspam spam
## nonspam 1852 1
## spam 0 1212
#training misclassification rate
MCR.boost.train<-
1-sum(diag(conf.mat.boost.train))/sum(conf.mat.boost.train)
MCR.boost.train
## [1] 0.0003262643
#holdout predicted probabilities
prob.boost<-predict(spam.boost,
newdata=spam.test.boost,
n.trees=5000,
type="response")
pred.boost<-ifelse(prob.boost>=0.5,
"spam",
"nonspam")
pred.boost<-factor(pred.boost,
levels=levels(spam.test$type))
conf.mat.boost<-table(spam.test$type,
pred.boost)
conf.mat.boost
## pred.boost
## nonspam spam
## nonspam 904 31
## spam 36 565
#holdout misclassification rate
MCR.boost<-
1-sum(diag(conf.mat.boost))/sum(conf.mat.boost)
MCR.boost
## [1] 0.04361979
#Number incorrectly classified
1536-sum(diag(conf.mat.boost))
## [1] 67
################################################################################
######################## Final Comparison ######################################
################################################################################
#Default tree
MCR.tree<-
1-sum(diag(conf.mat))/sum(conf.mat)
#Pruned default tree
MCR.prune<-
1-sum(diag(conf.mat.prune))/sum(conf.mat.prune)
#Large tree
MCR.ctrl<-
1-sum(diag(conf.mat.ctrl))/sum(conf.mat.ctrl)
#Pruned large tree
MCR.prune.ctrl<-
1-sum(diag(conf.mat.prune.ctrl))/sum(conf.mat.prune.ctrl)
results<-data.frame(
Model=c("Default tree",
"Pruned tree",
"Large tree",
"Pruned large tree",
"Gini tree",
"Bagging",
"Random forest",
"Boosting"),
Holdout.Misclassification=c(
MCR.tree,
MCR.prune,
MCR.ctrl,
MCR.prune.ctrl,
MCR.Gini,
MCR.bag,
MCR.rf,
MCR.boost
)
)
results
## Model Holdout.Misclassification
## 1 Default tree 0.09114583
## 2 Pruned tree 0.09114583
## 3 Large tree 0.08854167
## 4 Pruned large tree 0.07617188
## 5 Gini tree 0.24674479
## 6 Bagging 0.05598958
## 7 Random forest 0.04687500
## 8 Boosting 0.04361979
#Rank from best to worst
results[order(results$Holdout.Misclassification),]
## Model Holdout.Misclassification
## 8 Boosting 0.04361979
## 7 Random forest 0.04687500
## 6 Bagging 0.05598958
## 4 Pruned large tree 0.07617188
## 3 Large tree 0.08854167
## 1 Default tree 0.09114583
## 2 Pruned tree 0.09114583
## 5 Gini tree 0.24674479