3
$\begingroup$

I am using the emmeans package to run post-hoc analysis on linear mixed models. The results provide what I would expect except for the standard error. I run the example code provided in the documentation:

warp.lm <- lm(breaks ~ wool*tension, data = warpbreaks)
warp.emm <- emmeans(warp.lm, ~ tension | wool)
warp.emm
# wool = A:
#   tension   emmean       SE df lower.CL upper.CL
# L       44.55556 3.646761 48 37.22325 51.88786
# M       24.00000 3.646761 48 16.66769 31.33231
# H       24.55556 3.646761 48 17.22325 31.88786
# 
# wool = B:
#   tension   emmean       SE df lower.CL upper.CL
# L       28.22222 3.646761 48 20.88992 35.55453
# M       28.77778 3.646761 48 21.44547 36.11008
# H       18.77778 3.646761 48 11.44547 26.11008
# 
# Confidence level used: 0.95

I can't understand how the SE values are all the same, when the means (and everything else) are completely different? What do these SE values actually represent? Thank you

$\endgroup$
2
  • $\begingroup$ Look at vignette(“FAQs”). emmeans() summarizes am model, not its underlying data. If you fit a model based on an underlying assumption of equal variances, and the design is balanced, then the SEs will be equal because the model assumes that to be true. $\endgroup$
    – Russ Lenth
    Oct 1, 2018 at 20:38
  • $\begingroup$ @rvl But the question remains why this standard error equals the standard error of the regression model intercept: all.equal(summary(warp.emm)$SE[1], summary(warp.lm)$coef["(Intercept)", "Std. Error"]) $\endgroup$
    – Roland
    Oct 2, 2018 at 11:51

1 Answer 1

6
$\begingroup$

OK, let us dissect this model. First, the model itself:

> getOption("contrasts")
        unordered           ordered 
"contr.treatment"      "contr.poly" 

> warp.lm <- lm(breaks ~ wool*tension, data = warpbreaks)
> summary(warp.lm)
... (some output omitted) ...
Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)      44.556      3.647  12.218 2.43e-16
woolB           -16.333      5.157  -3.167 0.002677
tensionM        -20.556      5.157  -3.986 0.000228
tensionH        -20.000      5.157  -3.878 0.000320
woolB:tensionM   21.111      7.294   2.895 0.005698
woolB:tensionH   10.556      7.294   1.447 0.154327

Residual standard error: 10.94 on 48 degrees of freedom
Multiple R-squared:  0.3778,    Adjusted R-squared:  0.3129 
F-statistic: 5.828 on 5 and 48 DF,  p-value: 0.0002772

This model has an underlying assumption that the error SD is homogeneous, and its estimated value is 10.94. With the default contrast coding ("contr.treatment"),

  1. The intercept is an estimate of the cell mean when each factor is at its first level
  2. The main-effect coefficients (for a model with interaction) are estimates of certain comparisons between cell means. In particular, they are comparisons of cell means where one factor is held constant while the other one changes.
  3. The interaction coefficients are estimates of certain interaction contrasts (namely, differences of differences)

We can observe these results in the output from emmeans() and its relatives. For (1), note that the first result below matches the intercept, in both the estimate and the standard error:

> emmeans(warp.lm, ~ wool * tension)
 wool tension   emmean       SE df lower.CL upper.CL
 A    L       44.55556 3.646761 48 37.22325 51.88786
 B    L       28.22222 3.646761 48 20.88992 35.55453
 A    M       24.00000 3.646761 48 16.66769 31.33231
 B    M       28.77778 3.646761 48 21.44547 36.11008
 A    H       24.55556 3.646761 48 17.22325 31.88786
 B    H       18.77778 3.646761 48 11.44547 26.11008

For (2), the first, second, and fourth results below match the model summary in both estimate (with signs reversed) and standard error:

> pairs(emmeans(warp.lm, ~ wool*tension))
 contrast    estimate       SE df t.ratio p.value
 A,L - B,L 16.3333333 5.157299 48   3.167  0.0302
 A,L - A,M 20.5555556 5.157299 48   3.986  0.0030
 A,L - B,M 15.7777778 5.157299 48   3.059  0.0398
 A,L - A,H 20.0000000 5.157299 48   3.878  0.0041
 A,L - B,H 25.7777778 5.157299 48   4.998  0.0001
 B,L - A,M  4.2222222 5.157299 48   0.819  0.9627
 B,L - B,M -0.5555556 5.157299 48  -0.108  1.0000
 B,L - A,H  3.6666667 5.157299 48   0.711  0.9797
 B,L - B,H  9.4444444 5.157299 48   1.831  0.4561
 A,M - B,M -4.7777778 5.157299 48  -0.926  0.9377
 A,M - A,H -0.5555556 5.157299 48  -0.108  1.0000
 A,M - B,H  5.2222222 5.157299 48   1.013  0.9115
 B,M - A,H  4.2222222 5.157299 48   0.819  0.9627
 B,M - B,H 10.0000000 5.157299 48   1.939  0.3919
 A,H - B,H  5.7777778 5.157299 48   1.120  0.8706

For (3), the first two results below match the model summary in both estimate and standard error.

> contrast(emmeans(warp.lm, ~ wool*tension), interaction = "pairwise")
 wool_pairwise tension_pairwise  estimate       SE df t.ratio p.value
 A - B         L - M             21.11111 7.293523 48   2.895  0.0057
 A - B         L - H             10.55556 7.293523 48   1.447  0.1543
 A - B         M - H            -10.55556 7.293523 48  -1.447  0.1543

Notice that within each set of results above, the standard errors stay constant within a table. That is a consequence of the fact that the design is balanced (equal numers of observations in each cell) and the model assumption of a constant error SD. If you fit a different model (using, say, nlme::gls() that allows for nonhomogeneous error variances, then you will get unequal SEs in both the model summary and the emmeans results.

$\endgroup$

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.