Scaling behavior ov bVar from lmer models
Locked
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.@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.@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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html |
Re: Scaling behavior ov bVar from lmer models
|
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: [hidden email] [mailto:[hidden email]] On Behalf Of Alan Bergland Sent: Tuesday, March 21, 2006 6:31 AM To: [hidden email] 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.@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.@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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html |
Re: Scaling behavior ov bVar from lmer models
|
|
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: [hidden email] > [mailto:[hidden email]] On Behalf Of Alan Bergland > Sent: Tuesday, March 21, 2006 6:31 AM > To: [hidden email] > 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.@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.@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 > > ______________________________________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html |
«
Return to R help
|
1 view|%1 views
Loading...
| Powered by Nabble | Edit this page |