[R-meta] multilevel and multivariate model

[James Pustejovsky]

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 using 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
>
> _______________________________________________
> R-sig-meta-analysis mailing list @ R-sig-meta-analysis using r-project.org
> To manage your subscription to this mailing list, go to:
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>

	[[alternative HTML version deleted]]