# SE for estimated marginal means

subject <- factor(rep(c(1,2,3,4,5,6,7,8),each=4))
dep <- c(5,4,9,3,4,4,2,1,10,7,8,7,1,2,1,1,5,10,1,7,3,2,1,4,3,8,7,3,1,1,2,1,1,1,2,1,15,10,20,11,2,2,1,3,11,12,9,7,2,3,1,2,11,9,8,9,3,4,2,1,14,20,11,21)

f1 <- factor(rep(c("Female","Male"), each=32))
f2 <- factor(rep(c("day1","day2"),times=32))

data <- data.frame(sub=subject, dep=dep, f1=f1, f2=f2)

m <- lmer(dep ~ f1*f2 + (1|sub), data=data)


I've been trying without success to calculate standard errors as shown with

library(emmeans)
emmeans(m, ~ f1|f2)


Can somebody help me to understand how to get the standard errors without using the library?

EDIT (unbalanced)

subject <- factor(rep(c(1,2,3,4,5,6,7),each=4, times=2))
dep <- c(5,4,9,3,4,4,2,1,10,7,8,7,1,2,1,1,5,10,1,7,3,2,1,4,3,8,7,3,1,1,2,1,15,10,20,11,2,2,1,3,11,12,9,7,2,3,1,2,11,9,8,9,3,4,2,1)

f1 <- factor(rep(c(rep("Female",times=16),rep("Male",times=12)), times=2))
f2 <- factor(c("day1","day1","day1","day1", rep(c("day1","day2"),times=24), "day2","day2","day2","day2"))

data <- data.frame(sub=subject, dep=dep, f1=f1, f2=f2)

m <- lmer(dep ~ f1*f2 + (1|sub), data=data)


The following code illustrates how this computation is done:

grid <- with(data, expand.grid(f1 = levels(f1), f2 = levels(f2)))
X <- model.matrix(~ f1 * f2, data = grid)

V <- vcov(m)
betas <- fixef(m)
grid$$emmean <- c(X %*% betas) grid$$SE <- sqrt(diag(X %*% V %*% t(X)))
grid

• Thanks very much for your answer Dimitris, it's exactly what I was looking for! Just for my understanding, are you also able to calculate the SE with based on predicted values data\$predict <- predict(m)? Or is that not possible? Nov 2, 2018 at 16:39
• Yes, you could also calculate in the same manner SEs for predicted values. However, one thing you need to be aware of is the use of the model.matrix(), when in the formula of the model you have used functions that depend on the data (e.g., poly(), splines::ns(), splines::bs(), etc.). In this case, you need to ensure that X is correctly calculated. Nov 2, 2018 at 16:57
• Thanks Dimitris. Could you please show me how you would do it with the example data from above? I'm not really sure how I can calculate SE with predicted values. I tried using function(x) sd(x)/sqrt(length(x)) but the SE are not the same as the ones that emmeans provide... Nov 2, 2018 at 17:41

Look at the model summary:

> m
Linear mixed model fit by REML ['lmerMod']
Formula: dep ~ f1 * f2 + (1 | sub)
Data: data
REML criterion at convergence: 371.7578
Random effects:
Groups   Name        Std.Dev.
sub      (Intercept) 1.222
Residual             4.768
Number of obs: 64, groups:  sub, 8
Fixed Effects:
(Intercept)         f1Male         f2day2  f1Male:f2day2
3.9375         3.1250         0.1250         0.0625


This gives estimates of the subject SD (1.222) and the error SD (4.768). This is a balanced experiment, and each mean consists of 16 observations. However, there are only 8 different subject effects. The random component of each cell mean includes the average of all 8 subject effects, and the average of 16 of the residual effects. So its SD is estimated as:

> sqrt(1.222^2 / 8 + 4.768^2 / 16)
[1] 1.267882


This is the same (except for slight error due to roundoff) as the SE displayed in the emmeans results:

> emmeans(m, ~ f1*f2)
f1     f2   emmean       SE    df lower.CL upper.CL
Female day1 3.9375 1.267923 40.78 1.376463 6.498537
Male   day1 7.0625 1.267923 40.78 4.501463 9.623537
Female day2 4.0625 1.267923 40.78 1.501463 6.623537
Male   day2 7.2500 1.267923 40.78 4.688963 9.81103

• Thanks for your answer rvl, that makes sense. But if the data is unbalanced (as my real data set is), how do you calculate the SE for each of the cell means? emmeans gives me a different SE value for each mean. I edited my question with an example of an unbalanced data set. Nov 4, 2018 at 18:43
• If it’s unbalanced, then you don’t do it by hand, you use Dimitris’s answer. Nov 4, 2018 at 18:46
• BTW, what’s the problem with using the emmeans package? It already computes what you want. Why go to extra trouble? Nov 4, 2018 at 19:50
• there is nothing wrong with the package, I use it all the time. I just want to have a better understanding how the SE is calculated Nov 4, 2018 at 21:46
• Well, the thing to know then is that the V in Dimitris's answer is the inverse of a matrix, and for an unbalanced mixed model, it is a complicated one. That V is what stands in conceptually for sigma^2/n that comes up in the formula for the SE of a mean in a simple one-sample case. Nov 4, 2018 at 21:53