[R-meta] multilevel and multivariate model

Tue Jun 27 16:39:23 CEST 2023

Yes, exactly. In your code, your predictor x varies at the lowest level of
nesting, so I don't think it makes much sense to aggregate over that level.
If your predictors were all about sample- or study-level characteristics,
then aggregating might be more reasonable.

On Tue, Jun 27, 2023 at 8:59 AM Filippo Gambarota <
filippo.gambarota using gmail.com> wrote:

> 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 using 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 using 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 using 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
>>> > ________________________________
>>> > From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org>
>>> on behalf of Filippo Gambarota via R-sig-meta-analysis <
>>> r-sig-meta-analysis using r-project.org>
>>> > Sent: Tuesday, 27 June 2023 17:46
>>> > To: R Special Interest Group for Meta-Analysis <
>>> r-sig-meta-analysis using r-project.org>
>>> > Cc: Filippo Gambarota <filippo.gambarota using 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 using 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 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
>>> > > _______________________________________________
>>> > > 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
>>> >
>>> >
>>> >
>>> > --
>>> > 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
>>>
>>>
>>>
>>> --
>>> 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
>>>
>>
>
> --
> *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/>
>

	[[alternative HTML version deleted]]