[R-sig-ME] Problems with convergence

[Ben Bolker]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

  Hmm.  Have you tried with nlminb?  This works for me ...

glmer1a <- glmer(cbind(nEvents,total-nEvents) ~ -1 + trt +
     factor(id) +
      (0+trt12|id), data=thedata, family=binomial, nAGQ=7)

glmer1X <- update(glmer1a,
                  control=glmerControl(optimizer="nlminbwrap"))


  I don't know whether it works better overall (it does use a
different convergence criterion ...) or whether this is just the luck
of the draw on this particular case.

   cheers
    Ben Bolker

PS I don't remember when nlminbwrap was added to the code base, but
it's pretty simple:

function (par, fn, lower, upper, control = list(), ...)
{
    res <- nlminb(start = par, fn, gradient = NULL, hessian = NULL,
        scale = 1, lower = lower, upper = upper, control = control,
        ...)
    list(par = res$par, fval = res$objective, conv = res$convergence,
        message = res$message)
}


On 15-04-02 05:27 AM, Ken Beath wrote:
> It would be nice to have this work properly, as I need it for
> certain things and it seems that other people are having similar
> problems.
> 
> Getting it to work by increasing the quadrature points is a bit of
> an aberration, it is not what typically happens, and I've put an
> example at the end. At least in this one the profiling works which
> means the maximum must be fairly close to that obtained from the
> optimisation.
> 
> My feeling on this, is that possibly the problem is not with the
> optimiser, seeing that it fails with so many optimisers, but rather
> with the calculation of the marginal likelihood. These optimisers
> don't tend to stop with 0.001 gradients. When I have time I will
> find in the code how the node locations are calculated and see what
> is happening.
> 
> Anyway, here is one that fails irrespective of  nAGQ value.
> 
> thedata <- structure(list(nEvents = c(10L, 53L, 17L, 18L, 22L, 6L,
> 16L, 14L, 13L, 18L, 15L, 19L, 52L, 19L, 8L, 16L, 50L, 8L, 9L, 4L, 
> 26L, 45L, 18L, 20L, 5L, 16L, 18L, 7L, 3L, 19L, 30L, 26L, 66L, 23L,
> 29L, 18L, 72L, 25L, 9L, 2L), total = c(200, 200, 200, 200, 200,
> 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
> 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200,
> 200, 200, 200, 200, 200, 200, 200, 200, 200), trt = c(0, 0, 0, 0,
> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1,
> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), id = structure(c(1L, 2L,
> 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L,
> 18L, 19L, 20L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L,
> 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L), .Label = c("1", "2", "3",
> "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15",
> "16", "17", "18", "19", "20"), class = "factor"), trt12 = c(-0.5,
> -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
> -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, 0.5, 0.5, 0.5, 0.5,
> 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
> 0.5, 0.5, 0.5)), .Names = c("nEvents", "total", "trt", "id",
> "trt12"), row.names = c(NA, 40L), class = "data.frame")
> 
> glmer1a <- glmer(cbind(nEvents,total-nEvents) ~ -1 + trt +
> factor(id) + (0+trt12|id), data=thedata, family=binomial, nAGQ=7)
> 
> glmer1b <- glmer(cbind(nEvents,total-nEvents) ~ -1 + trt +
> factor(id) + (0+trt12|id), data=thedata, family=binomial, nAGQ=21)
> 
> 
> 
> 
> 
> On 1 April 2015 at 02:25, Viechtbauer Wolfgang (STAT) < 
> wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:
> 
>> This discussion piqued my interest. The model that Ken was
>> fitting is in essence one of the models that is fitted by the
>> rma.glmm() function in the metafor package. This is sometimes
>> called the unconditional model with fixed study effects. To
>> illustrate:
>> 
>> ### original data
>> 
>> thedata <- structure(list(nEvents=c(10L,53L,17L,18L,22L,6L,16L, 
>> 14L,13L,18L,15L,19L,52L,19L,8L,16L,50L,8L,9L,4L, 
>> 26L,45L,18L,20L,5L,16L,18L,7L,3L,19L,30L,26L,66L, 
>> 23L,29L,18L,72L,25L,9L,2L),total=c(200,200,200,200, 
>> 200,200,200,200,200,200,200,200,200,200,200,200,200, 
>> 200,200,200,200,200,200,200,200,200,200,200,200,200, 
>> 200,200,200,200,200,200,200,200,200,200),trt=c(0, 
>> 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1, 
>> 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),id=structure(c(1L, 
>> 2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,12L,13L,14L,15L, 
>> 16L,17L,18L,19L,20L,1L,2L,3L,4L,5L,6L,7L,8L,9L, 
>> 10L,11L,12L,13L,14L,15L,16L,17L,18L,19L,20L),.Label=c("1", 
>> "2","3","4","5","6","7","8","9","10","11","12","13", 
>> "14","15","16","17","18","19","20"),class="factor")),.Names=c("nEvents",
>>
>> 
"total","trt","id"),row.names=c(NA,40L),class="data.frame")
>> 
>> ### restructure data as needed for input into rma.glmm()
>> 
>> dat <- cbind(thedata[1:20,], thedata[21:40,]) dat$id <- dat$id <-
>> dat$trt <- dat$trt <- NULL colnames(dat) <- c("ci", "n2i", "ai",
>> "n1i")
>> 
>> library(metafor) library(lme4)
>> 
>> ### model fitted by Ken res1 <-
>> glmer(cbind(nEvents,total-nEvents) ~ trt + factor(id) + 
>> (0+trt|id), data=thedata, family=binomial)
>> 
>> ### fit unconditional model with fixed study effects via
>> rma.glmm() res2 <- rma.glmm(measure="OR", ai=ai, n1i=n1i, ci=ci,
>> n2i=n2i, data=dat, nAGQ=1)
>> 
>> ### to get exact equivalence, use +-1/2 coding for the random
>> effects thedata$trt12 <- thedata$trt - 1/2 res3 <-
>> glmer(cbind(nEvents,total-nEvents) ~ -1 + trt + factor(id) + 
>> (0+trt12|id), data=thedata, family=binomial)
>> 
>> summary(res1) summary(res2) summary(res3)
>> 
>> ### end example
>> 
>> A few notes:
>> 
>> 1) rma.glmm() uses nAGQ=7 by default, so I switched that to 1 for
>> the comparison.
>> 
>> 2) Some discussion of the 0/1 versus +-1/2 coding can be found in
>> Turner et al. (2000) and Higgins et al. (2001). I tend to prefer
>> the +-1/2 coding, so that is also what is currently implemented
>> in rma.glmm(), but I may add the 0/1 coding as an option.
>> 
>> 3) A nice discussion of the model is provided by Senn (2000). He
>> also discusses a variety of other modeling options, including a
>> model using random study effects.
>> 
>> 4) In fact, the unconditional model with random study effects can
>> be fitted with:
>> 
>> rma.glmm(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat,
>> model=" UM.RS")
>> 
>> (which makes use of glmer() underneath). As discussed by Senn,
>> this model may violate what he calls the 'concurrent control
>> principle', but his wording is cautious ('may violate', 'may be
>> regarded as undesirable'), which reflects the lack of a thorough
>> discussion in the literature comparing the various models.
>> 
>> 5) Yet another option is the (mixed-effects) conditional logistic
>> model. See, for example, Stijnen et al. (2010). This model is
>> obtained when conditioning on the total number of events within
>> each study and leads to non-central hypergeometric distributions
>> for the data within each study. This model can be fitted with:
>> 
>> rma.glmm(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, 
>> model="CM.EL")
>> 
>> Sorry, it's slow (I haven't found a clever way of speeding up
>> the integration over the non-central hypergeometric
>> distributions). Much faster, thanks to lme4, is:
>> 
>> rma.glmm(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat,
>> model=" CM.AL")
>> 
>> which uses an approximation to the exact conditional likelihood.
>> 
>> 6) And of course there are Bayesian implementations of such
>> models.
>> 
>> 7) With respect to the model fitted by Ken, it's maybe
>> interesting to note that NOT using the Laplace approximation, but
>> something like 7 quadrature points, does not cause any
>> convergence warnings:
>> 
>> glmer(cbind(nEvents,total-nEvents) ~ trt + factor(id) +
>> (0+trt|id), data=thedata, family=binomial, nAGQ=7)
>> 
>> Alright, I'll shut up now.
>> 
>> References mentioned above:
>> 
>> Higgins, J. P. T., Whitehead, A., Turner, R. M., Omar, R. Z., &
>> Thompson, S. G. (2001). Meta-analysis of continuous outcome data
>> from individual patients. Statistics in Medicine, 20(15),
>> 2219-2241.
>> 
>> Senn, S. (2000). The many modes of meta. Drug Information
>> Journal, 34, 535-549.
>> 
>> Stijnen, T., Hamza, T. H., & Ozdemir, P. (2010). Random effects 
>> meta-analysis of event outcome in the framework of the
>> generalized linear mixed model with applications in sparse data.
>> Statistics in Medicine, 29(29), 3046-3067.
>> 
>> Turner, R. M., Omar, R. Z., Yang, M., Goldstein, H., & Thompson,
>> S. G. (2000). A multilevel model framework for meta-analysis of
>> clinical trials with binary outcomes. Statistics in Medicine,
>> 19(24), 3417-3432.
>> 
>> Best, Wolfgang
>> 
>> -- Wolfgang Viechtbauer, Ph.D., Statistician Department of
>> Psychiatry and Neuropsychology School for Mental Health and
>> Neuroscience Faculty of Health, Medicine, and Life Sciences 
>> Maastricht University, P.O. Box 616 (VIJV1) 6200 MD Maastricht,
>> The Netherlands +31 (43) 388-4170 | http://www.wvbauer.com
>> 
>>> -----Original Message----- From: R-sig-mixed-models
>>> [mailto:r-sig-mixed-models-bounces at r- project.org] On Behalf Of
>>> Ben Bolker Sent: Monday, March 30, 2015 22:21 To:
>>> r-sig-mixed-models at r-project.org Subject: Re: [R-sig-ME]
>>> Problems with convergence
>>> 
>>> Ken Beath <ken.beath at ...> writes:
>>> 
>>>> Yes, I was demonstrating that it fails convergence and then
>>>> as a consequence fails to profile. I have my doubts about 
>>>> convergence for the bobyqa algorithm, I have other
>>>> applications where it doesn't converge
>>> properly.
>>>> For some of my own work I've used nlminb followed by
>>>> Nelder-Mead if
>>> there
>>>> is a convergence failure. Not optimal but it seems to work.
>>> 
>>> I'm still not sure whether you expect it to converge (I think
>>> you do), or whether you are just pointing out that the
>>> convergence warning in this case is probably justified (in the
>>> face of so many convergence warnings that turn out to be false
>>> positives, this is a useful piece of information).
>>> 
>>>> While it is fairly heavily parameterised it is a real model, 
>>>> a frequentist implementation of Smith, T. C., Spiegelhalter,
>>>> D. J., & Thomas, a.
>>> (1995).
>>>> Bayesian approaches to random-effects meta-analysis: a
>>>> comparative
>>> study.
>>>> Statistics in Medicine, 14(24), 2685–99. The reason for
>>>> having studies
>>> as
>>>> fixed effects is probably philosophical, the overall success
>>>> rates are
>>> not
>>>> likely to be given by normally distributed random effects,
>>>> and are in
>>> many
>>>> cases specifically chosen.
>>> 
>>> I can appreciate that, but I still think it's unrealistic to
>>> expect to be able to fit 22 parameters to 40 observations
>>> except under very special circumstances.  One point about
>>> switching from the Bayesian to the frequentist world is that
>>> the Bayesians (by definition) put priors on their parameters,
>>> which provides a degree of regularization that is not by
>>> default available to frequentist methods.  What priors did
>>> Smith et al. use?  It might be worth trying this in blme with 
>>> priors on the fixed effects ...
>>> 
>>>> I did find that one of the data sets that I have also failed,
>>>> but
>>> fitted
>>>> with a commercial program that is based on the EM algorithm.
>>>> For this
>>> type
>>>> of problem it is actually faster, as any type of quasi-Newton
>>>> needs to calculate lots of derivatives.
>>> 
>>> I could whine about the difficulty of finding globally robust, 
>>> reliable, and fast optimization algorithms, but I won't.  I can
>>> certainly appreciate that there are more reliable methods for
>>> particular sub-classes of problems.
>>> 
>>>> Anyway, I'm going to keep looking at the methods, and
>>>> eventually the
>>> code
>>>> for glmer and may eventually have some suggestions.
>>> 
>>> Would be happy to hear them.
>>> 
>>> It's worth pointing out that lme4 is using a preliminary
>>> "nAGQ=0" step, which ignores the terms contributed by the
>>> integrals over the distributions of the conditional modes and
>>> as a result is able to fit both the fixed-effect parameters and
>>> the conditional modes in a single linear-algebra step, reducing
>>> the dimensionality of the nonlinear optimization to the length
>>> of the variance-covariance parameter vector ...
>>> 
>>>> On 19 March 2015 at 14:45, Ben Bolker <bbolker <at>
>>>> gmail.com> wrote:
>>>> 
>>>>> Ken Beath <ken.beath <at> ...> writes:
>>>>> 
>>>>>> The following code shows that there are convergence
>>>>>> problem
>>> messages
>>>>>> where there is a problem with convergence. The profiling
>>>>>> shows that the maximum found is not the correct one. This
>>>>>> is simulated data
>>> for
>>>>>> a binary meta-analysis with fixed effect for study and
>>>>>> random
>>> effect
>>>>>> for treatment.
>>>>> 
>>> 
>>> [paragraph snipped to try to make Gmane happy]
>>> 
>>>>> However, may I comment that this is a slightly ridiculous
>>>>> scenario? The data set here has 40 observations, and the
>>>>> model tries to fit 22 parameters.  The model that treats id
>>>>> as a random effect works much better.  I can believe there
>>>>> are scenarios where you really do want study as a fixed
>>>>> effect, but did you expect it to be practical here?
>>>>> 
>>>>> But maybe you're just trying to show that this is a "true
>>>>> positive" case for the convergence warnings.
>>>>> 
>>>>> Some random code I wrote while diagnosing what was going
>>>>> on:
>>>>> 
>>>>> library(ggplot2); theme_set(theme_bw())
>>>>> 
>>>>> ## proportion + weights is a little easier to handle 
>>>>> thedata <- transform(thedata,prop=nEvents/total)
>>>>> 
>>>>> ggplot(thedata,aes(trt,prop))+geom_point(aes(size=total))+ 
>>>>> geom_line(aes(group=id),colour="gray") glmer1 <-
>>>>> glmer(prop~trt+factor(id)+(0+trt|id), 
>>>>> weights=total,data=thedata,family=binomial)
>>>>> 
>>>>> ## id as RE glmer2 <- glmer(prop~trt+(1|id)+(0+trt|id), 
>>>>> weights=total,data=thedata,family=binomial)
>>>>> 
>>>>> dd <- update(glmer1,devFunOnly=TRUE) pars <-
>>>>> unlist(getME(glmer1,c("theta","fixef"))) library("bbmle") 
>>>>> ss <- slice2D(pars,dd) library("lattice") plot(ss) ## too
>>>>> complex, but too much work to cut down significantly
>>>>> 
>>>>>> library(lme4)
>>>>>> 
>>>>>> thedata <-
>>>>>> structure(list(nEvents=c(10L,53L,17L,18L,22L,6L,16L, 
>>>>>> 14L,13L,18L,15L,19L,52L,19L,8L,16L,50L,8L,9L,4L, 
>>>>>> 26L,45L,18L,20L,5L,16L,18L,7L,3L,19L,30L,26L,66L, 
>>>>>> 23L,29L,18L,72L,25L,9L,2L),total=c(200,200,200,200, 
>>>>>> 200,200,200,200,200,200,200,200,200,200,200,200,200, 
>>>>>> 200,200,200,200,200,200,200,200,200,200,200,200,200, 
>>>>>> 200,200,200,200,200,200,200,200,200,200),trt=c(0, 
>>>>>> 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1, 
>>>>>> 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),id=structure(c(1L, 
>>>>>> 2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,12L,13L,14L,15L, 
>>>>>> 16L,17L,18L,19L,20L,1L,2L,3L,4L,5L,6L,7L,8L,9L, 
>>>>>> 10L,11L,12L,13L,14L,15L,16L,17L,18L,19L,20L),.Label=c("1",
>>>>>>
>>>>>> 
"2","3","4","5","6","7","8","9","10","11","12","13",
>>>>>> 
>>> "14","15","16","17","18","19","20"),class="factor")),.Names=c("nEvents",
>>>>>>
>>> 
"total","trt","id"),row.names=c(NA,40L),class="data.frame")
>>>>>> 
>>>>>> glmer1<-glmer(cbind(nEvents,total-nEvents)~trt+factor(id)+
>>>>>
>>>>>> 
##   (0+trt|id),data=thedata,family=binomial)
>>>>>> 
>>>>>> # while glmer has problems with component 9 it is 8 with
>>>>>> a problem
>>>>> profile
>>>>>> # I've use devtol so the discrepancy is printed 
>>>>>> prof.glmer1<-profile(glmer1,which=8,devtol=1.0e-3)
>> 
> 
> 
> 

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.11 (GNU/Linux)

iQEcBAEBAgAGBQJVHU/nAAoJEOCV5YRblxUHxPEIAKQ6kISS1aANugvV6wvq7/6C
tIkLUZIxxqCn1UzfbB4naqVQ/4B9ipNlpWV8FEVJs5H7IxoQY1iCojONqihmiHO6
SE/m44U4xEsX94q47d/vWCNqbwk9sGyXZj7R4lAOsPZiUqkvaKKV4Py8IF+O2nwM
RjuzSGGw3ZusE7gIX6UB79lpgVvE8RZUOzJwR8pjKcBhbBHSYslAmdd9o8HpXTkm
EWhZ3PpFEnWopXfeG2uN5HWGD8OjP7w8578ZHKLJhOW/9VEZQ9ykI/I8/yOlb7OJ
dfbzRpQe/cKFghMYvF/xR2y0tdp98OwQUE7tKJneOJLqV2Q1u9w/m4oCRXVYNmA=
=8saW
-----END PGP SIGNATURE-----