Suppose I have a dataset with a continuous outcome, a dichotomous variable representing treatment, and another dichotomous variable representing membership in some group. As an example we can use R's built-in CO2 dataset, taking uptake as the outcome, predicted by Treatment (chilled vs not chilled) and plant Type (Quebec vs. Mississippi). My research questions are:

  • Q1: Does treatment affect uptake for Quebec-type plants?
  • Q2: Does treatment affect uptake for Mississippi-type plants?

I could run two separate regressions estimating the effect of treatment within each subset of the data. But suppose that I expect the two types to have similar untreated outcomes, and similar responses to other covariates that might be in the model (not interacting with treatment or type, just soaking up additional variance). Therefore, I'd expect to have higher power to detect the two treatment effects if I could use all my data in a single model.

I could do a multiple regression interacting treatment with type:

summary(lm(uptake ~ Treatment * Type, data=CO2))

# Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                        35.333      1.747  20.225  < 2e-16 ***
# Treatmentchilled                   -3.581      2.471  -1.449 0.151141    
# TypeMississippi                    -9.381      2.471  -3.797 0.000284 ***
# Treatmentchilled:TypeMississippi   -6.557      3.494  -1.877 0.064213 .  
# F (3,80) = 23.82, p=4.1e-11

The coefficients and t-tests answer a different set of questions:

  • A1: Does treatment affect uptake for Quebec-type plants? (same as my question 1)
  • A2: Does treatment affect uptake for Mississippi-type plants any more or less than treatment affects Quebec-type plants? (not the same as my question 2)

Let's say I don't care about A2 at all and I only want to know if the effect on Mississippi-type plants is significantly different from zero. I could run the same model but using Mississippi as the reference class:

summary(lm(uptake ~ Treatment * (Type=="Quebec"), data=CO2))
# (Intercept)                             25.952      1.747  14.855  < 2e-16 ***
# Treatmentchilled                       -10.138      2.471  -4.103 9.74e-05 ***
# Type == "Quebec"TRUE                     9.381      2.471   3.797 0.000284 ***
# Treatmentchilled:Type == "Quebec"TRUE    6.557      3.494   1.877 0.064213 .  

The Treatmentchilled parameter is precisely what you'd have calculated for Mississippi-type plants using the previous model's Treatmentchilled + Treatmentchilled:TypeMississippi parameters, but now it has a simple t-test telling me whether it's statistically significant.

Is there anything wrong with running a model twice with different reference classes in order to test the treatment effect within each class? Specifically:

  1. Does it introduce any multiple-comparison or interpretation issues that are different from what I'd have to deal with if I ran only the first model and were interested in the treatment and interaction effects?
  2. Is there a more elegant way to accomplish what I want?

1 Answer 1


Question 1: The model is fundamentally the same regardless of your choice of reference level. Although the coefficient values reported by summary() will differ between the two choices of reference, predictions made for any particular scenario would be identical. There are no additional multiple-comparison issues.

Question 2: There are better ways to proceed. Trying to work solely from the model summary() is likely to lead to confusion when there are interactions. Use post-modeling tools instead.

For example, the Anova() function in the R car package can evaluate overall statistical significance of predictors in a way that incorporates all interaction terms appropriately, even if the design is imbalanced. The emmeans package provides tools for comparing scenarios of interest from many types of models, with appropriate corrections for multiple comparisons.

  • 1
    $\begingroup$ Thanks, very helpful! For future readers, the command that does what I specifically wanted is pairs(emmeans::emmeans(lm_model, ~ Treatment * Type), by="Type") where lm_model is the result from lm() $\endgroup$
    – octern
    May 12, 2023 at 18:56

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.