0
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ this has pulls out the relevant parts of a lmer model which may help $\endgroup$
    – user20650
    Jul 26, 2023 at 14:42

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.