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

[Andrzej Galecki]

Unfortunately, two lines of my code were merged  together. They are 
separated in this email.

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       # This line was  lumped  with previous one
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