[R-meta] Multilevel meta-analysis with a categorical moderator | subgroup analysis using meta-regression

3 messages · Wolfgang Viechtbauer, Katharina Agethen

Wed, May 14, 2025 7:07 AM #

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:

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 at 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 at gmail.com>
> Sent: Thursday, April 18, 2024 22:27
> To: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer at maastrichtuniversity.nl>
> Cc: Katharina Agethen <katharina.agethen at th-owl.de>; R Special Interest Group
> for Meta-Analysis <r-sig-meta-analysis at 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 at 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 at 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 at gmail.com>
> > Sent: Thursday, April 18, 2024 16:31
> > To: Viechtbauer, Wolfgang (NP)
> <mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>
> > Cc: Katharina Agethen <mailto:katharina.agethen at th-owl.de>; R Special Interest
> Group
> > for Meta-Analysis <mailto:r-sig-meta-analysis at 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 at 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 at th-owl.de>
> > > Sent: Thursday, April 18, 2024 16:22
> > > To: Viechtbauer, Wolfgang (NP)
> > <mailto:mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>;
> > > Reza Norouzian <mailto:mailto:rnorouzian at gmail.com>
> > > Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at 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 at maastrichtuniversity.nl>
> > > Gesendet: Thursday, April 18, 2024 12:09 PM
> > > An: Reza Norouzian <mailto:mailto:rnorouzian at gmail.com>; Katharina Agethen
> > > <mailto:mailto:katharina.agethen at th-owl.de>
> > > Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at 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

Wolfgang Viechtbauer

Wed, May 14, 2025 8:15 AM #

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 at th-owl.de>
Sent: Wednesday, May 14, 2025 16:07
To: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer at maastrichtuniversity.nl>;
Reza Norouzian <rnorouzian at gmail.com>
Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at 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 at 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 at gmail.com>
Sent: Thursday, April 18, 2024 22:27
To: Viechtbauer, Wolfgang (NP)
<mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>
Cc: Katharina Agethen <mailto:katharina.agethen at th-owl.de>; R Special Interest
Group
for Meta-Analysis <mailto:r-sig-meta-analysis at 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 at 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 at 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 at gmail.com>
Sent: Thursday, April 18, 2024 16:31
To: Viechtbauer, Wolfgang (NP)
<mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>
Cc: Katharina Agethen <mailto:katharina.agethen at th-owl.de>; R Special
Interest
Group
for Meta-Analysis <mailto:r-sig-meta-analysis at 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 at 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 at th-owl.de>
Sent: Thursday, April 18, 2024 16:22
To: Viechtbauer, Wolfgang (NP)
<mailto:mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>;
Reza Norouzian <mailto:mailto:rnorouzian at gmail.com>
Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at 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 at maastrichtuniversity.nl>
Gesendet: Thursday, April 18, 2024 12:09 PM
An: Reza Norouzian <mailto:mailto:rnorouzian at gmail.com>; Katharina Agethen
<mailto:mailto:katharina.agethen at th-owl.de>
Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at 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

Katharina Agethen

Wed, May 14, 2025 11:00 AM #

Thanks for the quick response and clarifying, Wolfgang!


Yes, you are right. gamma2.2 and gamma2.3 are estimated to be zero, whereas the estimate for gamma2.1 is 0.0121. I misunderstood the part of the "Interpreting profile likelihoods plots" regarding the peak at zero. So, thank you :)


Best,

Katharina

Von: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer at maastrichtuniversity.nl>

Gesendet: Mittwoch, 14. Mai 2025 17:15:33
An: Katharina Agethen
Cc: R Special Interest Group for Meta-Analysis
Betreff: RE: [R-meta] Multilevel meta-analysis with a categorical moderator | subgroup analysis using meta-regression

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 at th-owl.de>
> Sent: Wednesday, May 14, 2025 16:07
> To: Viechtbauer, Wolfgang (NP) <wolfgang.viechtbauer at maastrichtuniversity.nl>;
> Reza Norouzian <rnorouzian at gmail.com>
> Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at 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 at 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 at gmail.com>
> > Sent: Thursday, April 18, 2024 22:27
> > To: Viechtbauer, Wolfgang (NP)
> <mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>
> > Cc: Katharina Agethen <mailto:katharina.agethen at th-owl.de>; R Special Interest
> Group
> > for Meta-Analysis <mailto:r-sig-meta-analysis at 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 at 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 at 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 at gmail.com>
> > > Sent: Thursday, April 18, 2024 16:31
> > > To: Viechtbauer, Wolfgang (NP)
> > <mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>
> > > Cc: Katharina Agethen <mailto:katharina.agethen at th-owl.de>; R Special
> Interest
> > Group
> > > for Meta-Analysis <mailto:r-sig-meta-analysis at 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 at 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 at th-owl.de>
> > > > Sent: Thursday, April 18, 2024 16:22
> > > > To: Viechtbauer, Wolfgang (NP)
> > > <mailto:mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>;
> > > > Reza Norouzian <mailto:mailto:rnorouzian at gmail.com>
> > > > Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at 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 at maastrichtuniversity.nl>
> > > > Gesendet: Thursday, April 18, 2024 12:09 PM
> > > > An: Reza Norouzian <mailto:mailto:rnorouzian at gmail.com>; Katharina Agethen
> > > > <mailto:mailto:katharina.agethen at th-owl.de>
> > > > Cc: R Special Interest Group for Meta-Analysis <r-sig-meta-analysis at 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