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

[Etienne Laliberté]

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