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



Correlations of fixed effects in linear mixed models

Ever wondered what those correlations are in a linear mixed model? For example:

> (lm.full<-lmer(wear~material-1+(1|Subject), data = BHHshoes))
Linear mixed model fit by REML
Formula: wear ~ material - 1 + (1 | Subject)
Data: BHHshoes
AIC BIC logLik deviance REMLdev
62.9 66.9 -27.5 53.8 54.9
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 6.1009 2.470
Residual 0.0749 0.274
Number of obs: 20, groups: Subject, 10
Fixed effects:
Estimate Std. Error t value
materialA 10.630 0.786 13.5
materialB 11.040 0.786 14.1
Correlation of Fixed Effects:
matrlA
materialB 0.988

The estimated correlation between \hat{\beta}_1 and \hat{\beta}_2 is 0.988.  Note that

\hat{\beta}_1 = (Y_{1,1} + Y_{2,1} + \dots + Y_{10,1})/10=10.360

and 

\hat{\beta}_2 = (Y_{1,2} + Y_{2,2} + \dots + Y_{10,2})/10 = 11.040

From this we can recover the correlation 0.988 as follows:

> b1.vals<-subset(BHHshoes,material=="A")$wear
> b2.vals<-subset(BHHshoes,material=="B")$wear
>
> vcovmatrix<-var(cbind(b1.vals,b2.vals))
>
> covar<-vcovmatrix[1,2]
> sds<-sqrt(diag(vcovmatrix))
> covar/(sds[1]*sds[2])
b1.vals
0.98823

By comparison, in the linear model version of the above:

> summary(lm<-lm(wear~material-1,BHHshoes))
> X<-model.matrix(lm)
> 2.49^2*solve(t(X)%*%X)
materialA materialB
materialA 0.62001 0.00000
materialB 0.00000 0.62001

because Var(\hat{\beta}) = \hat{\sigma}^2 (X^T X)^{-1}.