With LM and GLM the predict function can return the standard error for the predicted values on either the observed data or on new data. This is then used to draw confidence or prediction intervals around the fitted regression lines. The confidence intervals (CI) focus on the regression lines and can be interpreted as (assuming that we draw 95% CI): “If we would repeat our sampling X times the regression line would fall between this interval 95% of the time”. On the other hand the prediction interval focus on single data point and could be interpreted as (again assuming that we draw 95% CI): “If we would sample X times at these particular value for the explanatory variables, the response value would fall between this interval 95% of the time”. The wikipedia page has some nice explanation about the meaning of confidence intervals.

For GLMM the predict function does not allow one to derive standard error, the reason being (from the help page of predict.merMod): “There is no option for computing standard errors of predictions because it is difficult to define an efficient method that incorporates uncertainty in the variance parameters”. This means there is for now no way to include in the computation of the standard error for predicted values the fact that the fitted random effect standard deviation are just estimates and may be more or less well estimated. We can however still derive confidence or prediction intervals keeping in mind that we might underestimate the uncertainty around the estimates.

library(lme4) #first case simple lmer, simulate 100 data points from 10 groups with one continuous fixed effect variable x<-runif(100,0,10) f<-gl(n = 10,k = 10) data<-data.frame(x=x,f=f) modmat<-model.matrix(~x,data) #the fixed effect coefficient fixed<-c(1,0.5) #the random effect rnd<-rnorm(10,0,0.7) #the simulated response values data$y<-rnorm(100,modmat%*%fixed+rnd[f],0.3) #model m<-lmer(y~x+(1|f),data) #first CI and PI using predict-like method, using code posted here: http://glmm.wikidot.com/faq newdat<-data.frame(x=seq(0,10,length=20)) mm<-model.matrix(~x,newdat) newdat$y<-mm%*%fixef(m) #predict(m,newdat,re.form=NA) would give the same results pvar1 <- diag(mm %*% tcrossprod(vcov(m),mm)) tvar1 <- pvar1+VarCorr(m)$f[1] # must be adapted for more complex models newdat <- data.frame( newdat , plo = newdat$y-1.96*sqrt(pvar1) , phi = newdat$y+1.96*sqrt(pvar1) , tlo = newdat$y-1.96*sqrt(tvar1) , thi = newdat$y+1.96*sqrt(tvar1) ) #second version with bootMer #we have to define a function that will be applied to the nsim simulations #here we basically get a merMod object and return the fitted values predFun<-function(.) mm%*%fixef(.) bb<-bootMer(m,FUN=predFun,nsim=200) #do this 200 times #as we did this 200 times the 95% CI will be bordered by the 5th and 195th value bb_se<-apply(bb$t,2,function(x) x[order(x)][c(5,195)]) newdat$blo<-bb_se[1,] newdat$bhi<-bb_se[2,] plot(y~x,data) lines(newdat$x,newdat$y,col="red",lty=2,lwd=3) lines(newdat$x,newdat$plo,col="blue",lty=2,lwd=2) lines(newdat$x,newdat$phi,col="blue",lty=2,lwd=2) lines(newdat$x,newdat$tlo,col="orange",lty=2,lwd=2) lines(newdat$x,newdat$thi,col="orange",lty=2,lwd=2) lines(newdat$x,newdat$bhi,col="darkgreen",lty=2,lwd=2) lines(newdat$x,newdat$blo,col="darkgreen",lty=2,lwd=2) legend("topleft",legend=c("Fitted line","Confidence interval","Prediction interval","Bootstrapped CI"),col=c("red","blue","orange","darkgreen"),lty=2,lwd=2,bty="n")

This looks pretty familiar, the prediction interval being always bigger than the confidence interval.

Now in the help page for the predict.merMod function the authors of the lme4 package wrote that bootMer should be the prefered method to derive confidence intervals from GLMM. The idea there is to simulate N times new data from the model and get some statistic of interest. In our case we are interested in deriving the bootstrapped fitted values to get confidence interval for the regression line. bb$t is a matrix with the observation in the column and the different bootstrapped samples in the rows. To get the 95% CI for the fitted line we then need to get the [0.025*N,0.975*N] values of the sorted bootstrapped values.

The bootstrapped CI falls pretty close to the “normal” CI, even if for each bootstrapped sample new random effect values were calculated (because use.u=FALSE per default in bootMer)

Now let’s turn to a more complex example, a Poisson GLMM with two crossed random effects:

#second case more complex design with two crossed RE and a poisson response x<-runif(100,0,10) f1<-gl(n = 10,k = 10) f2<-as.factor(rep(1:10,10)) data<-data.frame(x=x,f1=f1,f2=f2) modmat<-model.matrix(~x,data) fixed<-c(-0.12,0.35) rnd1<-rnorm(10,0,0.7) rnd2<-rnorm(10,0,0.2) mus<-modmat%*%fixed+rnd1[f1]+rnd2[f2] data$y<-rpois(100,exp(mus)) m<-glmer(y~x+(1|f1)+(1|f2),data,family="poisson") #for GLMMs we have to back-transform the prediction after adding/removing the SE newdat<-data.frame(x=seq(0,10,length=20)) mm<-model.matrix(~x,newdat) y<-mm%*%fixef(m) pvar1 <- diag(mm %*% tcrossprod(vcov(m),mm)) tvar1 <- pvar1+VarCorr(m)$f1[1]+VarCorr(m)$f2[1] ## must be adapted for more complex models newdat <- data.frame( x=newdat$x, y=exp(y), plo = exp(y-1.96*sqrt(pvar1)) , phi = exp(y+1.96*sqrt(pvar1)) , tlo = exp(y-1.96*sqrt(tvar1)) , thi = exp(y+1.96*sqrt(tvar1)) ) #second version with bootMer predFun<-function(.) exp(mm%*%fixef(.)) bb<-bootMer(m,FUN=predFun,nsim=200) bb_se<-apply(bb$t,2,function(x) x[order(x)][c(5,195)]) newdat$blo<-bb_se[1,] newdat$bhi<-bb_se[2,] #plot plot(y~x,data) lines(newdat$x,newdat$y,col="red",lty=2,lwd=3) lines(newdat$x,newdat$plo,col="blue",lty=2,lwd=2) lines(newdat$x,newdat$phi,col="blue",lty=2,lwd=2) lines(newdat$x,newdat$tlo,col="orange",lty=2,lwd=2) lines(newdat$x,newdat$thi,col="orange",lty=2,lwd=2) lines(newdat$x,newdat$bhi,col="darkgreen",lty=2,lwd=2) lines(newdat$x,newdat$blo,col="darkgreen",lty=2,lwd=2) legend("topleft",legend=c("Fitted line","Confidence interval","Prediction interval","Bootstrapped CI"),col=c("red","blue","orange","darkgreen"),lty=2,lwd=2,bty="n")

Again in this case the bootstrapped CI falled pretty close to the “normal” CI. We have here seen three different way to derive intervals representing the uncertainty around the regression lines (CI) and the response points (PI). Choosing among them would depend on what you want to see (what is the level of uncertainty around my fitted line vs if I sample new observations which value will they take), but also for complex model on computing power, as bootMer can take some time to run for GLMM with many observations and complex model structure.

very well explained!

lines 40 & 41 in the first code chunk and 36 & 37 in the second

newdat$blo<-bb_se[1,]

newdat$bhi<-bb_se[2,]

should probably be

newdat$blo<-bb_se[5,]

newdat$bhi<-bb_se[195,]

let me know if correct

Thanks for your comment and interest!

Actually bb$t is a matrix with B x n dimension were B is the number of bootstrapped replication (here 200) and n the number of observation (here 20). We then have the observations in the columns and the bootstrapped replication in the lines. We are here interested in the lower and upper limit of the 95% confidence interval so to get for each column the 5th and 195th value on the sorted vector of bootstrapped replication. This is what the apply call does. We then have 2 x n matrix (bb_se) with the first line being the lower limit and the second line the upper one.

As always in R there are 1000s of ways to do one thing!

Thanks for your great post. I have tried to replicate it using a model containing a spline term, but with the splines the matrix is non-conformable, so running line 22 gives an error. In contrast, at line 23 you mention that predict could substitute line 22, but predict works fine. I have reproduced an example using your data:

> library(splines)

> m mm newdat mm newdat$y<-mm%*%fixef(m)

Error in mm %*% fixef(m) : non-conformable arguments

Do you have any suggestion for how to overcome this problem?

Also, I echo Pankil's concern – when I run those two lines I get a prediction interval that is tiny and doesn't include the mean, so I think he is correct inasmuch as at least line 41 should perhaps be something different?

Thanks,

Sheena

I agree with the others, the apply command simply sorts the matrix rows, so the 1st and 2nd row would give the lowest bootstrapped values out of the 200 simulations. Pankil’s solution should work.

Cheers,

Steffen.

Hi Steffen,

Thanks for your comment, note that the command to create the bb_se matrix object is NOT:

bb_se<-apply(bb$t,2,function(x) x[order(x)])

Which would just sort the COLUMNS, it is rather

bb_se<-apply(bb$t,2,function(x) x[order(x)][c(5,195)])

Which in addition to sorting the columns take only the 5th and 195th values from each columns, representing a 95% bootstrapped confidence interval.

This is a SUPER helpful method – thank you! And thanks for the clarification in the bb_se matrix object specification – again, super helpful.

I’m trying to do this for a glmer with interactions that I have centered. I was able to specify it with the un-centered variables interacting, but can’t figure out how to work in the centered variables correctly when the un-centered variables are also still in the model.

When I get to line 21 I get the warning:

“Error in diag(mm %*% tcrossprod(vcov(f.sl.dea.m1.f.glmer), mm)) :

error in evaluating the argument ‘x’ in selecting a method for function ‘diag’: Error in .local(x, y, …) :

Dimensions of x and y are not compatible for tcrossprod”

I’m surprised by this, since it seemed to work for a model with an interaction in the example on http://glmm.wikidot.com/faq and when I just interacted the unaltered predictors. Adding ‘new’ centered variables seems to be the problem.

Any idea what the problem might be? Do you have an example with an interaction?

Thanks for your comment!

I must admit that I do not see why one would leave un-centered predictors alongside an interaction term with this predictor but then centered … I would rather either only use un-centered predictors (in both the main effect and the interaction) or only centered ones. But would be happy to be proven wrong on this.

Concerning the error you should check that dim(vcov(m) and dim(mm) both return the same number of columns.

Very helpful code indeed – thank you. I too am getting problems with lines 41 and 42 in the first chunk of code. Actually, the posted code itself does not quite reproduce the posted figure – in the figure the green lines for the bootstrapped CIs are close to each of their equivalent ‘normal’ CIs. But copying and pasting the code into R does not quite reproduce this – rather the two bootstrapped CI lines are low and close to each other, and both below the line of best fit (red line). So there does appear to be something amiss in the code with setting up newdat$blo and

newdat$bhi. Any thoughts on this would be much appreciated. Thank you, Lewis

Hi Lewis,

Thanks for your kind words.

Just re-ran the code on a newer machine with updated packages and I can reproduce the figure of the post. Could it be that you did not copy the whole line 39? I guess that you might have only got: x[order(x)], instead of x[order(x)][c(5,195)]. Let me know if your problem persist …

You are spot on! You quickly uncovered my deep ignorance, and I appreciate it! Great – with the updated code I’ve now got some impressive looking bootstrapped CIs.

Is it straightforward to adapt your code to handle a second continuous covariate? I.e.

m<-lmer(y~x+z+(1|f),data)

I am interested in plotting, with CIs, y against x controlling for z, and y against z controlling for x.

Great! I was starting to doubt about my own code … Adding more covariates is fairly straightforward, you have to make sure that in lines 19-20 you properly define the newdat object with some code like: newdat <- expand.grid(x=seq(0,10,length=20),z=mean(z)), and mm <- model.matrix(~x+z,newdat), after that the bootstrap will take care of itself! The basic idea is to vary one continuous covariates at a time (x) while maintaining all the others at their mean values (z), and this for each covariates. I actually implemented a function some time ago to automatize this, have alook here: http://wp.me/p2SacB-7V, maybe it'll work for you …

Thank you – this has worked like a dream. Sometimes the bootMer function grumbles (‘In checkConv(attr(opt, “derivs”), opt$par, ctrl = control$checkConv, … :

unable to evaluate scaled gradient’).

In such instances, is it best to stick to ‘standard’, non-boot-strapped CIs, as per your first chunk of code?

On a second issue – is it also straightforward to adapt your code for the following:

m<-lmer(y~x+z+p+(1|f),data), where x and z are covariates, but p is a fixed factor (i.e. rather than a covariate like x and z)? Or would you recommend running the model separately for each level of p, and in turn providing a plot with several regression lines (and CIs) – one for each level?

Thanks Lewis for your kind words. The error that you are getting also appeared in some of my analysis, I guess that as long as you have enough bootstrapped replicates and the warning does not appear that often (say less than 5% of the case) you should be fine. In the end it is a trade-off between how much confidence you have that your model log-likelihood function behave nicely (see ?profile for instance) vs the computational costs of running potentially many bootstrapped replicates so that replicates that failed to converge do not affect your estimate … If you are using lmer I would tend more towards using standard CIs, I tend to think that glmer models are a bit more tricky so in these case I would tend more towards bootstrapped ones, note that this is just gut-feeling!

For the second issue, if p is a fixed factor with more than one level you would include it as a factor variable in the model (check that str(data) indicate that p is factor) and then when defining the newdat object you would do something like: newdat <- expand.grid(x=seq(0,1,length=10),y=0,p=c("Group1","Group2")), the rest of the code stay the same, at the end the first 10 values would correspond to Group1 and the 10 next to Group2. You could browse through some other of my posts on plotting regression lines where I give some more infos.