Convergence testing

library(lattice)
library(lme4)
library(blme)
library(reshape2)
library(ggplot2); theme_set(theme_bw())
library(gridExtra)  ## for grid.arrange
library(bbmle) ## for slice2D; requires *latest* r-forge version (r121)
source("allFit.R")

Load data:

dList <- load("data.RData")

Spaghetti plot: don’t see much pattern other than (1) general increasing trend; (2) quantized response values (table(dth$Estimate) or unique(dth$Estimate) also show this); (3) skewed residuals

sort(unique(dth$Estimate))
##  [1]  0.0  0.5  1.0  1.5  2.0  2.5  3.0  4.0  5.0  6.0  7.0  7.5  8.0 10.0
## [15] 12.0 15.0 17.0 20.0 25.0 30.0 35.0 40.0 50.0 60.0 70.0 75.0 90.0
(p0 <- ggplot(dth,aes(Actual,Estimate))+geom_point()+
    geom_line(aes(group=factor(pid)))+
    stat_summary(fun.y=mean,geom="line",colour="red",lwd=2))

plot of chunk spaghett1

getcor <- function(object) attr(VarCorr(object)[[1]],"correlation")[2,1]
f1 <- lmer(Estimate~Actual+(Actual|pid),data=dth,REML=FALSE)
getcor(f1)

logging (there are zero values: it’s a bit sloppy, but I’m just going to use log(0.5+x) as the response – log(1+x) may look even better …) helps a great deal with the fit (no longer singular), and with the Q-Q

f1L <- lmer(log10(0.5+Estimate)~Actual+(Actual|pid),data=dth,REML=FALSE)
VarCorr(f1L)
##  Groups   Name        Std.Dev. Corr 
##  pid      (Intercept) 0.28539       
##           Actual      0.00633  -0.20
##  Residual             0.27445

Q-Q plot shows the problem …

grid.arrange(qqmath(f1),qqmath(f1L),nrow=1)

plot of chunk qq

and helps the general residual plot

grid.arrange(plot(f1,type=c("p","smooth")),
             plot(f1L,type=c("p","smooth")),nrow=1)

plot of chunk diag

Double-check that the same thing happens with lme4.0

f2 <- lme4.0::lmer(Estimate~Actual+(Actual|pid),data=dth,REML=FALSE)
attr(lme4.0::VarCorr(f2)[[1]],"correlation")[2,1]
## [1] 1
aa <- allFit(f1)
## 
## Attaching package: 'minqa'
## 
## The following objects are masked from 'package:nloptr':
## 
##     bobyqa, newuoa
ss <- summary.allfit(aa)

results: warning for L-BFGS-B (which does indeed fail); all other results nearly identical (nlminb gives cor slightly <1) (only sd/cor estimates shown; log-likelihoods and fixed effect estimates are consistent with this picture) (Facet and axis tick labels below are slightly wonky/need tweaking.)

names(which(!sapply(ss$msgs,is.null)))
## [1] "optimx.L-BFGS-B"
(g0 <- ggplot(melt(ss$sdcor),aes(Var2,value,colour=Var1))+
    geom_point(position=position_dodge(width=0.5))+
    facet_wrap(~Var2,scale="free"))

plot of chunk sdcor1

blmer

L-BFGS-B fails completely (unusual error, missing value pops up in internal boundary test – should investigate) none of the other cases gives convergence warnings, which is slightly alarming since both versions of BOBYQA (minqa and nloptr) clearly fail, with NLL 50-100 log-likelihood units worse than the other fits … I think the problem here is actually something in blmer where the warnings fail to get stored properly …

## hack to refit lmer as glmer
ff <- getCall(f1)
ff[[1]] <- quote(blmer)
f3 <- eval(ff)
## Warning: Model failed to converge with max|grad| = 0.446135 (tol = 0.002, component 3)
## Warning: Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
aa3 <- allFit(f3)
## bobyqa. :
## Warning: Model failed to converge with max|grad| = 0.446135 (tol = 0.002, component 3)
## Warning: Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
## [OK]
## Nelder_Mead. : [OK]
## optimx.nlminb : [OK]
## optimx.L-BFGS-B :
## Warning: convergence code 9999 from optimx
## [OK]
## nloptWrap.NLOPT_LN_NELDERMEAD : [OK]
## nloptWrap.NLOPT_LN_BOBYQA :
## Warning: Model failed to converge with max|grad| = 2.06062 (tol = 0.002,
## component 1)
## [OK]
ss3 <- summary.allfit(aa3)

Despite the failure of blme to store the warnings, we can see from the output that the warnings do correspond to the incorrect results (i.e., both BOBYQA implementations). It’s interesting that BOBYQA fails here, since in most other cases we have found it to be the most robust algorithm …

g0 %+% melt(cbind(ss3$sdcor,NLL=-2*(ss3$ll-max(ss3$ll))))

plot of chunk blmeplot

Visualizing likelihood surface

rmat <- matrix(c(0,1,-0.2,0.2,0,1),byrow=TRUE,ncol=2,
               dimnames=list(paste0("theta",1:3),c("lwr","upr")))
dd1 <- update(f1,devFunOnly=TRUE)
pp1 <- getME(f1,"theta")
s2D.1 <- bbmle:::slice2D(pp1,dd1,
                       tranges=rmat)
##
dd3 <- update(f3,devFunOnly=TRUE)
pp3 <- getME(aa3[[2]],"theta")  ## use a version that actually worked
s2D.3 <- bbmle:::slice2D(pp3,dd3,
                       tranges=rmat)

Left plot, lmer; right, blmer.

grid.arrange(splom(s2D.1),splom(s2D.3),nrow=1)

plot of chunk sploms

Zoom in a bit:

rmat2 <- matrix(c(0,1,-0.1,0.1,0,0.1),byrow=TRUE,ncol=2,
               dimnames=list(paste0("theta",1:3),c("lwr","upr")))
s2D.1b <- bbmle:::slice2D(pp1,dd1,
                       tranges=rmat2)
s2D.3b <- bbmle:::slice2D(pp3,dd3,
                       tranges=rmat2)
grid.arrange(splom(s2D.1b),splom(s2D.3b),nrow=1)

plot of chunk sploms2

Visualizations could use a lot of work (scale bar and/or contours; allow for more points, e.g. visualize where non-converging fits ended up …)

Stuff left out (so far):