1.6.1 Example: partially linear model

\[ y_i = \theta_0 d_i + f_0(x_i) + \epsilon_i \]

• Compare the estimates from

  1. \(\En[d_i(y_i - \tilde{\theta} d_i - \hat{f}(x_i)) ] = 0\)

  and

  2. \(\En[(d_i - \hat{m}(x_i))(y_i - \hat{\mu}(x_i) - \theta (d_i - \hat{m}(x_i)))] = 0\)

  where \(m(x) = \Er[d|x]\) and \(\mu(y) = \Er[y|x]\)

1.6.2 Lessons from the example

• Need an extra condition on moments – Neyman orthogonality \[ \partial \eta \Er[\psi(W;\theta_0,\eta_0)](\eta-\eta_0) = 0 \]

• Want estimators faster than \(n^{-1/4}\) in the prediction norm, \[ \sqrt{\En[(\hat{\eta}(x_i) - \eta(x_i))^2]} \lesssim_P n^{-1/4} \]

• Also want estimators that satisfy something like \[ \sqrt{n} \En[(\eta(x_i)-\hat{\eta}(x_i))\epsilon_i] = o_p(1) \]
  • Sample splitting will make this easier

2.1.1 Predicting house prices

• Example from Mullainathan and Spiess (2017)
• Training on 10000 observations from AHS
• Predict log house price using 150 variables
• Holdout sample of 41808

2.1.2 Predicting pipeline revenues

• Data on US natural gas pipelines
  • Combination of FERC Form 2, EIA Form 176, and other sources, compiled by me
  • 1996-2016, 236 pipeline companies, 1219 company-year observations
• Predict: \(y =\) profits from transmission of natural gas
• Covariates: year, capital, discovered gas reserves, well head gas price, city gate gas price, heating degree days, state(s) that each pipeline operates in

load("~/natural-gas-pipelines/dataAndCode/pipelines.Rdata")
# data has problems before 1996 due to format change
data <- subset(data,report_yr>=1996)
# replace NA state weights with 0's
data[,59:107][is.na(data[,59:107])] <- 0
# spaces in variable names will create problems later
names(data) <- gsub(" ",".",names(data))
summary(data[,c("transProfit","transPlant_bal_beg_yr","cityPrice","wellPrice")])

##   transProfit         transPlant_bal_beg_yr   cityPrice      
##  Min.   : -31622547   Min.   :0.000e+00     Min.   : 0.4068  
##  1st Qu.:   2586031   1st Qu.:2.404e+07     1st Qu.: 3.8666  
##  Median :  23733170   Median :1.957e+08     Median : 5.1297  
##  Mean   :  93517513   Mean   :7.772e+08     Mean   : 5.3469  
##  3rd Qu.: 129629013   3rd Qu.:1.016e+09     3rd Qu.: 6.5600  
##  Max.   :1165050214   Max.   :1.439e+10     Max.   :12.4646  
##  NA's   :2817         NA's   :2692          NA's   :1340     
##    wellPrice     
##  Min.   :0.0008  
##  1st Qu.:2.1230  
##  Median :3.4370  
##  Mean   :3.7856  
##  3rd Qu.:5.1795  
##  Max.   :9.6500  
##  NA's   :2637

library(GGally)
ggpairs(data[,c("transProfit","transPlant_bal_beg_yr","cityPrice","wellPrice")],
        lower=list(continuous="smooth"))  + theme_minimal()

2.1.3 Predicting pipeline revenues : methods

• OLS : 67 covariates (year dummies and state(s) create a lot)
• Lasso
• Random forests
• Randomly choose 75% of sample to fit the models, then look at prediction accuracy in remaining 25%

## Create X matrix for OLS and random forests
xnames <-c("transPlant_bal_beg_yr",   "reserve",   "wellPrice",  "cityPrice",
           "plantArea",  "heatDegDays",
           names(data)[59:107] )
yname <- "transProfit"
fmla <- paste(yname,"~",paste(xnames,collapse=" + "),"+ as.factor(report_yr)")
ols <- lm(fmla,data=data,x=TRUE,y=TRUE)
X <- ols$x[,!(colnames(ols$x) %in% c("(Intercept)")) &
            !is.na(ols$coefficients)]
y <- ols$y
train <- runif(nrow(X))<0.75

# OLS prediction on training set
y.t <- y[train]
X.t <- X[train,]
ols <- lm(y.t ~ X.t)
y.hat.ols <- ols$coefficients[1] + X %*% ols$coefficients[2:(length(ols$coef))]
df <- data.frame(y=y, y.hat=y.hat.ols, train=train, method="ols")

## Lasso
library(glmnet)
# Create larger X matrix for lasso
fmla.l <- paste(yname,"~ (",
                paste(xnames,collapse=" + "),")*(report_yr + transPlant_bal_beg_yr +
    reserve + wellPrice + cityPrice + plantArea + heatDegDays) + ",
                paste(sprintf("I(%s^2)",xnames[1:6],collapse=" + "))
                )
reg <- lm(fmla.l, data=data, x=TRUE,y=TRUE)
Xl <- reg$x[,!(colnames(reg$x) %in% c("(Intercept)")) &
             !is.na(reg$coefficients)]
lasso <- cv.glmnet(Xl[train,],y[train],alpha=1,parallel=FALSE,
                   standardize=TRUE, intercept=TRUE, nfolds = 50)
y.hat.lasso <- predict(lasso, Xl, s=lasso$lambda.min, type="response")
df <- rbind(df,  data.frame(y=y, y.hat=as.vector(y.hat.lasso),
                            train=train, method="lasso"))

## Random forest
library(grf)
rf <- regression_forest(X[train,],y[train],tune.parameters = TRUE)
y.hat.rf  <-  predict(rf, X)$predictions
df <- rbind(df, data.frame(y=y, y.hat=y.hat.rf, train=train,
                           method="random forest"))

# Neural network
library(RSNNS)
n <- nrow(X[train,])
p <- ncol(X)
rn <- floor(n^(1/(2*(1+1/(1+p))))/2)
xn <- normalizeData(X)
yn <- normalizeData(y)
nn <- mlp(x=xn[train,], y=yn[train], linOut=TRUE, size=rn)
yn.hat.nn <- predict(nn, xn)
y.hat.nn <- denormalizeData(yn.hat.nn, getNormParameters(yn))
df <- rbind(df, data.frame(y=y, y.hat=y.hat.nn, train=train,
                           method="neural network"))

ggplot(data=df,aes(x=y,y=y.hat,colour=method,shape=train)) +  geom_point(alpha=0.5) +
  geom_line(aes(y=y)) + theme_minimal()

df$trainf <- factor(df$train, levels=c("TRUE", "FALSE"))
df$error <- df$y - df$y.hat
ggplot(data=df,aes(x=error,colour=method)) +
  geom_density() + theme_minimal() +
  xlim(quantile(df$error,c(0.01,0.99))) +
  facet_grid(trainf ~ .,labeller=label_both)

library(kableExtra)
fn <- function(df)  with(df,c(mean((y.hat - y)^2)/var(y),
                              mean(abs(y.hat - y))/mean(abs(y-mean(y)))))
tab1 <-unlist(by(subset(df,train), df$method[train],  FUN=fn))
tab1 <- (matrix(tab1,nrow=2))
rownames(tab1) <- c("relative MSE","relative MAE")
colnames(tab1) <- c("OLS","Lasso","Random forest","Neural Network")
tab2 <- unlist(by(subset(df,!train), df$method[!train],  FUN=fn))
tab2 <- (matrix(tab2,nrow=2))
rownames(tab2) <- c("relative MSE","relative MAE")
colnames(tab2) <- c("OLS","Lasso","Random forest","Neural Network")

kable_styling(kable(tab1, caption="Training sample", format="html", digits=3),
              bootstrap_options = c("striped", "hover", "condensed",
              "responsive"), full_width=F)

kable_styling(kable(tab2, caption="Hold-out sample", format="html", digits=3),
              bootstrap_options = c("striped", "hover", "condensed",
              "responsive"), full_width=F)

2.2.1 Statistical properties of Lasso

• Model : \[ y_i = x_i'\beta_0 + \epsilon_i \]
  • \(\Er[x_i \epsilon_i] = 0\)
  • \(\beta_0 \in \R^n\)
  • \(p\), \(\beta_0\), \(x_i\), and \(s\) implicitly depend on \(n\)
  • \(\log p = o(n^{1/3})\)
    • \(p\) may increase with \(n\) and can have \(p>n\)
• Sparsity \(s\)
  • Exact : \(\norm{\beta_0}_0 = s = o(n)\)
  • Approximate : \(|\beta_{0,j}| < Aj^{-a}\), \(a > 1/2\), \(s \propto n^{1/(2a)}\)

2.2.2 Rate of convergence

• With \(\lambda = 2c \sqrt{n} \Phi^{-1}(1-\gamma/(2p))\) \[ \sqrt{\En[(x_i'(\hat{\beta}^{lasso} - \beta_0))^2 ] } \lesssim_P \sqrt{ (s/n) \log (p) }, \]

\[ \norm{\hat{\beta}^{lasso} - \beta_0}_2 \lesssim_P \sqrt{ (s/n) \log (p) }, \]

and

\[ \norm{\hat{\beta}^{lasso} - \beta_0}_1 \lesssim_P \sqrt{ (s^2/n) \log (p) } \]

• Constant \(c>1\)

  • Small \(\gamma \to 0\) with \(n\), and \(\log(1/\gamma) \lesssim \log(p)\)

  • Rank like condition on \(x_i\)

• near-oracle rate

2.2.3 Other statistical properties

• Inference on \(\beta\): not the goal in our motivating examples
  • Difficult, but some recent results
  • See Lee et al. (2016), Taylor and Tibshirani (2017), Caner and Kock (2018)
• Model selection: not the goal in our motivating examples
  • Under stronger conditions, Lasso correctly selects the nonzero components of \(\beta_0\)
  • See Belloni and Chernozhukov (2011)

2.2.4 Post-Lasso

• Two steps :

  1. Estimate \(\hat{\beta}^{lasso}\)

  2. \({\hat{\beta}}^{post} =\) OLS regression of \(y\) on components of \(x\) with nonzero \(\hat{\beta}^{lasso}\)

• Same rates of convergence as Lasso
• Under some conditions post-Lasso has lower bias
  • If Lasso selects correct model, post-Lasso converges at the oracle rate

2.3.1 Regression trees

• \(y_i \in R\) on \(x_i \in \R^p\)
• Want to estimate \(\Er[y | x]\)
• Locally constant estimate \[ \hat{t}(x) = \sum_m^M c_m 1\{x \in R_m \} \]
• Rectangular regions \(R_m\) determined by tree

2.3.2 Simulated data

n <- 1000
x <- runif(n)
y <- runif(n)
f <- function(x,z) {
  1/3*(sin(5*x)*sqrt(z)*exp(-(z-0.5)^2))
}
f0 <- f(x,y)
z <- f0 + rnorm(n)*0.1
tree.df <- data.frame(x=x,y=y,z=z)
# plot true function and data
x.g <- seq(0,1,length.out=100)
y.g <- seq(0,1,length.out=100)
f0.g <- t(outer(x.g,y.g,f))

library(plotly)
fig <- plot_ly( colors="YlOrBr")
fig <- add_markers(fig,x=x,y=y,z=z, size=0.3, opacity=0.2)
fig <- add_surface(fig, x=x.g, y=y.g, z=f0.g, opacity=1)
fig

2.3.3 Estimated tree

# fit regression tree
tree <- ctree(z ~ x + y, data=tree.df)

# plot estimate
x.g <- seq(0,1,length.out=100)
y.g <- seq(0,1,length.out=100)
df <- expand.grid(x.g,y.g)
names(df) <- c("x","y")
fhat.g <- matrix(predict(tree, newdata=df),nrow=length(x.g), byrow=TRUE)
library(plotly)
fig <- plot_ly(colors="YlOrBr")
fig <- add_markers(fig,x=x,y=y,z=z, size=0.3, opacity=0.2)
fig <- add_surface(fig, x=x.g, y=y.g, z=fhat.g, opacity=1)
fig

2.3.4 Estimated tree

plot(tree)

2.3.5 Tree algorithm

• For each region, solve \[ \min_{j,s} \left[ \min_{c_1} \sum_{i: x_{i,j} \leq s, x_i \in R} (y_i - c_1)^2 + \min_{c_2} \sum_{i: x_{i,j} > s, x_i \in R} (y_i - c_2)^2 \right] \]
• Repeat with \(R = \{x:x_{i,j} \leq s^*\} \cap R\) and \(R = \{x:x_{i,j} \leq s^*\} \cap R\)
• Stop when \(|R| =\) some chosen minimum size
• Prune tree \[ \min_{tree \subset T} \sum (\hat{f}(x)-y)^2 + \alpha|\text{terminal nodes in tree}| \]

2.3.6 Random forests

• Average randomized regression trees
• Trees randomized by
  • Bootstrap or subsampling
  • Randomize branches: \[ \min_{j \in S,s} \left[ \min_{c_1} \sum_{i: x_{i,j} \leq s, x_i \in R} (y_i - c_1)^2 + \min_{c_2} \sum_{i: x_{i,j} > s, x_i \in R} (y_i - c_2)^2 \right] \] where \(S\) is random subset of \(\{1, ..., p\}\)
• Variance reduction

2.3.7 Rate of convergence: regression tree

• \(x \in [0,1]^p\), \(\Er[y|x]\) Lipschitz in \(x\)
• Crude calculation for single tree, let denote \(R_i\) node that contains \(x_i\) \[ \begin{align*} \Er(\hat{t}(x_i) - \Er[y|x_i])^2 = & \overbrace{\Er(\hat{t}(x_i) - \Er[y|x\in R_i])^2}^{variance} + \overbrace{(\Er[y|x \in R_i] - \Er[y|x])^2}^{bias^2} \\ = & O_p(1/m) + O\left(L^2 \left(\frac{m}{n}\right)^{2/p}\right) \end{align*} \] optimal \(m = O(n^{2/(2+p)})\) gives \[ \Er[(\hat{t}(x_i) - \Er[y|x_i])^2] = O_p(n^{\frac{-2}{2+p}}) \]

2.3.8 Rate of convergence: random forest

• Result from Biau (2012)
• Assume \(\Er[y|x]=\Er[y|x_{(s)}]\), \(x_{(s)}\) subset of \(s\) variables, then \[ \Er[(\hat{r}(x_i) - \Er[y|x_i])^2] = O_p\left(\frac{1}{m\log(n/m)^{s/2p}}\right) + O_p\left(\left(\frac{m}{n}\right)^{\frac{0.75}{s\log 2}} \right) \] or with optimal \(m\) \[ \Er[(\hat{t}(x_i) - \Er[y|x_i])^2] = O_p(n^{\frac{-0.75}{s\log 2+0.75}}) \]

2.3.9 Other statistical properties

• Pointwise asymptotic normality : Wager and Athey (2018)

2.3.10 Simulation study

• Partially linear model
• DGP :
  • \(x_i \in \R^p\) with \(x_{ij} \sim U(0,1)\)
  • \(d_i = m(x_i) + v_i\)
  • \(y_i = d_i\theta + f(x_i) + \epsilon_i\)
  • \(m()\), \(f()\) either linear or step functions
• Estimate by OLS, Lasso, and random forest
  • Lasso & random forest use orthogonal moments \[ \En[(d_i - \hat{m}(x_i))(y_i - \hat{\mu}(x_i) - \theta (d_i - \hat{m}(x_i)))] = 0 \]

rm(list=ls())
p <- 4 # number of x's
mu.linear <- function(x) x%*%rep(1,p)
m.linear <- function(x) x%*%rep(2,p)
mu.step <- function(x) (x>0.5)%*%rep(1,p)
m.step <- function(x) (x>0.5)%*%rep(2,p)
theta <- 1

simulate <- function(n,p,mu,m) {
  theta <- 1
  x <- matrix(runif(n*p), ncol=p)
  d <- m(x) + rnorm(n)
  y <- theta*d + mu(x) + rnorm(n)
  data.frame(y=y,d=d,x=x)
}

library(grf)
library(hdm)
df <- simulate(100,p,mu.linear,m.linear)
mrfparams <- NULL
murfparams <- NULL
n.save <- NULL
partial.linear.rf <- function(df) {
  x.names <- names(df)[grep("x.",names(df))]
  if (is.null(mrfparams) || n.save!=nrow(df)) {
    # to save time, we only tune once per cluster worker and data set
    # size
    cat("tuning")
    m.rf  <- regression_forest(df[,x.names], df$d, num.trees=1000,
                               tune.parameters=TRUE)
    mrfparams <<- m.rf$tunable.params
    mu.rf  <- regression_forest(df[,x.names], df$y, num.trees=1000,
                                tune.parameters=TRUE)
    n.save <<- nrow(df)
    murfparams <<- mu.rf$tunable.params
  } else {
    cat("not tuning")
    m.rf  <- regression_forest(df[,x.names], df$d, num.trees=200,
                               tune.parameters=FALSE,
                               min.node.size =
                                 as.numeric(mrfparams["min.node.size"]),
                               alpha = as.numeric(mrfparams["alpha"]),
                               imbalance.penalty=as.numeric(mrfparams["imbalance.penalty"]),
                               sample.fraction = as.numeric(mrfparams["sample.fraction"]),
                               mtry=as.numeric(mrfparams["mtry"]))
    mu.rf  <- regression_forest(df[,x.names], df$y, num.trees=200,
                                tune.parameters=FALSE,
                               min.node.size =
                                 as.numeric(murfparams["min.node.size"]),
                               alpha = as.numeric(murfparams["alpha"]),
                               imbalance.penalty=as.numeric(murfparams["imbalance.penalty"]),
                               sample.fraction = as.numeric(murfparams["sample.fraction"]),
                               mtry=as.numeric(murfparams["mtry"]))

  }
  vhat <- df$d - predict(m.rf)$predictions
  ehat <- df$y - predict(mu.rf)$predictions
  lm(ehat ~ vhat)
}

## Manual sample splitting --- this turns out to be unneccessary. The
## default behavior of predict.regression_forest is to return
## predictions on the training data using only trees that were not fit
## on each observation. In other words, regression_forest already does
## the sample splitting for us.
##
## rf.tuneOnce <- function(x.names, y.name) {
##   parms <- NULL
##   function(df) {
##     if (is.null(parms)) {
##       rf  <- regression_forest(df[,x.names], df[,y.name], num.trees=500,
##                                  tune.parameters=TRUE)
##       parms <<- rf$tunable.params
##       rf
##     } else {
##       rf <- regression_forest(df[,x.names], df[,y.name], num.trees=200,
##                               tune.parameters=FALSE,
##                               honesty=FALSE,
##                               min.node.size =
##                                 as.numeric(parms["min.node.size"]),
##                               alpha = as.numeric(parms["alpha"]),
##                               imbalance.penalty=as.numeric(parms["imbalance.penalty"]),
##                               sample.fraction = as.numeric(parms["sample.fraction"]),
##                               mtry=as.numeric(parms["mtry"]))
##     }
##   }
## }
## n.save.split <- NULL
## m.hat.rf <- NULL
## mu.hat.rf  <- NULL
## partial.linear.split.rf <- function(df , splits=3) {
##   x.names <- names(df)[grep("x.",names(df))]
##   if (is.null(n.save.split) || n.save.split != nrow(df)) {
##     n.save.split <<- nrow(df)
##     m.hat.rf <<- rf.tuneOnce(x.names,"d")
##     mu.hat.rf <<- rf.tuneOnce(x.names,"y")
##   }
##   df$group <- sample(1:splits, nrow(df), replace=TRUE)
##   vhat <- df$d
##   ehat <- df$y
##   for(g in 1:splits) {
##     sdf <- subset(df, group!=g)
##     m <- m.hat.rf(sdf)
##     mu <- mu.hat.rf(sdf)
##     vhat[df$group==g] <- df$d[df$group==g] -
##       predict(m, newx=df[df$group==g,x.names])$predictions
##     ehat[df$group==g] <- df$y[df$group==g] -
##       predict(mu, newx=df[df$group==g,x.names])$predictions
##   }
##   lm(ehat ~ vhat)
## }


partial.linear.lasso <- function(df) {
  x.names <- names(df)[grep("x.",names(df))]
  fmla <- as.formula(paste(c("y ~ d",x.names), collapse=" + "))
  rlassoEffects(fmla, data=df, I = ~ d)
}

#summary(partial.linear.lasso(df))

# simulate a bunch of times in parallel
simulations <- 500 # number of simulations
library(parallel)
cl <- makeCluster(detectCores()/2)  # change as you see fit
clusterEvalQ(cl,library(hdm))
clusterEvalQ(cl,library(grf))

# R Socket cluster spawns new R sessions with empty environments, we
# need to make sure they load any needed libraries and have access to
# things from the main environment that they use
design <- c("linear") #,"step")
sim.df <- data.frame()
start.time <- Sys.time()
for (d in design) {
  if (d=="linear") {
    m <- m.linear
    mu <- mu.linear
  } else {
    m <- m.step
    mu <- mu.step
  }
  for (p in c(2,4,6,8)) {
    for (n in c(100, 200, 400, 800, 1600)) {
      clusterExport(cl,c("simulate","partial.linear.lasso",
                         "partial.linear.rf","p","mu","m",
                         "mrfparams","murfparams", "n.save"))
#                         "partial.linear.split.rf", "n.save.split",
#                         "m.hat.rf","mu.hat.rf","rf.tuneOnce"))

      thetas <- parSapply(cl, rep(n,simulations), function(n)
      {
        df <- simulate(n, p, mu, m)
        x.names <- names(df)[grep("x.",names(df))]
        fmla <- as.formula(paste(c("y ~ d",x.names), collapse=" + "))
        c(lm(fmla,data=df)$coef[2],
          partial.linear.rf(df)$coef[2],
          #partial.linear.split.rf(df)$coef[2],
          partial.linear.lasso(df)$coefficients)
      }
      )
      tmp <- (data.frame(t(thetas)) - 1)*sqrt(n)
      names(tmp) <- c("OLS","Random.Forest","Lasso")
      tmp$n <- n
      tmp$p <- p
      tmp$design <- d
      sim.df <- rbind(sim.df, tmp)
      cat("finished sample size ",n,"\n")
      cat("Elapsed time ", Sys.time()-start.time,"\n")
    }
    cat("finished p = ",p,"\n")
  }
  cat("finished design = ", d,"\n")
}
stopCluster(cl)
save(sim.df, file="partialLinearSim.Rdata")

library(ggplot2)
library(reshape2)
library(latex2exp)
TeX <- latex2exp::TeX
load("partialLinearSim.Rdata") # see partialLinearSim.R for simulation
                               # code
df <- melt(sim.df, measure.vars=c("OLS","Random.Forest","Lasso"))
ggplot(subset(df,p==2), aes(x=value, colour=variable)) +
  facet_grid(n ~ design) + geom_density() + theme_minimal() +
  xlab(TeX('$\\sqrt{n}(\\hat{\\theta}-\\theta_0)$')) +
  ggtitle("p=2")

ggplot(subset(df,p==4), aes(x=value, colour=variable)) +
  facet_grid(n ~ design) + geom_density() + theme_minimal() +
  xlab(unname(TeX("$\\sqrt{n}(\\hat{\\theta} - \\theta_0)$"))) +
  ggtitle("p=4")

ggplot(subset(df,p==6), aes(x=value, colour=variable)) +
  facet_grid(n ~ design) + geom_density() + theme_minimal() +
  xlab(unname(TeX("$\\sqrt{n}(\\hat{\\theta} - \\theta_0)$"))) +
  ggtitle("p=6")

ggplot(subset(df,p==8), aes(x=value, colour=variable)) +
  facet_grid(n ~ design) + geom_density() + theme_minimal() +
  xlab(unname(TeX("$\\sqrt{n}(\\hat{\\theta} - \\theta_0)$"))) +
  ggtitle("p=8")

2.4.1 Deep Neural Networks

• Many hidden layers
  • \(x^{(0)} = x\)
  • \(x^{(\ell)}_j = \psi(a_j^{(\ell)} x^{(\ell-1)} + b_j^{(\ell)})\)

nn <- mlp(x=xn, y=yn, linOut=TRUE, size=c(5,10,3,5, 6))
plot.nnet(nn, x.lab=colnames(mod$x)[2:ncol(mod$x)], y.lab="transProfit")

2.4.2 Rate of convergence

• Chen and White (1999)
• \(f(x) = \Er[y|x]\) with Fourier representation \[ f(x) = \int e^{i a'x} d\sigma_f(a) \] where \(\int (\sqrt{a'a} \vee 1) d|\sigma_f|(a) < \infty\)
• Network sieve : \[ \begin{align*} \mathcal{G}_n = \{ & g: g(x) = \sum_{j=1}^{r_n} \beta_j (a_j'a_j \vee 1)^{-1} \psi(a_j'x + b_j), \\ & \norm{\beta}_1 \leq B_n \} \end{align*} \]

2.4.3 Simulation Study

• Same setup as for random forests earlier
• Partially linear model
• DGP :
  • \(x_i \in \R^p\) with \(x_{ij} \sim U(0,1)\)
  • \(d_i = m(x_i) + v_i\)
  • \(y_i = d_i\theta + f(x_i) + \epsilon_i\)
  • \(m()\), \(f()\) either linear or step functions
• Estimate by OLS, Neural network with & without cross-fitting
  • Using orthogonal moments \[ \En[(d_i - \hat{m}(x_i))(y_i - \hat{\mu}(x_i) - \theta (d_i - \hat{m}(x_i)))] = 0 \]

rm(list=ls())
p <- 4 # number of x's
mu.linear <- function(x) x%*%rep(1,p)
m.linear <- function(x) x%*%rep(2,p)
mu.step <- function(x) (x>0.5)%*%rep(1,p)
m.step <- function(x) (x>0.5)%*%rep(2,p)
theta <- 1

simulate <- function(n,p,mu,m) {
  theta <- 1
  x <- matrix(runif(n*p), ncol=p)
  d <- m(x) + rnorm(n)
  y <- theta*d + mu(x) + rnorm(n)
  data.frame(y=y,d=d,x=x)
}

# There are many R packages for neural networks. I first tried neuralnet,
# but found it slow, then I tried RSNNS, and it seems fast
# enough. There might be a better option though
library(RSNNS)
df <- simulate(1000,p,mu.linear,m.linear)


partial.linear.net <- function(df, crossfits=5) {
  x.names <- names(df)[grep("x.",names(df))]
  fmla <- as.formula(paste(c("y ~ 1",x.names), collapse=" + "))
  n <- nrow(df)
  p <- ncol(df)
  n <- round(n*(1-1/pmax(1,crossfits)))
  rn <- floor(n^(1/(2*(1+1/(1+p)))))
  n <- nrow(df)
  ehat <- rep(NA, n)
  vhat <- rep(NA, n)
  if (crossfits>1) {
    ki <- sample(1:crossfits, n, replace=TRUE)
    for (k in 1:crossfits) {
      ynn <- mlp(df[ki!=k,x.names], df$y[ki!=k], size=c(rn),
                 linOut=TRUE, learnFunc="SCG")
      yhat <- predict(ynn, df[ki==k, x.names])
      ehat[ki==k] <- df$y[ki==k]-yhat
      dnn <-  mlp(df[ki!=k,x.names], df$d[ki!=k], size=c(rn),
                  linOut=TRUE, learnFunc="SCG")
      dhat <- predict(dnn, df[ki==k,x.names])
      vhat[ki==k] <- df$d[ki==k]-dhat
    }
  } else {
    ynn <- mlp(df[,x.names], df$y, size=c(rn), linOut=TRUE,
               learnFunc="SCG")
    yhat <- predict(ynn, )
    ehat <- df$y-yhat
    dnn <-  mlp(df[,x.names], df$d, size=c(rn), linOut=TRUE,
                learnFunc="SCG")
    dhat <- predict(dnn)
    vhat <- df$d-dhat
  }
  lm(ehat ~ vhat)
}

partial.linear.lasso <- function(df) {
  x.names <- names(df)[grep("x.",names(df))]
  fmla <- as.formula(paste(c("y ~ d",x.names), collapse=" + "))
  rlassoEffects(fmla, data=df, I = ~ d)
}

#summary(partial.linear.lasso(df))

# simulate a bunch of times in parallel
simulations <- 500 # number of simulations
library(parallel)
cl <- makeCluster(pmax(detectCores()-2,1))  # change as you see fit
clusterEvalQ(cl,library(hdm))
clusterEvalQ(cl,library(RSNNS))

# R Socket cluster spawns new R sessions with empty environments, we
# need to make sure they load any needed libraries and have access to
# things from the main environment that they use
design <- c("linear","step")
sim.df <- data.frame()
start.time <- Sys.time()
for (d in design) {
  if (d=="linear") {
    m <- m.linear
    mu <- mu.linear
  } else {
    m <- m.step
    mu <- mu.step
  }
  for (p in c(2,4,6,8)) {
    for (n in c(100, 200, 400, 800, 1600)) {
      clusterExport(cl,c("simulate",
                         "partial.linear.net","p","mu","m"))
#                         "partial.linear.split.rf", "n.save.split",
#                         "m.hat.rf","mu.hat.rf","rf.tuneOnce"))

      thetas <- parSapply(cl, rep(n,simulations), function(n)
      {
        df <- simulate(n, p, mu, m)
        x.names <- names(df)[grep("x.",names(df))]
        fmla <- as.formula(paste(c("y ~ d",x.names), collapse=" + "))
        c(lm(fmla,data=df)$coef[2],
          partial.linear.net(df,crossfits=5)$coef[2],
          partial.linear.net(df,crossfits=0)$coef[2])
          #partial.linear.split.rf(df)$coef[2],
          #partial.linear.lasso(df)$coefficients)
      }
      )
      tmp <- (data.frame(t(thetas)) - 1)*sqrt(n)
      names(tmp) <- c("OLS","Neural Network",
                      "Neural Network (no cross fitting)")
      tmp$n <- n
      tmp$p <- p
      tmp$design <- d
      sim.df <- rbind(sim.df, tmp)
      cat("finished sample size ",n,"\n")
      cat("Elapsed time ", Sys.time()-start.time,"\n")
    }
    cat("finished p = ",p,"\n")
  }
  cat("finished design = ", d,"\n")
}
stopCluster(cl)
save(sim.df, file="partialLinearSimNet.Rdata")

library(ggplot2)
library(reshape2)
library(latex2exp)
TeX <- latex2exp::TeX
load("partialLinearSimNet.Rdata") # see partialLinearSim.R for simulation
                               # code
df <- melt(sim.df, measure.vars=names(sim.df)[1:3])
ggplot(subset(df,p==2), aes(x=value, colour=variable)) +
  facet_grid(n ~ design) + geom_density() + theme_minimal() +
  xlab(TeX('$\\sqrt{n}(\\hat{\\theta}-\\theta_0)$')) +
  ggtitle("p=2")

ggplot(subset(df,p==4), aes(x=value, colour=variable)) +
  facet_grid(n ~ design) + geom_density() + theme_minimal() +
  xlab(unname(TeX("$\\sqrt{n}(\\hat{\\theta} - \\theta_0)$"))) +
  ggtitle("p=4")

ggplot(subset(df,p==6), aes(x=value, colour=variable)) +
  facet_grid(n ~ design) + geom_density() + theme_minimal() +
  xlab(unname(TeX("$\\sqrt{n}(\\hat{\\theta} - \\theta_0)$"))) +
  ggtitle("p=6")

ggplot(subset(df,p==8), aes(x=value, colour=variable)) +
  facet_grid(n ~ design) + geom_density() + theme_minimal() +
  xlab(unname(TeX("$\\sqrt{n}(\\hat{\\theta} - \\theta_0)$"))) +
  ggtitle("p=8")

3.1.1 Cross-fitting

• Randomly partition into \(K\) subsets \((I_k)_{k=1}^K\)
• \(I^c_k = \{1, ..., n\} \setminus I_k\)
• \(\hat{\eta}_k =\) estimate of \(\eta\) using \(I^c_k\)
• Estimator: \[ \begin{align*} 0 = & \frac{1}{K} \sum_{k=1}^K \frac{K}{n} \sum_{i \in I_k} \psi(w_i;\hat{\theta},\hat{\eta}_k) \\ 0 = & \frac{1}{K} \sum_{k=1}^K \En_k[ \psi(w_i;\hat{\theta},\hat{\eta}_k)] \end{align*} \]

3.1.2 Assumptions

• Linear score \[ \psi(w;\theta,\eta) = \psi^a(w;\eta) \theta + \psi^b(w;\eta) \]
• Near Neyman orthogonality: \[ \lambda_n := \sup_{\eta \in \mathcal{T}_n} \norm{\partial \eta \Er\left[\psi(W;\theta_0,\eta_0)[\eta-\eta_0] \right] } \leq \delta_n n^{-1/2} \]

3.1.3 Proof outline:

• Let \(\hat{J} = \frac{1}{K} \sum_{k=1}^K \En_k [\psi^a(w_i;\hat{\eta}_k)]\), \(J_0 = \Er[\psi^a(w_i;\eta_0)]\), \(R_{n,1} = \hat{J}-J_0\)

• Show: \[ \small \begin{align*} \sqrt{n}(\hat{\theta} - \theta_0) = & -\sqrt{n} J_0^{-1} \En[\psi(w_i;\theta_0,\eta_0)] + \\ & + (J_0^{-1} - \hat{J}^{-1}) \left(\sqrt{n} \En[\psi(w_i;\theta_0,\eta_0)] + \sqrt{n}R_{n,2}\right) + \\ & + \sqrt{n}J_0^{-1}\underbrace{\left(\frac{1}{K} \sum_{k=1}^K \En_k[ \psi(w_i;\theta_0,\hat{\eta}_k)] - \En[\psi(w_i;\theta_0,\eta_0)]\right)}_{R_{n,2}} \end{align*} \]

• Show \(\norm{R_{n,1}} = O_p(n^{-1/2} + r_n)\)

• Show \(\norm{R_{n,2}}= O_p(n^{-1/2} r_n' + \lambda_n + \lambda_n')\)

3.1.4 Proof outline: Lemma 6.1

Lemma 6.1

1. If \(\Pr(\norm{X_m} > \epsilon_m | Y_m) \to_p 0\), then \(\Pr(\norm{X_m}>\epsilon_m) \to 0\).

2. If \(\Er[\norm{X_m}^q/\epsilon_m^q | Y_m] \to_p 0\) for \(q\geq 1\), then \(\Pr(\norm{X_m}>\epsilon_m) \to 0\).

3. If \(\norm{X_m} = O_p(A_m)\) conditional on \(Y_m\) (i.e. for any \(\ell_m \to \infty\), \(\Pr(\norm{X_m} > \ell_m A_m | Y_m) \to_p 0\)), then \(\norm{X_m} = O_p(A_m)\) unconditionally

3.1.5 Proof outline: \(R_{n,1}\)

\[ R_{n,1} = \hat{J}-J_0 = \frac{1}{K} \sum_k \left( \En_k[\psi^a(w_i;\hat{\eta}_k)] - \Er[\psi^a(W;\eta_0)] \right) \]

• \(\norm{\En_k[\psi^a(w_i;\hat{\eta}_k)] - \Er[\psi^a(W;\eta_0)]} \leq U_{1,k} + U_{2,k}\) where \[ \begin{align*} U_{1,k} = & \norm{\En_k[\psi^a(w_i;\hat{\eta}_k)] - \Er[\psi^a(W;\hat{\eta}_k)| I^c_k]} \\ U_{2,k} = & \norm{ \Er[\psi^a(W;\hat{\eta}_k)| I^c_k] - \Er[\psi^a(W;\eta_0)]} \end{align*} \]

3.1.6 Proof outline: \(R_{n,2}\)

• \(R_{n,2} = \frac{1}{K} \sum_{k=1}^K \En_k\left[ \psi(w_i;\theta_0,\hat{\eta}_k) - \psi(w_i;\theta_0,\eta_0) \right]\)
• \(\sqrt{n} \norm{\En_k\left[ \psi(w_i;\theta_0,\hat{\eta}_k) - \psi(w_i;\theta_0,\eta_0) \right]} \leq U_{3,k} + U_{4,k}\) where

\[ \small \begin{align*} U_{3,k} = & \norm{ \frac{1}{\sqrt{n}} \sum_{i \in I_k} \left( \psi(w_i;\theta_0, \hat{\eta}_k) - \psi(w_i;\theta_0,\eta_0) - \Er[ \psi(w_i;\theta_0, \hat{\eta}_k) - \psi(w_i;\theta_0,\eta_0)] \right) } \\ U_{4,k} = & \sqrt{n} \norm{ \Er[ \psi(w_i;\theta_0, \hat{\eta}_k) | I_k^c] - \Er[\psi(w_i;\theta_0,\eta_0)]} \end{align*} \]

• \(U_{4,k} = \sqrt{n} \norm{f_k(1)}\) where

\[ f_k(r) = \Er[\psi(W;\theta_0,\eta_0 + r(\hat{\eta}_k - \eta_0)) | I^c_k] - \Er[\psi(W;\theta_0,\eta_0)] \]

3.1.7 Asymptotic normality

\[ \sqrt{n} \sigma^{-1} (\hat{\theta} - \theta_0) = \frac{1}{\sqrt{n}} \sum_{i=1}^n \bar{\psi}(w_i) + O_p(\rho_n) \leadsto N(0,I) \]

• \(\rho_n := n^{-1/2} + r_n + r_n' + n^{1/2} (\lambda_n +\lambda_n') \lesssim \delta_n\)

• Influence function \[\bar{\psi}(w) = -\sigma^{-1} J_0^{-1} \psi(w;\theta_0,\eta_0)\]

• \(\sigma^2 := J_0^{-1} \Er[\psi(w;\theta_0,\eta_0) \psi(w;\theta_0,\eta_0)'](J_0^{-1})'\)

3.1.8 Creating orthogonal moments

• Need \[ \partial \eta\Er\left[\psi(W;\theta_0,\eta_0)[\eta-\eta_0] \right] \approx 0 \]

• Given an some model, how do we find a suitable \(\psi\)?

3.1.9 Orthogonal scores via concentrating-out

• Original model: \[ (\theta_0, \beta_0) = \argmax_{\theta, \beta} \Er[\ell(W;\theta,\beta)] \]
• Define \[ \eta(\theta) = \beta(\theta) = \argmax_\beta \Er[\ell(W;\theta,\beta)] \]
• First order condition from \(\max_\theta \Er[\ell(W;\theta,\beta(\theta))]\) is \[ 0 = \Er\left[ \underbrace{\frac{\partial \ell}{\partial \theta} + \frac{\partial \ell}{\partial \beta} \frac{d \beta}{d \theta}}_{\psi(W;\theta,\beta(\theta))} \right] \]

3.1.10 Orthogonal scores via projection

• Original model: \(m: \mathcal{W} \times \R^{d_\theta} \times \R^{d_h} \to \R^{d_m}\) \[ \Er[m(W;\theta_0,h_0(Z))|R] = 0 \]
• Let \(A(R)\) be \(d_\theta \times d_m\) moment selection matrix, \(\Omega(R)\) \(d_m \times d_m\) weighting matrix, and \[ \begin{align*} \Gamma(R) = & \partial_{v'} \Er[m(W;\theta_0,v)|R]|_{v=h_0(Z)} \\ G(Z) = & \Er[A(R)'\Omega(R)^{-1} \Gamma(R)|Z] \Er[\Gamma(R)'\Omega(R)^{-1} \Gamma(R) |Z]^{-1} \\ \mu_0(R) = & A(R)'\Omega(R)^{-1} - G(Z) \Gamma(R)'\Omega(R)^{-1} \end{align*} \]
• \(\eta = (\mu, h)\) and \[ \psi(W;\theta, \eta) = \mu(R) m(W;\theta, h(Z)) \]

3.1.11 Example: average derivative

• \(x,y \in \R^1\), \(\Er[y|x] = f_0(x)\), \(p(x) =\) density of \(x\)

• \(\theta_0 = \Er[f_0'(x)]\)

• Joint objective \[ \min_{\theta,f} \Er\left[ (y - f(x))^2 + (\theta - f'(x)^2) \right] \]

• \(f_\theta(x) = \Er[y|x] + \theta \partial_x \log p(x) - f''(x) - f'(x) \partial_x \log p(x)\)

• Concentrated objective: \[ \min_\theta \Er\left[ (y - f_\theta(x))^2 + (\theta - f_\theta'(x)^2) \right] \]

• First order condition at \(f_\theta = f_0\) gives \[ 0 = \Er\left[ (y - f_0(x))\partial_x \log p(x) + (\theta - f_0'(x)) \right] \]

3.1.12 Example : average derivative with endogeneity

• \(x,y \in \R^1\), \(p(x) =\) density of \(x\)

• Model : \(\Er[y - f(x) | z] = 0\) \(\theta_0 = \Er[f_0'(x)]\)

• Joint objective: \[ \min_{\theta,f} \Er\left[ \Er[y - f(x)|z]^2 + (\theta - f'(x))^2 \right] \]

• \(f_\theta(x) = (T^* T)^{-1}\left((T^*\Er[y|z])(x) - \theta \partial_x \log p(x)\right)\)
  • where \(T:\mathcal{L}^2_{p} \to \mathcal{L}^2_{\mu_z}\) with \((T f)(z) = \Er[f(x) |z]\)
  • and \(T^*:\mathcal{L}^2_{\mu_z} \to \mathcal{L}^2_{p}\) with \((T^* g)(z) = \Er[g(z) |x]\)

• Orthogonal moment condition : \[ 0 = \Er\left[ \Er[y - f(x) | z] (T (T^* T)^{-1} \partial_x \log p)(z) + (\theta - f'(x)) \right] \]

3.1.13 Example: average elasticity

• Demand \(D(p)\), quantities \(q\), instruments \(z\) \[\Er[q-D(p) |z] = 0\]

• Average elasticity \(\theta \Er[D'(p)/D(p) | z ]\)

• Joint objective : \[ \min_{\theta,D} \Er\left[ \Er[q - D(p)|z]^2 + (\theta - D'(p)/D(p))^2 \right] \]

3.1.14 Example: control function

\[ \begin{align*} 0 = & \Er[d - p(x,z) | x,z] \\ 0 = & \Er[y - x\beta - g(p(x,z)) | x,z] \end{align*} \]

3.2.1 Fixed, finite groups

• \(G_1, ..., G_K\) finite partition of support \((x)\)

• Estimate \(\Er[y(1) - y(0) | x \in G_k]\) as above

• pros: easy inference, reveals some heterogeneity

• cons: poorly chosen partition hides some heterogeneity, searching partitions violates inference

3.2.2 Generic Machine Learning Inference on Heterogenous Treatment Effects in Randomized Experiments

• Chernozhukov, Demirer, et al. (2018)

• Use machine learning to find partition with sample splitting to allow easy inference

• Randomly partition sample into auxillary and main samples

• Use any method on auxillary sample to estimate \[S(x) = \widehat{\Er[y(1) - y(0) | x]}\] and \[B(x) = \widehat{\Er[y(0)|x]}\]

3.2.3 Best linear projection of CATE

• Randomly partition sample into auxillary and main samples

• Use any method on auxillary sample to estimate \[S(x) = \widehat{\Er[y(1) - y(0) | x]}\] and \[B(x) = \widehat{\Er[y(0)|x]}\]

• Use main sample to regress with weights \((P(x)(1-P(X)))^{-1}\) \[ y = \alpha_0 + \alpha_1 B(x) + \beta_0 (d-P(x)) + \beta_1 (d-P(x))(S(x) - \Er[S(x)]) + \epsilon \]

• \(\hat{\beta}_0, \hat{\beta}_1 \to_p \argmin_{b_0,b_1} \Er[(s_0(x) - b_0 - b_1 (S(x)-E[S(x)]))^2]\)

3.2.4 Inference on CATE

• Inference on \(\Er[y(1) - y(0) | x] = s_0(x)\) challenging when \(x\) high dimensional and/or few restrictions on \(s_0\)

• Pointwise results for random forests : Wager and Athey (2018), Athey, Tibshirani, and Wager (2016)

• Recent review of high dimensional inference : Alexandre Belloni, Chernozhukov, Chetverikov, Hansen, et al. (2018)

3.2.5 Random forest asymptotic normality

• Wager and Athey (2018)

• \(\mu(x) = \Er[y|x]\)

• \(\hat{\mu}(x)\) estimate from honest random forest

  • honest \(=\) trees independent of outcomes being averaged

  • sample-splitting or trees formed using another outcome

• Then \[ \frac{\hat{\mu}(x) - \mu(x)}{\hat{\sigma}_n(x)} \leadsto N(0,1) \]
  • \(\hat{\sigma}_n(x) \to 0\) slower than \(n^{-1/2}\)

3.2.6 Random forest asymptotic normality

• Pointwise result, how to do inference on:
  • \(H_0: \mu(x_1) = \mu(x_2)\)
  • \(\{x: \mu(x) \geq 0 \}\)
  • \(\Pr(\mu(x) \leq 0)\)

3.2.7 Uniform inference

• Alexandre Belloni, Chernozhukov, Chetverikov, Hansen, et al. (2018)
• Alexandre Belloni, Chernozhukov, Chetverikov, and Wei (2018)