[R] Scaling behavior ov bVar from lmer models

[Doran, Harold]

Alan:

I think you need to multiply the values in the bVar slot by th residual
variance to get the correct estimates.  

-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Alan Bergland
Sent: Tuesday, March 21, 2006 6:31 AM
To: r-help at stat.math.ethz.ch
Subject: [R] Scaling behavior ov bVar from lmer models

Hi all,

    To follow up on an older thread, it was suggested that the following
would produce confidence intervals for the estimated BLUPs from a linear
mixed effect model:


OrthoFem<-Orthodont[Orthodont$Sex=="Female",]
fm1OrthF. <- lmer(distance~age+(age|Subject), data=OrthoFem)

fm1.s <- coef(fm1OrthF.)$Subject
fm1.s.var <- fm1OrthF. at bVar$Subject
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))

        [,1]     [,2]     [,3]
 [1,] 11.08578 14.48469 17.88359
 [2,] 13.86631 17.26521 20.66412
 [3,] 13.37444 16.77334 20.17224
 [4,] 13.55727 16.95617 20.35508
 [5,] 14.96331 18.36221 21.76112
 [6,] 13.88490 17.28381 20.68271
 [7,] 12.65522 16.05412 19.45303
 [8,] 16.00368 19.40258 22.80149
 [9,] 12.95778 16.35669 19.75559
[10,] 15.62506 19.02396 22.42287
[11,] 15.73831 19.13721 22.53612

 > fm1.s[,2]+outer(fm1.s0.a, c(-2,0,2))
            [,1]      [,2]      [,3]
 [1,] 0.07354181 0.3758789 0.6782160
 [2,] 0.05062219 0.3529593 0.6552964
 [3,] 0.09632606 0.3986632 0.7010003
 [4,] 0.10176055 0.4040977 0.7064348
 [5,] 0.08322913 0.3855662 0.6879033
 [6,] 0.21706649 0.5194036 0.8217407
 [7,] 0.33132614 0.6336632 0.9360004
 [8,] 0.05382874 0.3561658 0.6585029
 [9,] 0.37048154 0.6728186 0.9751558
[10,] 0.22354726 0.5258844 0.8282215
[11,] 0.34756193 0.6498990 0.9522361
 >

These matrices contain in the middle column the BLUPs for each female
individual in the study, and on the left and right columns the upper and
lower confidence intervals.


However, when the response variable is scaled up by a factor of 100, the
CIs do not scale accordingly.  Recall that the variance, SE, and CI
scale with the magnitude of the variable at hand.  I.e.,

x<-c(1,2,3,4,5)
var(x)
 >2.5

x2<-x*10
var(x2)
 >250


So, to rerun the above example with "distance" scaled up by a factor of
100:

OrthFem$distance2<-OrthoFen$distance*100
fm1OrthF. <- lmer(distance2~age+(age|Subject), data=OrthoFem)

fm1.s <- coef(fm1OrthF.)$Subject
fm1.s.var <- fm1OrthF. at bVar$Subject
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[,1]+outer(fm1.s0.s, c(-2,0,2))
          [,1]     [,2]     [,3]
 [1,] 1445.070 1448.469 1451.868
 [2,] 1723.122 1726.521 1729.920
 [3,] 1673.935 1677.334 1680.733
 [4,] 1692.218 1695.617 1699.016
 [5,] 1832.822 1836.221 1839.620
 [6,] 1724.982 1728.381 1731.780
 [7,] 1602.014 1605.412 1608.811
 [8,] 1936.859 1940.258 1943.657
 [9,] 1632.270 1635.669 1639.068
[10,] 1898.997 1902.396 1905.795
[11,] 1910.322 1913.721 1917.120
 > fm1.s[,2]+outer(fm1.s0.a, c(-2,0,2))
          [,1]     [,2]     [,3]
 [1,] 37.28555 37.58789 37.89023
 [2,] 34.99359 35.29593 35.59827
 [3,] 39.56398 39.86632 40.16865
 [4,] 40.10743 40.40977 40.71210
 [5,] 38.25429 38.55662 38.85896
 [6,] 51.63802 51.94036 52.24270
 [7,] 63.06399 63.36632 63.66866
 [8,] 35.31425 35.61658 35.91892
 [9,] 66.97953 67.28186 67.58420
[10,] 52.28610 52.58844 52.89077
[11,] 64.68757 64.98990 65.29224

Note that the BLUPs scaled accordingly, but the CIs do not.  For
example, take fm1.s[,2][11,] from both examples:

##unscaled
[11,] 0.34756193 0.6498990 0.9522361


##scaled
[11,] 64.68757 64.98990 65.29224

In both cases the upper CI is ~0.3 greater than the BLUP.


This seems very strange to me.  Am I totally misusing the bVar slot?  If
so, is there a way to obtain variances, or SE's associated with each
BLUP?

Thanks,
Alan

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html