[R] split-plot analysis with lme()

[Spencer Graves]

      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
>