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

[Douglas Bates]

On Tue, Apr 19, 2011 at 4:41 AM, S Ellison <S.Ellison at lgc.co.uk> wrote:
> 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?

The problems are due to the optimizer used in the lme4 package, which
is the one from the nlminb function in R.  It has a habit of getting
stuck at the lower bound.  The lme4a package uses a different
optimizer from the minqa package and is more successful in this case.

> lmer(MG~(1|Run)+(1|unit), data=mgb, subset=unit!="20", verbose=2L)
npt = 4 , n =  2
rhobeg =  0.2 , rhoend =  2e-07
   0.020:   9:     -66.6254;0.298266 0.332244
  0.0020:  14:     -66.6367;0.289025 0.365102
 0.00020:  17:     -66.6367;0.288926 0.363631
 2.0e-05:  19:     -66.6367;0.289092 0.363130
 2.0e-06:  22:     -66.6367;0.289109 0.363093
 2.0e-07:  25:     -66.6367;0.289110 0.363092
At return
 27:    -66.636729: 0.289110 0.363092

Linear mixed model fit by REML ['merMod']
Formula: MG ~ (1 | Run) + (1 | unit)
   Data: mgb
 Subset: unit != "20"
REML criterion at convergence: -66.6367

Random effects:
 Groups   Name        Variance  Std.Dev.
 unit     (Intercept) 0.0004821 0.02196
 Run      (Intercept) 0.0007604 0.02758
 Residual             0.0057678 0.07595
Number of obs: 33, groups: unit, 11; Run, 3

Fixed effects:
            Estimate Std. Error t value
(Intercept)  2.77470    0.02173   127.7


>
> 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}}
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>