lmer errors and predict function in R v 3.6.1 / obtain confidence intervals from lmer model for a given independent variable

Question

This is two questions in one, I hope this is OK.

First, I am trying to obtain confidence interval values from an lmer object from the lme4 package. I've previously used R v 3.4.4 and the models run just fine, and I can find and display the confidence intervals for the mean fit on the plot. I've recently upgraded to R v 3.6.1, and am now getting an error message when using the predict function. I've shown this below as part of the code

My main question is, how do I calculate the upper and lower confidence limits for a given value? For traditional lm I would use:

new.dat <- data.frame(variable = ##)
predict(lm_object, newdata = new.dat, interval = 'confidence')

But this doesn't work for lmer objects.

Here is the data:

TYPE <- c(rep("A", 100), rep("B", 31), rep("C", 18))
MAX<-c(NA,32.6,19.5,23.5,0,17.3,0,31,35.3,23.9,20.8,18.3,10.6,19.4,0,9,14.5,
       27.1,0,27.5,21,0,14.7,23.7,17.4,13.7,30.7,25.3,NA,0,16.5,0,NA,18.5,23.9,
       8.6,11.9,21.5,0,20.3,10.1,0,20.2,33.6,40.6,21.9,16.6,18.3,0,28.3,36.4,0,
       29.4,25.7,24.8,25,0,36.9,19,19.3,27.8,20.4,19.2,0,25.5,26.3,30.6,0,27.8,
       5.7,0,21,19.7,15.3,0,16.5,14.5,17.2,31.7,13,21.5,20,32.5,0,6.8,26.2,0,
       24.6,21.2,0,32.3,17.3,29.2,43.1,26.2,0,29.5,26.1,36.8,10.9,45.76,17.41,
       62.475,0,11.82,57.12,0,41.35,52.935,13.01,0,56.095,60.345,56.645,78.775,
       69.565,47.98,15.28,16.46,12.91,0,14.76,29.185,29.26,0,77.72,78.25,0,45.875,
       0,40.27,16.43,27.065,45.44,71.38,21.875,0,33.625,45.825,51.79,39.705,27.46,
       36.61,44.21,62.38,0,120.295,26.61,0)
STEM_PERCENT<-c(-8.708496564,-8.708496564,-8.708496564,-8.708496564,0,-8.708496564,
                0,-8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,0,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,0,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,0,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,0,-8.708496564,-8.708496564,0,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,0,-8.708496564,-8.708496564,-8.708496564,0,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,0,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,-8.708496564,-8.708496564,
                -8.708496564,0,-8.708496564,-8.708496564,0,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-8.708496564,0,-8.708496564,-8.708496564,
                -8.708496564,-8.708496564,-7.541459043,0,-7.541459043,0,0,0,0,
                -7.541459043,-7.541459043,0,0,-7.541459043,-7.541459043,-7.541459043,
                -7.541459043,-7.541459043,0,-7.541459043,8.156584524,-7.541459043,0,
                0,8.156584524,-7.541459043,8.156584524,0,-7.541459043,0,-7.541459043,
                0,0,0,0,-7.541459043,-15.08291809,0,0,0,-7.541459043,-7.541459043,
                -8.156584524,8.156584524,0,0,-16.31316905,0,-16.31316905,0,8.156584524)
val<- data.frame(TYPE, MAX, STEM_PERCENT)
na.strings=c("",NA)
val<-subset(val,MAX!= "NA")

And here is the code which I use to run the mixed effects analysis, which works just fine in R v 3.4.4

library(lme4)
library(merTools)
mod3<-lmer(STEM_PERCENT~1+MAX+(1|TYPE)+(0+MAX|TYPE),data=val)

Here, I get the following error message, which I never received in R v 3.4.4

Warning message:
In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.107451 (tol = 0.002, component 1)

code continues

mod4<-lmer(STEM_PERCENT~1+MAX+(1+MAX|TYPE),data=val) # also produces error not seen before


val2<-expand.grid(MAX=seq(0,120,length=1000),
                  TYPE=levels(val$TYPE))

val2$STEM_PERCENT<-predict(mod3,newdata=val2)

val3<-data.frame(MAX=seq(0,120,0.1))
val3$STEM_PERCENT<-predict(mod3,newdata=val3,re.form=~0)

CI <- cbind(val2, predictInterval(mod4, val2))


plot(val$MAX,val$STEM_PERCENT,pch=19,col=(1+as.integer((val$TYPE))),
     bty="l",cex=0.5,xaxt="n",las=1,ylim=c(-30,10),
     xlim=c(0,120),yaxt="n",xlab=NA,ylab=NA)

axis(side=2,tck=-0.02,at=seq(-30,10,10),cex.axis=0.5,
     font.axis=1,las=2,mgp=c(0,.5,0),labels=T)
axis(side=1,tck=-0.02,at=seq(0,120,20),
     cex.axis=0.5,labels=T,mgp=c(0,-0.1,0))


xv<-seq(0,120,0.01)
typea<-rep("A",length(xv))
yv<-predict(mod3,list(MAX=xv,TYPE=typea),type="response")

Here is where I get the error message in R v 3.6.1, which I don't get in v 3.4.4. The same is true of each instance of yv

Error in rep(0, nobs) : invalid 'times' argument

code continues

lines(xv,yv,col="red",lwd=1.5)

typeb<-rep("B",length(xv))
yv<-predict(mod3,list(MAX=xv,TYPE=typeb),type="response")
lines(xv,yv,col="green",lwd=1.5)

typec<-rep("C",length(xv))
yv<-predict(mod3,list(MAX=xv,TYPE=typec),type="response")
lines(xv,yv,col="blue",lwd=1.5)

lines(val3$MAX,val3$STEM_PERCENT,lwd=2)

ssmin<-smooth.spline(CI$MAX,CI$lwr,df=4)
lines(ssmin,lty=2,col="black",lwd=1)

ssmax<-smooth.spline(CI$MAX,CI$upr,df=4)
lines(ssmax,lty=2,col="black",lwd=1)

Which makes the following plot in R v 3.4.4 with no problems.

I'd like to be able to find upper and low confidence interval values around the mean at a given value of MAX, e.g. upr and lwr limits of the fit when MAX = 50. It looks like it should be somehwere around -2 and -12 for upr and lwr limits respectively. I'd also like to know what's going on with the error messages, which I've not previously received.

Any help and advice gratefully received, thanks

Relevant are most likely not the different R versions but the different versions of lme4. — Roland
The newdata parameter is documented to expect a data.frame and, indeed, yv<-predict(mod3, data.frame(MAX=xv,TYPE=typea), type="response") works just fine. If you haven't already, you should probably study help("convergence") regarding the "no convergence" warnings (which are not errors). — Roland
Thanks, @Roland. I solved the convergence problem by restarting the models with the parameter estimates before the algorithm failed to converge, thus: ss <- getME(mod3,c("theta","fixef")) new_mod <- update(mod3,start=ss) For others, there is more help on the convergence problem here Can you offer advice on the CI query? — A.Benson

A.Benson A.Benson · Accepted Answer · 2019-07-18T20:46:27

I figured out how to determine CI values. Hopefully others may benefit from the solution. As usual, it was quite simple.

First, round the CIs to 0 d.p

CI$MAX<-as.numeric(round(CI$MAX,0))

Then, determine mean upper and lower values at MAX = 50

low<-mean(CI$lwr[which(CI$MAX==50)])
high<-mean(CI$upr[which(CI$MAX==50)])

Add lines to plot to check

abline(v = 50)
abline(h = low)
abline(h = high)

lmer errors and predict function in R v 3.6.1 / obtain confidence intervals from lmer model for a given independent variable

1 Answers