I'm trying to understand how standard errors for the parameter estimates are calculated in linear mixed models, and why I don't get the same output with different methods. I've made the following example for a simple linear mixed model using package lme4:
library("lme4")
library("lmerTest")
library("effect")
library("emmeans")
response <- c(33,85,77,43,93,87,24,81,65,56,74,96,47,57,94)
ind <- c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5)
treatment <- c("A","B","C","A","B","C","A","B","C","A","B","C","A","B","C")
df <- data.frame(response, ind, treatment)
mod <- lmer(response ~ treatment + (1 | ind), data = df)
summary(mod)
as.data.frame(effect("treatment", mod))
emmeans(mod, spec = c("treatment"))
summary(mod)
produces the following output, where we get the standard errors (for the fixed effects):
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: response ~ treatment + (1 | ind)
Data: df
REML criterion at convergence: 100.4
Scaled residuals:
Min 1Q Median 3Q Max
-1.6290 -0.5492 0.2168 0.6793 1.1625
Random effects:
Groups Name Variance Std.Dev.
ind (Intercept) 3.551 1.884
Residual 164.783 12.837
Number of obs: 15, groups: ind, 5
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 40.600 5.802 11.989 6.997 1.45e-05 ***
treatmentB 37.400 8.119 8.000 4.607 0.00174 **
treatmentC 43.200 8.119 8.000 5.321 0.00071 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) trtmnB
treatmentB -0.700
treatmentC -0.700 0.500
We can also get standard errors (and confidence intervals) from e.g. the effects and emmeans packages (which produce the same output), and for as.data.frame(effect("treatment", mod))
it looks like this:
treatment fit se lower upper
1 A 40.6 5.802299 27.95788 53.24212
2 B 78.0 5.802299 65.35788 90.64212
3 C 83.8 5.802299 71.15788 96.44212
The Estimate/fit produces identical values (just with the difference that they are already summed in the effect("treatment", mod)
output). For the standard errors, we get the same value for intercept/treatment A (5.80), but different values for treatment B and C (8.12 and 5.80). I'm not too familiar with the details of mixed models, and I might miss something obvious here, but I don't understand why this is the case. My questions are (1) how are standard errors for the parameters calculated in linear mixed models, and (2) why does summary(mod)
and effect("treatment", mod)
give different values, and (3) which one would be more "correct" to report?