Let's say I have a 2x2 design where participants are either in condition A or condition B and, within each condition, either get exposed to exposure C or exposure D.

First, I want to test whether exposure D > exposure C, within each condition. I can do one of two things.

First, I can subset the data to the appropriate condition and check the coefficient on a dummy variable set to 1 for exposure_d.

lm(outcome ~ exposure_d, data = subset(data, condition == "A")
lm(outcome ~ exposure_d, data = subset(data, condition == "B")

Or I can specify a full model and run a contrast

model <- lm(outcome ~ exposure_d * condition_a, data = data)
means = emmeans(model, ~ condition_a | exposure_d)
contrast_table = contrast(means, method = "pairwise")

Are these methods formally equivalent? I tested the two on simulated data and they seem to produce identical coefficients and standard errors, but I worry things might change when I start adding in more control variables (especially since the contrast() function throws the warning that the cell means and standard errors may not be accurate in the presence of interactions). Which of these methods is preferred?


They are not equivalent, because when you leave everything together, you use all the data to estimate the error SD.

Another way to look at it is that the EMMs are based on the model you fitted, and when you fit a different model, you get different results. When you keep all the data together in the one model, the one you show has an underlying assumption that the variance is homogeneous. If you separate the data according to one factor, those two models together could be viewed as one model that allows separate error variances for the two levels of that factor. Yet another model would be one where the interaction is excluded. With these different models, you will get the same estimates (if the data are balanced), but the SEs will differ.

It is typically not a good idea to subset the data. You get more power (and/or lower SEs) by using more data to form estimates. The exception might be models that allow various types of heterogeneity. Such models have more parameters and may not perform better.

  • $\begingroup$ Gotcha, thanks Russ! I was just concerned about the warning I often get that "NOTE: Results may be misleading due to involvement in interactions." Should I be worried about this? $\endgroup$ – Parseltongue Apr 14 at 18:19
  • $\begingroup$ Absolutely you should! It calls attention to the fact that if your two factors interact very much, you have no business computing marginal means because they don't make sense. So do emmip(model, exposure_d ~ condition_a). Subjectively, does it make sense to compute marginal means? $\endgroup$ – Russ Lenth Apr 14 at 18:29
  • $\begingroup$ I see, so it's just a warning to be cautious in the presence of interactions. In my case, I'm experimentally manipulating both exposure and condition, so I think it makes sense to compute marginal means for each cell in the 2x2. One more follow-up -- I just tested on simulated data and discovered that both the subset approach and the full model approach yield identical results with the only difference being degrees of freedom. This could perhaps be due to the precision outputted by emmeans. Is there a way to output more precise results after contrast? $\endgroup$ – Parseltongue Apr 14 at 18:47
  • $\begingroup$ Do the plot. Just because you experimentally manipulated both factors does not mean that their effects are not interactive. "Interaction" does not refer to physical interaction, it's the relationship of their effects pn the response variable. $\endgroup$ – Russ Lenth Apr 14 at 18:52
  • $\begingroup$ Apologies for being obtuse, and I appreciate your patience. What do you mean by "do the plot" -- are you referring to a basic interaction plot, or plots within the emmeans package? $\endgroup$ – Parseltongue Apr 15 at 5:10

This may be an example showing why subsets of the same data give different results and why the assumption of equal variance across covariates matters.

Here, group A has much higher dispersion compared to the other groups:



reps <- c(10, 10, 10, 10, 10, 10)
data <- data.table(
    x= rep(LETTERS[1:length(reps)], times= reps),
    y= rnorm(n= sum(reps), mean= rep(20, times= sum(reps)), 
        sd= rep(c(5, 1, 1, 1, 1, 1), times= reps))

ggplot(data= data, aes(x= x, y= y)) +
    geom_quasirandom(width = 0.25)

enter image description here

When we fit the model including all groups, the mean of group A has an unrealistically small error:

fit_all <- summary(lm(y ~  0 + x, data= data))

   Estimate Std. Error t value Pr(>|t|)    
xA  18.0842     0.6976   25.92   <2e-16 ***
xB  19.8818     0.6976   28.50   <2e-16 ***
xC  19.6121     0.6976   28.11   <2e-16 ***
xD  19.2338     0.6976   27.57   <2e-16 ***
xE  19.3902     0.6976   27.80   <2e-16 ***
xF  19.7211     0.6976   28.27   <2e-16 ***

This is because the variance of the estimates $\beta$ is:

$$ var(\hat{\beta}) = \hat{\sigma}^{2}(\mathbf{X'X})^{-1} $$

where $\mathbf{X}$ is the design matrix and $\hat{\sigma}$ is the mean squared error, so it is common to all factor levels

Fitting the same model using only group A and B makes the errors larger but using all but group A the errors are smaller:

fit_ab <- summary(lm(y ~ 0 + x, data= data[x %in% c('A', 'B')]))

   Estimate Std. Error t value Pr(>|t|)    
xA   18.084      1.139   15.88 4.93e-12 ***
xB   19.882      1.139   17.46 9.89e-13 ***

fit_woa <- summary(lm(y ~ 0 + x, data= data[x != 'A']))

   Estimate Std. Error t value Pr(>|t|)    
xB  19.8818     0.2969   66.96   <2e-16 ***
xC  19.6121     0.2969   66.05   <2e-16 ***
xD  19.2338     0.2969   64.77   <2e-16 ***
xE  19.3902     0.2969   65.30   <2e-16 ***
xF  19.7211     0.2969   66.42   <2e-16 ***
  • $\begingroup$ Thanks! This is helpful. $\endgroup$ – Parseltongue Apr 15 at 5:13
  • $\begingroup$ In the case of unequal variances across groups (as in your simulated example), would you run two separate models (subsets) or a single model that you run contrasts on? $\endgroup$ – Parseltongue Apr 16 at 19:49

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.