[R-sig-ME] lme with varIdent in a split-plot design

[Etienne Laliberté]

Thanks Luca. Yes, you're right that in that case the model is overly
complicated as only "fert" is important.

What puzzled me (and still does) is that I had used varIdent(form = ~ 1
| fert) before with other response variables from this same experimental
design and using all fixed terms (i.e. fert * graz), without problems.

Thanks for your help.

Etienne

Le lundi 19 juillet 2010 à 01:28 -0400, Luca Borger a écrit :
> Hello,
> 
> > Error in MEestimate(lmeSt, grps) :
> >  Singularity in backsolve at level 0, block 1
> >
> > Can someone please explain me what this error means?
> 
> 
> This message usually comes when the fixed effects parts of the model ("level 
> 0") is too complicated, given the data. In fact, as you say, your dataset is 
> "minimally replicated" and the message disappears when you take the 
> interaction out. You can 'force' an estimate by slightly changing the 
> variance structure (if this is adequate to do so, I'm not entirely sure):
> 
> mod1b <- update(mod1, weights = varPower(form = ~ as.numeric(fert)))
> mod1c <- update(mod1, weights = varPower())
> 
> However, given your data, the variance structure does not appear to be 
> supported, and neither is the effect of grazing:
> 
> mod3b <- update(mod1, .~. - fert:graz)
> mod4 <- update(mod3b, .~. - graz)
> mod5 <- update(mod4, .~. - fert)
> 
> 
> > AIC(mod1,mod1b,mod1c,mod3,mod3b,mod4,mod5)
>       df        AIC
> mod1  18  25.365540
> mod1b 19  26.318237
> mod1c 19  27.022812
> mod3  14   8.683218
> mod3b 10   2.501522
> mod4   8  -9.972859
> mod5   4 -10.700400
> >
> 
> 
> HTH and probably others can provide better suggestions.
> 
> 
> Cheers,
> 
> Luca
> 
> 
> 
> 
> ----- Original Message ----- 
> From: "Etienne Laliberté" <etiennelaliberte at gmail.com>
> To: "r-sig-mixed-models" <r-sig-mixed-models at r-project.org>
> Sent: Sunday, July 18, 2010 7:34 PM
> Subject: [R-sig-ME] lme with varIdent in a split-plot design
> 
> 
> > I'm using lme() to analyse soil pH data from a simple, balanced,
> > minimally replicated agricultural split-plot experiment with 2 block
> > (replicates), each with five whole plots (corresponding to 5 fertilizer
> > levels). Within each of the 10 whole plots are 3 subplots corresponding
> > to 3 grazing intensity treatments, for a total of 30 subplots.
> >
> > There was evidence for heterogeneity when I plotted residuals against
> > fertilizer levels, so I tried using varIdent(form = ~1 | fert). However,
> > I get this error message:
> >
> > Error in MEestimate(lmeSt, grps) :
> >  Singularity in backsolve at level 0, block 1
> >
> > Can someone please explain me what this error means?
> >
> > To make my example reproducible, I'm attaching the R code (with the real
> > data and design):
> >
> > # get pH dataset
> > pH <- c(5.0, 4.9, 4.8, 5.0, 4.9, 4.9, 5.5, 5.5, 5.6, 5.2, 5.0, 5.2, 4.9,
> > 4.9, 5.0, 5.1, 5.2, 5.0, 4.9, 5.1, 4.9, 5.1, 5.0, 5.2, 5.6, 5.4, 5.2,
> > 4.6, 4.7, 4.8)
> > block <- gl(2, 15)
> > fert <- factor(c(500, 500, 500, 250, 250, 250, 0, 0, 0, 50, 50, 50, 100,
> > 100, 100, 500, 500, 500, 100, 100, 100, 50, 50, 50, 0, 0, 0, 250, 250,
> > 250 ))
> > graz <- factor(c("mod", "lax", "hard", "lax", "mod", "hard", "hard",
> > "mod", "lax", "lax", "mod", "hard", "hard", "mod", "lax", "lax", "mod",
> > "hard", "lax", "hard", "mod", "lax", "mod", "hard", "mod", "lax",
> > "hard", "lax", "mod", "hard"))
> > pHdata <- data.frame(pH, block, fert, graz)
> >
> > # run model
> > library(nlme)
> > mod1 <- lme(pH ~ fert * graz, random = ~ 1 | block / fert, data =
> > pHdata)
> > anova(mod1)
> >
> > # get residuals and fitted values
> > mod1.res <- residuals(mod1)
> > mod1.fit <- fitted(mod1)
> >
> > # check residuals against fitted
> > plot(mod1.res ~ mod1.fit)
> >
> > # check residuals against fert
> > plot(mod1.res ~ fert) # heterogeneity
> >
> > # update model with different variances for fert
> > mod2 <- update(mod1, weights = varIdent(form = ~ 1 | fert) ) # returns
> > error message I don't understand
> >
> > # try something else: remove interaction and use different variances for
> > fert
> > mod3 <- update(mod1, .~. - fert:graz, weights = varIdent(form = ~ 1 |
> > fert) ) # this works
> >
> >
> >
> > Many thanks
> >
> > Etienne
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> > 
> 


-- 
Etienne Laliberté
================================
School of Forestry
University of Canterbury
Private Bag 4800
Christchurch 8140, New Zealand
Phone: +64 3 366 7001 ext. 8365
Fax: +64 3 364 2124
www.elaliberte.info