Re: Strip-split-plot analysis with lme: problem nearly solved

["Paul Quataert" <paul.quataert@freebel.net>]

Dear Splussers,

 

The suggestion of Kevin Wright possibly solved the problem. However there is something wrong with the degrees of freedom. So I still have some doubt.

 

From: "Kevin Wright" <kwright@eskimo.com>

 

> To have A and B crossed within Block you need to have something like
> random=list(Block=pdBlocked(list(~1,pdIdent(~A-1),pdIdent(~B-1))))
 
> I believe this will give you a random intercept for each level of block
> and then random effects for A nested within Block and B nested within Block.

 

In fact, as suggested by Joris De Wolf, the answer on how to analyse a strip-split-plot with lme, can be found in "Mixed-Effect Models in S and S-plus". Specifically, their example Cell Culture Bioassay with Crossed Random Effects in Section 4.2.2 (Patterned Variance-Covariance Matrices for the Random Effects) is a strip-split-plot design.

 

However, according to "classical old-fashioned (?) theory" the denominator degrees of freedom of the F-tests are not OK. Further exploration showed that this generally the case with lme if random effect structures are specified in a non-standard way (as is the case in section 4.2.2). 

 

As a simple illustration the Oats example (the first example of 4.2.2) which is a split-plot experiment.

 

The "classical" analysis with aov gives:

 

Efit <- aov(yield~nitro*Variety+Error(Block/Variety),data="">summary(Efit)

 Error: Block
           Df Sum of Sq Mean Sq F Value Pr(F)
 Residuals  5     15875    3175             
 
 Error: Variety %in% Block
           Df Sum of Sq Mean Sq F Value Pr(F)
   Variety  2      1786     893    1.49 0.272
 Residuals 10      6013     601             
 
 Error: Within
               Df Sum of Sq Mean Sq F Value Pr(F)
         nitro  1     19536   19536     116  0.00
 nitro:Variety  2       168      84       0  0.61
     Residuals 51      8606     169             

The same result (in another format) is obtained by lme (df of tests equal with above)

 

Mfit <- lme(yield~nitro*Variety,data="">anova(Mfit)
               numDF denDF F-value p-value
   (Intercept)     1    51     245  <.0001
         nitro     1    51     116  <.0001
       Variety     2    10       1   0.272
 nitro:Variety     2    51       0   0.610

This is not true with alternative (but equivalent !) models. The degrees of freedom of the denominator of the F-test are different (all other estimates seem to be equal, only the p-values differ).

 

# Compound symmetry model

 

MfitA <- lme(yield~nitro*Variety,data="">

    random=list(Block=pdCompSymm(~Variety-1)))
anova(MfitA)
               numDF denDF F-value p-value
   (Intercept)     1    61     245  <.0001
         nitro     1    61     116  <.0001
       Variety     2    61       1   0.234
 nitro:Variety     2    61       0   0.610

#  pdBlocked (as necessary for a strip-split-plot analysis)

 

MfitB <- lme(yield~nitro*Variety,data="">

    random=list(Block=pdBlocked(list(pdIdent(~1),pdIdent(~Variety-1)))))
anova(MfitB)
               numDF denDF F-value p-value
   (Intercept)     1    61     245  <.0001
         nitro     1    61     116  <.0001
       Variety     2    61       1   0.234
 nitro:Variety     2    61       0   0.610

Can anybody explain why their is a difference or is this simply a bug?

 

Thanks,

Paul Quataert

Institute for Forestry and Game Management

Belgium