First set the pre-ambles for the model.
Then we generate the data for the simulation. \[ \tau(x) = \left(1 + \frac{1}{1 + e^{-20(x_1 - \frac{1}{3})}}\right) \left(1 + \frac{1}{1 + e^{-20(x_2 - \frac{1}{3})}}\right) \]
effect = function(x) {
(1 + 1 / (1 + exp(-20 * (x[1] - 1/3)))) * (1 + 1 / (1 + exp(-20 * (x[2] - 1/3))))
}
d = 2
X = matrix(runif(n * d, 0, 1), n, d) # features
W = rbinom(n, 1, 0.5) #treatment condition
Y = (W - 0.5) * apply(X, 1, effect) + sigma * rnorm(n)
X.test = matrix(runif(n.test * d, 0, 1), n.test, d)
true.eff = apply(X.test, 1, effect)Now run the random forest using the random_forest
command from grf package. Along with the options as
specified. For full documentation of using the package, check https://grf-labs.github.io/grf/reference/causal_forest.html.
Also, we construct the 95% coverage from the variance estimated using the infitestimal jacknife estimation.
forest = causal_forest(X, Y, W, num.trees = ntree, sample.fraction = sample.fraction, min.node.size = 1)
predictions = predict(forest, X.test,estimate.variance = TRUE)
se.hat <- sqrt(predictions$variance.estimates)
rf.cov = abs(predictions$predictions - true.eff) <= 1.96 * se.hat
rf.covered = mean(rf.cov)
rf.mse = mean((predictions$predictions - true.eff)^2)Nest, we construct the KNN estimator using knn.reg from
the package FNN. The variance is given by
\[ \widehat{Var}(\hat{\tau}(x) ) = \widehat{\tau^2}(x) - \big(\hat{\tau}(x)^2 \big) / (k - 1). \]
We consider both k = 7 and k = 50.
k.small = 7
knn.0.mu = knn.reg(X[W==0,], X.test, Y[W==0], k = k.small)$pred
knn.1.mu = knn.reg(X[W==1,], X.test, Y[W==1], k = k.small)$pred
knn.0.mu2 = knn.reg(X[W==0,], X.test, Y[W==0]^2, k = k.small)$pred
knn.1.mu2 = knn.reg(X[W==1,], X.test, Y[W==1]^2, k = k.small)$pred
knn.0.var = (knn.0.mu2 - knn.0.mu^2) / (k.small - 1)
knn.1.var = (knn.1.mu2 - knn.1.mu^2) / (k.small - 1)
knn.tau = knn.1.mu - knn.0.mu
knn.se = sqrt(knn.0.var + knn.1.var)
knn.cov = abs(knn.tau - true.eff) <= 1.96 * knn.se
knn.covered = mean(knn.cov)
knn.mse = mean((knn.tau - true.eff)^2)
k.big = 50
knnbig.0.mu = knn.reg(X[W==0,], X.test, Y[W==0], k = k.big)$pred
knnbig.1.mu = knn.reg(X[W==1,], X.test, Y[W==1], k = k.big)$pred
knnbig.0.mu2 = knn.reg(X[W==0,], X.test, Y[W==0]^2, k = k.big)$pred
knnbig.1.mu2 = knn.reg(X[W==1,], X.test, Y[W==1]^2, k = k.big)$pred
knnbig.0.var = (knnbig.0.mu2 - knnbig.0.mu^2) / (k.big - 1)
knnbig.1.var = (knnbig.1.mu2 - knnbig.1.mu^2) / (k.big - 1)
knnbig.tau = knnbig.1.mu - knnbig.0.mu
knnbig.se = sqrt(knnbig.0.var + knnbig.1.var)
knnbig.cov = abs(knnbig.tau - true.eff) <= 1.96 * knnbig.se
knnbig.covered = mean(knnbig.cov)
knnbig.mse = mean((knnbig.tau - true.eff)^2)simu.fun = function(d) {
X = matrix(runif(n * d, 0, 1), n, d) # features
W = rbinom(n, 1, 0.5) #treatment condition
Y = (W - 0.5) * apply(X, 1, effect) + sigma * rnorm(n)
X.test = matrix(runif(n.test * d, 0, 1), n.test, d)
true.eff = apply(X.test, 1, effect)
#
# random forest
#
forest = causal_forest(X, Y, W, num.trees = ntree, sample.fraction = sample.fraction, min.node.size = 1)
predictions = predict(forest, X.test,estimate.variance = TRUE)
se.hat <- sqrt(predictions$variance.estimates)
rf.cov = abs(predictions$predictions - true.eff) <= 1.96 * se.hat
rf.covered = mean(rf.cov)
rf.mse = mean((predictions$predictions - true.eff)^2)
k.small = 7
knn.0.mu = knn.reg(X[W==0,], X.test, Y[W==0], k = k.small)$pred
knn.1.mu = knn.reg(X[W==1,], X.test, Y[W==1], k = k.small)$pred
knn.0.mu2 = knn.reg(X[W==0,], X.test, Y[W==0]^2, k = k.small)$pred
knn.1.mu2 = knn.reg(X[W==1,], X.test, Y[W==1]^2, k = k.small)$pred
knn.0.var = (knn.0.mu2 - knn.0.mu^2) / (k.small - 1)
knn.1.var = (knn.1.mu2 - knn.1.mu^2) / (k.small - 1)
knn.tau = knn.1.mu - knn.0.mu
knn.se = sqrt(knn.0.var + knn.1.var)
knn.cov = abs(knn.tau - true.eff) <= 1.96 * knn.se
knn.covered = mean(knn.cov)
knn.mse = mean((knn.tau - true.eff)^2)
k.big = 50
knnbig.0.mu = knn.reg(X[W==0,], X.test, Y[W==0], k = k.big)$pred
knnbig.1.mu = knn.reg(X[W==1,], X.test, Y[W==1], k = k.big)$pred
knnbig.0.mu2 = knn.reg(X[W==0,], X.test, Y[W==0]^2, k = k.big)$pred
knnbig.1.mu2 = knn.reg(X[W==1,], X.test, Y[W==1]^2, k = k.big)$pred
knnbig.0.var = (knnbig.0.mu2 - knnbig.0.mu^2) / (k.big - 1)
knnbig.1.var = (knnbig.1.mu2 - knnbig.1.mu^2) / (k.big - 1)
knnbig.tau = knnbig.1.mu - knnbig.0.mu
knnbig.se = sqrt(knnbig.0.var + knnbig.1.var)
knnbig.cov = abs(knnbig.tau - true.eff) <= 1.96 * knnbig.se
knnbig.covered = mean(knnbig.cov)
knnbig.mse = mean((knnbig.tau - true.eff)^2)
c(rf.covered = rf.covered,
rf.mse = rf.mse,
knn.covered = knn.covered,
knn.mse = knn.mse,
knnbig.covered = knnbig.covered,
knnbig.mse = knnbig.mse)
}
results.raw = lapply(dvals, function(d) {
print(paste("NOW RUNNING:", d))
res.d = sapply(1:simu.reps, function(iter) simu.fun(d))
res.fixed = data.frame(t(res.d))
print(paste("RESULT AT", d, "IS", colMeans(res.fixed)))
res.fixed
})## [1] "NOW RUNNING: 2"
## [1] "RESULT AT 2 IS 0.9267" "RESULT AT 2 IS 0.0733130342931873"
## [3] "RESULT AT 2 IS 0.926075" "RESULT AT 2 IS 0.284735781192886"
## [5] "RESULT AT 2 IS 0.94535" "RESULT AT 2 IS 0.0412793870690497"
## [1] "NOW RUNNING: 3"
## [1] "RESULT AT 3 IS 0.933625" "RESULT AT 3 IS 0.0492377323825507"
## [3] "RESULT AT 3 IS 0.9245" "RESULT AT 3 IS 0.29537483991694"
## [5] "RESULT AT 3 IS 0.915525" "RESULT AT 3 IS 0.0529459567121625"
## [1] "NOW RUNNING: 4"
## [1] "RESULT AT 4 IS 0.940475" "RESULT AT 4 IS 0.0408055644982495"
## [3] "RESULT AT 4 IS 0.922775" "RESULT AT 4 IS 0.300349708474548"
## [5] "RESULT AT 4 IS 0.84505" "RESULT AT 4 IS 0.07857191037481"
## [1] "NOW RUNNING: 5"
## [1] "RESULT AT 5 IS 0.949875" "RESULT AT 5 IS 0.0347339669693442"
## [3] "RESULT AT 5 IS 0.920225" "RESULT AT 5 IS 0.314260952972744"
## [5] "RESULT AT 5 IS 0.76905" "RESULT AT 5 IS 0.11035786902589"
## [1] "NOW RUNNING: 6"
## [1] "RESULT AT 6 IS 0.9534" "RESULT AT 6 IS 0.0315981733251053"
## [3] "RESULT AT 6 IS 0.91" "RESULT AT 6 IS 0.334370383389004"
## [5] "RESULT AT 6 IS 0.680725" "RESULT AT 6 IS 0.149320778550552"
## [1] "NOW RUNNING: 8"
## [1] "RESULT AT 8 IS 0.9579" "RESULT AT 8 IS 0.0285608214236088"
## [3] "RESULT AT 8 IS 0.893725" "RESULT AT 8 IS 0.381910923441542"
## [5] "RESULT AT 8 IS 0.561525" "RESULT AT 8 IS 0.217390592079868"results.condensed = lapply(results.raw, function(RR) {
RR.mu = colMeans(RR)
RR.var = sapply(RR, var) / (nrow(RR) - 1)
rbind("mu"=RR.mu, "se"=sqrt(RR.var))
})
results.condensed## [[1]]
## rf.covered rf.mse knn.covered knn.mse knnbig.covered knnbig.mse
## mu 0.926700000 0.073313034 0.926075000 0.28473578 0.945350000 0.041279387
## se 0.001984224 0.001204711 0.001866831 0.00333079 0.003319211 0.001083068
##
## [[2]]
## rf.covered rf.mse knn.covered knn.mse knnbig.covered knnbig.mse
## mu 0.933625000 0.049237732 0.92450000 0.295374840 0.915525000 0.052945957
## se 0.002736677 0.001002018 0.00117893 0.002268658 0.004043808 0.001117326
##
## [[3]]
## rf.covered rf.mse knn.covered knn.mse knnbig.covered knnbig.mse
## mu 0.940475000 0.040805564 0.922775000 0.300349708 0.845050000 0.078571910
## se 0.002003937 0.000844419 0.001522778 0.002558538 0.004267201 0.001406867
##
## [[4]]
## rf.covered rf.mse knn.covered knn.mse knnbig.covered knnbig.mse
## mu 0.94987500 0.0347339670 0.920225000 0.314260953 0.769050000 0.110357869
## se 0.00225763 0.0007970887 0.001670208 0.002724653 0.004518928 0.001768862
##
## [[5]]
## rf.covered rf.mse knn.covered knn.mse knnbig.covered knnbig.mse
## mu 0.953400000 0.0315981733 0.910000000 0.334370383 0.680725000 0.149320779
## se 0.002149813 0.0008507164 0.001782743 0.002560929 0.005020535 0.002157183
##
## [[6]]
## rf.covered rf.mse knn.covered knn.mse knnbig.covered knnbig.mse
## mu 0.957900000 0.0285608214 0.893725000 0.381910923 0.561525000 0.217390592
## se 0.002481437 0.0008453122 0.001670405 0.003758819 0.004464925 0.002404839results.parsed = lapply(results.condensed, function(RR) {
apply(RR, 2, function(arg) {
paste0(round(arg[1], 2), " (", round(100 * arg[2], 0), ")")
})
})
results.table = data.frame(cbind(d=dvals, Reduce(rbind, results.parsed)))
results.table## d rf.covered rf.mse knn.covered knn.mse knnbig.covered knnbig.mse
## init 2 0.93 (0) 0.07 (0) 0.93 (0) 0.28 (0) 0.95 (0) 0.04 (0)
## X 3 0.93 (0) 0.05 (0) 0.92 (0) 0.3 (0) 0.92 (0) 0.05 (0)
## X.1 4 0.94 (0) 0.04 (0) 0.92 (0) 0.3 (0) 0.85 (0) 0.08 (0)
## X.2 5 0.95 (0) 0.03 (0) 0.92 (0) 0.31 (0) 0.77 (0) 0.11 (0)
## X.3 6 0.95 (0) 0.03 (0) 0.91 (0) 0.33 (0) 0.68 (1) 0.15 (0)
## X.4 8 0.96 (0) 0.03 (0) 0.89 (0) 0.38 (0) 0.56 (0) 0.22 (0)results.table = results.table[,c(1, 3, 5, 7, 2, 4, 6)]
xtab = xtable(results.table)
print(xtab, include.rownames = FALSE)## % latex table generated in R 4.3.3 by xtable 1.8-4 package
## % Sun Apr 21 16:01:57 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{lllllll}
## \hline
## d & rf.mse & knn.mse & knnbig.mse & rf.covered & knn.covered & knnbig.covered \\
## \hline
## 2 & 0.07 (0) & 0.28 (0) & 0.04 (0) & 0.93 (0) & 0.93 (0) & 0.95 (0) \\
## 3 & 0.05 (0) & 0.3 (0) & 0.05 (0) & 0.93 (0) & 0.92 (0) & 0.92 (0) \\
## 4 & 0.04 (0) & 0.3 (0) & 0.08 (0) & 0.94 (0) & 0.92 (0) & 0.85 (0) \\
## 5 & 0.03 (0) & 0.31 (0) & 0.11 (0) & 0.95 (0) & 0.92 (0) & 0.77 (0) \\
## 6 & 0.03 (0) & 0.33 (0) & 0.15 (0) & 0.95 (0) & 0.91 (0) & 0.68 (1) \\
## 8 & 0.03 (0) & 0.38 (0) & 0.22 (0) & 0.96 (0) & 0.89 (0) & 0.56 (0) \\
## \hline
## \end{tabular}
## \end{table}