[R-meta] simulating multivariate random effect meta-analysis

[James Pustejovsky]

Hi Filippo,

This looks reasonable to me. Two suggestions:

* Your approach to simulating the sampling errors (what you call e)
involves a degree of approximation since 1) multivariate d statistics
aren't exactly multivariate normally distributed (although they're quite
close, especially when the primary study sample size is large) and 2) the
formula you used to calculate the variance of d is an approximation to the
actual sampling variance. These approximations seem pretty reasonable to me
given how large your per-group sample sizes are. If you had per-group
sample sizes < 20, I might be more concerned about using approximations.

* Just in terms of computation, you might find it useful to write a little
function that generates the vector of effect size estimates for a single
study, then apply it (using tapply or pmap or similar) rather than
simulating each set of errors for the full set of studies. The sampling
variance matrix V is (k p) X (k p), so it gets quite big in your example,
and most of the entries are zeros. Writing the function and then applying
it across study design parameters would let you avoid the giant matrix,
since the relevant vcov matrices would only ever be as large as p.

James

On Sat, Jan 14, 2023 at 6:41 AM Filippo Gambarota <
filippo.gambarota using gmail.com> wrote:

> Hi,
> I'm trying to simulate a very simple multivariate meta-analysis model.
> I'm following section 5.3 of Cheung (2015, meta-analysis, a structural
> equation modeling approach) and the blog post by James
> (https://www.jepusto.com/simulating-correlated-smds/). I would like to
> have some suggestions about the general approach. In particular, I'm
> not sure if the way I use the d + random effect for calculating the
> sampling variance and then sampling the error term is correct.
> This is my code:
>
> # Packages
>
> library(metafor)
> library(tidyr)
> library(MASS)
>
> # Parameters
>
> set.seed(2023)
>
> k <- 500 # number of studies
> p <- 3 # number of outcomes
> d <- 0.5 # outcomes effect size
> minn <- 50 # minimum sample size
> meann <- 100 # average sample size
> rs <- 0.6 # sampling error correlation
> r <- 0.6 # pre-post correlation (for the vi)
> rho <- 0.7 # the between outcomes correlation
> tau <- 0.4 # the outcomes tau
>
> ```
> # Setup
>
> sim <- tidyr::expand_grid(
>   paper = 1:k,
>   outcome = paste0("y", 1:p),
>   d = d
> )
>
> n <- minn + rpois(k, meann - minn) # generate sample size from a poisson
> sim$n <- rep(n, each = p) # append to the dataset
>
> # between-outcomes variance covariance matrix
> Cmat <- rho + diag(1 - rho, nrow = p)
> Tmat <- diag(tau, nrow = p) %*% Cmat %*% diag(tau, nrow = p)
>
> # study-specific adjustment to the average effect
> sim$u <- c(t(MASS::mvrnorm(k, rep(0, p), Tmat)))
>
> # known sampling variance computed using d + u
> sim$vi <- with(sim, 2*(1-rpp) * (1/n + 1/n) + (d + u)^2 / (2*(n + n)))
>
> # block variance covariance matrix of sampling errors
> V <- metafor::vcalc(vi = vi, cluster = paper, obs = outcome, rho = rs,
> data = sim)
>
> # errors
> sim$e <- MASS::mvrnorm(1, rep(0, nrow(V)), Sigma = V)
>
> # observed effect size yi
> sim$yi <- with(sim, d + u + e)
>
> # fitting the correct model according to the parameters
> fit <- rma.mv(yi, V,
>               mods = ~ 0 + outcome,
>               random = ~outcome|paper,
>               struct = "CS",
>               data = sim,
>               sparse = TRUE)
>
> summary(fit)
> ```
> What do you think? I'm aware that the simulation is very simple.
> Thanks
> --
> Filippo Gambarota
> PhD Student - 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]]