[R-sig-ME] Problems with convergence

[Ben Bolker]

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.
 
  I'm not 100% sure what the point is here (that's a real question,
not a criticism).  I *think* it's that we can get a convergence
failure with a computed gradient size that seems fairly innocuous
(only 0.003, which is perhaps below what I where I would have set
the new cutoff to eliminate a lot of the probably-false-positives
that people have been getting with larger data sets), which is 
nevertheless a real convergence failure.  This underscores the
need to set a gradient tolerance that is data-set-size-dependent
(which the developers have discussed a bit, but maybe not in a 
public forum).

  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)
>