Created
March 15, 2013 20:46
BLUPs
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
> # Calculate the predicted random effects by hand for the ergoStool data | |
> (fm1<-lmer(effort~Type-1 + (1|Subject),ergoStool)) | |
Linear mixed model fit by REML | |
Formula: effort ~ Type - 1 + (1 | Subject) | |
Data: ergoStool | |
AIC BIC logLik deviance REMLdev | |
133 143 -60.6 122 121 | |
Random effects: | |
Groups Name Variance Std.Dev. | |
Subject (Intercept) 1.78 1.33 | |
Residual 1.21 1.10 | |
Number of obs: 36, groups: Subject, 9 | |
Fixed effects: | |
Estimate Std. Error t value | |
TypeT1 8.556 0.576 14.8 | |
TypeT2 12.444 0.576 21.6 | |
TypeT3 10.778 0.576 18.7 | |
TypeT4 9.222 0.576 16.0 | |
Correlation of Fixed Effects: | |
TypeT1 TypeT2 TypeT3 | |
TypeT2 0.595 | |
TypeT3 0.595 0.595 | |
TypeT4 0.595 0.595 0.595 | |
> | |
> ## Here are the BLUPs we will estimate by hand: | |
> ranef(fm1) | |
$Subject | |
(Intercept) | |
1 1.7088e+00 | |
2 1.7088e+00 | |
3 4.2720e-01 | |
4 -8.5439e-01 | |
5 -1.4952e+00 | |
6 -1.3546e-14 | |
7 4.2720e-01 | |
8 -1.7088e+00 | |
9 -2.1360e-01 | |
> | |
> ## this gives us all the variance components: | |
> VarCorr(fm1) | |
$Subject | |
(Intercept) | |
(Intercept) 1.7755 | |
attr(,"stddev") | |
(Intercept) | |
1.3325 | |
attr(,"correlation") | |
(Intercept) | |
(Intercept) 1 | |
attr(,"sc") | |
[1] 1.1003 | |
> | |
> # First, calculate the predicted random effect for subject 1: | |
> | |
> ## The variance for the random effect subject is the term C_{u,Y}: | |
> covar.u.y<-VarCorr(fm1)$Subject[1] | |
> | |
> | |
> # Estimated covariance between u_1 and Y_1 | |
> ## make up a var-covar matrix from this: | |
> (cov.u.Y<-matrix(covar.u.y,1,4)) | |
[,1] [,2] [,3] [,4] | |
[1,] 1.7755 1.7755 1.7755 1.7755 | |
> | |
> | |
> # Estimated variance matrix for Y_1 | |
> (V.Y<-matrix(1.7755,4,4)+diag(1.2106,4,4)) | |
[,1] [,2] [,3] [,4] | |
[1,] 2.9861 1.7755 1.7755 1.7755 | |
[2,] 1.7755 2.9861 1.7755 1.7755 | |
[3,] 1.7755 1.7755 2.9861 1.7755 | |
[4,] 1.7755 1.7755 1.7755 2.9861 | |
> | |
> # Extract observations for subject 1 | |
> (Y<-matrix(ergoStool$effort[1:4],4,1)) | |
[,1] | |
[1,] 12 | |
[2,] 15 | |
[3,] 12 | |
[4,] 10 | |
> | |
> # Estimated fixed effects | |
> (beta.hat<-matrix(fixef(fm1),4,1)) | |
[,1] | |
[1,] 8.5556 | |
[2,] 12.4444 | |
[3,] 10.7778 | |
[4,] 9.2222 | |
> | |
> # Predicted random effect | |
> cov.u.Y %*% solve(V.Y)%*%(Y-beta.hat) | |
[,1] | |
[1,] 1.7087 | |
> | |
> # Compare with ranef command | |
> ranef(fm1)$Subject[1,1] | |
[1] 1.7088 | |
> | |
> # Calculate predicted random effects for all subjects | |
> t(cov.u.Y %*% solve(V.Y)%*%(matrix(ergoStool$effort,4,9)-matrix(fixef(fm1),4,9))) | |
[,1] | |
[1,] 1.7087e+00 | |
[2,] 1.7087e+00 | |
[3,] 4.2717e-01 | |
[4,] -8.5435e-01 | |
[5,] -1.4951e+00 | |
[6,] -1.3906e-14 | |
[7,] 4.2717e-01 | |
[8,] -1.7087e+00 | |
[9,] -2.1359e-01 | |
> ranef(fm1) | |
$Subject | |
(Intercept) | |
1 1.7088e+00 | |
2 1.7088e+00 | |
3 4.2720e-01 | |
4 -8.5439e-01 | |
5 -1.4952e+00 | |
6 -1.3546e-14 | |
7 4.2720e-01 | |
8 -1.7088e+00 | |
9 -2.1360e-01 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment