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

You can just do:

sds <- svd(chol(fm3$G))$d / sqrt(fm3$sigma2)
sds^2 / sum(sds^2)

to get those 'Proportion of Variance' values. In fact, the scaling by sqrt(fm3$sigma2) is irrelevant then, so this is the same as:

sds <- svd(chol(fm3$G))$d
sds^2 / sum(sds^2)

Or you can do the same as you did in those helper functions you wrote.

I find it a bit strange to talk of the contribution of the intercept and slope variances to their sum (since they are on different scales), but maybe this makes sense when using a singular value decomposition of the Cholesky factorization of G instead of just doing summary(princomp(fm$G)) directly.

Best,
Wolfgang

-----Original Message-----
From: Luke Martinez [mailto:martinezlukerm at gmail.com]
Sent: Tuesday, 14 September, 2021 2:12
To: Viechtbauer, Wolfgang (SP)
Cc: R meta
Subject: Re: [R-meta] Clarification on ranef.rma.mv()

Hi Wolfgang,

You'll see the meaningful results when you run the following:

lapply(rePCA_lme(fm2), summary)

The first column in the output (below) is the amount of total
between-subjects variation explained by the random intercepts. The
second column is the amount of total between-subjects variation
cumulatively explained by the random intercepts AND random slopes.

The idea is to find out whether some overfitting of random-effects
structure has taken place or not. For example, in this case, it seems
that random slopes are only responsible for a tiny amount of total
between-subjects variation as 99.41% of such variation is explained by
the random intercepts alone. So, maybe the model's random structure
can be simplified as:

lme(distance ~ age, random = ~1 | Subject , data = Orthodont)

This is a useful feature possibly for metafor as well, so I wonder how
to add some names and "prcomplist" class to do the same thing to
rma.mv just like I extended that from lmer to lme?

I don't think the fact that yi values have their own variance would
interfere with the central concept here.

Thanks,
Luke

$Subject
Importance of components:
                                        [,1]    [,2]
Standard deviation     1.7794 0.13681
Proportion of Variance 0.9941 0.00588
Cumulative Proportion  0.9941 1.00000

On Mon, Sep 13, 2021 at 1:17 PM Viechtbauer, Wolfgang (SP)
<wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:

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 at 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 at 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 at 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 at 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 at 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 at 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 at 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