In linear mixed models, we fit models like these (the Ware-Laird formulation--see Pinheiro and Bates 2000, for example):
\begin{equation} Y = X\beta + Zu + \epsilon \end{equation}
Let u\sim N(0,\sigma_u^2), and this is independent from \epsilon\sim N(0,\sigma^2).
Given Y, the ``minimum mean square error predictor'' of u is the conditional expectation:
\begin{equation} \hat{u} = E(u\mid Y) \end{equation}
We can find E(u\mid Y) as follows. We write the joint distribution of Y and u as:
\begin{equation} \begin{pmatrix} Y \\ u \end{pmatrix} = N\left( \begin{pmatrix} X\beta\\ 0 \end{pmatrix}, \begin{pmatrix} V_Y & C_{Y,u}\\ C_{u,Y} & V_u \\ \end{pmatrix} \right) \end{equation}
V_Y, C_{Y,u}, C_{u,Y}, V_u are the various variance-covariance matrices.
It is a fact (need to track this down) that
\begin{equation} u\mid Y \sim N(C_{u,Y}V_Y^{-1}(Y-X\beta)), Y_u - C_{u,Y} V_Y^{-1} C_{Y,u}) \end{equation}
This apparently allows you to derive the BLUPs:
\begin{equation} \hat{u}= C_{u,Y}V_Y^{-1}(Y-X\beta)) \end{equation}
Substituting \hat{\beta} for \beta, we get:
\begin{equation} BLUP(u)= \hat{u}(\hat{\beta})C_{u,Y}V_Y^{-1}(Y-X\hat{\beta})) \end{equation}
Here is a working example:
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 |
1 comment:
I winder how will u fit repeat measurements on subjects in the model? for example u have three repeats fro each kin for each subject. Thanks
Post a Comment