Using the glmer() function in the LME4-library in R I computed logistic models, of the form: $Y ~ cat1 * cont1 + (1|Subject)$ where, obviously, Y is the binomial outcome variable (0 or 1), cat1 is a categorial variable (0,1,2) and cont1 is a continuous variable). Then, using confint(model, method = "boot") I computed confidence interval on the variables.

Now I would like to plot a graph of the chance P(Y==1), I want to plot P against cont1 for every cat1.

So you'd say: $X = B(0) + B(cont1) * cont1 + B(cat1:1) * (cat1==1) + .... +$ etc And then: $P(Y==1) = 1/(1-exp(-X))$

Which does exactly what I expect. But now, I want to incorporate the bootstrapped confidence intervals (so not std. error * 1.96!!) in the graph. I have the numbers, I do not know how to interpret them, what would be the formula for e.g. the 97.5 % line and the 2.5 % line?

Thanks in advance!

EDIT Is this the correct way? Basically taking 10000 samples with replacement, same size as original data, creating the model, computing the predictions, and taking the 97.5th and 2.5th intervals of the predictions.

prediction_pars = expand.grid(cont1= seq(-4,4,.05), cat1= as.factor(c(1,2,3)));

predictions = array(dim = c(10000, dim(prediction_pars)[1]));

for (i in 1:10000){
  new_sample = data[sample(nrow(data), samplesize, replace = T) , ];
  new_model  = glmer (Y ~ cont1 * cat1 + (1|Subject), dat=newdat, family="binomial");
  predictions[i , ] = predict(new_model, newdat = new_sample, re.form = NA);

hi = lo = array(dim = dim(prediction_pars)[1]);
for (i in (1, dim(prediction_pars)[1])){
  hi[i] = sort(predictions[,1]) [9750];
  lo[i] = sort(predictions[,1]) [ 250];
  • $\begingroup$ the natural thing to do (although it wouldn't account for bias correction) would be to do something like dd <- expand.grid(cont1=seq(...),cat1=unique(data$cat1)); confint(fitted_model,FUN=function(x) predict(x,newdata=dd)) $\endgroup$
    – Ben Bolker
    Commented Jun 7, 2016 at 0:22

1 Answer 1


This is an old problem without an efficient solution. As you know, confidence intervals and prediction intervals are very different things. Bootstrapped variance estimates for parameters will not give you robust prediction intervals. A theoretically correct approach would require you to iteratively bootstrap the data by hand, fit mixed models and obtain predictions, then you must calculate bias, and using a theoretically correct approach to obtain bootstrapped prediction interval estimates. This would probably take days to do.

Furthermore, the question of whom exactly you're predicting effects in needs to be addressed. Predictions can be obtained at any heirarchical level of correlation. Mixed models capture individual level effects. So you might be giving predictions for a particular individual (in which case the random effect is a fixed, non-zero value), or it could be marginalized to a population averaged result. If you consult the default predict.merMod method, you'll see the option re.form which requires you to tell R among exactly whom you're prediction.

Also this is an illuminating read:


  • $\begingroup$ Thank you, sorry for the late reaction. I did as follows: take 10000 samples with replacement, same size as original data, then for each sample: compute the model, predict the points on the line (with re.form=NA, as I want general predictions), store the points, and take the 97.5% and 2.5% intervals of the points. It takes around 30 hours. As for prediction v. confidence intervals, I want to visualise the area in the plot in which the model graph would lie with 95% certainty if I were to do the experiment again, I care not for individual measurements. Is this the correct way? $\endgroup$
    – Coen
    Commented Jun 15, 2016 at 11:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.