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
>
______________________________________________
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.
