[R-meta] multilevel and multivariate model

Tue Jun 27 16:45:07 CEST 2023

great! thank you now it is very clear!

On Tue, 27 Jun 2023 at 16:39, James Pustejovsky <jepusto using gmail.com> wrote:

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

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