Supervised ML: Tree-based methods

Goals:

• Illustrate tree-based methods

• and the difference between global and local methods

Data generating process

Consider a model with only one covariate: \(Y = \underbrace{2X - X^3}_{CEF} + \epsilon\), where \(X \sim uniform(0,1)\) and \(\epsilon \sim N(0,1/3)\).

Plot a random draw with \(N=500\) and the true CEF.

# Load the packages required for later
if (!require("tidyverse")) install.packages("tidyverse", dependencies = TRUE); library(tidyverse)
if (!require("rpart")) install.packages("rpart", dependencies = TRUE); library(rpart)
if (!require("partykit")) install.packages("partykit", dependencies = TRUE); library(partykit)
if (!require("grf")) install.packages("grf", dependencies = TRUE); library(grf)
if (!require("rpart.plot")) install.packages("rpart.plot", dependencies = TRUE); library(rpart.plot)

set.seed(1234) # For replicability

cef = function(x){-x^3 + 2*x}

n = 500

x = runif(n)
y = cef(x) + rnorm(n,0,1/3)

df = data.frame(x=x,y=y)
ggplot(df) + stat_function(fun=cef,size=1) + 
            geom_point(aes(x=x,y=y),color="black",alpha = 0.4)

Different ways to fit the data

OLS

Plain OLS without any polynomials fits a straight line through the data cloud missing the non-linear shape but providing the best linear approximation:

df$y_hat_ols = predict(lm(y~x))
ggplot(df,aes(x=x,y=y)) + stat_function(fun=cef,size=1) + 
  geom_point(color="black",alpha = 0.4) + 
  geom_point(aes(x=x,y=y_hat_ols),shape="square",color="blue") + 
  geom_smooth(formula="y~x", method='lm') 

Regression tree

The regression tree provides a step function approximation of the non-linear shape:

tree = rpart(y~x,data = df)
rpart.plot(tree)

df$y_hat_tree = predict(tree)

ggplot(df) + stat_function(fun=cef,size=1) + 
  geom_point(aes(x=x,y=y),color="black",alpha = 0.4) +
  geom_point(aes(x=x,y=y_hat_tree),shape="square",color="blue") 

Notably the leaves are more narrow in the left part region where the curve is steeper, while more observation are pooled on the right where the functions is relatively flat. This makes intuitive sense.

Random Forest

Random Forest combines many of such stepwise approximation and provides a relatively good and smooth approximation of non-linear CEF:

rf = regression_forest(as.matrix(x),y,tune.parameters = "all",num.trees = 2000)
df$y_hat_rf = predict(rf)$predictions

ggplot(df) + stat_function(fun=cef,size=1) + 
  geom_point(aes(x=x,y=y),color="black",alpha = 0.4) +
  geom_point(aes(x=x,y=y_hat_rf),shape="square",color="blue")

Global vs. local predictors

To understand the difference between global and local methods, consider that we want to predict the outcome value at \(x_0 = 0.2\). Each of the three estimators uses just a weighted average of observed outcomes to form this prediction.

test_point = 0.2

OLS

Recall that in matrix notation the parameter estimates are \(\beta = (X'X)^{-1}X'Y\), where \(X\) contains a constant. Now the fitted value at the test point is \(\hat{y}_0 = [1~x_0] \beta = \underbrace{[1~x_0] (X'X)^{-1}X'}_{w_{ols}} Y\), where \(w_{ols}\) contains the weights that each outcome receives when forming the fitted value.

We graph the scatterplot of the observed data with larger points indicating a larger absolute outcome weight, and the color indicating whether the outcome weight is positive or negative:

X = cbind(rep(1,n),x)
predict_ols = as.numeric(c(1,test_point) %*% solve(t(X) %*% X) %*% t(X) %*% y)
df$w_ols = t(c(1,test_point) %*% solve(t(X) %*% X) %*% t(X))

w_sign_ols = rep("negative",n)
w_sign_ols[df$w_ols>0] = "positive"
w_sign_ols = factor(w_sign_ols,level=c("negative","positive"))

ggplot(df) + stat_function(fun=cef,size=1) + 
  geom_point(aes(x=x,y=y,size=abs(w_ols),color=w_sign_ols),alpha = 0.2) +
  geom_point(x=test_point,y=predict_ols,shape="cross",size = 3, stroke = 2,color="yellow") + 
  scale_color_manual(values=c("red","blue"))

A global method like OLS uses all available data points to form the prediction and outcomes that are far away of the point of interest receive negative weights.

Regression trees

The predictions of regression trees are formed as an average within the leaf of the point of interest. This means that all outcomes in that leaf receive the same weight and all other outcomes a weight of zero:

predict_tree = predict(tree,newdata=data.frame(x=test_point))

tree2 = as.party(tree)
nodes = predict(tree2,df, type = "node")
node_test = predict(tree2, newdata=data.frame(x=test_point), type = "node")

df$w_tree = (nodes==node_test)  / sum(nodes==node_test)

w_sign_tree = rep("zero",n)
w_sign_tree[df$w_tree>0] = "positive"
w_sign_tree = factor(w_sign_tree,level=c("zero","positive"))

ggplot(df) + stat_function(fun=cef,size=1) + 
  geom_point(aes(x=x,y=y,size=w_tree,color=w_sign_tree),alpha = 0.2) +
  geom_point(x=test_point,y=predict_tree,shape="cross",size = 3, stroke = 2,color="yellow") + 
  scale_color_manual(values=c("darkgray","blue"))

The problem is that predictions at the boundary discard outcomes that are very close.

Random Forest

The Random Forest creates a smoother weighting function and gives larger weights to closer outcomes and ignores outcomes that are further away.

predict_rf = predict(rf,newdata=as.matrix(test_point))$predictions
df$w_rf = t( as.matrix( get_forest_weights(rf,newdata=as.matrix(test_point)) ) )

w_sign_rf = rep("zero",n)
w_sign_rf[df$w_rf>0] = "positive"
w_sign_rf = factor(w_sign_rf,level=c("zero","positive"))

ggplot(df) + stat_function(fun=cef,size=1) + 
  geom_point(aes(x=x,y=y,size=w_rf,color=w_sign_rf),alpha = 0.2) +
  geom_point(x=test_point,y=predict_rf,shape="cross",size = 3, stroke = 2,color="yellow") + 
  scale_color_manual(values=c("darkgray","blue"))

Take-away:

• Tree based methods use only outcomes close to the prediction point, while global methods like OLS use a (at least for me) less intuitive weighting with negative weights

Suggestions to play with the toy model:

Feel free to play around with the code. This is useful to sharpen and challenge your understanding of the methods. Think about the consequences of a modifications before you run it and check whether the results are in line with your expectation. Some suggestions:

• Modify DGP (functional form of CEF, noise term, …)

• Decrease and increase the number of observations

• Decrease and increase the number of trees for the random forest

• Add squared or cubic covariates to OLS specification