After reading a few of the papers on causal forests, it seems like the idea is the trees will decide to split more often on variables that cause the treatment effect to vary (that is, moderators of the experimental effect). However, the object returned by causal_forest seems to work like any other machine learning algorithm focused on prediction—not causal statistical inference—in the way the predict function associated with it works.

How can I use the causal_forest function to find where the treatment varies?

For example, I simulate a randomized experiment with a binary outcome where there is a positive effect for women and a negative effect for black men. There are also nuisance variables in here. How would I use causal_forest to tell me that, "The trees are tending to split most on gender, and then when it leads to male, it the trees also are more likely to split on race." This would help show me where the treatment effects are occurring.

Here are the data:

set.seed(1839)
n <- 5000
X <- data.frame(
  condition = factor(sample(c("control", "treatment"), n, TRUE)),
  gender = factor(sample(c("male", "female"), n, TRUE)),
  race = factor(sample(c("black", "white", "hispanic", "asian"), n, TRUE)),
  generation = factor(sample(c("millennial", "x", "babyboomer"), n, TRUE)),
  has_kids = factor(sample(c("no", "yes"), n, TRUE))
)
y <- factor(ifelse(
  X$gender == "female" & X$condition == "treatment",
  sample(c("positive", "negative"), n, TRUE, c(.50, .50)),
  ifelse(X$race == "black" & X$condition == "treatment",
    sample(c("positive", "negative"), n, TRUE, c(.12, .88)),
    sample(c("positive", "negative"), n, TRUE, c(.40, .60))
  )
))

And then I fit a model using grf::causal_forest:

library(grf)
mod <- causal_forest(
  model.matrix(~ -1 + ., X[, -1]), # convert factors to dummies
  (as.numeric(y) - 1), # convert outcome to 0 or 1
  (as.numeric(X[, 1]) - 1), # convert treatment to 0 or 1
  num.trees = 4000
)

Looking at mod, I see some variable importance issues that hints at where the splits are occurring most:

Variable importance: 
    1     2     3     4     5     6     7     8 
0.233 0.203 0.267 0.061 0.055 0.046 0.081 0.054 

However, this doesn't capture how multiple X variables may depend on one another in their interaction with W on Y (what may be called a three-way interaction in the general linear model world).

If I use the predict() function, I can only look at individual-level conditional treatment effects and their associated variance estimates. For example, just the first row of my data:

predict(mod, model.matrix(~ -1 + ., X[1, -1]), estimate.variance = TRUE)

However, how could I estimate the treatment effect and variance for the category women, collapsing across all other variables? Or the combination of male and Black, collapsing across all other variables? Is this as simple as averaging the treatment effects on the training set within the groups of interest? And if so, how are confidence intervals calculated from those?

Additionally, can the causal_forest function tell me where these variations are most likely to occur? It seems like this is possible from the papers I have read on causal forests (as well as earlier papers on causal trees and transformed outcome trees, etc.), but if I may very well be mistaken.

Description of the bug I am estimating causal forest on a large dataset (around 2-3 GB with 1.6 million observations and 8 independent variables, as well as 1 dependent variable) with 4000 trees. When I use a binary treatment W, the forest runs fine. However, when I switch to a continuous treatment W, R crashes.

I am running it on a HPC with 32 CPUs and 12GB RAM per CPU.

The error message is (I masked some of my memory addresses with *):

*** Error in `/PATH/lib/R/bin/exec/R': break adjusted to free malloc space: 0x000******** ***
======= Backtrace: =========
/lib64/libc.so.6(+0x82257)[0x7f5*********]
/lib64/libc.so.6(+0x82cea)[0x7f5*********]
/lib64/libc.so.6(__libc_malloc+0xc7)[0x7f5*********]
/usr/local/Cluster-Apps/dmtcp/dmtcp-2.6.0-intel-17.0.4/lib/dmtcp/libdmtcp_alloc.so(malloc+0x22)[0x7f5*********]
/PATH/lib/R/bin/exec/../../../libstdc++.so.6(_Znwm+0x15)[0x7f5*********]
/PATH/lib/R/library/grf/libs/grf.so(_ZNSt6vectorImSaImEE17_M_default_appendEm+0xc3)[0x7f5*********]
/PATH/lib/R/library/grf/libs/grf.so(_ZNK3grf4Tree15find_leaf_nodesERKNS_4DataERKSt6vectorImSaImEE+0xb7)[0x7f5*********]
/PATH/lib/R/library/grf/libs/grf.so(_ZNK3grf11TreeTrainer21repopulate_leaf_nodesERKSt10unique_ptrINS_4TreeESt14default_deleteIS2_EERKNS_4DataERKSt6vectorImSaImEEb+0xe0)[0x7f5*********]
/PATH/lib/R/library/grf/libs/grf.so(_ZNK3grf11TreeTrainer5trainERKNS_4DataERNS_13RandomSamplerERKSt6vectorImSaImEERKNS_11TreeOptionsE+0x408)[0x7f5*********]
/PATH/lib/R/library/grf/libs/grf.so(_ZNK3grf13ForestTrainer14train_ci_groupERKNS_4DataERNS_13RandomSamplerERKNS_13ForestOptionsE+0x135)[0x7f5*********]
/PATH/lib/R/library/grf/libs/grf.so(_ZNK3grf13ForestTrainer11train_batchEmmRKNS_4DataERKNS_13ForestOptionsE+0x1a9)[0x7f5*********]
/PATH/lib/R/library/grf/libs/grf.so(_ZNSt17_Function_handlerIFSt10unique_ptrINSt13__future_base12_Result_baseENS2_8_DeleterEEvENS1_12_Task_setterIS0_INS1_7_ResultISt6vectorIS0_IN3grf4TreeESt14default_deleteISA_EESaISD_EEEES3_ENSt6thread8_InvokerISt5tupleIJMNS9_13ForestTrainerEKFSF_mmRKNS9_4DataERKNS9_13ForestOptionsEEPKSL_mmSO_SP_EEEESF_EEE9_M_invokeERKSt9_Any_data+0x67)[0x7f5*********]
/PATH/lib/R/library/grf/libs/grf.so(+0x526e9)[0x7f5*********]
/lib64/libpthread.so.0(+0x620b)[0x7f5*********]
/PATH/lib/R/library/grf/libs/grf.so(_ZNSt6thread11_State_implINS_8_InvokerISt5tupleIJZNSt13__future_base17_Async_state_implINS1_IS2_IJMN3grf13ForestTrainerEKFSt6vectorISt10unique_ptrINS5_4TreeESt14default_deleteIS9_EESaISC_EEmmRKNS5_4DataERKNS5_13ForestOptionsEEPKS6_mmSH_SI_EEEESE_EC4EOSQ_EUlvE_EEEEE6_M_runEv+0x116)[0x7f5*********]
/PATH/lib/R/bin/exec/../../../libstdc++.so.6(+0xc819d)[0x7f5*********]
/usr/local/Cluster-Apps/dmtcp/dmtcp-2.6.0-intel-17.0.4/lib/dmtcp/libdmtcp.so(+0x35eeb)[0x7f5*********]
/lib64/libpthread.so.0(+0x7ea5)[0x7f5*********]
/usr/local/Cluster-Apps/dmtcp/dmtcp-2.6.0-intel-17.0.4/lib/dmtcp/libdmtcp.so(+0x35cad)[0x7f5*********]
/lib64/libc.so.6(clone+0x6d)[0x7f5*********]

Does continuous treatment consume more memory than binary treatment?

GRF version development

Description of the bug Download the dataset at this link. This is a dataset by LaLonde (1986), Dehejia and Wahba (1999) on the NSW training program. I am not sure if there is copyright - if so it belongs to them.

The data set is very imbalanced - there are around 10000 untreated example and less than 300 treated example.

When I run the code below in the development version, the average_treatment_effect(cps.forest, target.sample = "all") returns NaN. However, the code returns numbers in the releasedversion 0.10.2. It only returns NaN in the development branch.

Is this supposed to be a feature (i.e. returns NaN when data is very imbalanced), or is it actually a bug?

The average treatment effect for the treated and the overlap weighted treatment effects are all normal. Although the numbers differ slightly from the released branch as well.

Steps to reproduce

#compile the grf development code and include it into the library
library(dplyr)
select <- dplyr::select
TREES_CPS=1000
SEED=6000
set.seed(SEED)
dat=read.csv("bug.csv")
forest=causal_forest(seed=SEED,num.trees = TREES_CPS, X=as.matrix(select(dat, re74, re75, age, education, black, hispanic, married, nodegree)), Y=dat$re78, W=dat$treat)
average_treatment_effect(forest, target.sample = "all")

Output in Development Version

> average_treatment_effect(forest, target.sample = "all")
estimate  std.err 
     NaN      NaN 

Output in Release Version

> average_treatment_effect(forest, target.sample = "all")
 estimate   std.err 
-2058.140  1401.162 

GRF version development

When honesty is enabled, the training subsample is further split in half before performing splitting. With small datasets, this may not leave enough information for the algorithm to determine high-quality splits.

This issue is still pending a concrete proposal on how it should be addressed.

I understand that the default option of the causal_forest function is for it to use orthogonalization; that is, it incorporates both the actual values of the outcome (Y) and treatment (W), as well as the predicted values of the outcome (Y.hat) and treatment (W.hat) to estimate CATEs as well as to make splits throughout the forest. However, I have a question related to this:

Q: Is there any way to turn "off" the orthogonalization? That is, for demonstration purposes (not practical purposes), I wanted to simulate how performance of the causal_forest degrades if it is trained on only the actual values of the outcome (Y) and treatment (W), and not the residual values (i.e., actual - predicted) ?

Hello, I am trying to use the Causal survival forest prediction function. I have some questions about the "causal_survival_forest" function as follows:

1. I cannot use the arguments "target" and "horizon ";
2. In the details about the function as follows,why "'D[Y >= Y.max] <- 1' and 'Y[Y >= Y.max] <- Y.max'." and not "'D[Y >= Y.max] <- 0' and 'Y[Y >= Y.max] <- Y.max'." Details An important assumption for identifying the conditional average treatment effect tau(X) is that there exists a fixed positive constant M such that the probability of observing an event time past the maximum follow-up time Y.max is at least M. This may be an issue with data where most endpoint observations are censored. The suggested resolution is to re-define the estimand as the treatment effect up to some suitable maximum follow-up time Y.max. One can do this in practice by thresholding Y before running causal_survival_forest: 'D[Y >= Y.max] <- 1' and 'Y[Y >= Y.max] <- Y.max'. For details see Cui et al. (2020). The computational complexity of this estimator scales with the cardinality of the event times Y. If the number of samples is large and the Y grid dense, consider rounding the event times (or supply a coarser grid with the 'failure.times' argument).
3. When I used the "causal_survival_forest" function to analysis my RCT (random control trial) data, it's expected that tau(X) = E[Y(1) - Y(0) | X = x] > 0 (1:treatment arm; 0: control arm) where E[Y] is the expected survival time = integral(survival function), but the results of this analysis was tau(X) = E[Y(1) - Y(0) | X = x] < 0 contrary to what I thought, which were very unreasonable.

Could you give me an answer.

Thanks very much! Xin

A start to providing the functionality requested in https://github.com/grf-labs/grf/issues/403#issuecomment-486915981:

Finally, it would be nice to provide a short function which performs the above steps, or at least add this information to our algorithm reference.

This does not include the general calculation if multiple, or non {0,1} outcomes exist, but happy to consider that if helpful.

Hi, I am writing this in an effort to better understand the predict() and average_treatment_effect() functions, particularly in regards to when one should be utilized over the other to estimate conditional average treatment effects (CATEs) in an honest causal forest. Additional details related to this query are listed below:

(1) Research Goal: After building an honest causal_forest on my dataset, I would now like to calculate Conditional Average Treatment Effects (CATEs) within specific strata of covariates on the same dataset.

(2) Initial Plan: Based on this application paper by Athey & Wager, it seems that I should be using the predict() function in order to estimate these CATEs using an "honest" approach. Based on the documentation on predict(), it appears that by default, this function estimates the treatment effects such that these effects are estimated for every observation using only the trees in the forest which did not use that particular observation when it was modeled--so, out-of-bag estimation.

(3) Question: However, I understand that there is additionally an average_treatment_effect function, which can also supposedly estimate Conditional Average Treatment Effects in a causal forest in a doubly robust fashion. I would like to better understand the differences between these two functions (predict() and average_treatment_effect()), and the different circumstances in which one function should be used over the other to estimate CATEs on data within an honest causal forest. Evidently, the math behind each function differs, as shown in my attached code that I ran on a mock dataset. This attached R code can be used to reproduce the very different conditional average treatment effects that I calculated for a specific subset of individuals when I used the predict() approach vs. average_treatment_effect().

In addition to explaining which function is better for calculating conditional average treatment effects in various circumstances, could someone also please explain in layman's terms what each function is doing behind the scenes? For example, is the average_treatment_effect function using all of the trees to estimate the treatment effects (and not out-of-bag?)? Also - how are propensity scores utilized for the average_treatment_effect function?

Thanks!

Steps to reproduce Please find the attached code. cates_cf.txt

GRF version 2.0.2

For propensity score matching, there are some ways to test whether the samples are indeed matched well.

Are there similar methods to exam the quality of matching for the causal forest?

Thanks a lot for your time!

Right now GRF is based on an IID assumption. It would be nice to be able to use GRF on data with a hierarchical structure. This is especially relevant in RCTs where studies occur across administrative units like provinces, towns, school districts, etc.

I'm trying to build a single causal tree using the following code:

model <- causal_forest(X_train, y_train, z_train, num.trees = 1)

However, I noticed that the causal_forest method has the parameter sample.fraction, which defines the fraction of the data that is used to build each tree (and is 0.5 by default). Because I want to use the entire data set to build the causal tree, I want to set this to 1, but when I run the following code:

model <- causal_forest(X_train, y_train, z_train, num.trees = 1, sample.fraction=1)

I get the following error message:

"Error in causal_train(data$default, data$sparse, outcome.index, treatment.index, : When confidence intervals are enabled, the sampling fraction must be less than 0.5."

Could you please tell me how to disable confidence intervals in order to build a tree using the entire sample? Thanks in advance!

Hi grf team,

I am currently trying to follow how instrumental forests select balanced splits. I have read the corresponding section in the algorithm reference for causal forests, but I suppose that this does not fully apply to instrumental forests.

In particular, I'm interested in:

1. how min.node.size is determined
2. what the node size measure is which is used together with alpha and imbalance.penalty
3. what changes in 1. and 2. if stabilize.splits is set to FALSE

Thank you for your support & best regards, Jens

I was wondering if you might be able to help me with using the causal_survival_forest function to estimate individual treatment effect in R environment. I am having some difficulty understanding how to assess model calibration and check discrimination performance when using this function. Would you happen to have any recommendations or resources that might be helpful in this regard? I would really appreciate any guidance you might be able to provide.

I am trying to understand Causal Forest's behavior in the tail of the distribution of the covariates. I run simulations of the type below and often find that CF estimates are constant towards the tails (see plot) which means that CFs are biased there. In the simulation below I compare this to the naïve approach using random forests for each of Y0 and Y1 to predict the outcomes and then obtain tau by the difference in predictions \hat{Y1}-\hat{Y0}. That approach seemingly does not suffer from the bias problem in the tails (but has larger variance throughout). In larger samples (e.g. n=1e4) the problem seems to persist. Is there anything I can do to get a better fit in the tails? Thanks.

library(grf)
library(ranger)

## Simulate non-linear ps-scores and y-models
e.x = function(x) 1/ (1 + exp(- ( 3 * x) ))
mu.0 = function(x) sin(-1/2- 4*x)
mu.1 = function(x) sin(1/2 + 4*x)
tau  = function(x) (mu.1(x) - mu.0(x))
m.x = function(x) mu.0(x) + e.x(x) * tau(x)

## Sample training data
set.seed(2023)
n = 1000
sd.y = 0.6
Z = rnorm(n, mean=0, sd = 0.3)
y0 = mu.0(Z) + rnorm(n, sd = sd.y)
y1 = mu.1(Z) + rnorm(n, sd = sd.y)
e  = e.x(Z)
W  = rbinom(n, 1, e)
y  = W*y1 + (1-W) * y0
df = data.frame(Z = Z, y, y0, y1, e, W = factor(W), W.num = W)

## Estimate causal forest with tuning on
cf = causal_forest(X = cbind(df$Z), Y = df$y, W = df$W.num, num.trees = 2e3, tune.parameters = 'all')

## Estimate heterogeineity using conditional means random forests
rf0 = ranger(y ~ Z, data = df[df$W.num==0,], num.trees = 5e3 )
rf1 = ranger(y ~ Z, data = df[df$W.num==1,], num.trees = 5e3 )

## Create test data
Z.test  = seq(-1,1,0.01)
df.test = data.frame(Z = Z.test)
tau.test = tau(Z.test)
tau.test.cf = predict(cf, newdata = df.test)[,1]
tau.test.cdmrf = predict(rf1, data= df.test)$predictions - predict(rf0, data= df.test)$predictions

## Compare fits
par(mfrow=c(1,2))
plot(Z.test,tau.test, ylim=c(-3,3),ty='l', main='Predictions vs true effect')
lines(Z.test,tau.test.cf,col=2)
lines(Z.test,tau.test.cdmrf,col=4)
plot(Z.test,tau.test.cf-tau.test, ylim=c(-3,3),col=2, ty='l', main= 'Error')
lines(Z.test,tau.test.cdmrf-tau.test,col=4)
abline(h=0)
legend('bottomright', legend = c('Causal Forest', 'Standard Forests'),lty=c(1,1),col=c(2,4))

Hi, everyone. Thank you for posting this fancy package. I want to test whether a heterogeneous treatment effect exists in my instrumental forest model. For your information, the sample size is around 2,300 and the number of covariates is around 40 in my data.

There are 2 major challenges that I am facing: 1) test_calibration function does not support instrumental_forest. While the function supports some forest models, it is impossible to find any treatment heterogeneity in my instrumental forest model using the function. Is there any technical difficulty in supporting instrumental_forest in the test_calibration function? I wonder if there would be any plans for updates as best_linear_projection started to support instrumental_forest recently.

2) Rank average treatment effect is unstable. Since I can not use the test_calibration function, I have tried to use the rank_average_treatment function. However, I found that the p-values vary significantly based on the parameters I am using. For example, if I change tune.parameters from ‘all’ to c('sample.fraction', 'mtry', 'min.node.size', 'alpha', 'imbalance.penalty'), the p-value increases from 0.06 to 0.97, or decreases from 0.59 to 0.20 depending on different Ys. Moreover, the p-value also varies a lot if I change the seed of the instrumental forest model(e.g., seed=123 to seed=119, and so on). The following is the code that I'm using:

set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion", sample.kind = "Rejection") 
cf.priority <- instrumental_forest( X[train, ], Y[train], W[train], Z[train], 
num.trees = 50000, 
#tune.parameters = 'all', 
tune.parameters = c('sample.fraction', 'mtry', 'min.node.size', 'alpha', 'imbalance.penalty'), 
tune.num.trees = 4000, tune.num.reps = 250, tune.num.draws = 4500)
set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion", sample.kind = "Rejection")

# Estimate AUTOC on held-out data.
cf.eval <- instrumental_forest( X[-train, ], Y[-train], W[-train], Z[-train], 
num.trees = 50000, 
#tune.parameters = 'all', 
tune.parameters = c('sample.fraction', 'mtry', 'min.node.size', 'alpha', 'imbalance.penalty'),
tune.num.trees = 4000, tune.num.reps = 250, tune.num.draws = 4500)

Can we say that if a certain hyperparameter set yields a low p-value, then the tuning is valid? I wonder if there is any rule of thumb in tuning these hyperparameters. Thank you for your time and all the work!

Best, Minje

Description of the bug Regression forest allows the X matrix to have some incomplete cases, but local linear forest returns an error. Not sure if there's some technical reason why ll forests can't handle incomplete cases?

Steps to reproduce

#Toy data 
Y    <- as.vector(rnorm(100))
X    <- data.frame(x1 = rnorm(100), x2 = rnorm(100))
#
#Add NAs 
X$x1 <- ifelse(X$x1 > 0,X$x1,NA)
#
#Let's try an r forest 
regression_forest(Y = Y,
                  X = X)
#R forest runs fine 
#
# Now let's try ll forest 
ll_regression_forest(Y = Y,
                     X = X)
#
#ll forest returns: Error in validate_X(X) : The feature matrix X contains at least one NA.

GRF version GRF version 2.2.0

Hi everyone, Thank you for posting this package. I am trying to understand the intuition of the calibration test. I see from the source code that the calibration test does the following for a causal forest:

preds <- predict(forest)$predictions mean.pred <- mean(preds) DF <- data.frame( target = unname(forest$Y.orig - forest$Y.hat), mean.forest.prediction = unname(forest$W.orig - forest$W.hat) * mean.pred,differential.forest.prediction = unname(forest$W.orig - forest$W.hat) *(preds - mean.pred))

summary(lm(target~ mean.forest.prediction + differential.forest.prediction +0, data=DF))

The target are the orthogonalized outcomes. But then, the target is not regressed on the mean forest prediction and the differential forest prediction alone, but on the product of those two and the orthogonalized treatments....

I understand that those orthogonalized outcomes are the outcome variable of the forest. But i don't understand why the mean forest prediction needs to be multiplied by the orthogonalized treatments for the test to work.

I am just so curious why the description of the function says that the test computes the best linear predictor of the target estimand using the forest prediction as well as the mean forest prediction as the sole two regressors. It seems to me that the test uses the forest prediction and the mean forest prediction multiplied by the orthogonalized treatment status as the sole two regressors.

Or is this clarification redundant?

Any guidance on this would be greatly appreciated.

Lucy