Is it preferable to subset data to test specific hypotheses or specify a full model and run contrasts?

Question

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

library(emmeans)
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?

Answer 1

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.

Answer 2

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:

library(ggplot2)
library(ggbeeswarm)
library(data.table)

set.seed(1234)

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)

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))

Coefficients:
   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')]))

Coefficients:
   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']))

Coefficients:
   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 ***