I'm running a mixed model on some data. I want to calculate confidence intervals for my model.

For this I have adapted the following code section from Predictions and/or confidence (or prediction) intervals on predictions (lme4). The problem is that I have a hard time understanding what the code actually does, but so far it's the only way I have found to calculate some reasonable confidence intervals for my result.

So my question is: what are pvar1 and tvar1?

I'm guessing it is something like $\text{SE}^2$, but I can't quite see it, and I can't describe their difference, or their normal statistical term. I'm guessing that the difference is one is CI while the other is PI?

Further when the code then calculates the confidence intervals it says prediction +/- 2*sqrt(pvar1). Shouldn't it be the t.crit value or 1.96 for a normal distribution?

library(ggplot2) # Plotting
fm1 <- lmer(
    formula = distance ~ age*Sex + (age|Subject)
    , data = Orthodont
newdat <- expand.grid(
    , Sex=c("Male","Female")
    , distance = 0
mm <- model.matrix(terms(fm1),newdat)
newdat$distance <- mm %*% fixef(fm1)
pvar1 <- diag(mm %*% tcrossprod(vcov(fm1),mm))
tvar1 <- pvar1+VarCorr(fm1)$Subject[1]  ## must be adapted for more complex models
newdat <- data.frame(
    , plo = newdat$distance-2*sqrt(pvar1)
    , phi = newdat$distance+2*sqrt(pvar1)
    , tlo = newdat$distance-2*sqrt(tvar1)
    , thi = newdat$distance+2*sqrt(tvar1)
#plot confidence
g0 <- ggplot(newdat, aes(x=age, y=distance, colour=Sex))+geom_point()
g0 + geom_errorbar(aes(ymin = plo, ymax = phi))+
    labs(title="CI based on fixed-effects uncertainty ONLY")
#plot prediction
g0 + geom_errorbar(aes(ymin = tlo, ymax = thi))+
    labs(title="CI based on FE uncertainty + RE variance")

Thank you for your help.

  • 3
    $\begingroup$ The value '2' is a rough t critical value; I'm inclined to agree, I'd tend to use either an actual $t$ critical value or 1.96. At first glance, the first variance estimates seem to relate to the contribution of the fixed effect parameter uncertainty to the variance of the predicted value and the second looks to be sum of that and the contribution the random effect variance. $\endgroup$
    – Glen_b
    May 30 '13 at 15:31
  • 1
    $\begingroup$ You're correct on all counts. The justification for using 2 instead of a more precise value is to emphasize that these are approximate values -- e.g. they don't account at all for uncertainty in the random effects parameters. $\endgroup$
    – Ben Bolker
    May 30 '13 at 21:42
  • $\begingroup$ OK. Thank you. Because my next question was going to be on how many degrees of freedom should one find the t.crit-value. I wonder if it is the number of data-points or the number of subjects. $\endgroup$ Jun 2 '13 at 15:21
  • $\begingroup$ Hi,did you modified the above code to get predictions ci for a glmer model? I fitted a glmer model (logistic) and I used this the code, but what I'm getting doesn't make sense Predictions and ci are far over 1 and under 0. I can't figure out where in the code I have to specify that my model is a logistic one. Any advice. Also, I'm wondering if you modify the code and you are willing to share it. tanks! $\endgroup$
    – Rafael
    Jun 4 '14 at 8:15
  • $\begingroup$ would such an approach work with glmer as well as lmer? $\endgroup$
    – user553480
    Jun 22 '20 at 15:01

the code then calculates the confidence intervals it says prediction +/- 2*sqrt(pvar1). Shouldn't it be the t.crit value or 1.96 for a normal distribution?

This is correct. 1.96 is a more accurate critical value, however, statistical inference in mixed models is plagued with problems due to the presence of the random effects. It is difficult to quantify the uncertainty in the random effects.

2 is often use, as Ben Bolker points out in the comments, in order to empasize that these are approximate values.

Accordingly, any p values or confidence intervals should be considered approximate.

  • $\begingroup$ Does this answer your question ? If so, please consider marking it as the accepted answer, or if not please let us know why so that it can be improved. $\endgroup$ Nov 10 '20 at 18:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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