[R-meta] multilevel and multivariate model

[Viechtbauer, Wolfgang (NP)]

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