[R] Interpretation of hypothesis tests for mixed models

[Douglas Bates]

Olof Leimar <Olof.Leimar at zoologi.su.se> writes:

> My question concerns the logic behind hypothesis tests for fixed-effect
> terms in models fitted with lme. Suppose the levels of Subj indicate a
> grouping structure (k subjects) and Trt is a two-level factor (two
> treatments) for which there are several (n) responses y from each
> treatment and subject combination. If one suspects a subject by
> treatment interaction, either of the following models seem natural
>   
> > fm1 <- lme(y ~ Trt, random = list(Subj = pdDiag(~ Trt)))
> > fm2 <- lme(y ~ trt, random = ~ 1 | Subj/Trt)

I don't think those models are equivalent.  To get equivalent models I
think you need random = list(Subj = pdCompSymm(~ Trt - 1)) in the first
model.

For example

> library(nlme)
Loading required package: nls 
Loading required package: lattice 
> data(Machines)
> fm1 <- lme(score ~ Machine, data = Machines, random = ~ 1 | Worker/Machine)
> logLik(fm1)
`log Lik.' -107.8438 (df=6)
> fm2 <- lme(score ~ Machine, data = Machines, random=list(Worker=pdCompSymm(~Machine - 1)))
> logLik(fm2)
`log Lik.' -107.8438 (df=6)
> summary(fm1)
Linear mixed-effects model fit by REML
 Data: Machines 
       AIC      BIC    logLik
  227.6876 239.2785 -107.8438

Random effects:
 Formula: ~1 | Worker
        (Intercept)
StdDev:    4.781049

 Formula: ~1 | Machine %in% Worker
        (Intercept)  Residual
StdDev:    3.729536 0.9615768

Fixed effects: score ~ Machine 
               Value Std.Error DF   t-value p-value
(Intercept) 52.35556  2.485829 36 21.061606  <.0001
MachineB     7.96667  2.176974 10  3.659514  0.0044
MachineC    13.91667  2.176974 10  6.392665  0.0001
 Correlation: 
         (Intr) MachnB
MachineB -0.438       
MachineC -0.438  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.26958756 -0.54846582 -0.01070588  0.43936575  2.54005852 

Number of Observations: 54
Number of Groups: 
             Worker Machine %in% Worker 
                  6                  18 
> summary(fm2)
Linear mixed-effects model fit by REML
 Data: Machines 
       AIC      BIC    logLik
  227.6876 239.2785 -107.8438

Random effects:
 Formula: ~Machine - 1 | Worker
 Structure: Compound Symmetry
         StdDev    Corr       
MachineA 6.0626364            
MachineB 6.0626364 0.621      
MachineC 6.0626364 0.621 0.621
Residual 0.9615618            

Fixed effects: score ~ Machine 
               Value Std.Error DF   t-value p-value
(Intercept) 52.35556  2.485416 46 21.065106  <.0001
MachineB     7.96667  2.177416 46  3.658770   7e-04
MachineC    13.91667  2.177416 46  6.391367  <.0001
 Correlation: 
         (Intr) MachnB
MachineB -0.438       
MachineC -0.438  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.26962147 -0.54844560 -0.01068651  0.43936587  2.54004419 

Number of Observations: 54
Number of Groups: 6 

This still doesn't get around the problem of the different definitions
of the degrees of freedom in the denominator of the F tests.  That
occurs because the definition of the degrees of freedom is not
invariant with respect to the different formulation of the models.

I prefer to formulate such a model as nested random effects rather
than trying to decide how the model matrices should be defined in a
model with a compound symmetry structure for the variance-covariance
term.  I think the definition of the degrees of freedom is cleaner
in that case too.