[R-meta] Multilevel meta-analysis with a categorical moderator | subgroup analysis using meta-regression
Viechtbauer, Wolfgang (NP)
wo||g@ng@v|echtb@uer @end|ng |rom m@@@tr|chtun|ver@|ty@n|
Wed May 14 17:15:33 CEST 2025
Dear Katharina,
If I understand you correctly, gamma2.1 > 0, but gamma2.2 and gamma2.3 are estimated to be zero. That is not an inherent problem. If a variance component is estimated to be zero, then the peak of its profile log-likelihood will also be at zero. See also:
https://wviechtb.github.io/metafor/reference/profile.rma.html#interpreting-profile-likelihood-plots
So, no, I would not switch to struct=("DIAG","ID"), especially since the model with struct=c("DIAG","DIAG")) seems to fit better according to the LRT.
Best,
Wolfgang
> -----Original Message-----
> From: Katharina Agethen <katharina.agethen using th-owl.de>
> Sent: Wednesday, May 14, 2025 16:07
> To: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer using maastrichtuniversity.nl>;
> Reza Norouzian <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
>
> 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)
> <mailto: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 <mailto:rnorouzian using gmail.com>
> > Sent: Thursday, April 18, 2024 22:27
> > 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
> >
> > 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
More information about the R-sig-meta-analysis
mailing list