Let's say I run the following regression specification:

lm(outcome ~ treatment) 

Where treatment is a factor variable that can take four values A, B, C, or D . I suspect that there is this predictor, X that moderates the relationship between the treatment and the outcome. However, I don't really care about comparisons across treatment groups -- what I want to be able to say is that, within a given treatment (e.g. A), people with high values of X have better outcomes then individuals with low values of X.

What comes to mind is some kind of interaction effect like:

lm(outcome ~ treatment*X)

But that will only tell me how the treatment varies at different levels of X, in comparison to some left out reference group (rather than in comparison to the same treatment but low values of X).

What is the right way of going about doing this?

Edit: Ran the nested specification suggested by @cdtip and it resulted in the following output:

             Estimate Std. Error  t value  Pr(>|t|) CI Lower CI Upper  DF
(Intercept)   -1.33113   4.96200 -0.26826 0.7889642 -11.1572  8.49498 118
treatmentC    -5.90437   6.94422 -0.85026 0.3969036 -19.6558  7.84706 118
treatmentB    24.92750   9.27155  2.68860 0.0082136   6.5673 43.28769 118
treatmentA    7.26947    5.73785  1.26693 0.2076739  -4.0930 18.63197 118
treatmentD:X  0.65245    1.27021  0.51366 0.6084530  -1.8629  3.16782 118
treatmentC:X  5.15284    1.64776  3.12717 0.0022231   1.8898  8.41586 118
treatmentB:X -2.65812    1.08513 -2.44959 0.0157709  -4.8070 -0.50927 118
treatmentA:X -0.69822    0.67031 -1.04163 0.2997141  -2.0256  0.62919 118

I'm somewhat unclear how to interpret this. The coefficient on treatmentC:X is 5.15 and significant. But it seems like it's significant relative to the left out condition which is treatmentD, rather than making a within treatment comparison. How am I supposed to interpret these results?

It seems like the / operator is just shorthand for:

lm(outcome ~ treatment + treatment:x)

(so removing the main effect of X).

I have found references to "subgroup analysis" and "heterogeneous treatment effect", but these seem to refer to a singular treatment condition, in which the size of the effect varies by some moderator.

EDIT: In order to test the overall moderating effect of the moderate X I've conducted a Wald Test like so:

m1 <- lm_robust(outcome ~ treatment / X - 1, data = data)
m2 <- lm_robust(outcome ~ treatment - 1, data = data)
waldtest(m1, m2)

Wald test

Model 1: outcome ~ treatment/X - 
Model 2: outcome ~ treatment - 1
  Res.Df Df Chisq Pr(>Chisq)    
1    640                        
2    644 -4  18.9    0.00084 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

How do I interpret this? I assume this is evidence that of an overall moderating effect of X, but how do I interpret 18.9?


It seems like what you want is to avoid reference cell coding. You want a slope and intercept for the effect of X on outcome at each level of treatment. You need to modify @cdtip's code to exclude an intercept, and you should get what you need:

lm(outcome ~ treatment/X - 1)

Now, each coefficient treatmentA:X, etc., is the effect of X on the outcome for those in treatment A. By removing the intercept, you get a slope for each treatment level. To test whether there is an overall moderating effect of X, you can use anova() to compare the above model with a model that excludes X (but is otherwise the same); a significant p-value indicates the presence of moderation.

In response to question edits:

This parameterization essentially creates 4 separate regression models of outcome on X, one for each treatment group, but estimated together in a single model. The "main effects" are the intercepts of each treatment group-specific regression. They are not treatment effects. Treatment effects would be comparisons between these intercepts, which represent the difference in the treatment group means when X is 0. The interaction terms are interpreted simply as the slopes of each treatment group-specific regression. So from this coefficient table, you can describe the slope and intercept of the regression of outcome on X in each treatment group. This parameterization does not allow you to make any comparisons between treatment groups, but it is possible to reparameterize to get the contrasts you want, if any.

| cite | improve this answer | |
  • $\begingroup$ Just to be clear. With this specification, if I see a significant effect on an interaction term that implies that, within that treatment, the moderator was significantly predictive of the outcome? $\endgroup$ – Parseltongue Aug 19 '19 at 19:02
  • $\begingroup$ Also do you have a sense of how to test the overall effect of X with cluster-robust standard errors? I'm using lm_robust for clustering, and it doesn't look like anova is compatible $\endgroup$ – Parseltongue Aug 19 '19 at 19:03
  • 1
    $\begingroup$ Your clarification statement is correct. If you are willing to migrate to the survey package, the functions there are compatible with anova(). $\endgroup$ – Noah Aug 19 '19 at 19:13
  • 1
    $\begingroup$ I edited my answer. $\endgroup$ – Noah Aug 19 '19 at 19:50
  • 1
    $\begingroup$ No. With only categorical predictors, a fully saturated model is just estimating subgroup means. As long as as many coefficients are being estimated as there are means, the fit of the model will be identical, regardless of the parameterization. In R, when you remove an intercept in the presence of a categorical predictor, an extra coefficient is estimated to account for what would have been the reference category, so the same number of coefficients is estimated either way. With continuous predictors, however, your worry is valid. $\endgroup$ – Noah Aug 22 '19 at 16:56

If I'm understanding right, you may be wanting to do an analysis of covariance (ANCOVA). There is a helpful tutorial here (http://r-eco-evo.blogspot.com/2011/08/comparing-two-regression-slopes-by.html).

Specifically, read down to the bottom where comparing slopes and intercepts is discussed. This was achieved through a slightly different version of your code.

lm(outcome ~ treatment/x) # assuming x is a continuous predictor

Because of the nested design, you are essentially only comparing the values of x within their respective treatments.

Hope that helps!

| cite | improve this answer | |
  • $\begingroup$ Thanks so much! This certainly seems like the right direction. Could you take a look at my edit? I'm not sure how to make sense of the resulting values. $\endgroup$ – Parseltongue Aug 13 '19 at 23:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.