[R-meta] Testing for subgroup differences - meta (byvar) vs metafor

Zelniker, Thomas TZELNIKER at BWH.HARVARD.EDU
Wed Dec 6 20:10:55 CET 2017 

Thank you both so much. This was very insightful! I truly appreciate your support.



From: James Pustejovsky [mailto:jepusto at gmail.com]
Sent: Wednesday, December 06, 2017 12:18 PM
To: Viechtbauer Wolfgang (SP) <wolfgang.viechtbauer at maastrichtuniversity.nl>
Cc: Zelniker, Thomas <TZELNIKER at BWH.HARVARD.EDU>; r-sig-meta-analysis at r-project.org
Subject: Re: [R-meta] Testing for subgroup differences - meta (byvar) vs metafor

I'll add one further difference between metagen and rma. In metagen, the test of the moderator is based on the usual large-sample chi-squared approximation, even if the CIs for the sub-group average effects use the Knapp-Hartung adjustment:

metagen_fit <- metagen(MD, seMD, byvar = region, data = Fleiss93cont, tau.common = TRUE, method.tau = "REML", sm = "MD", hakn = TRUE)
metagen_fit$Q.b.random
pchisq(metagen_fit$Q.b.random, df = metagen_fit$df.Q.b, lower.tail = FALSE)

If you fit the same model in metafor, the test of the moderator appears to use the meta-regression version of the Knapp-Hartung adjustment and an F-test with (k - p) degrees of freedom, where k is the number of studies and p is the number of coefficients in the meta-regression:

rma_fit <- rma(y = MD, sei = seMD, mods = ~ region, data = Fleiss93cont, method = "REML", test = "knha")
pf(rma_fit$QM, df1 = 1, df2 = 3, lower.tail = FALSE) # as an F test
2 * pt(abs(rma_fit$zval[2]), df = 3, lower.tail = FALSE) # as a t test

There is a similar discrepancy if we instead use rma.mv<http://rma.mv> and test = "t". In that case (with two variance components), I doubt either version of the test is that great and it would probably be better to use some sort of Satterthwaite approximation for the degrees of freedom. Or get yourself more studies.

James

On Wed, Dec 6, 2017 at 3:23 AM, Viechtbauer Wolfgang (SP) <wolfgang.viechtbauer at maastrichtuniversity.nl<mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>> wrote:
Hi Thomas,

The results are the same, but you just have to go about things in different ways to get the same results. Start with:

metagen(MD, seMD, byvar = region, data = Fleiss93cont, tau.common = TRUE,
        method.tau = "REML", sm = "MD")
rma(y = MD, sei = seMD, mods = ~ region - 1, data = Fleiss93cont,
    method = "REML")

So, assuming a common tau^2 value across subgroups and not using the K&H adjustment. By using 'mods = ~ region - 1' you get the estimate for each level of 'region'. If you leave out the - 1, then one level (in this case 'Asia') becomes the reference level and the coefficient for 'Europe' is the difference between the two levels.

Now let's use the K&H / H&K adjustment:

metagen(MD, seMD, byvar = region, data = Fleiss93cont, tau.common = TRUE,
        method.tau = "REML", sm = "MD", hakn = TRUE)
rma(y = MD, sei = seMD, mods = ~ region - 1, data = Fleiss93cont,
    method = "REML", test = "knha")

You will again get the same estimates, but the CIs are different. metagen() estimates a common tau^2 value and then applies the K&H adjustment separately within each subset. rma() also estimates a common tau^2 value but then applies the meta-regression version of the K&H adjustment. This will give different results. To get the same results, use:

res <- rma(y = MD, sei = seMD, mods = ~ region - 1, data = Fleiss93cont,
           method = "REML", test = "knha")
rma(y = MD, sei = seMD, data = Fleiss93cont, method = "REML", test = "knha",
    subset=region=="Asia", tau2 = res$tau2)
rma(y = MD, sei = seMD, data = Fleiss93cont, method = "REML", test = "knha",
    subset=region=="Europe", tau2 = res$tau2)

So, first estimate the common tau^2 value and then fit separate RE models in each subset with the K&H adjustment and specify the value of the tau^2 from the first step.

Now let's do this without assuming a common tau^2 value:

metagen(MD, seMD, byvar = region, data = Fleiss93cont, tau.common = FALSE,
        method.tau = "REML", sm = "MD")
rma.mv<http://rma.mv>(MD, V= seMD^2, mods = ~ region - 1, method = "REML", random = ~region|study,
       struct="DIAG", data=Fleiss93cont)

As you show, here one can use rma.mv<http://rma.mv>() to accomplish the same thing and the results are the same. And now with K&H adjustment:

metagen(MD, seMD, byvar = region, data = Fleiss93cont, tau.common = FALSE,
        method.tau = "REML", sm = "MD", hakn = TRUE)
rma(y = MD, sei = seMD, data = Fleiss93cont, method = "REML", test = "knha",
    subset=region=="Asia")
rma(y = MD, sei = seMD, data = Fleiss93cont, method = "REML", test = "knha",
    subset=region=="Europe")

So, for rma(), just fit RE models in each subset using the K&H adjustment and then the results are the same. And tau^2 is then of course estimated separately within each subset.

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces at r-project.org<mailto:r-sig-meta-analysis-bounces at r-project.org>] On Behalf Of Zelniker, Thomas
Sent: Wednesday, 06 December, 2017 3:08
To: r-sig-meta-analysis at r-project.org<mailto:r-sig-meta-analysis at r-project.org>
Subject: [R-meta] Testing for subgroup differences - meta (byvar) vs metafor

Dear all,

I am performing a meta-analysis of e.g. 4 trials. Each of these trials reports results for 2 subgroups (group A and group B). My aim is to meta-analyze and test for interaction between these two subgroups.
The R package meta allows doing this very easily by including the term byvar. However, I am not able to replicate this with the package metafor. Could anyone explain to me how the byvar statement works?

What would be the best approach?

Please find below an example indicating the differences between metafor and meta:

library(meta)
library(metafor)
rm(list = ls())
data(Fleiss93cont)
# Add some (fictious) grouping variables:
Fleiss93cont$age <- c(55, 65, 55, 65, 55)
Fleiss93cont$region <- c("Europe", "Europe", "Asia", "Asia", "Europe")
Fleiss93cont$MD <- with(Fleiss93cont, mean.e - mean.c)
Fleiss93cont$seMD <- with(Fleiss93cont, sqrt(sd.e^2/n.e + sd.c^2/n.c))

meta_mod <- metagen(MD, seMD,
                    byvar = region,
                    data = Fleiss93cont,
                    tau.common = F,
                    method.tau = "REML",
                    hakn = T,
                    sm = "MD")

# Here the results:
meta_mod
rma(y = MD, sei = seMD, mods = ~ region, data = Fleiss93cont, method = "REML", test = "knha", intercept = T)
rma.mv<http://rma.mv>(MD, V= seMD^2, mods = ~ region, method = "REML",  test="t", random = ~region|study, struct="DIAG", data=Fleiss93cont)

Best regards,
Thomas

________________________________
Dr. Thomas Zelniker, MD, MSc
Research Fellow
TIMI Study Group
Brigham and Women's Hospital | Harvard Medical School
60 Fenwood Rd. | Suite 7022-7024W | Boston, MA 02115 | Ph: 1-617-278-0326<tel:1-617-278-0326>
tzelniker at bwh.harvard.edu<mailto:tzelniker at bwh.harvard.edu><mailto:tzelniker at partners.org<mailto:tzelniker at partners.org>>  | www.timi.org<http://www.timi.org><http://www.timi.org/>
_______________________________________________
R-sig-meta-analysis mailing list
R-sig-meta-analysis at r-project.org<mailto:R-sig-meta-analysis at r-project.org>
https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis



The information in this e-mail is intended only for the person to whom it is
addressed. If you believe this e-mail was sent to you in error and the e-mail
contains patient information, please contact the Partners Compliance HelpLine at
http://www.partners.org/complianceline . If the e-mail was sent to you in error
but does not contain patient information, please contact the sender and properly
dispose of the e-mail.

	[[alternative HTML version deleted]]

More information about the R-sig-meta-analysis mailing list