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

[Viechtbauer, Wolfgang (SP)]

Going back to the earlier example: 

library(nlme)
fm2 <- lme(distance ~ age, random = ~ age | Subject, data = Orthodont)
rePCA_lme(fm2)

This is essentially:

svd(chol(as.matrix(getVarCov(fm2))))$d / sigma(fm2)
svd(chol(as.matrix(getVarCov(fm2))))$u

Doing the same with rma.mv():

library(metafor)
Orthodont$id <- 1:nrow(Orthodont)
fm3 <- rma.mv(distance ~ age, 0, random = list(~ age | Subject, ~ 1 | id), struct="GEN", data = Orthodont)
svd(chol(fm3$G))$d / sqrt(fm3$sigma2)
svd(chol(fm3$G))$u

No idea how to interpret this and whether the idea underlying this generalizes to the case where each yi value has its own variance (instead of estimating a single residual variance).

Best,
Wolfgang

>-----Original Message-----
>From: Luke Martinez [mailto:martinezlukerm using gmail.com]
>Sent: Monday, 13 September, 2021 19:49
>To: Viechtbauer, Wolfgang (SP)
>Cc: R meta
>Subject: Re: [R-meta] Clarification on ranef.rma.mv()
>
>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