For research purposes, I need or would like to understand the EBLUP calculation of random effect.
Marginally, y follows :
$$y_i \sim N(X_i\beta + Z_iu_i, Z_iDZ_i^\top + \sigma^2_eI_{n_i})$$ and the random effects,
$$\textbf{u}_i \sim N(\boldsymbol{0}, D)$$ and they have a multivariate normal distribution with $DZ^\top$ being the covariance term between $b$ and $y$.
The conditional mode/mean in this case is then: $$E({b_i}|y_i) = DZ^\top(ZDZ^\top + \sigma^2_e I_{n_i})^{-1}(y-X\beta)$$
However, after writing up code, I realize I am not getting the same from the call ranef to the lmeObject.
#Get back Random effects manual calculation
library(JM)
fitLMEFULL <- lme(log(serBilir) ~ drug * year, random = ~ year | id, data = pbc2)
longitudinalDFFULL <-pbc2
N = length(unique(longitudinalDFFULL$id))
D <- getVarCov(fitLMEFULL)
formYz <- formula(fitLMEFULL$modelStruct$reStruct[[1]])
mfZ <- model.frame(terms(formYz), data = longitudinalDFFULL)
TermsZ <- attr(mfZ, "terms")
Z <- model.matrix(formYz, mfZ)
s <- fitLMEFULL$sigma
residualsM<- residuals(fitLMEFULL)
id <- fitLMEFULL$groups[[1]]
calculatedRANEFF <- ranef(fitLMEFULL)
ni <- as.vector(tapply(id, id, length))
calRanEff = dim(calculatedRANEFF)
checkcalRanEff <- matrix(0, ncol=calRanEff[2], nrow=calRanEff[1])
for(i in 1:N){
if(ni[i]==1){
#Inverse of (ZDZ^T + \sigma^2_e*I_n)
V_inv = solve(t(Z[id==i])%*%D%*%(Z[id==i]) + diag(s^2, ni[i]))
checkcalRanEff[i,] <- t(D%*%(Z[id==i])%*%V_inv%*%residualsM[id==i])
} else {
V_inv = solve(Z[id==i,]%*%D%*%t(Z[id==i,]) + diag(s^2, ni[i]))
checkcalRanEff[i,] <- t(D%*%t(Z[id==i,])%*%V_inv%*%residualsM[id==i])
}
}
checkcalRanEff
ranef(fitLMEFULL)
For some reason, mathjax is not rendering matrices as multiline so sorry for the poor math formatting on my part. Would appreciate any help on where my code is wrong is possible.
