[R-meta] multilevel and multivariate model

Tue Jun 27 15:58:32 CEST 2023

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