Message-ID: <CAFUVuJzBGksoXQU4n4v=x72v-zwQhgwmPTvANB96Qx2gW8+6SA@mail.gmail.com>
Date: 2023-06-27T13:53:26Z
From: James Pustejovsky
Subject: [R-meta] multilevel and multivariate model
In-Reply-To: <CAKWaqjbDVHQX_CgRWVE1UMJziwvuu6yESeKiK0Foehc=jQFcDw@mail.gmail.com>
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 at 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 at 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 at r-project.org>
> on behalf of Filippo Gambarota via R-sig-meta-analysis <
> r-sig-meta-analysis at r-project.org>
> > Sent: Tuesday, 27 June 2023 17:46
> > To: R Special Interest Group for Meta-Analysis <
> r-sig-meta-analysis at r-project.org>
> > Cc: Filippo Gambarota <filippo.gambarota at 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 at 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 at 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 at 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 at 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 at 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 at 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]]