Lecture 16
2025-04-03
Outline of This week
- Decision Trees
- Classification Trees
- Regression Trees
- Pruning
library(tidyverse)
library(cluster)
library(factoextra)
library(ggplot2)
library(pdxTrees)
library(gridExtra)
library(gglm)
library(rpart)
library(rpart.plot)
library(parttree)
library(yardstick)Trees
Why use trees
- Advantages
- Easy to interpret
- Even compared to OLS
- Similar to human decison making
- Easy to display graphically
- Easy to interpret even by non experts
- Don’t need dummy vairables
- Easy to interpret
- Disadvantages
- Less accurate
- non-robust
- small changes in input data can lead to very different trees
- Solutions to disadvantages: Next week
Regression Trees
Much of this is taken from Professor Wells
Create a Tree (set of rules) to partition data
Step 1: Divide the predictor space into J distinct and non overlapping regions
Step 2: For every observation in region \(R_j\) predict the mean (\(\bar{y}\)) for the training observations in region \(R_j\)
RBS
In theory, these regions could be any shape. In practice, we will create high dimensional boxes.
Goal: Find boxes to minimize RSS \(\sum_{j=1}^J\sum_{i\in R_j}(y_i-\hat{y_{R_j}})^2\)
Computationally infeasible to consider every possible partition, so instead we will use recursive binary splitting.
- Top down approach: begins at the top
- Binary: each node splits into two branches
- Greedy: at each step, make best local choice
Step 1: Find the predictor \(X_j\) and the cutpoint \(s\) such that splitting the predictor space into the two regions leads to the greates possible reduction in RSS. In other words given \[R_1(j,s)=\{X|X_j<s\}\text{ and } R_2(j,s)=\{X|X_j\geq s\}\] find j that minimizes \[\sum_{i:x_i\in R_1(j,s)}(y_i-\hat{y_{R_1}})^2+\sum_{i:x_i\in R_2(j,s)}(y_i-\hat{y_{R_2}})^2\]
Step 2 onwards: Repeat
Tree Pruning
Much of this is taken from Professor Wells
This should give good predictions on the training set but will likely overfit the data. Smaller trees are likely to have lower variance and easier interpretations, however they may have more bias.
Options:
- Stop tree when RSS reduction slows down
- Pruning
Cost Complexity Pruning
Also called weakest link pruning. We can consider a sequence of trees indexed by a tuning parameter \(\alpha\)
- For each \(\alpha\) there exists a
unique subtree T that minimizes RSS+\(\alpha\)|T| where |T| is the number of
terminal nodes
- In other words, trees are penalized based on number of terminal nodes
- Analagous to penalty parameter in Penalized Regression (e.g. Lasso)
- As \(\alpha\) increases from 0 (i.e. the full tree), branches get pruned in a predictable way
- To find the optimal value of \(\alpha\): cross validation, train/test, etc
- Choosing the best subtree
- Smallest rMSE
- Smallest tree with rMSE within 1 sd of lowest rMSE
Creating Trees by hand
Base Data
set.seed(2020)
X=rep(1:4,4)
Y=c(1,2,3,4,2,3,4,1,3,4,1,2,4,1,2,3)
Z=sample(1:5,16,TRUE)
toyData=data.frame(X,Y,Z)
ggplot(data=toyData,aes(x=X,y=Y,color=Z))+
geom_point()
Useful Functions
groupSSE<-function(dataFrame){
Zmean=mean(dataFrame$Z)
SSE=sum((dataFrame$Z-Zmean)^2)
return(SSE)
}
splitChoices<-function(dataFrame,otherSSE){
### Split X
for (i in c(1.5,2.5,3.5)){
group1<-dataFrame%>%filter(X>i)
group2<-dataFrame%>%filter(X<i)
print(groupSSE(group1)+groupSSE(group2)+otherSSE)
}
### Split Y
for (i in c(1.5,2.5,3.5)){
group1<-dataFrame%>%filter(Y>i)
group2<-dataFrame%>%filter(Y<i)
print(groupSSE(group1)+groupSSE(group2)+otherSSE)
}
}
splitChoicesGeneral<-function(DataFrames){
otherSSE=0
for(i in 1:length(DataFrames)){
otherSSE=otherSSE+groupSSE(DataFrames[[i]])
}
for(i in 1:length(DataFrames)){
print(paste("i =",i))
tempotherSSE=otherSSE-groupSSE(DataFrames[[i]])
splitChoices(DataFrames[[i]],tempotherSSE)
}
}Step 1
DataFrames<-list()
DataFrames[[1]]<-toyData
splitChoicesGeneral(DataFrames)## [1] "i = 1"
## [1] 32.41667
## [1] 32.875
## [1] 32.91667
## [1] 26.91667
## [1] 18.875
## [1] 25.41667First Split is Y=2.5
toyDataL<-toyData%>%filter(Y<2.5)
toyDataR<-toyData%>%filter(Y>2.5)
ggplot(data=toyData,aes(x=X,y=Y,color=Z))+
geom_point()+
geom_hline(yintercept=2.5)
Step 2
DataFrames<-list()
DataFrames[[1]]<-toyDataL
DataFrames[[2]]<-toyDataR
splitChoicesGeneral(DataFrames)## [1] "i = 1"
## [1] 18.83333
## [1] 17.75
## [1] 18.5
## [1] 18.75
## [1] 18.875
## [1] 18.875
## [1] "i = 2"
## [1] 18.20833
## [1] 18.375
## [1] 18.20833
## [1] 18.875
## [1] 18.875
## [1] 18.375- Split L by X=2.5 is the next best
- Question 1 Why are there two 18.875 SSE values in each of the possible splits?
toyDataLL<-toyDataL%>%filter(X<2.5)
toyDataLR<-toyDataL%>%filter(X>2.5)
ggplot(data=toyData,aes(x=X,y=Y,color=Z))+
geom_point()+
geom_hline(yintercept=2.5)+
geom_vline(xintercept=2.5)
- Question 2 Draw the tree and graph that this corresponds to, what are the roots/branches/leafs?
The below graph is purely for spacing, ignore it
toyDataLL<-toyDataL%>%filter(X<2.5)
toyDataLR<-toyDataL%>%filter(X>2.5)
ggplot(data=toyData,aes(x=X,y=Y,color=Z))+
geom_point()+
geom_hline(yintercept=2.5)+
geom_vline(xintercept=2.5)
This is an example of how to do it in R
Question 3 What does flip=TRUE do, and why do we need it?
DataFrames<-list()
DataFrames[[1]]<-toyDataLL
DataFrames[[2]]<-toyDataLR
DataFrames[[3]]<-toyDataR
mean(DataFrames[[1]]$Z)## [1] 3.5mean(DataFrames[[2]]$Z)## [1] 4.25mean(DataFrames[[3]]$Z)## [1] 2my_tree <- rpart(Z ~
X+Y,
control = rpart.control(
minsplit = 1, xval = 1, maxdepth = 2, cp = 0.03),
data = toyData)
rpart.plot(my_tree)
ggplot(data=toyData,aes(x=X,y=Y,color=Z))+
geom_point()+
geom_parttree(data =my_tree, aes(fill=Z),flip= TRUE, alpha = 0.2 )+scale_colour_viridis_c(aesthetics = c('fill'), option = "magma", name = "Pred", begin = 0, end = 1 )
- Question 4 Repeat the above process for seed=1 going
1 level deeper (not necessarily 1 layer deeper)
Creating Trees in R
Walking through an example using trees from Professor Wells
my_pdxTrees <- get_pdxTrees_parks(park = c("Kenilworth Park", "Westmoreland Park",
"Woodstock Park","Berkeley Park", "Powell Park"))
## Create average crown width (East West vs North South)
my_pdxTrees <- my_pdxTrees %>% mutate(Crown_Width = (Crown_Width_EW + Crown_Width_NS)/2)
dim(my_pdxTrees)## [1] 1039 35names(my_pdxTrees)## [1] "Longitude" "Latitude"
## [3] "UserID" "Genus"
## [5] "Family" "DBH"
## [7] "Inventory_Date" "Species"
## [9] "Common_Name" "Condition"
## [11] "Tree_Height" "Crown_Width_NS"
## [13] "Crown_Width_EW" "Crown_Base_Height"
## [15] "Collected_By" "Park"
## [17] "Scientific_Name" "Functional_Type"
## [19] "Mature_Size" "Native"
## [21] "Edible" "Nuisance"
## [23] "Structural_Value" "Carbon_Storage_lb"
## [25] "Carbon_Storage_value" "Carbon_Sequestration_lb"
## [27] "Carbon_Sequestration_value" "Stormwater_ft"
## [29] "Stormwater_value" "Pollution_Removal_value"
## [31] "Pollution_Removal_oz" "Total_Annual_Services"
## [33] "Origin" "Species_Factoid"
## [35] "Crown_Width"## Goal: Predict Carbon Sequestration
my_pdxTrees %>% select(Carbon_Sequestration_lb, Crown_Width, Tree_Height) %>% drop_na() %>% cor()## Carbon_Sequestration_lb Crown_Width Tree_Height
## Carbon_Sequestration_lb 1.0000000 0.6126951 0.4362020
## Crown_Width 0.6126951 1.0000000 0.5980118
## Tree_Height 0.4362020 0.5980118 1.0000000Predicting Carbon Sequestration
- Can we predict carbon sequestration based on other tree features?
g1 <- ggplot(my_pdxTrees, aes(x = Crown_Width, y = Carbon_Sequestration_lb ))+geom_point(alpha =.5, shape = 16)+theme_bw()
g2<-ggplot(my_pdxTrees, aes(x = Tree_Height, y = Carbon_Sequestration_lb ))+geom_point(alpha =.5, shape = 16)+theme_bw()
g3<-ggplot(my_pdxTrees, aes( x = Carbon_Sequestration_lb ))+geom_histogram(color = "white", bins = 20)+theme_bw()
g4<-ggplot(my_pdxTrees, aes( y = Tree_Height, x = Crown_Width, color = Carbon_Sequestration_lb ))+geom_point(alpha =.75, shape = 16)+theme_bw()+
scale_colour_viridis_c(begin = .15, end = .8, option = "magma", name = "Carbon")
grid.arrange(g1,g2,g3,g4,ncol =2)
Linear Regression
tree_lm<-lm(Carbon_Sequestration_lb ~Crown_Width + Tree_Height, data=my_pdxTrees)
summary(tree_lm)##
## Call:
## lm(formula = Carbon_Sequestration_lb ~ Crown_Width + Tree_Height,
## data = my_pdxTrees)
##
## Residuals:
## Min 1Q Median 3Q Max
## -87.395 -13.283 -4.912 10.982 121.950
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.08819 2.03721 -1.516 0.129853
## Crown_Width 0.88769 0.04947 17.944 < 2e-16 ***
## Tree_Height 0.10140 0.02848 3.560 0.000388 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.46 on 1031 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.383, Adjusted R-squared: 0.3818
## F-statistic: 320 on 2 and 1031 DF, p-value: < 2.2e-16## Diagnostic Plots
gglm(tree_lm)
Regression Tree
my_tree <- rpart(Carbon_Sequestration_lb ~ Tree_Height + Crown_Width, data = my_pdxTrees,
control = rpart.control(cp = .005))
rpart.plot(my_tree,extra = 0, box.palette = "-Pu")
g4+
geom_parttree(data =my_tree, aes(fill=Carbon_Sequestration_lb,flipaxes = TRUE), alpha = 0.2 )+scale_colour_viridis_c(aesthetics = c('fill'), option = "magma", name = "Pred", begin = 0, end = 1 )
Interpretation
Crown_Widthis the most important predictor ofCarbon_Sequestration_lbAfter accounting for width,
Tree_Heighthas some impact onCarbon_Sequestration_lbVery narrow and very wide trees tend to have low
Carbon_Sequestration_lbTrees of moderate width and height have largest
Carbon_Sequestration_lb
Tree Accuracy
- Let’s create a test set consisting of two other parks:
my_pdxTrees_test <- get_pdxTrees_parks(park = c("Mt Scott Park", "Glenwood Park"))my_pdxTrees_test<- my_pdxTrees_test %>% mutate(Crown_Width = (Crown_Width_EW + Crown_Width_NS)/2)- We’ll measure the accuracy of our model using root Mean Square Error:
\[\textrm{rMSE} = \sqrt{\frac{1}{n} \sum_{i = 1}^n (y_i - \hat{y}_i)^2}\]
tree_preds<-predict(my_tree, newdata = my_pdxTrees_test)Tree_rMSE <- sqrt(mean((tree_preds - my_pdxTrees_test$Carbon_Sequestration_lb)^2, na.rm = T))
data.frame(Tree_rMSE)## Tree_rMSE
## 1 15.37258- And compared to the linear model:
lm_preds<-predict(tree_lm, my_pdxTrees_test)
lm_rMSE<-sqrt(mean((lm_preds - my_pdxTrees_test$Carbon_Sequestration_lb)^2, na.rm = T))
data.frame(lm_rMSE)## lm_rMSE
## 1 16.87355Why did the tree model outperform the linear model?
- Nevertheless, what are some downsides to the tree model?
Pruning
plotcp(my_tree, upper = "size")
Comparison
results <- data.frame(model = "full", obs = my_pdxTrees_test$Carbon_Sequestration_lb, preds = predict(my_tree, my_pdxTrees_test))
pruned_tree <- prune(my_tree, cp = .01)
results <- results %>% rbind(data.frame(model = "pruned", obs = my_pdxTrees_test$Carbon_Sequestration_lb, preds = predict(pruned_tree, my_pdxTrees_test)))
very_pruned_tree <- prune(my_tree, cp = .021)
results <- results %>% rbind(data.frame(model = "very pruned", obs = my_pdxTrees_test$Carbon_Sequestration_lb, preds = predict(very_pruned_tree, my_pdxTrees_test)))
results <- results %>% rbind(data.frame(model = "linear", obs = my_pdxTrees_test$Carbon_Sequestration_lb, preds = lm_preds ))
results %>% group_by(model) %>% rmse(truth = obs, estimate = preds) %>% arrange(.estimate)## # A tibble: 4 × 4
## model .metric .estimator .estimate
## <chr> <chr> <chr> <dbl>
## 1 pruned rmse standard 14.3
## 2 full rmse standard 15.4
## 3 linear rmse standard 16.9
## 4 very pruned rmse standard 16.9Horizontal axis gives values of complexity parameter (cp)
Upper scale indicates number of terminal nodes for given tree
Vertical axis gives the cross-validated relative root mean squared error
Dotted horizontal line has height equal to 1 standard error above smallest rMSE
par(mfrow=c(1,3))
rpart.plot(my_tree, main = "Full Tree", extra = 0, box.palette = "-Pu", cex = 1)
rpart.plot(pruned_tree, main = "Pruned Tree", extra = 0, box.palette = "-Pu", cex = 1 )
rpart.plot(very_pruned_tree, main = "Very Pruned Tree", extra = 0, box.palette = "-Pu", cex = 1)
Trees using RPart
set.seed(1)
tree_model1 <- rpart(Carbon_Sequestration_lb ~
Tree_Height + Crown_Width,
data = my_pdxTrees)
set.seed(1)
tree_model2 <- rpart(Carbon_Sequestration_lb ~
Tree_Height + Crown_Width,
control = rpart.control(
minsplit = 20, xval = 10, maxdepth = 10, cp = 0.005),
data = my_pdxTrees)- minsplit: minimum observations per node
- xval: number of cross validation folds
- max depth: max depth of a node in the tree
- cp: min reduction in RSS to attempt a split
Plotting in R
plot(tree_model2)
text(tree_model2, pretty = 0, cex = .5)
rpart.plot(tree_model2)
plotcp(tree_model2)
- Based on CV plot 10 leaves, CP=0.0077 gives lowest error
- 7 leaves, CP=0.012 smallest tree within 1 SE of best
pruned_tree <- prune(tree_model2, cp = 0.0077)par(mfrow=c(1,2))
rpart.plot(tree_model2)
rpart.plot(pruned_tree)
Test Error Rates
results <- data.frame(model = "full",
obs = my_pdxTrees_test$Carbon_Sequestration_lb,
preds = predict(tree_model2, my_pdxTrees_test))
results <- rbind(results,
data.frame(model = "pruned",
obs = my_pdxTrees_test$Carbon_Sequestration_lb,
preds = predict(pruned_tree, my_pdxTrees_test)))- And use
rmsefromyardstickto assess:
results %>% group_by(model) %>%
rmse(truth = obs, estimate = preds) %>% arrange(.estimate)## # A tibble: 2 × 4
## model .metric .estimator .estimate
## <chr> <chr> <chr> <dbl>
## 1 pruned rmse standard 14.2
## 2 full rmse standard 15.4Text Problems
- Question 5 8.4.1
- Question 6 Work through 8.3.2
Wrap-up
Lab goals:
In this week we covered Tree based models.
Course Schedule:
- This week: Trees
- Next week: Ensemble Models
Reminders for next class:
- Labs due Friday
- Please do readings
- Project 2 has been released, it is due next week
- Final Project proposal is due this Friday
- You have a homework on Clustering due this Friday