[R-meta] multilevel and multivariate model

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