# R, glmer(), plot bootstrapped CI in graph

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?

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

• 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)) – Ben Bolker Jun 7 '16 at 0:22

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.