3
$\begingroup$

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:

library(estimatr)
library(lmtest)
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 - 
    1
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?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
+50
$\begingroup$

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.

Improve this answer
$\endgroup$
15
  • 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
    Commented Aug 19, 2019 at 19:13
  • 1
    $\begingroup$ I edited my answer. $\endgroup$
    – Noah
    Commented Aug 19, 2019 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
    Commented Aug 22, 2019 at 16:56
  • 1
    $\begingroup$ waldtest does the same thing as anova for linear models as long as you set test = "F", so that is a valid method. waldtest is compatible with lm_robust models, too, if you feel comfortable using them. Otherwise, with lm models, you can supply waldtest with a function to extract robust standard errors. waldtest(fit_1, fit_2, type = "F", vcov = function(x) sandwich::vcovHC(x, type = "HC2") with lm objects will give the same results as waldtest on lm_robust objects. $\endgroup$
    – Noah
    Commented Jun 11, 2020 at 19:02
  • 1
    $\begingroup$ The p-value on the test indicates whether there is evidence that the coefficients that differ between the full and reduced models are jointly equal to zero, indicating that that the moderator has neither a main effect nor an interaction with the treatment on the outcome. Rejecting the null hypothesis means that the moderator does play some role. You can also think of it as testing whether the change in $R^2$ from the reduced to full model is different from zero. $\endgroup$
    – Noah
    Commented Jun 11, 2020 at 19:07
4
$\begingroup$

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!

Improve this answer
$\endgroup$
1
  • $\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
    Commented Aug 13, 2019 at 23:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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