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

[Luca Borger]

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
>