# How to "trick" R into producing treatment effects in a treatment * dummy model?

I'm running some basic regressions which can be specified compactly as the formula: y~treatment*dummy.

Say there are several ($m$) treatments, $T_1,\ldots,T_m$ (with $T_1$ being the reference/control); the dummy is also multifaceted (categorical), taking $n$ values $D_1,\ldots,D_n$

Then (suppressing the observation index and error) the above formula specification basically returns the specified formula as

$$y=T^TBD$$

Where $T$ is the $m\times 1$ vector $[1, T_2, \ldots, T_m]$ of treatment indicators (excluding the reference treatment), $B$ is the $m\times n$ matrix of coefficients $\{\beta_{i,j}\}_{i=1,j=1}^{m\quad n}$, and $D$ is the $n\times 1$ vector $[1, D_2,\ldots, D_n]$ of dummy indicators (excluding the reference category).

This is all well and good, but the resulting coefficients in $B$ don't really have any clean interpretation, especially for my application. In particular, I'm looking for significant treatment effects--consider trying to answer the following: was Treatment 5 significantly better among individuals in category 3?

In the above specification, we'd be examining $\mathbb{E}[y|T_5,D_3]-\mathbb{E}[y|T_1,D_3]=\beta_{5,0}+\beta_{5,3}$, so we could add the coefficients we get out and use, e.g., a Wald test to determine significance.

However, consider the equivalent specification (I'm 100% sure someone besides has written it this way before since it took me all of 20 minutes to come up with):

$$y = \delta_0 + \sum_{j=2}^n\beta_j D_j + \sum_{i=2}^m \sum_{j=1}^n \gamma_{i,j}T_iD_j$$

Now the treatment effect is $\mathbb{E}[y|T_5,D_3]-\mathbb{E}[y|T_1,D_3]=\gamma_{5,3}$.

So this latter formulation has the convenient property that we can read our treatment effects right off our regression summary (especially including standard errors); its major drawback is that there's no way to supply this as an R formula parsimoniously, or at least I can't see a way to.

Does anyone have any experience with some secret formula or package to deal with this (I imagine exceedingly common) specification?

I have to struggle a bit with the mathematical notation, if I but in the following formula that you provided:

$$y = \delta_0 + \sum_{j=1}^n\beta_j D_j + \sum_{i=1}^m \sum_{j=0}^n \gamma_{i,j}T_iD_j$$

Isn't that essentially saying that there should be one coefficient per dummy variable (excluding the first), and then one coefficient per combination of T and D? In that case, you could simply write the formula as:

y ~ group + treatment:group


But when including interaction terms, it is often recommended that you should include all the main terms that are included in the interaction terms:

y ~ group * treatment

• Nevertheless, my answer remains the same. Using R notation, you can do it like group:treatment or group + group:treatment or group*treatment, you will still get the same number of categories. When including interaction terms, the usual recommendation in textbooks is to also include the main effects, which in this case corresponds to y ~ group * treatment.
– JonB
Commented Sep 16, 2015 at 12:15
• actually y ~ group + treatment:group does not exclude the reference level of group, so your answer is correct. I just get the feeling you're not understanding why I disagree with the "textbook recommendation" Commented Sep 16, 2015 at 15:35
• I've been thinking about it some more, and I think that the group + group:treatment approach makes sense, so go for it!
– JonB
Commented Sep 16, 2015 at 17:40
• By the way, check here for a simpler way to implement this in R ;-) Commented Oct 3, 2015 at 21:12