[R-sig-ME] Eqn. (5.22)

[Andrzej Galecki]

Dear Mixed List:

I am trying to follow eqn (5.22) in Chapter 5 of  "lme4: Mixed-effects 
modeling with R" book by Douglas Bates
and cross-check the results from a linear mixed effects model  fit.

Specifically, in the code below I extract pieces of information from 
object "fm2" and compute
expressions on the left and right hand side of (5.22). For some reason 
left side is (slightly) different from right hand side.
What am I doing wrong? Same with object "fm2a".

Thank you

Andrzej Galecki


# http://lme4.r-forge.r-project.org/book/Ch5.pdf
# Eqn. (5.22) on p.89

library(lme4)
# sessionInfo()    # see postscript later in this email

fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy)
# fm2a <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) 

fm <- fm2

###### Right hand side of (5.20)
Zt <- slot(fm,"Zt")      # Z transposed
# image(Zt)
STs  <- expand(fm)       # Expand ST slot

#summary(STs)

P    <- STs$P  # permutation matrix P

S    <- STs$S  # Diagonal scale matrix S
# summary(S)
T    <- STs$T  # Unit lower-triangular matrix T
#summary(T)     # All off-diagonal equal to 0      

Lmbda <- T %*% S     #  Formula w/out number on p.84      
U <- t(Zt) %*% Lmbda
U1 <-   (t(U) %*% U + Diagonal(ncol(U)))
U1P <- P %*% U1 %*% t(P)
#image(U1P)

######  Left hand side of (5.20)
Ls  <- slot(fm,"L") 
#image(Ls)
L  <- as(Ls, "sparseMatrix")
LtL <- tcrossprod(L)

# Compare left and right hand sides of (5.20)

max(abs(LtL -U1P ))    # LtL not equal to U1P
# [1] 0.003046723

LtL[1:2,1:2]                 # Block extracted
#2 x 2 sparse Matrix of class "dsCMatrix"
#                     
# [1,] 10.60194 10.32746
# [2,] 10.32746 16.63319

        
U1P[1:2,1:2]             # Block
# 2 x 2 sparse Matrix of class "dgCMatrix"
                     
# [1,] 10.60194 10.32846
# [2,] 10.32846 16.63624



PS.

######## sessionInfo ######  

 > sessionInfo()
R version 2.10.1 (2009-12-14)
i386-pc-mingw32

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252  
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                         
[5] LC_TIME=English_United States.1252   

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base    

other attached packages:
[1] lme4_0.999375-32   Matrix_0.999375-33 lattice_0.17-26  

loaded via a namespace (and not attached):
[1] grid_2.10.1 nlme_3.1-96