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

[Luke Martinez]

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