Can you include interactions with the treatment variable as a covariate in linear regression?

I'm quite new to regression modelling. I'm trying to predict the average treatment effect of a weight loss intervention by comparing outcomes in control and treated groups. The model currently includes gender, age, age-squared, and treatment as covariates. But being male predicts better outcomes in the treated group and worse outcomes in controls. Adding an interaction term between gender and treatment improves model fit (about 20 AIC points) and the term is significant. But I'm not sure if it makes sense to do that, if the goal is to assess the effect of the treatment. I can't immediately think why it would be a problem, but I haven't seen it done before either.

• Can you tell us more about the study design? Were participants randomized to intervention/control arm? Why use a quadratic effect of age instead of something less biased like a spline? Commented Apr 23, 2022 at 14:46
• It is an observational (non-randomised) study. I hadn't heard of splines; only quadratics were covered in my stats classes. Commented Apr 23, 2022 at 15:24
• See mkspline for how to generate the predictors for restricted cubic splines and incorporate them into your model for any continuous covariate (age, BMI,...). They generally work much better than quadratic polynomials.
– EdM
Commented Apr 23, 2022 at 20:53
• Are you developing this model to make causal statements about the treatment effect? Since your study is observational, not randomized, and the only covariates you have are age and gender, you won't be able to make a strong argument that the estimated effect is close to the causal effect of the treatment. Commented Apr 24, 2022 at 23:06
• @dipetkov - yes, that was my conclusion! Commented Apr 26, 2022 at 17:39

If it makes a better prediction, then of course it makes sense to include the interaction in the model. In assessing the treatment effect, one has to recognize that by including the interaction with gender, you are saying that the treatment effect depends on gender. So what I suggest doing is estimating the means and comparisons separately at each gender. Something like this, if using R:

library(emmeans)
EMM <- emmeans(model, ~ treat | gender)
EMM ## display the estimated marginal means
contrast(EMM, "pairwise")  ## pairwise comparisons


The marginal means are computed at the average age. Since you have age squared in the model, it is a little tricky because if age^2 is a separate variable, it will also by default use the average of (age^2), which makes no sense because that's not the same as (average age)^2. Your best bet is to fit the model with those variables computed, e.g., poly(age, degree = 2).

Stata help?

Can a Stata expert please show us how to do this in Stata?

• Thanks. I'm using Stata. I am also doing subgroup analyses by age and gender, but I need to produce a combined model as well. Commented Apr 23, 2022 at 15:25
• I also need to do regression adjustment, for which I won't be able to include the interactions with treatment, so I guess the combined model with those interactions might give a better estimate of average treatment effect than RA? Commented Apr 23, 2022 at 15:28
• I would recommend AGAINST using some average treatment effect. It would give you the treatment effect in a case that cannot occur! Surely there's a way to do this in Stata. I think there is something called margins in Stata that's something like emmeans. Commented Apr 23, 2022 at 15:31