[R-sig-ME] Silent conversion problem in lmer - should have said '_convergence_ problem'

[S Ellison]

Subject says it all - my only excuse is the midnight timestamp.
S!

>>> "S Ellison" <S.Ellison at lgc.co.uk> 19/04/2011 00:49 >>>
Using the randomised block design data and code below, I'm seeing a
consistent convergence (problem|feature) in lmer. Essentially, a
variance estimate is set to zero at iteration 1, but the optimisation
then can't dig itself out of the hole. This happens whether Run is set
to random or fixed.

This can apparently be avoided with a better choice of start value
(see
below), but obviously one doesn't always have the luxury of a decent
starting guess .

Now I know that one doesn't willingly use complicated optimisation on
data with two degrees of freedom in a variance estimate, and I also
know
that these variances can be extracted from this data by 2-way anova
and
a bit of arithmetic. But I'd kind of hoped lmer would manage anyway -
lme does, on setting Run to a fixed effect - and it's unsettling (me)
to
find a well-developed piece of code that does the 'set this to zero if
negative' thing that kills further optimisation with no warning.

Having peeked at the lme4 code, I can't see an easy fix to suggest, so
apologies for the helpless plea: if the 'set to zero' can't easily be
avoided, would it be possible in the mean time to at least warn that
something has been set to zero and consequently that convergence may
not
be optimal?

Steve Ellison


Example code follows.
Platform is win 7, 32-bit, R version 2.12.2 (2011-02-25) with lme4
updated yesterday.

require(lme4)

mgb <-
structure(list(unit = structure(c(10L, 5L, 1L, 12L, 6L, 4L, 9L, 
11L, 8L, 2L, 7L, 3L, 12L, 5L, 1L, 10L, 11L, 4L, 9L, 6L, 7L, 2L, 
3L, 8L, 1L, 5L, 12L, 10L, 4L, 11L, 6L, 9L, 7L, 2L, 3L, 8L), .Label =
c("10", 
"14", "2", "20", "23", "34", "37", "43", "51", "56", "60", "65"
), class = "factor"), MG = c(2.8048, 2.6143, 2.8601, 2.8125, 
2.67788333333333, 2.87218333333333, 2.60828333333333, 2.77158333333333,

2.86959166666667, 2.83259166666667, 2.90769166666667, 2.80179166666667,

2.76881666666667, 2.82161666666667, 2.83231666666667, 2.88741666666667,

2.8035, 2.8723, 2.6975, 2.7232, 2.81370833333333, 2.84940833333333, 
2.84570833333333, 2.85160833333333, 2.72210833333333, 2.86660833333333,

2.84610833333333, 2.75790833333333, 3.47419166666667, 2.67299166666667,

2.74289166666667, 2.67809166666667, 2.6723, 2.6619, 2.7912, 2.6971
), Run = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("Run 1", 
"Run 2", "Run 3"), class = "factor")), .Names = c("unit", "MG", 
"Run"), row.names = c(1L, 2L, 10L, 12L, 13L, 22L, 23L, 24L, 25L, 
34L, 35L, 36L, 3L, 4L, 5L, 11L, 14L, 15L, 16L, 18L, 26L, 28L, 
29L, 31L, 6L, 7L, 8L, 9L, 17L, 19L, 20L, 21L, 27L, 30L, 32L, 
33L), class = "data.frame")

st <-
list(structure(0.1, .Dim = c(1L, 1L), .Dimnames = list("(Intercept)", 
    "(Intercept)")), structure(0.3, .Dim = c(1L, 1L), 
    .Dimnames = list("(Intercept)", "(Intercept)")))

lmer(MG~(1|Run)+(1|unit), data=mgb, subset=unit!="20", verbose=TRUE)
	#Note the zero intermediate value at iteration 1

#compare with a rather better answer:
lmer(MG~(1|Run)+(1|unit), data=mgb, subset=unit!="20", start=st,
verbose=TRUE)



*******************************************************************
This email and any attachments are confidential. Any\ us...{{dropped:16}}