Re: [R] split-plot analysis with lme()

From: i.m.s.white <i.m.s.white_at_ed.ac.uk>
Date: Tue 17 Oct 2006 - 14:26:46 GMT

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@ed.ac.uk              *

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Wed Oct 18 01:00:35 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 17 Oct 2006 - 15:30:11 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.