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

[Viechtbauer, Wolfgang (SP)]

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