[R] split-plot analysis with lme()

[i.m.s.white]

Thanks, that clarifies things. And this gets all 5 interaction degrees
of freedom:

oats <- read.table("testlme.dat", head=T)
# This is a subset of the standard data set with
# the combination Variety=Golden Rain, nitro=0 deleted
oats$nitro <- factor(oats$nitro)
attach(oats)
library(nlme)
M <- model.matrix(~Variety*nitro)
fit <- lme(yield ~ Variety+nitro+M[,7:11], random=~1|Block/Variety)
anova(fit)


On Sun, Oct 15, 2006 at 08:28:43AM -0700, Spencer Graves wrote:
>      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
> >  

-- 
************************************************
*    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              *