[R] Scaling behavior ov bVar from lmer models

Tue Mar 21 16:15:58 CET 2006

Hi Harold,

Ahh, yes, that makes things scale appropriately.

 > lmer(distance~age+(age|Subject), data=OrthoFem)@bVar$Subject[2,2,]* 
(attr(VarCorr(lmer(distance~age+(age|Subject), data=OrthoFem)), "sc")^2)
[1] 0.01020546 0.01020546 0.01020546 0.01020546 0.01020546 0.01020546  
0.01020546 0.01020546 0.01020546 0.01020546 0.01020546

OrthoFem$distance2<-OrthoFem$Distance*100
 > lmer(distance2~age+(age|Subject), data=OrthoFem)@bVar$Subject[2,2,] 
*(attr(VarCorr(lmer(distance2~age+(age|Subject), data=OrthoFem)),  
"sc")^2)
[1] 102.0546 102.0546 102.0546 102.0546 102.0546 102.0546 102.0546  
102.0546 102.0546 102.0546 102.0546

Thanks,
Alan

On Mar 21, 2006, at 8:14 AM, Doran, Harold wrote:

> 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
>