[R-meta] multilevel and multivariate model

Thank you James, your blog posts are always a great resource :), just for
being sure... what do you mean by "predictors that vary within the sample"?
Using a predictor referring to the lowest nesting level (e.g., outcome in
my case)?

On Tue, 27 Jun 2023 at 15:53, James Pustejovsky <jepusto at gmail.com> wrote:

Yefeng's observations are right on target.

Regarding aggregating over the lowest level, if the model does not
include predictors that vary within the sample, then aggregating is
equivalent to estimating a model with the disaggregated effect sizes but
dropping the lowest-level random effects. More detailed explanation here:
https://www.jepusto.com/sometimes-aggregating-effect-sizes-is-fine/

If the model does include predictors that vary within sample, then
aggregating is not a good idea because you're removing relevant variation
in the outcome.

James

On Tue, Jun 27, 2023 at 5:24?AM Filippo Gambarota via R-sig-meta-analysis
<r-sig-meta-analysis at r-project.org> wrote:

thank you Yefeng. You are right, also robust variance is a good
option. My only concern is about omitting the 4th level as presented
by Wolfgang here
(https://www.metafor-project.org/doku.php/analyses:konstantopoulos2011)
about including or not the nested effect. In fact, my data structure
in the lowest level has only 2 (sometimes 1) effect within each
experiment thus maybe that level can be problematic to estimate. The
only way to reduce model complexity should be to aggregate the lowest
level reducing the model to a 3-level.

On Tue, 27 Jun 2023 at 10:11, Yefeng Yang <yefeng.yang1 at unsw.edu.au>
wrote:

Hi Filippo,

Besides what James and Wolfgang said, you can try different random

structures and choose the so-called "best" (3 levels vs. 4 levels; nested
vs. cross-classified) from the perspective of information criteria like
(c)AIC. But it should be noted that some researchers think the structure
should be informed by your understanding of the data generation mechanism
rather than relying on some sort of measures. The point estimate is usually
not biased. The important thing that needs to take care of is the sampling
variance of the point estimate. In this regard, either a 3- or 4-level
model can work well generally. Robust variance estimation also works well.

Best,
Yefeng

on behalf of Filippo Gambarota via R-sig-meta-analysis <
r-sig-meta-analysis at r-project.org>

Sent: Tuesday, 27 June 2023 17:46
To: R Special Interest Group for Meta-Analysis <

r-sig-meta-analysis at r-project.org>

Cc: Filippo Gambarota <filippo.gambarota at gmail.com>
Subject: Re: [R-meta] multilevel and multivariate model

Thank you James, Yes my parameters seem to be correctly recovered. My
only doubt was about the 4-level model that looks a little bit
complicated but in fact if I have effects nested within experiments
nested within papers seems to be the only solution. Beyond
aggregating, I'm not aware of any other solution for this situation.
What do you think?
Thank you!
Filippo

On Tue, 27 Jun 2023 at 08:52, Viechtbauer, Wolfgang (NP) via
R-sig-meta-analysis <r-sig-meta-analysis at r-project.org> wrote:

Hi James,

No particular reason -- added to my to-do list (probably will make

it an option, like 'sparse=TRUE' as in rma.mv()).

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis [mailto:

r-sig-meta-analysis-bounces at r-project.org] On

Behalf Of James Pustejovsky via R-sig-meta-analysis
Sent: Monday, 26 June, 2023 19:25
To: R Special Interest Group for Meta-Analysis
Cc: James Pustejovsky
Subject: Re: [R-meta] multilevel and multivariate model

Filippo, Your code looks right to me. To double check it for

yourself, try

increasing K to something much larger and verify that rma.mv()

returns

estimates that are very close to the specified parameter values.

(This

might take some time to compute because of the vcalc() call, which

returns

a square matrix with the same number of rows as dat. If nrow(dat)

is large,

then V can get very big.)

Wolfgang, is there a reason that V returns a dense k x k matrix? As

far as

I can see, V should always have block diagonal form since the

cluster

argument is required. Given that, would it be feasible for vcalc()

to

return a Matrix::bdiag (or the metafor equivalent)?

James

On Mon, Jun 26, 2023 at 10:22?AM Filippo Gambarota via

R-sig-meta-analysis <

r-sig-meta-analysis at r-project.org> wrote:

Hi,
I'm working on a meta-analysis with a multilevel (effects within

the

same paper with independent groups) and multivariate (effects

within

the same paper/experiment measured on the same participants i.e.
different outcomes). I'm not sure if my workflow is at least

capturing

the effect size dependency appropriately.

I have some simulated data with the same structure:

```
library(metafor)

seqw <- function(x){
  unlist(sapply(x, function(i) 1:i))
}

set.seed(2023)

K <- 50 # number of papers
J <- sample(1:3, K, replace = TRUE) # number of studies, within

each paper

Z <- sample(1:2, sum(J), replace = TRUE) # number of outcomes per
study/paper

dat <- data.frame(paper = rep(rep(1:K, J), Z),
                  exp = rep(seqw(J), Z),
                  effect = seqw(Z))
head(dat)
```
the `paper` variable is the paper, the `exp` is the experiment
(different experiments have different subjects) and `effect` is

the

outcome within each experiment (1 and/or 2).

Then I simulate a 4-level model:

```
set.seed(2023)
# residual variance components
tau2 <- 0.3
omega2 <- 0.1
zeta2 <- 0.1

# random effects
b0_i <- rnorm(K, 0, sqrt(tau2))
b0_ij <- rnorm(sum(J), 0, sqrt(omega2))
b0_ijz <- rnorm(nrow(dat), 0, sqrt(zeta2))

# add to dataframe
dat$b0_i <- b0_i[dat$paper]
dat$b0_ij <- rep(b0_ij, Z)
dat$b0_ijz <- b0_ijz
dat$vi <- runif(nrow(dat), 0.05, 0.1)
```
Now I create the block variance-covariance matrix where sampling
errors are correlated within each experiment and independent

across

experiments and papers:

```
set.seed(2023)
# create block-matrix
V <- vcalc(vi, cluster = paper, subgroup = exp, obs = effect, rho

=

0.7, data = dat)
# sampling errors
e_ij <- MASS::mvrnorm(1, mu = rep(0, nrow(V)), Sigma = V)
```
Finally I add a dummy variable for the outcome 1 or 2 and

simulate the

observed effects:

```
b0 <- 0.1
b1 <- 0.1
# moderator
dat$x <- ifelse(dat$effect == 1, 1, 0)

# simulate effect
dat$yi <- with(dat, (b0 + b0_i + b0_ij + b0_ijz) + b1*x + e_ij)
dat$x <- factor(dat$x)
dat$exp <- factor(dat$exp)
```
Finally my model should be written as:

```
fit <- rma.mv(yi, V, mods = ~0 + x, random =

~1|paper/exp/effect, data

= dat, sparse = TRUE)
```
My question regards if the simulated data structure is correctly
captured by the proposed model.
Thank you!

Filippo

--
Filippo Gambarota, PhD
Postdoctoral Researcher - University of Padova
Department of Developmental and Social Psychology
Website: filippogambarota.xyz
Research Groups: Colab   Psicostat

--
Filippo Gambarota, PhD
Postdoctoral Researcher - University of Padova
Department of Developmental and Social Psychology
Website: filippogambarota.xyz
Research Groups: Colab   Psicostat

--
Filippo Gambarota, PhD
Postdoctoral Researcher - University of Padova
Department of Developmental and Social Psychology
Website: filippogambarota.xyz
Research Groups: Colab   Psicostat

--
*Filippo Gambarota, PhD*
Postdoctoral Researcher - University of Padova
Department of Developmental and Social Psychology
Website: filippogambarota.xyz
Research Groups: Colab <http://colab.psy.unipd.it/>   Psicostat
<https://psicostat.dpss.psy.unipd.it/>