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?
Best Answer
By default in R, treatment contrasts are used for factors. This means that what you get in the output from
summary(mod)
are the differences from the reference level fortreatment
. E.g., 37.4 is the difference between treatment B and treatment A.If you want to get the mean for treatment B, you will need to add the coefficients. For the standard errors, you also need to account for the covariance between the estimates of the fixed effects. The following code illustrates how this is done (which essentially what effects and emmeans do under the hood):