Transform Z matrix

My goal is to incorporate correlation structure between the random effects which is by default not supported by rstanarm/lme4. The idea is to transform the Z matrix in the linear predictor X\beta + Zb such that it accounts for the correlation structure.

I have slightly adapted the stan_glmer() source code by adding one line of code at line 47 (here you find full function code)

glmod$reTrms$Ztlist[[1]] <- Lt %*% glmod$reTrms$Ztlist[[1]]

Lt is the matrix which transforms Z.

This method works perfectly fine with lme4. However, with stan_glmer() the modes of the random effect b are completely wrong.

Interestingly, already a very small change in Z leads to completely different random effects. Below a small example which shows this unexpected behavior

library(rstanarm)
library(lme4)
# load adapted stan_glmer() function
source("https://raw.githubusercontent.com/retodomax/cowfit/master/R/stan_adapted.R")

Lt <- diag(length(levels(lme4::sleepstudy$Subject)))
Lt[1,2] <- 0.01    # The transformed Z matrix will
                   # be almost the same as the original

fit_lme4 <- lmer(Reaction ~ (1|Subject), data = lme4::sleepstudy)
fit_stan <- stan_glmer(Reaction ~ (1|Subject), data = lme4::sleepstudy)
fit_adapt <- stan_adapt(Reaction ~ (1|Subject), data = lme4::sleepstudy)

re <- data.frame(lme4 = ranef(fit_lme4)$Subject[[1]],
                 stan = ranef(fit_stan)$Subject[[1]],
                 adapt = ranef(fit_adapt)$Subject[[1]])
pairs(re)

Do you have an idea what could be the problem? Thank you

Sorry, it looks like your question fell through… Overall @jonah knows much more about internals of rstanarm than I do.

I think that what you aim for should be possible without hacks in brms, where you can use (1 | gr(Subject, cov = your_covariance_matrix) to force a covariance matrix (this was introduced in brms to model phylogenetic covariance, see https://cran.r-project.org/web/packages/brms/vignettes/brms_phylogenetics.html for some more info).

2 Likes

Thanks for your reply! Yes, in the meantime I also found brms which mostly serves what I need.

Maybe a small suggestion to the brms developer (@paul.buerkner). For some applications in breeding, it is computationally expensive to calculate the covariance matrix but much simpler to calculate the Cholesky factor. Maybe an additional argument to gr() would be useful in order to directly specify the Cholesky factor e.g. brm(y ~ (1|Subject, chol = your_Cholesky_factor)).

I have implemented some hackes such that this is possible (find code here) but a solid implementation in brms would be much better. However, I am not sure if there are many people which would use such a feature, so maybe it is not top priority.

Thanks @martinmodrak, I missed this question.

@retodomax I’m not sure why modifying stan_glmer() in this way doesn’t work as expected. Unfortunately I haven’t had time to look into it yet, but if I can find some time I will. It’s possible @bgoodri will know immediately, although I know he’s super busy with the RStan release at the moment.

I just cli

This link is giving me 404: Not Found. Is there an updated link?

I’m glad you found brms. And I agree it often makes sense to calculate the Cholesky factor of the covariance matrix. I would think that’s the case in many applications, so it could be worth opening an issue on Github with a feature request for this.

2 Likes

In stan_glmer, there isn’t really a \mathbf{Z} \mathbf{b}; instead there is a \mathbf{Z}\mathbf{L} \mathbf{u} with the elements of \mathbf{u} having iid standard normal priors. Aside from \mathbf{L} being unknown, I think stan_glmer is already doing what you are trying to accomplish.

@jonah sorry I broke the link, now it should work again…

After a while, I tracked the problem down to the first line in the file ‘model/eta_add_Zb.stan’

if (special_case) for (i in 1:t) eta += b[V[i]];

It seems that if there is a special_case (only one random term of the form (1|factor), the Z matrix is completely ignored and random effects are directly assigned to observations according to V. V is created with the function make_V() which does not take any information about values in the Z matrix (usually it does not have to if we only fit a random intercept for each factor level).

rstanarm seems to do the sampling very efficiently but brms seems to be easier to configure in order to transform the Z matrix in a special user defined way .

Anyway, thanks for your help solving this problem!

1 Like