Standard errors in linear mixed models

Question

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?

Answer 1

By default in R, treatment contrasts are used for factor. This means that what you get in the output from summary(mod) are the differences from the reference level for treatment. 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):

coefs <- fixef(mod)
V <- vcov(mod)

# mean and std. error for treatment B
DF <- data.frame(treatment = factor("B", levels = LETTERS[1:3]))
X <- model.matrix(~ treatment, data = DF)
c(X %*% coefs)
sqrt(diag(X %*% V %*% t(X)))


# mean and std. error for treatment C
DF <- data.frame(treatment = factor("C", levels = LETTERS[1:3]))
X <- model.matrix(~ treatment, data = DF)
c(X %*% coefs)
sqrt(diag(X %*% V %*% t(X)))