[R] Confidence Intervals for Random Effect BLUP's

Rick Bilonick rab at nauticom.net
Fri Nov 9 18:46:35 CET 2007


I want to compute confidence intervals for the random effect estimates
for each subject. From checking on postings, this is what I cobbled
together using Orthodont data.frame as an example. There was some
discussion of how to properly access lmer slots and bVar, but I'm not
sure I understood. Is the approach shown below correct?

Rick B.

# Orthodont is from nlme (can't have both nlme and lme4 loaded at same
time!)
# OrthoFem<-Orthodont[Orthodont$Sex=="Female",] 
# http://tolstoy.newcastle.edu.au/R/help/06/03/23758.html

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



More information about the R-help mailing list