Loading web-font TeX/Math/Italic

Search

Friday, March 15, 2013

How are the random effects (BLUPs) `predicted' in linear mixed models?




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:




> # 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
view raw blups.R hosted with ❤ by GitHub



1 comment:

Fariba Izadi said...

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