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"
),
- The intercept is an estimate of the cell mean when each factor is at its first level
- 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.
- 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.
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$all.equal(summary(warp.emm)$SE[1], summary(warp.lm)$coef["(Intercept)", "Std. Error"])
$\endgroup$