5
$\begingroup$ 
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)
Improve this question
$\endgroup$

2 Answers 2

Reset to default
7
$\begingroup$

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
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – locus
    Commented Nov 2, 2018 at 16:39
  • 2
    $\begingroup$ 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. $\endgroup$
    – Dimitris Rizopoulos
    Commented Nov 2, 2018 at 16:57
  • $\begingroup$ 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... $\endgroup$
    – locus
    Commented Nov 2, 2018 at 17:41
3
$\begingroup$

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
Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – locus
    Commented Nov 4, 2018 at 18:43
  • $\begingroup$ If it’s unbalanced, then you don’t do it by hand, you use Dimitris’s answer. $\endgroup$
    – Russ Lenth
    Commented Nov 4, 2018 at 18:46
  • $\begingroup$ BTW, what’s the problem with using the emmeans package? It already computes what you want. Why go to extra trouble? $\endgroup$
    – Russ Lenth
    Commented Nov 4, 2018 at 19:50
  • $\begingroup$ 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 $\endgroup$
    – locus
    Commented Nov 4, 2018 at 21:46
  • $\begingroup$ 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. $\endgroup$
    – Russ Lenth
    Commented Nov 4, 2018 at 21:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.