I have some experimental data in which treatments are nested within sessions which are nested within groups. It's tricky, because it seems like a simple between subjects comparison (TX1 vs TX2 vs TX3), but it's not, because they were grouped into nonindependent sessions within nonindependent groups. For example, subjects 1,2,3 interacted with each other, and subjects 1-6 (group 1) came from the same subject "pool." I'd like to run this as a nested (TX within session within group) model using lme4/lmerTest, check the omnibus F tests sing the anova() command, and then do specific comparisons between treatments and controls. However, when I use the difflsmeans() function in lmerTest, the degrees of freedom for the individual comparisons seems to reflect the df for the entire model, not for the specific (smaller) comparisons. So I have a situation where the df is much higher than the number of individuals in both groups. I've made some randomized dummy data for reproducibility, and below are my commands and the output
Data
SubjectID Group Session TX Sex SessionSize Outcome
1 1 AM TX1 F 2 75
2 1 AM TX2 F 5 61
3 1 AM TX3 M 3 53
4 1 PM TX1 F 5 61
5 1 PM TX2 F 3 92
6 1 PM TX3 F 3 74
7 2 PM TX1 F 3 95
8 2 PM TX2 F 3 86
9 2 AM TX3 M 4 83
10 2 AM TX1 M 2 97
11 2 PM TX2 F 4 79
12 2 PM TX3 M 3 69
13 3 AM TX1 F 2 96
14 3 AM TX2 M 2 93
15 3 PM TX3 M 4 87
16 3 PM TX1 F 5 67
17 3 AM TX2 F 5 99
18 3 AM TX3 F 4 86
19 4 PM TX1 M 5 87
20 4 PM TX2 F 4 54
21 4 AM TX3 F 2 89
22 4 AM TX1 F 4 90
23 4 PM TX2 F 4 88
24 4 PM TX3 F 4 76
25 5 AM TX1 F 5 71
26 5 AM TX2 M 2 91
27 5 PM TX3 M 4 51
28 5 PM TX1 F 3 77
29 5 AM TX2 F 4 53
30 5 AM TX3 M 2 92
31 6 PM TX1 F 2 53
32 6 PM TX2 F 4 83
33 6 AM TX3 F 3 90
34 6 AM TX1 F 3 78
35 6 PM TX2 F 2 92
36 6 PM TX3 F 4 88
37 7 AM TX1 F 3 76
38 7 AM TX2 M 4 50
39 7 PM TX3 M 5 83
40 7 PM TX1 F 4 97
41 7 PM TX2 M 5 87
42 7 PM TX3 F 4 54
Commands/Output
> library(lmerTest)
> fit<-lmer(Outcome~Sex*TX+Session+SessionSize+(1+SessionSize|Group/Session),data=temp,REML=F)
> anova(fit)
Analysis of Variance Table of type III with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
Sex 57.67 57.67 1 38.770 0.32850 0.5699
TX 344.58 172.29 2 33.307 0.98146 0.3853
Session 15.56 15.56 1 12.521 0.08862 0.7708
SessionSize 446.26 446.26 1 18.654 2.54215 0.1276
Sex:TX 381.06 190.53 2 40.311 1.08537 0.3474
...
So far nothing out of the ordinary. From here let's imagine that the Sex*TX interaction is significant, and I'd like to know which treatments are different from each other, so I use difflsmeans()
> difflsmeans(fit,test.effs="Sex:TX")
Differences of LSMEANS:
Estimate Standard Error DF t-value Lower CI Upper CI p-value
...
...
Sex:TX F TX1 - F TX2 -2.2 5.963 36.8 -0.38 -14.32 9.85 0.7
Sex:TX F TX1 - M TX2 -0.8 7.972 34.7 -0.10 -16.97 15.41 0.9
Sex:TX F TX1 - F TX3 -0.5 6.430 34.2 -0.07 -13.55 12.58 0.9
Sex:TX F TX1 - M TX3 3.8 6.419 33.9 0.59 -9.28 16.81 0.6
...
...
Here is where I am having problems with the estimated DF. There are only 12 females in TX1, and 10 females in TX2, so shouldn't the degrees of freedom be closer to 20? Or am I misunderstanding what these DF are? I have the same problem when looking at the individual F tests generated using summary(fit), since TX1 is my comparison group.
> summary(fit)
...
Random effects:
Groups Name Variance Std.Dev. Corr
Session:Group (Intercept) 6.301 2.510
SessionSize 2.614 1.617 -1.00
Group (Intercept) 4.204 2.050
SessionSize 1.395 1.181 -1.00
Residual 175.543 13.249
Number of obs: 42, groups: Session:Group, 14; Group, 7
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 91.1585 8.1845 35.3100 11.138 4.17e-13 ***
SexM 14.4439 10.5401 39.6400 1.370 0.178
TXTX2 2.2383 5.9627 36.8300 0.375 0.710
TXTX3 0.4819 6.4296 34.1800 0.075 0.941
SessionPM -1.4267 4.7926 12.5200 -0.298 0.771
SessionSize -3.5824 2.2468 18.6500 -1.594 0.128
SexM:TXTX2 -15.9039 13.7195 40.6400 -1.159 0.253
SexM:TXTX3 -18.6945 12.9455 40.1600 -1.444 0.156
...
