# A few questions regarding the practice of heterogeneous treatment effect analysis (a.k.a, interaction detection or subgroup analysis) methods

Imagine I am looking at a randomized experiment between a control and one or more treatment conditions. For example, I have a treatment that aims to get people out of debt. I randomize people to condition and measure whether or not they did (positive outcome) or did not (negative outcome) get out of debt.

Getting an overall treatment effect is easy: I perform a logistic regression, outcome ~ condition. I can see if (any of) my treatment(s) outperform the control.

However, imagine I also have a wealth of "big data" on these people, which are often bucketed, categorical variables. Let's imagine I also know someone's gender, marital status, whether or not they have kids, if their kids are in college, their education, their race, their age, etc.

I would like to know:

First, is there variation in treatment effects for different subgroups? Second, if there is variation, where are these differences at? These can be marginal—just looking at one variable and collapsing across the others: Is the treatment effect for men significant? These can also be specific combinations of many variables: Is there a significant treatment effect for millennial Black males without a college education?

There are two ways people go about answering these questions:

1. Exhaustive subgroup analysis. This involves just analyzing a ton of stuff, generally reporting only what is significant. This has the problem of multiple comparisons and inflating Type I error: We could potentially be finding patterns in noise.

2. Pre-registration. This involves coming up with specific contrasts to test a priori and only testing those. This has the problem of Type II error: We could potentially be missing out on unexpected and important treatment effects, and we often do not have great priors a priori about where we should be seeing the effect.

I would like to do some type of regularization to balance Type I and Type II errors (since p-value adjustments here will be a nightmare with all the possible tests one could do in a "big data" scenario).

What I would like is to give my outcome vector, my condition vector, and a matrix of covariates to an algorithm and have it return the subgroups where there is a significant effect (defining significance as a p-value below .05, a confidence interval that doesn't include zero, or a Bayesian credible interval that doesn't include zero, for example). There are many papers on this currently: causal forests (now subsumed into generalized random forests, using the grf package), Bayesian additive regression trees (using one of the bartMachine or BayesTree or BART packages), support vector machine with LASSO constraints (using the FindIt package), etc.

However, I have found that many of the packages do not have tutorials or user-friendly walkthroughs/documentation of how to do this in practice.

My question is: Can someone provide an explanation of how to use one of these methods that are appropriate for my specific use case? Papers abound that explain the technicalities of them, but I am having a hard time finding out how to translate the details of those papers into how to use them in a valid fashion in practice using a package.

R code for a sample dataset is below. In this data set, there should be a significant positive treatment effect for women and a significant negative treatment effect for millennial Black men.

# generate data ----------------------------------------------------------------
set.seed(1839)
n <- 5000
dat <- 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))
)
dat$outcome <- factor(ifelse( dat$gender == "female" & dat$condition == "treatment", sample(c("positive", "negative"), n, TRUE, c(.50, .50)), ifelse( dat$generation == "millennial" &
dat$race == "black" & dat$condition == "treatment",
sample(c("positive", "negative"), n, TRUE, c(.12, .88)),
sample(c("positive", "negative"), n, TRUE, c(.40, .60))
)
))

# checking women effect --------------------------------------------------------
with(
dat[dat$gender == "female", ], prop.table(table(outcome, condition), 2) ) with( dat[dat$gender == "female", ],
summary(glm(outcome ~ condition, family = binomial()))
)

# checking black millennial men effect -----------------------------------------
with(
dat[dat$gender == "male" & dat$race == "black"
& dat$generation == "millennial", ], prop.table(table(outcome, condition), 2) ) with( dat[dat$gender == "male"
& dat$race == "black" & dat$generation == "millennial", ],
summary(glm(outcome ~ condition, family = binomial()))
)

# but note marginal effect of men not there ------------------------------------
summary(glm(
outcome ~ condition,
dat[dat\$gender == "male", ],
family = binomial()
))

# and overall effect not present -----------------------------------------------
summary(glm(outcome ~ condition, dat, family = binomial()))


Assume I did not know those two subgroups were significant a priori, how could I best find them while still keeping Type I error in check?

One of the issues I am running into is that I hardly ever want to go too granular with subgroup analyses. I would look at treatment effects by one variable (what would be a two-way interaction) or two (what would be a three-way interaction). But hardly ever more. The predict functions in grf and BayesTree require you to have test data that have values for all variables, obviously, which means it makes me do the analyses at the most granular level, which does not let me look at the average treatment effect for women as a whole, for example, or just Gen-X men, for example.