[R-meta] Clarification on ranef.rma.mv()

[Luke Martinez]

Hi Wolfgang,

My goal was to follow lme4:::rePCA.merMod assuming that there is an
underlying reason that lme4 developers chose not to directly apply PCA
to the G matrix.

But if you think for rma.mv() models applying the PCA directly to the
G matrix gives the same result, I will very much appreciate a quick
demonstration?

Thank you again,
Luke

On Mon, Sep 13, 2021 at 11:42 AM Viechtbauer, Wolfgang (SP)
<wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>
> I am not familiar with lme4::rePCA, but if you want to do PCA on the var-cov matrix of the random effects, why not just apply PCA directly to G?
>
> Best,
> Wolfgang
>
> >-----Original Message-----
> >From: Luke Martinez [mailto:martinezlukerm using gmail.com]
> >Sent: Monday, 13 September, 2021 18:25
> >To: Viechtbauer, Wolfgang (SP)
> >Cc: R meta
> >Subject: Re: [R-meta] Clarification on ranef.rma.mv()
> >
> >Sure, I recently wrote a function similar to ?lme4::rePCA for lme
> >models. I keep benefiting from using it. So, I was wondering if I
> >could extend that to rma.mv() models?
> >
> >Here is how it works for lme models:
> >
> >#-- Some helper functions:
> >
> >Tlist_lme <- function(fit) rev(pdMatrix(fit$modelStruct$reStruct,
> >factor = TRUE))
> >theta_lme <- function(fit) sapply(Tlist_lme(fit), function(i)
> >i[lower.tri(i, diag = TRUE)])
> >
> >lowerbd <- function(x){
> >dd <- diag(0, nrow=nrow(x))
> >dd[lower.tri(dd)] <- -Inf
> >dd[lower.tri(dd, diag=TRUE)]
> >}
> >
> >lwr_lme <- function(fit) sapply(Tlist_lme(fit), lowerbd)
> >
> >#-- main function:
> >
> >rePCA_lme <- function(x){
> >
> >  chfs <- Tlist_lme(x)
> >  nms <- names(chfs)
> >  unms <- unique(nms)
> >  names(unms) <- unms
> >
> >  svals <- function(m) { # this is applied each of the RE matrices
> >    vv <- svd(m, nv = 0L)
> >    names(vv) <- c("sdev", "rotation")
> >    vv$center <- FALSE
> >    vv$scale <- FALSE
> >    class(vv) <- "prcomp"
> >    vv
> >  }
> >
> >  structure(lapply(unms, function(m)
> >    svals(Matrix::bdiag(chfs[which(nms ==
> >                                     m)]))), class = "prcomplist")
> >}
> >
> >#-- Testing:
> >
> >fm2 <- lme(distance ~ age,random=~age | Subject , data = Orthodont)
> >
> >rePCA_lme(fm2)
> >
> >On Mon, Sep 13, 2021 at 11:10 AM Viechtbauer, Wolfgang (SP)
> ><wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> >>
> >> This gives the scaled (by the residual variance) Cholesky decomposition of the
> >var-cov matrix of the random effects.
> >>
> >> We can see this if we compare:
> >>
> >> t(Tlist_lme(fm1)$Subject) %*% Tlist_lme(fm1)$Subject * sigma(fm1)^2
> >>
> >> with:
> >>
> >> VarCorr(fm1)
> >>
> >> rma.mv() also uses a Cholesky decomposition to ensure that the var-cov matrix of
> >the random effects is (semi)positive definite. However, it does not scale things
> >by the residual variance, since those variances are heteroscedastic.
> >>
> >> I don't know why you would want to obtain just Tlist_lme(fm1), but if you are
> >actually after the var-cov matrix itself, it is stored as "G" in the model object.
> >To use your earlier example:
> >>
> >> library(metafor)
> >> dat <- dat.berkey1998
> >> res.mv <- rma.mv(yi~ outcome - 1, vi, data = dat, random = ~ outcome | trial,
> >struct = "UN")
> >> res.mv$G
> >>
> >> If you really want the Cholesky decomposition, then just do:
> >>
> >> C <- chol(res.mv$G)
> >> C
> >>
> >> and, as expected, t(C) %*% C is the same as res.mv$G.
> >>
> >> Best,
> >> Wolfgang
> >>
> >> >-----Original Message-----
> >> >From: Luke Martinez [mailto:martinezlukerm using gmail.com]
> >> >Sent: Monday, 13 September, 2021 17:27
> >> >To: Viechtbauer, Wolfgang (SP)
> >> >Cc: R meta
> >> >Subject: Re: [R-meta] Clarification on ranef.rma.mv()
> >> >
> >> >Dear Wolfgang,
> >> >
> >> >Thank you very much. For correlated random-effect structures in
> >> >rma.mv(), is it possible to extract the equivalent of the following
> >> >lme() information:
> >> >
> >> >Tlist_lme <- function(fit) rev(pdMatrix(fit$modelStruct$reStruct,
> >> >factor = TRUE))
> >> >
> >> >#----- Testing:
> >> >library(nlme)
> >> >fm1 <- lme(distance ~ age, data = Orthodont)
> >> >
> >> >Tlist_lme(fm1)
> >> >
> >> >On Mon, Sep 13, 2021 at 2:38 AM Viechtbauer, Wolfgang (SP)
> >> ><wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> >> >>
> >> >> Dear Luke,
> >> >>
> >> >> Sure:
> >> >>
> >> >> library(metafor)
> >> >>
> >> >> dat <- dat.berkey1998
> >> >> res.mv <- rma.mv(yi~ outcome - 1, vi, data = dat, random = ~ outcome | trial,
> >> >struct = "UN")
> >> >> res.mv$rho
> >> >>
> >> >> ran.mv <- ranef.rma.mv(res.mv)
> >> >> cor(matrix(ran.mv[[1]]$intrcpt, ncol=2, byrow=FALSE))[1,2]
> >> >>
> >> >> The two correlations are not the same for the reasons explained at the link
> >you
> >> >provided.
> >> >>
> >> >> If the dataset is not so nicely balanced, one can do something similar, but
> >the
> >> >restructuring of the output from ranef() into a wide format gets a bit more
> >> >tedious.
> >> >>
> >> >> For example:
> >> >>
> >> >> dat <- dat[-4,]
> >> >> res.mv <- rma.mv(yi~ outcome - 1, vi, data = dat, random = ~ outcome | trial,
> >> >struct = "UN")
> >> >> res.mv$rho
> >> >>
> >> >> ran.mv <- ranef.rma.mv(res.mv)
> >> >> ran.mv
> >> >>
> >> >> ran.mv <- ran.mv[[1]][1]
> >> >> ran.mv$study <- sapply(strsplit(rownames(ran.mv), "|", fixed=TRUE), tail, 1)
> >> >> ran.mv$arm   <- sapply(strsplit(rownames(ran.mv), "|", fixed=TRUE), head, 1)
> >> >> ran.mv
> >> >>
> >> >> wide <- reshape(ran.mv, direction="wide", idvar="study", v.names="intrcpt",
> >> >timevar="arm")
> >> >> rownames(wide) <- NULL
> >> >> wide
> >> >> cor(wide[2:3], use="pairwise.complete.obs")[1,2]
> >> >>
> >> >> There might be a more elegant way to do this, but this gets the job done.
> >> >>
> >> >> Best,
> >> >> Wolfgang
> >> >>
> >> >> >-----Original Message-----
> >> >> >From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org]
> >On
> >> >> >Behalf Of Luke Martinez
> >> >> >Sent: Monday, 13 September, 2021 8:38
> >> >> >To: R meta
> >> >> >Subject: [R-meta] Clarification on ranef.rma.mv()
> >> >> >
> >> >> >Dear Wolfgang and List Members,
> >> >> >
> >> >> >In ordinary multilevel models (lmer), one can use ranef() of a model
> >> >> >to get the correlations from the conditional modes of the random
> >> >> >effects (https://stats.stackexchange.com/q/153253/124093) which may
> >> >> >reveal how much random-effects for slopes and intercepts are roughly
> >> >> >correlated.
> >> >> >
> >> >> >Similarly, for a struct = "UN" model, I was wondering if
> >> >> >"ranef.rma.mv(fitted_model)" could reveal how much random-effects for
> >> >> >outcome levels are roughly correlated?
> >> >> >
> >> >> >For example, for the "res.mv" model below where the REML rho estimate
> >> >> >is "0.775", can ranef.rma.mv(res.mv) values give a rough estimate of
> >> >> >this correlation conditional on the observed data?
> >> >> >
> >> >> >dat <- dat.berkey1998
> >> >> >res.mv <- rma.mv(yi~ outcome - 1, vi, data = dat, random = ~ outcome |
> >> >> >trial, struct = "UN")
> >> >> >
> >> >> >ran.mv <- ranef.rma.mv(res.mv)
> >> >> >
> >> >> >Thank you very much for your help,
> >> >> >Luke