[R] Confidence Intervals for Random Effect BLUP's

[Douglas Bates]

On Nov 9, 2007 1:15 PM, Rick Bilonick <rab at nauticom.net> wrote:
> On Fri, 2007-11-09 at 18:55 +0000, Prof Brian Ripley wrote:
>
> > I think Bert's point is important: I picked up a student on it in a case
> > study presentation on this week because I could think of three
> > interpretations, none strictly confidence intervals.  I think 'tolerance
> > interval' is fairly standard for prediction of a random quantity: see
> > ?predict.lm.
> >
>
> I think prediction interval is what is usually used. Regardless, I'm not
> sure how "predict.lm" will be of much help because I asked specifically
> about BLUP's for random effects and the last time I checked lm did not
> handle mixed effects models. Neither predict.lme and predict.lmer
> provide intervals. Here is the code that I included in my original
> e-mail. My simple question is, will this code correctly compute a
> prediction interval for each subjects random effect? In particular, will
> the code handle the bVar slot correctly? Some postings warned about
> inappropriate access to slots. Here is the code that I asked about in my
> original e-mail:

Regarding the terminology, I prefer to call the quantities that are
returned by the ranef extractor "the conditional modes of the random
effects". If you want to be precise, these are the conditional modes
(for a linear mixed model they are also the conditional means) of the
random effects B given Y = y, evaluated at the parameter estimates.
One can also evaluate the condtional variance-covariance of B given Y
= y and hence obtain a prediction interval.  These are returned by
ranef when the optional argument "postVar" is TRUE.  I think the
intervals you want are what are shown in the "caterpillar" plots.  Try

library(lme4)
qqmath(ranef(fm1 <- lmer(Reaction ~ Days + (Days|Subject),
sleepstudy), postVar = TRUE))

BTW, the reason that I say "conditional modes", rather than
"conditional means", is so the term can apply to generalized linear
mixed models, nonlinear mixed models and generalized nonlinear mixed
models.  The conditional distribution of B given Y is a continuous
multivariate distribution that can be evaluated explicitly for a
linear mixed model but is known only up to a constant for the
generalized forms of the model.  We can still optimize the conditional
density of B given Y = y for particular values of the parameters but
doing so provides the conditional modes, not the conditional means.



>
> # OrthoFem has all the females from Orthodont from the nlme package
>
>
> library(lme4)
> fm1OrthF. <- lmer(distance~age+(age|Subject), data=OrthoFem)
>
> lmer(distance~age+(age|Subject), data=OrthoFem)@bVar$Subject[2,2,]*
>   (attr(VarCorr(lmer(distance~age+(age|
> Subject),data=OrthoFem)),"sc")^2)[1]
>
> (attr(VarCorr(fm1OrthF.),"sc")^2)[1]
>
> fm1.s <- coef(fm1OrthF.)$Subject
> fm1.s.var <- fm1OrthF. at bVar$Subject*(attr(VarCorr(fm1OrthF.),"sc")^2)[1]
> fm1.s0.s <- sqrt(fm1.s.var[1,1,])
> fm1.s0.a <- sqrt(fm1.s.var[2,2,])
> fm1.s[,1]+outer(fm1.s0.s, c(-2,0,2))
> fm1.s[,2]+outer(fm1.s0.a, c(-2,0,2))
>
> > fm1.s
>     (Intercept)       age
> F10    14.48493 0.3758608
> F09    17.26499 0.3529804
> F06    16.77328 0.3986699
> F01    16.95609 0.4041058
> F05    18.36188 0.3855955
> F07    17.28390 0.5193954
> F02    16.05461 0.6336191
> F08    19.40204 0.3562135
> F03    16.35720 0.6727714
> F04    19.02380 0.5258971
> F11    19.13726 0.6498911
>
> > fm1.s[,1]+outer(fm1.s0.s, c(-2,0,2))
>           [,1]     [,2]     [,3]
>  [1,] 12.21371 14.48493 16.75616
>  [2,] 14.99377 17.26499 19.53622
>  [3,] 14.50205 16.77328 19.04450
>  [4,] 14.68487 16.95609 19.22732
>  [5,] 16.09066 18.36188 20.63311
>  [6,] 15.01267 17.28390 19.55512
>  [7,] 13.78339 16.05461 18.32584
>  [8,] 17.13082 19.40204 21.67327
>  [9,] 14.08598 16.35720 18.62843
> [10,] 16.75257 19.02380 21.29502
> [11,] 16.86604 19.13726 21.40849
>
> > fm1.s[,2]+outer(fm1.s0.a, c(-2,0,2))
>            [,1]      [,2]      [,3]
>  [1,] 0.1738325 0.3758608 0.5778890
>  [2,] 0.1509522 0.3529804 0.5550087
>  [3,] 0.1966417 0.3986699 0.6006982
>  [4,] 0.2020775 0.4041058 0.6061340
>  [5,] 0.1835672 0.3855955 0.5876237
>  [6,] 0.3173671 0.5193954 0.7214236
>  [7,] 0.4315909 0.6336191 0.8356474
>  [8,] 0.1541852 0.3562135 0.5582417
>  [9,] 0.4707432 0.6727714 0.8747997
> [10,] 0.3238688 0.5258971 0.7279253
> [11,] 0.4478629 0.6498911 0.8519194
>
>
> This web page describes "bVar":
>
> http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/lme4/html/lmer-class.html
>
> bVar:
>         A list of the diagonal inner blocks (upper triangles only) of
>         the positive-definite matrices on the diagonal of the inverse of
>         ZtZ+Omega. With the appropriate scale factor (and conversion to
>         a symmetric matrix) these are the conditional
>         variance-covariance matrices of the random effects.
>
> Rick B.
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>