# Is it okay to include intervention usage in the model of a randomised-controlled trial?

Suppose we have a randomised-controlled trial where, for the $i$-th individual, we are given

• an outcome $Y_i$
• an indicator $Z_i$ ($Z_i = 1$ for intervention group, $Z_i = 0$ for control group)
• a vector of baseline covariates $X_i$

We decide to fit a linear model $Y_i = \alpha + \beta Z_i + \gamma X_i$.

We are interested to estimate $\mathbb E[Y|Z=1]-\mathbb E[Y|Z=0]$ which should be $\hat\beta$ in this model. Note that $X$ and $Z$ are independent.

Now suppose that we also have data $G_i$ about the extent to which the intervention was actually made use of. I'm aware that including a summand $\delta G_i$ in the model can lead to bias since $G$ and $Z$ are not independent. However, I was wondering if there are other ways to include $G$? How about an interaction term $\beta G_i Z_i$? Or a completely different model? Are there standard ways of doing this?

• Yes, just adding the interaction term is fine. If the $G_i$ are zero when there's no treatment you can just add a $\delta G_i$ term. Then $\beta$ is the treatment/no treatment effect and $\delta$ is the effect of dose (or whatever it is) within the treatment group. – Will Jun 10 '17 at 17:28
• Are you sure that, after adding $\delta G_i$, we still have $\hat\beta=\mathbb E[Y|Z=1]-E[Y|Z=0]$? Because $\mathbb E[G|Z=1]\neq\mathbb E[G|Z=0]=0$ so the term $\delta G_i$ wouldn't just be a constant in the expression for $\mathbb E[Y|Z]$. – akkarin Jun 10 '17 at 20:25
• Well you need to condition on $G_i$ too -- you have to treat at some dose. I'll write up an answer maybe tomorrow, but basically $\beta$ is the average treatment effect. I made a mistake above sorry, you do need the $Z_i$ interaction with the $G_i$ centred if you want to interpret the coefficients like this. – Will Jun 10 '17 at 22:17

I'd suggest a mixed-effects/multilevel model. Essentially you combine two levels: treatment vs. non-treatment, and within the treatment group.

So the model would be something like:

\begin{align} \mathbb{E}(Y_i) = \beta_0 + X_i\beta + \gamma_i Z_i \end{align}

but we additionally assume that $\gamma_i$ depends on the extent of treatment (which I'll call the dose) -- each individual who was treated now has their own parameter for the treatment effect, which we expect will depend on the dose.

So we have an additional linear model for the $\gamma_i$:

\begin{align} \gamma_i = \varphi_{0} + (G_i-\bar{G_i})\varphi_i + \xi_i. \end{align} with $\xi_i$ some error (probably zero-mean Gaussian).

Now, $\varphi_0$ is the average treatment effect -- i.e. the average benefit you get from the treatment when you're given an average dose. This is probably the most sensible way to compare the treatment vs. non-treatment groups.

If you're happy to assume that $\xi_i = 0$, i.e. there's no variability in the treatment effect at a given dose for a given patient, then you can just sub back in to get:

\begin{align} \mathbb{E}(Y_i) = \beta_0 + X_i\beta + (\varphi_0 + (G_i-\bar{G_i})\varphi_i) Z_i \end{align} which you can fit with any old linear regression routine.

If you want to do the full mixed-effects model, there are lots of packages in R that will do it for you, the most common is lme4.

Little advert for going Bayesian (ignore if you want):

I tend to do Bayesian stats personally and Bayesian methods are quite well suited to multilevel/hierarchical models -- the R package brms is absolutely brilliant, intuitive, very well documented and has very fast support from the developer in a Google group if you're struggling. One of the benefits of the Bayesian approach is that it's very easy to 'pool' your estimates -- you encourage the $\varphi_i$ towards zero if there's not sufficient data to estimate them, i.e. you encourage the individual effects towards their mean values.