R, glmer(), plot bootstrapped CI in graph

Question

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];
}

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:

https://cran.r-project.org/web/packages/merTools/README.html