1
$\begingroup$

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   
...
$\endgroup$
  • $\begingroup$ After some more digging, it seems to me that the difflsmeans procedure is roughly equivalent to performing a Fisher's LSD test in a typical ANOVA, which uses pooled variance and the total model degrees of freedom (rather than the N's of two comparison groups) to determine p values. $\endgroup$ – JT12 Oct 5 '16 at 13:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.