BLUPs

blups.R

> # 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