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

view raw corr_fixedeff_example1.R hosted with ❤ by GitHub

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

view raw corr_fixedeff_example2.R hosted with ❤ by GitHub

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

view raw corr_fixedef_lm.R hosted with ❤ by GitHub

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