R help - Oct 2006 - split-plot analysis with lme()

Dear R-help,

Why can't lme cope with an incomplete whole plot when analysing a split-plot
experiment? For example:

R : Copyright 2006, The R Foundation for Statistical Computing
Version 2.3.1 (2006-06-01)
> library(nlme)
> attach(Oats)
> nitro <- ordered(nitro)
> fit <- lme(yield ~ Variety*nitro, random=~1|Block/Variety)
> anova(fit)
              numDF denDF   F-value p-value
(Intercept)       1    45 245.14333  <.0001
Variety           2    10   1.48534  0.2724
nitro             3    45  37.68560  <.0001
Variety:nitro     6    45   0.30282  0.9322

# Excellent! However ---
> fit2 <- lme(yield ~ Variety*nitro, random=~1|Block/Variety, subset+
!(Variety == "Golden Rain" & nitro == "0"))
Error in MEEM(object, conLin, control$niterEM) : 
	Singularity in backsolve at level 0, block 1
-- 
************************************************
*    I.White                                   *
*    University of Edinburgh                   *
*    Ashworth Laboratories, West Mains Road    *
*    Edinburgh EH9 3JT                         *
*    Fax: 0131 650 6564   Tel: 0131 650 5490   *
*    E-mail: i.m.s.white at ed.ac.uk              *
************************************************

----- End forwarded message -----

-- 
************************************************
*    I.White                                   *
*    University of Edinburgh                   *
*    Ashworth Laboratories, West Mains Road    *
*    Edinburgh EH9 3JT                         *
*    Fax: 0131 650 6564   Tel: 0131 650 5490   *
*    E-mail: i.m.s.white at ed.ac.uk              *

The problem in your example is that 'lme' doesn't know how to 
handle the Variety*nitro interaction when all 12 combinations are not 
present.  The error message "singularity in backsolve" means that with
data for only 11 combinations, which is what you have in your example, 
you can only estimate 11 linearly independent fixed-effect coefficients, 
not the 12 required by this model:  1 for intercept + (3-1) for Variety 
+ (4-1) for nitro + (3-1)*(4-1) for Variety*nitro = 12. 

      Since 'nitro' is a fixed effect only, you can get what you want by
keeping it as a numeric factor and manually specifying the (at most 5, 
not 6) interaction contrasts  you want, something like the following: 

fit2. <- lme(yield ~ Variety+nitro+I(nitro^2)+I(nitro^3)
             +Variety:(nitro+I(nitro^2)), data=Oats,
             random=~1|Block/Variety,
            subset=!(Variety == "Golden Rain" & nitro ==
"0"))

      NOTE:  This gives us 4 degrees of freedom for the interaction.  
With all the data, we can estimate 6.  Therefore, there should be some 
way to get 5, but so far I haven't figured out an easy way to do that.  
Perhaps someone else will enlighten us both. 

      Even without a method for estimating an interaction term with 5 
degrees of freedom, I hope I've at least answered your basic question. 

      Best Wishes,
      Spencer Graves          

i.m.s.white wrote:
> Dear R-help,
>
> Why can't lme cope with an incomplete whole plot when analysing a
split-plot
> experiment? For example:
>
> R : Copyright 2006, The R Foundation for Statistical Computing
> Version 2.3.1 (2006-06-01)
>
>   
>> library(nlme)
>> attach(Oats)
>> nitro <- ordered(nitro)
>> fit <- lme(yield ~ Variety*nitro, random=~1|Block/Variety)
>> anova(fit)
>>     
>               numDF denDF   F-value p-value
> (Intercept)       1    45 245.14333  <.0001
> Variety           2    10   1.48534  0.2724
> nitro             3    45  37.68560  <.0001
> Variety:nitro     6    45   0.30282  0.9322
>
> # Excellent! However ---
>
>   
>> fit2 <- lme(yield ~ Variety*nitro, random=~1|Block/Variety,
subset>>
> + !(Variety == "Golden Rain" & nitro == "0"))
> Error in MEEM(object, conLin, control$niterEM) : 
> 	Singularity in backsolve at level 0, block 1
>