Tutorial 4: Subset Selection & LASSO

Ebeer

# set working directory using however you want to folder where data is stored.  I'll use 
ebeer <- read_csv("ebeer.csv")

# load ebeer, remove account number column
ebeer<-ebeer[-c(1)]

Training-Test Samples:

# drop the ID column, select customers that received a mailing only
ebeer_test<-subset(ebeer, mailing ==1)

# create ebeer rollout data
ebeer_rollout<-subset(ebeer, mailing ==0)

# rename ebeer_test ebeer
ebeer<-ebeer_test

Telco

# load telco
telco <- read_csv("telco.csv")
telco <- strings2factors(telco)

## 
## The following character variables were converted to factors
##  gender Partner Dependents PhoneService MultipleLines InternetService OnlineSecurity OnlineBackup DeviceProtection TechSupport StreamingTV StreamingMovies Contract PaperlessBilling PaymentMethod Churn

# drop ID column, divide Total charges by 1000
telco<-subset(telco, select=-customerID)
telco$TotalCharges<-telco$TotalCharges/1000

Training-Test:

# create 70% test and 30% holdout sample
set.seed(19103)
n <- nrow(telco)
sample <- sample(c(TRUE, FALSE), n, replace=TRUE, prob=c(0.7, 0.3))
telco.test <- telco[sample, ]
telco.holdout <- telco[!sample, ]

#call test telco, and full data set telco.all
telco.all<-telco
telco<-telco.test

Decision Trees

We’ll use ebeer for the trees. I first run a tree and graph it so we can talk about what it is. Later on I’ll explain the parameters

# DV needs to be factor variable so it knows to use a classification tree
tree<-tree(as.factor(respmail) ~ ., data=subset(ebeer, select = c(respmail, F, student)),mindev=.005)
par(mfrow=c(1,1))
plot(tree, col=8, lwd=2)
# cex controls the size of the type, 1 is the default.  
# label="yprob" gives the probability
text(tree, label = "yprob", cex=.75, font=2, digits = 2, pretty=0)

tree$frame

##        var    n  dev yval splits.cutleft splits.cutright yprob.0 yprob.1
## 1  student 4952 3720    0           <0.5            >0.5   0.876   0.124
## 2   <leaf> 2334    0    0                                  1.000   0.000
## 3        F 2618 2857    0           <2.5            >2.5   0.765   0.235
## 6        F  986  500    0           <1.5            >1.5   0.930   0.070
## 12  <leaf>  354    0    0                                  1.000   0.000
## 13  <leaf>  632  436    0                                  0.891   0.109
## 7        F 1632 2082    0           <4.5            >4.5   0.665   0.335
## 14  <leaf>  243  255    0                                  0.782   0.218
## 15  <leaf> 1389 1808    0                                  0.644   0.356

Here’s how to make predictions with trees, using the same data.

pred_tree <- predict(tree,ebeer, type = "vector")
head(pred_tree) %>% 
  kbl() %>%
  kable_styling()

0	1
0.644	0.356
1.000	0.000
1.000	0.000
0.891	0.109
0.782	0.218
0.891	0.109

Let’s look at the confusion matrix:

# take probability that responds
prob_resp <- pred_tree[,2]

# there is no predicted probability over 0.5
sum(prob_resp > 0.5)

## [1] 0

confusion_matrix <- (table(ebeer$respmail, prob_resp > 0.5))
confusion_matrix <- as.data.frame.matrix(confusion_matrix)
colnames(confusion_matrix) <- c("No")

confusion_matrix %>% 
  kbl() %>%
  kable_styling()

	No
0	4336
1	616

The misclassification rate is the proportion of responders, since all probabilities were under 0.5 so the model predicted no one responded. This is the same missclass. (This is a good example where using the default threshold of 0.5 would be not so good.)

confusion_matrix[1]/sum(confusion_matrix)

##      No
## 0 0.876
## 1 0.124

Residual mean deviance is the total residual deviance divided by the number of observations - number of terminal nodes, \(n-df\). The deviance is 2500.

Trees in R

We’ll use the tree package in R. Here’s how the tree fits the data:

# mean two graphs side-by-side
par(mfrow=c(1,2), oma=c(0,0,2,0))
# same model as above

tree <- tree(as.factor(respmail) ~ ., data=subset(ebeer, select = c(respmail, F, student)),mindev=.005)

plot(tree, col=8, lwd=2)

# cex controls the size of the type, 1 is the default.  
# label="yprob" gives the probability
text(tree, cex=.75, label="yprob", font=2, digits = 2, pretty = 0)



par(mai=c(.8,.8,.2,.2))

# create an aggregate table of response by frequency and student
tbl<- ebeer %>% group_by(student, F) %>% summarise(mean=mean(respmail)) %>% data.frame()

## `summarise()` has grouped output by 'student'. You can override using the `.groups` argument.

pred<-predict(tree,tbl, type = "vector")[,2]

tbl<-tbl %>% mutate(pred = pred)

# plot it
par(mai=c(.8,.8,.2,.2))
plot(tbl$F[1:12],tbl$mean[1:12], col = "red", xlab="Frequency", ylab="mean response",ylim=c(-.05,0.5), pch=20)
points(tbl$F[13:24],tbl$mean[13:24], col = "blue", pch=20)
legend(7.5, 0.5, legend=c("Student = no", "Student= yes"), col=c("red", "blue"), pch=20, cex=0.8)

# create predictions from tree for every F x student combo
newF <- seq(1,12,length=12)
lines(tbl$F[1:12], tbl$pred[1:12], col=2, lwd=2)
lines(tbl$F[1:12], tbl$pred[13:24], col=4, lwd=2)
mtext("A simple tree",outer=TRUE,cex=1.5)

Nonparametric

This method is nonparametric because it doesn’t make an assumption about the relationships between the independent and dependent variables.

By setting the mindev and mincut to zero, we can make the tree very complex and fit the data arbitrarily close.

# mean two graphs side-by-side
par(mfrow=c(1,2), oma = c(0, 0, 2, 0))
# same model as above
tree<-tree(as.factor(respmail) ~ ., data=subset(ebeer, select = c(respmail, F, student)),mindev=0, mincut=0)
plot(tree, col=8, lwd=2)
# cex controls the size of the type, 1 is the default.  
# label="yprob" gives the probability
text(tree, cex=.5, label="yprob", font=2, digits = 2, pretty = 0)

par(mai=c(.8,.8,.2,.2))

# create an aggregate table of response by frequency and student
tbl<- ebeer %>% group_by(student, F) %>% summarise(mean=mean(respmail)) %>% data.frame()

## `summarise()` has grouped output by 'student'. You can override using the `.groups` argument.

pred<-predict(tree,tbl, type = "vector")[,2]

tbl<-tbl %>% mutate(pred = pred)

# plot it
par(mai=c(.8,.8,.2,.2))
plot(tbl$F[1:12],tbl$mean[1:12], col = "red", xlab="Frequency", ylab="mean response",ylim=c(-.05,0.5), pch=20)
points(tbl$F[13:24],tbl$mean[13:24], col = "blue", pch=20)
legend(7.5, 0.5, legend=c("Student = no", "Student= yes"), col=c("red", "blue"), pch=20, cex=0.8)

# create predictions from tree for every F x student combo
newF <- seq(1,12,length=12)
lines(tbl$F[1:12], tbl$pred[1:12], col=2, lwd=2)
lines(tbl$F[1:12], tbl$pred[13:24], col=4, lwd=2)
mtext("A simple tree",outer=TRUE,cex=1.5)

Overfitting & K-fold cross validation

We begin by fitting a complicated tree by setting the mindev=0 and mincut to some low number or zero.

tree_complex<-tree(as.factor(respmail) ~ . , data=ebeer, mindev=0, mincut=100)

par(mfrow=c(1,1))
par(mai=c(.8,.8,.2,.2))
plot(tree_complex, col=10, lwd=2)
text(tree_complex, cex=.5, label="yprob", font=2, digits = 2, pretty = 0)
title(main="Classification Tree: complex")

I specify K=10 cross fold validation. The size is the resulting number of leaves. The out-of-sample error measure is the deviance. We want that as low as possible.

cv.tree_complex<-cv.tree(tree_complex, K=10)
cv.tree_complex$size

##  [1] 19 18 17 16 15 14 13 11 10  7  6  5  4  3  2  1

round(cv.tree_complex$dev)

##  [1] 2531 2531 2531 2531 2531 2531 2531 2531 2531 2531 2531 2531 2531 2589 2860 3724

par(mfrow=c(1,1))
plot(cv.tree_complex$size, cv.tree_complex$dev, xlab="tree size (complexity)", ylab="Out-of-sample deviance (error)", pch=20)

Choose the tree with the minimum out-of-sample error. Here the error remains the same after 4. So I choose the simplest model with the lowest OOS error, a tree with 4 leaves.

par(mfrow=c(1,2), oma = c(0, 0, 2, 0))
tree_cut<-prune.tree(tree_complex, best=4)
plot(tree_cut, col=10, lwd=2)
text(tree_cut, cex=1, label="yprob", font=2, digits = 2, pretty = 0)
title(main="A pruned tree")
summary(tree_cut)

## 
## Classification tree:
## snip.tree(tree = tree_complex, nodes = c(13L, 7L))
## Variables actually used in tree construction:
## [1] "student" "F"      
## Number of terminal nodes:  4 
## Residual mean deviance:  0.509 = 2520 / 4950 
## Misclassification error rate: 0.124 = 616 / 4952

pred_tree_ebeer<-predict(tree_cut, data=ebeer)[,2]

plot(roc(ebeer$respmail, pred_tree_ebeer), print.auc=TRUE,
     col="black", lwd=1, main="ROC curve", xlab="Specificity: true negative rate", ylab="Sensitivity: true positive rate", xlim=c(1,0))

Telco: In other data sets, the “optimal” size may be different. Here it is with telco.

# make a somewhat big tree
par(mfrow=c(1,2),oma = c(0, 0, 2, 0))
tree_telco<-tree(Churn ~ ., data=telco, mindev=0.005, mincut=0)
plot(tree_telco, col=10, lwd=2)
text(tree_telco, cex=.4, font=1, digits = 2, pretty = 0)

cv.tree_telco<-cv.tree(tree_telco, K=10)
cv.tree_telco

## $size
## [1] 9 8 7 6 5 4 3 2 1
## 
## $dev
## [1] 4363 4363 4363 4363 4476 4481 4579 4760 5698
## 
## $k
## [1]  -Inf  31.9  39.2  39.7  71.5  78.6 113.0 182.2 941.2
## 
## $method
## [1] "deviance"
## 
## attr(,"class")
## [1] "prune"         "tree.sequence"

plot(cv.tree_telco$size, cv.tree_telco$dev, xlab="tree size (complexity)", ylab="Out-of-sample deviance (error)", pch=20)
mtext("Another example: telco",outer=TRUE,cex=1.5)

The deviance doesn’t decrease any more once we get to 6 leaves. So I’ll set 6 as the best.

par(mfrow=c(1,1))
tree_cut<-prune.tree(tree_telco, best=6)
plot(tree_cut, col=10, lwd=2)
text(tree_cut, cex=1, font=1, digits = 2, pretty = 0, label="yprob")

With the telco data it looks like any number of leaves after 6 gives you the same OOS performance. So, going forward, the simplest model with the best OOS performance is 6.

Random Forests

We will fit random forests in R. We specify the number of trees and the minimum number of observations in a leaf (25), as well as the importance .

ebeer_rf<- ranger(respmail ~ ., data=ebeer, write.forest=TRUE, num.trees = 1000, min.node.size = 25, importance = "impurity", probability=TRUE, seed = 19103)

Variable importance:

par(mfrow=c(1,1))
par(mai=c(.9,.8,.2,.2))
sort(ebeer_rf$variable.importance, decreasing = TRUE)

##         F         M  firstpur   student         R age_class    single    gender   mailing 
##     124.0     117.1      93.7      83.1      50.0      31.9      12.5      10.3       0.0

barplot(sort(ebeer_rf$variable.importance, decreasing = TRUE), ylab = "variable importance")

Predictions on new data using out-of-the-bag (OOB) samples.

head(ebeer_rf$predictions) %>% 
  kbl() %>%
  kable_styling()

0	1
0.596	0.404
0.991	0.009
1.000	0.000
0.618	0.382
0.847	0.153
0.873	0.127

pred <- ebeer_rf$predictions[,2]

Confusion matrix

confusion_matrix <- (table(ebeer$respmail, pred > 0.5))
confusion_matrix <- as.data.frame.matrix(confusion_matrix)
colnames(confusion_matrix) <- c("No", "Yes")
confusion_matrix$Percentage_Correct <- confusion_matrix[1,]$No/(confusion_matrix[1,]$No+confusion_matrix[1,]$Yes)*100
confusion_matrix[2,]$Percentage_Correct <- confusion_matrix[2,]$Yes/(confusion_matrix[2,]$No+confusion_matrix[2,]$Yes)*100
print(confusion_matrix) %>% 
  kbl() %>%
  kable_styling()

##     No Yes Percentage_Correct
## 0 4290  46              98.94
## 1  586  30               4.87

	No	Yes	Percentage_Correct
0	4290	46	98.94
1	586	30	4.87

cat('Overall Percentage:', (confusion_matrix[1,1]+confusion_matrix[2,2])/nrow(ebeer)*100)

## Overall Percentage: 87.2

ROC

par(mfrow=c(1,1))
par(mai=c(.9,.8,.2,.2))
plot(roc(as.numeric(ebeer$respmail)-1, pred), print.auc=TRUE,
     col="black", lwd=1, main="ROC curve", xlab="Specificity: true negative rate", ylab="Sensitivity: true positive rate", xlim=c(1,0))
lines(roc(as.numeric(ebeer$respmail)-1, pred_tree_ebeer), print.auc=TRUE,  col="red", lwd=1)
legend('bottomright',legend=c("random forest", "decision tree"),col=c("black","red"), lwd=1)

Predictions on rollout data.

ebeer_rollout$p_rf <- predict(ebeer_rf, ebeer_rollout)$predictions[,2]