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

[Viechtbauer, Wolfgang (SP)]

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