[R-meta] Multilevel meta-analysis with a categorical moderator | subgroup analysis using meta-regression
Katharina Agethen
k@th@r|n@@@gethen @end|ng |rom th-ow|@de
Wed May 14 16:07:02 CEST 2025
Dear all,
I am coming back to this thread because I am struggling with the following issue after removing effect sizes from my dataset:
I fitted the following multi-level meta-regression model, allowing heterogeneity to differ between subgroups.
No. of samples/effect sizes per subgroup
sub1: 117/315
sub2: 24/28
sub3: 30/45
mod.pertypetau <- rma.mv(yi = z,
V = vz,
slab = samid,
data = dfa,
random = list(~pertype | samid, ~pertype | esid),
method = "REML",
dfs="contain",
mods = ~ pertype-1,
struct = c("DIAG", "DIAG"))
mod.pertypetau
When profiling the variance components, the peak for gamma2.2 and gamma 2.3 is at zero (also checked the CIs, which include zero). All other components are fine.
Is it advisable to use struct=("DIAG", "ID") instead?
mod.pertypetau.ID <- rma.mv(yi = z,
V = vz,
slab = samid,
data = dfa,
random = list(~pertype | samid, ~pertype | esid),
method = "REML",
dfs="contain",
mods = ~ pertype-1,
struct = c("DIAG", "ID"))
The profile plot for gamma2 seems to be okay, as it does not have its peak at zero (only a small gap on the far right side).
I compared the two models:
> anova(mod.pertypetau, mod.pertypetau.ID)
df AIC BIC AICc logLik LRT pval QE
Full 9 -283.8467 -248.2675 -283.3667 150.9233 3423.3270
Reduced 7 -276.1990 -248.5263 -275.9019 145.0995 11.6477 0.0030 3423.3270
In my opinion, this creates a bit of a dilemma. The statistical test suggests keeping the more complex model, but the more complex model had profiling issues.
Do you have any recommendations how to best handle this issue?
Thanks a lot,
Katharina
________________________________
Von: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>
Gesendet: Montag, 29. April 2024 13:08:54
An: Reza Norouzian
Cc: Katharina Agethen; R Special Interest Group for Meta-Analysis
Betreff: RE: [R-meta] Multilevel meta-analysis with a categorical moderator | subgroup analysis using meta-regression
I am coming back to this thread, because I just made a change in metafor that pertains to the example below. The example made me realize that the way rma.mv() counts levels wasn't quite right. Now, this example:
library(metafor)
dat <- dat.assink2016
# allow the mean effect and the between- and within-study variances to differ across 'deltype'
res <- rma.mv(yi, vi, mods = ~ 0 + deltype,
random = list(~ deltype | study, ~ deltype | id),
struct=c("DIAG","DIAG"), data=dat)
res
yields the following results:
Multivariate Meta-Analysis Model (k = 100; method: REML)
Variance Components:
outer factor: study (nlvls = 17)
inner factor: deltype (nlvls = 3)
estim sqrt k.lvl fixed level
tau^2.1 0.0000 0.0000 1 yes covert
tau^2.2 0.1961 0.4429 17 no general
tau^2.3 0.2140 0.4626 3 no overt
outer factor: id (nlvls = 100)
inner factor: deltype (nlvls = 3)
estim sqrt k.lvl fixed level
gamma^2.1 0.2585 0.5084 9 no covert
gamma^2.2 0.0813 0.2851 78 no general
gamma^2.3 0.0000 0.0000 13 no overt
Test for Residual Heterogeneity:
QE(df = 97) = 761.8270, p-val < .0001
Test of Moderators (coefficients 1:3):
QM(df = 3) = 19.5280, p-val = 0.0002
Model Results:
estimate se zval pval ci.lb ci.ub
deltypecovert -0.1468 0.1883 -0.7796 0.4356 -0.5158 0.2222
deltypegeneral 0.4647 0.1198 3.8795 0.0001 0.2299 0.6995 ***
deltypeovert 0.5428 0.2759 1.9672 0.0492 0.0020 1.0835 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
As can be seen, the fact that there is only one study with deltype 'covert' is correctly recognized and the corresponding variance component is fixed to 0. The results obtained are therefore now identical to what we get from:
# standard multilevel random-effects model within each deltype subset
res1 <- rma.mv(yi, vi, random = ~ 1 | study/id, data=dat, subset=deltype=="covert")
res2 <- rma.mv(yi, vi, random = ~ 1 | study/id, data=dat, subset=deltype=="general")
res3 <- rma.mv(yi, vi, random = ~ 1 | study/id, data=dat, subset=deltype=="overt")
as we can see by comparing the estimates of the fixed effects and variance components:
# compare the estimates of the fixed effects
round(coef(summary(res)), 4)
round(coef(summary(res1)), 4)
round(coef(summary(res2)), 4)
round(coef(summary(res3)), 4)
# compare the estimates of the variances
round(rbind(res$tau2, res$gamma2), 4)
round(sapply(list(res1, res2, res3), \(x) x$sigma2), 4)
Cluster-robust variance estimation for the 'full' model yields:
round(coef(summary(robust(res, cluster=study, clubSandwich=TRUE))), 4)
estimate se tval df pval ci.lb ci.ub
deltypecovert -0.1468 0.0183 -8.0333 1.0000 0.0788 -0.3790 0.0854
deltypegeneral 0.4647 0.1201 3.8710 15.5846 0.0014 0.2097 0.7198
deltypeovert 0.5428 0.2763 1.9646 1.9919 0.1889 -0.6506 1.7361
but the df=1 for 'deltypecovert' should be an indication that we should not trust the SE for this level. For the subset models, we can only get (the same) cluster-robust results for the other two levels:
#round(coef(summary(robust(res1, cluster=study, clubSandwich=TRUE))), 4)
round(coef(summary(robust(res2, cluster=study, clubSandwich=TRUE))), 4)
round(coef(summary(robust(res3, cluster=study, clubSandwich=TRUE))), 4)
Best,
Wolfgang
> -----Original Message-----
> From: Reza Norouzian <rnorouzian using gmail.com>
> Sent: Thursday, April 18, 2024 22:27
> To: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>
> Cc: Katharina Agethen <katharina.agethen using th-owl.de>; R Special Interest Group
> for Meta-Analysis <r-sig-meta-analysis using r-project.org>
> Subject: Re: [R-meta] Multilevel meta-analysis with a categorical moderator |
> subgroup analysis using meta-regression
>
> Forgive the typo in dat2: dat2 <- subset(dat1, !is.na(deltype))
>
> On Thu, Apr 18, 2024 at 2:05 PM Reza Norouzian <mailto:rnorouzian using gmail.com>
> wrote:
> Thanks, Wolfgang. If NAs are to be removed for the single subgroup model
> to work, then, I would probably prefer to use the "res123" models than to remove
> the NAs to get a single subgroup model;-)
>
> dat <- dat.assink2016
> dat1 <- transform(dat, deltype = ifelse(year %in% c(-7.5,-11.5), NA, deltype))
> dat2 <- subset(dat, !is.na(deltype))
> res <- rma.mv(yi, vi, mods = ~ 0 + deltype,
> random = list(~ deltype | study, ~ deltype | id),
> struct="DIAG", data=dat2)
>
> res123 <- lapply(na.omit(unique(dat1$deltype)), \(i) rma.mv(yi, vi,
> random=~1|study/id, data=dat1, subset=deltype==i))
>
> On Thu, Apr 18, 2024 at 11:57 AM Viechtbauer, Wolfgang (NP)
> <mailto:wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> There is a good chance that I might have misremembered this, but I just double-
> checked -- no NAs allowed in the random part.
>
> > -----Original Message-----
> > From: Reza Norouzian <mailto:rnorouzian using gmail.com>
> > Sent: Thursday, April 18, 2024 16:31
> > To: Viechtbauer, Wolfgang (NP)
> <mailto:wolfgang.viechtbauer using maastrichtuniversity.nl>
> > Cc: Katharina Agethen <mailto:katharina.agethen using th-owl.de>; R Special Interest
> Group
> > for Meta-Analysis <mailto:r-sig-meta-analysis using r-project.org>
> > Subject: Re: [R-meta] Multilevel meta-analysis with a categorical moderator |
> > subgroup analysis using meta-regression
> >
> > I would take the software developer's advice on this one;-)
> >
> > On Thu, Apr 18, 2024 at 9:28 AM Viechtbauer, Wolfgang (NP)
> > <mailto:mailto:wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> > Yes, rma.mv() is a bit picky with missing values in variables specified
> > via the 'random' argument.
> >
> > Best,
> > Wolfgang
> >
> > > -----Original Message-----
> > > From: Katharina Agethen <mailto:mailto:katharina.agethen using th-owl.de>
> > > Sent: Thursday, April 18, 2024 16:22
> > > To: Viechtbauer, Wolfgang (NP)
> > <mailto:mailto:wolfgang.viechtbauer using maastrichtuniversity.nl>;
> > > Reza Norouzian <mailto:mailto:rnorouzian using gmail.com>
> > > Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
> project.org>
> > > Subject: AW: [R-meta] Multilevel meta-analysis with a categorical moderator
> |
> > > subgroup analysis using meta-regression
> > >
> > > Thanks for the explanations and especially for pointing out the importance
> of
> > > profiling, Wolfgang and Reza!
> > >
> > > I have one more question:
> > > In my meta-analysis I have a different moderator consisting of two
> categories
> > > (generality = global vs. contextual). For some studies I have missing values
> > > (NAs) for the category.
> > >
> > > I have to delete those studies before specifying the single model, right?
> > >
> > > df.gen <- df[!is.na(df$generality), ]
> > >
> > > mod.generality <- rma.mv(yi = z,
> > > V = vz,
> > > slab = auth,
> > > data = df.gen,
> > > random = list(~generality | samid, ~generality |
> esid),
> > > test = "t",
> > > method = "REML",
> > > dfs="contain",
> > > mods = ~ generality-1,
> > > struct = c("DIAG", "DIAG"))
> > > profile(mod.generality)
> > > mod.generality.robust <- robust(mod.generality, cluster=df.gen$samid,
> clubSandwich = TRUE)
> > > summary(mod.generality.robust )
> > >
> > > Thanks,
> > > Katharina
> > >
> > > -----Ursprüngliche Nachricht-----
> > > Von: Viechtbauer, Wolfgang (NP)
> > <mailto:mailto:wolfgang.viechtbauer using maastrichtuniversity.nl>
> > > Gesendet: Thursday, April 18, 2024 12:09 PM
> > > An: Reza Norouzian <mailto:mailto:rnorouzian using gmail.com>; Katharina Agethen
> > > <mailto:mailto:katharina.agethen using th-owl.de>
> > > Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis using r-
> project.org>
> > > Betreff: RE: [R-meta] Multilevel meta-analysis with a categorical moderator
> |
> > > subgroup analysis using meta-regression
> > >
> > > I wanted to double-check that the equivalence also holds when using cluster-
> > > robust variance estimation, so I ran the following example.
> > >
> > > ###################################
> > >
> > > library(metafor)
> > >
> > > dat <- dat.assink2016
> > >
> > > # standard multilevel random-effects model
> > > res <- rma.mv(yi, vi, random = ~ 1 | study/id, data=dat) res
> > >
> > > # allow the mean effect and the between- and within-study variances to
> differ across 'deltype'
> > > res <- rma.mv(yi, vi, mods = ~ 0 + deltype,
> > > random = list(~ deltype | study, ~ deltype | id),
> > > struct="DIAG", data=dat)
> > > res
> > >
> > > # profile the variance components
> > > profile(res)
> > >
> > > # standard multilevel random-effects model within each deltype subset
> > > res1 <- rma.mv(yi, vi, random = ~ 1 | study/id, data=dat,
> subset=deltype=="covert")
> > > res2 <- rma.mv(yi, vi, random = ~ 1 | study/id, data=dat,
> subset=deltype=="general")
> > > res3 <- rma.mv(yi, vi, random = ~ 1 | study/id, data=dat,
> subset=deltype=="overt")
> > >
> > > # compare the estimates of the fixed effects round(coef(summary(res)), 4)
> > > round(coef(summary(res1)), 4) round(coef(summary(res2)), 4)
> > > round(coef(summary(res3)), 4)
> > >
> > > # compare the estimates of the variances round(rbind(res$tau2, res$gamma2),
> 4)
> > > round(sapply(list(res1, res2, res3), \(x) x$sigma2), 4)
> > >
> > > # cluster-robust variance estimation
> > > round(coef(summary(robust(res, cluster=study, clubSandwich=TRUE))), 4)
> > > #round(coef(summary(robust(res1, cluster=study, clubSandwich=TRUE))), 4) #
> cannot be run
> > > round(coef(summary(robust(res2, cluster=study, clubSandwich=TRUE))), 4)
> > > round(coef(summary(robust(res3, cluster=study, clubSandwich=TRUE))), 4)
> > >
> > > ###################################
> > >
> > > The profiling reveals that tau^2_1 is not identifiable -- the profile is
> > > entirely flat -- because there is a single study for deltype 'covert'. Note
> > > that when fitting model 'res1', rma.mv() detects this issue and sets the
> > > study-level sigma^2 equal to 0. However, this check isn't done when fitting
> > > model 'res'. So while the fixed effect estimate for deltype 'covert' is
> fine,
> > > the corresponding SE doesn't make sense. This goes to show that profiling is
> > > important.
> > >
> > > The example also shows that the rest of the estimates match with the two
> > > approaches. This is also true when using cluster-robust inference methods
> > > (except that this method cannot be used on 'res1' since there is a single
> > > cluster.
> > >
> > > Best,
> > > Wolfgang
[[alternative HTML version deleted]]
More information about the R-sig-meta-analysis
mailing list