[R-meta] test different levels of a random effect

[Viechtbauer, Wolfgang (NP)]

Dear Yefeng,

(1) Why would you want to add coef(res1) if you then test (coef(res1) + a) - (coef(res1) + b) = a - b, which doesn't involve coef(res1)?

(2) There are some references in the docs:

https://wviechtb.github.io/metafor/reference/ranef.html

(3) Not sure how exactly you want to do this. Maybe this?

res <- rma.mv(yi, V, random = list(~ 1 | study.id, ~ 1 | effect.size.id, ~ 1 | species.id), data=dat, sparse = T)
agg <- aggregate(dat, cluster=species.id, V=res$M)

But the same species show up in different studies, so there are non-zero covariances with estimates belonging to different clusters (species) and aggregate() cannot account for this.

Best,
Wolfgang

> -----Original Message-----
> From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org> On Behalf
> Of Yefeng Yang via R-sig-meta-analysis
> Sent: Thursday, November 21, 2024 05:18
> To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
> project.org>
> Cc: Yefeng Yang <yefeng.yang1 using unsw.edu.au>
> Subject: Re: [R-meta] test different levels of a random effect
>
>  Dear Wolfgang, James and other interested parties,
>
> Thank you very much for your solutions and suggestions. They help a lot. Blups-
> based and likelihood solutions make sense.
>
> It is amazing to see now metafor provides the full var-cov matrices of the BLUPs
> of the random effects. I  still have three additional questions.
>
> (1) If I am keen to test the difference between the "fixed effect plus blup" of
> one species (e.g., coef(res1) + ranef(res1)$species.id$intrcpt[1] ) and "fixed
> effect plus blup" of the other species (e.g., coef(res1) +
> ranef(res1)$species.id$intrcpt[2] ), how to do it properly?
>
> To be more precise, it is not very clear to me how to add coef(res1) and its
> error se(res1) to the following code.
> sub <- ranef(res1, vcov=TRUE)
> sub$pred <- sub$pred$species.id[1:2,]
> sub$vcov <- sub$vcov$species.id[1:2,1:2]
> rma.mv(intrcpt, V=sub$vcov, mods = c(0,1), data=sub$pred)
>
> (2) Would you like to point me the literature providing the details and formulas
> underpinning the estimation of blups and standard error of rma.mv()?
>
> (3) I came across a blog from Wolfgang showing how to incorporate the marginal
> variance-covariance matrix of the effect estimates into the aggregate() to the
> cluster-level mean effect. The so-called within-study averaging.
> https://www.metafor-project.org/doku.php/tips:forest_plot_with_aggregated_values
>
> I am wondering whether we can use this aggregation approaches to get the
> species-level estimates first and then use Wald-type method to test their
> difference.
>
> Very much appreciated the effort and time you put into my question.
>
> Thanks,
> Yefeng
>
> ________________________________
> From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org> on behalf
> of James Pustejovsky via R-sig-meta-analysis <r-sig-meta-analysis using r-project.org>
> Sent: Thursday, November 21, 2024 6:43
> To: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
> project.org>
> Cc: James Pustejovsky <jepusto using gmail.com>
> Subject: Re: [R-meta] test different levels of a random effect
>
> Hi Yefeng,
>
> Yefeng, for your aim (1), whether an estimate associated with a certain
> species is significantly different from 0 (in the sense of P value), I like
> what Wolfgang suggested of using the BLUPs to determine the species for
> which there is evidence of non-zero effect. A challenge is determining
> whether/how to do some sort of familywise error or false discovery rate
> control, since you'd be running 345 hypothesis tests.
>
> For your aim (2), whether any two of such species-specific estimates of
> interest are different, if you are looking for a way to test the *joint*
> null hypothesis that all species-specific estimates are equal (i.e., that
> there are no differences between the effects for any pair of species), then
> it seems like the most direct route would be to test that the species.id
> variance component is equal to zero. You could do this with a likelihood
> ratio test comparing the model with species.id random effects to the model
> without such random effects.
>
> James
>
> James
>
> On Wed, Nov 20, 2024 at 11:59 AM Viechtbauer, Wolfgang (NP) via
> R-sig-meta-analysis <r-sig-meta-analysis using r-project.org> wrote:
>
> > Dear Yefeng,
> >
> > Trying to estimate 345 different variance components in approach 2. seems
> > quite ambitious to me.
> >
> > As for 1., it is worth pointing out that using fixed versus random effects
> > for species is somewhat similar, except that the estimates you get from
> > using random effects exhibit a natural shrinkage effect. This can be nicely
> > visualized as follows:
> >
> > res1 <- rma.mv(yi, V, random = list(~ 1 | study.id, ~ 1 | effect.size.id,
> > ~ 1 | species.id), data=dat, sparse = TRUE)
> > res2 <- rma.mv(yi, V, mods = ~ species.id - 1, random = list(~ 1 |
> > study.id, ~ 1 | effect.size.id), data=dat, sparse = TRUE)
> >
> > plot(coef(res2), coef(res1) + ranef(res1)$species.id$intrcpt, pch=21,
> > bg="gray",
> >      xlab="Fixed Species Effects", ylab="Random Species Effects",
> >      xlim=c(-1,3), ylim=c(-1,3),
> >      panel.first={abline(0,1);
> >                   abline(h=coef(res1), lty="dotted");
> >                   abline(v=coef(res1), lty="dotted")})
> >
> > If you install the latest 'devel' version of metafor, you can use:
> >
> > ranef(res1, vcov=TRUE)
> >
> > to obtain also the full var-cov matrices of the BLUPs of the random
> > effects. Then one could test pairs of species against each other while
> > accounting for the covariance:
> >
> > sub <- ranef(res1, vcov=TRUE)
> > sub$pred <- sub$pred$species.id[1:2,]
> > sub$vcov <- sub$vcov$species.id[1:2,1:2]
> > sub
> >
> > rma.mv(intrcpt, V=sub$vcov, mods = c(0,1), data=sub$pred)
> >
> > Best,
> > Wolfgang