[R-meta] Comparing estimates of independent subgroups

[Viechtbauer, Wolfgang (SP)]

Hi Roger,

This will test the difference between the slopes using the subgrouping approach:

library(metafor)
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
dat$alloc <- ifelse(dat$alloc == "random", "random", "other")

res1 <- rma(yi, vi, mods = ~ year, data=dat, subset=alloc=="random")
res2 <- rma(yi, vi, mods = ~ year, data=dat, subset=alloc=="other")

dat.comp <- data.frame(estimate = c(coef(res1)[2], coef(res2)[2]), stderror = c(res1$se[2], res2$se[2]),
                       meta = c("random","other"), tau2 = round(c(res1$tau2, res2$tau2),3))
dat.comp
rma(estimate, sei=stderror, mods = ~ meta, method="FE", data=dat.comp, digits=3)
with(dat.comp, round(c(zval = (estimate[1] - estimate[2])/sqrt(stderror[1]^2 + stderror[2]^2)), 3))

This will test the difference between the slopes using a single model, allowing for different amounts of heterogeneity in the two subgroups:

rma.mv(yi, vi, mods = ~ year * alloc, random = ~ alloc | trial, struct="DIAG", data=dat, digits=3)
rma.mv(yi, vi, mods = ~ alloc-1 + year : alloc, random = ~ alloc | trial, struct="DIAG", data=dat, digits=3)

The z-value (i.e., -0.595; and hence the p-value) is identical.

As for your last question:

Use "list(~ inner | outer1, ~ innter | outer1.outer2)" where you have to create "outer1.outer2" yourself. For example:

outer1 outer2 outer1.outer2
---------------------------
1      A      1.A
1      B      1.B
2      A      2.A
3      A      3.A
3      B      3.B
3      C      3.C
4      B      4.B

"list(~ inner | outer1, ~ innter | interaction(outer1, outer2))" might also work, but I think it is safer to create the 'interaction' outside of the model fitting function.

Best,
Wolfgang

>-----Original Message-----
>From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-
>project.org] On Behalf Of Martineau, Roger (AAFC/AAC)
>Sent: Tuesday, 15 January, 2019 14:54
>To: r-sig-meta-analysis using r-project.org
>Subject: [R-meta] Comparing estimates of independent subgroups
>
>Dear list members,
>
>I apologize for the length of this request but this is an issue raised by
>a reviewer for which I have no definite answer.
>
>In http://www.metafor-
>project.org/doku.php/tips:comp_two_independent_estimates, Dr. Viechtbauer
>describes how to compare summary estimates obtained from 3 independent
>subgroups of studies. In the example, the categorical variable “alloc” is
>dichotomized into “random” and “other” and RR from each category are
>compared using 1) a procedure similar to conducting a Wald test, and the
>rma.mv function allowing the residual heterogeneity to differ in random
>and other studies.
>
>In animal science, the data set often include studies in which the levels
>of a categorical variable of interest are not evaluated in the same
>study. For example, a subgroup of studies could have been conducted in
>Ayrshire cows and another subgroup of studies could have used Holstein
>cows. Ayrshire and Holstein cows are 2 different breeds of dairy cows
>which may respond differently to a dietary treatment: the infusion of
>casein may increase milk yield (MY) to a greater extent in Holstein cows
>versus Ayrshire cows. My first question is about the validity of
>comparing the slopes of casein between breed (MY = casein * breed), even
>if no study used the 2 breeds to test the effect of casein; and how to do
>it?.
>
>Can I use the approach in http://www.metafor-
>project.org/doku.php/tips:comp_two_independent_estimates , and compare
>“year” according to “alloc” (i.e., year*alloc)?.
>
>library(metafor)
>dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg,
>data=dat.bcg)
>dat$alloc <- ifelse(dat$alloc == "random", "random", "other")
>dat
>
>### I added moderator year in the model:
>res1 <- rma(yi, vi, mods = ~ year, data=dat, subset=alloc=="random")
>res2 <- rma(yi, vi, mods = ~ year, data=dat, subset=alloc=="other")
>
>### The following script needs to be fixed in order to compare the
>estimates of year:
>dat.comp <- data.frame(estimate = c(coef(res1), coef(res2)), stderror =
>c(res1$se, res2$se),
>                       meta = c("random","other"), tau2 =
>round(c(res1$tau2, res2$tau2),3))
>dat.comp
>rma(estimate, sei=stderror, mods = ~ meta, method="FE", data=dat.comp,
>digits=3)
>with(dat.comp, round(c(zval = (estimate[1] -
>estimate[2])/sqrt(stderror[1]^2 + stderror[2]^2)), 3))
>
>### The final model could be:
>rma.mv(yi, vi, mods = ~ year * alloc, random = ~ alloc | trial,
>struct="DIAG", data=dat, digits=3)
>rma.mv(yi, vi, mods = ~ alloc-1 + year : alloc, random = ~ alloc | trial,
>struct="DIAG", data=dat, digits=3)
>### The slopes of year do not differ (P = 0.552) between random and
>other.
>
>Finally, I cannot use this approach for a multilevel model because we get
>this error message:
>Cannot use '~ inner | outer1/outer2' type terms in the 'random' argument.
>Can this be fixed ?
>
>Best regards,
>
>Roger ☺
>
>Roger Martineau, mv Ph.D.
>Nutrition et Métabolisme des ruminants
>Centre de recherche et de développement
>sur le bovin laitier et le porc
>Agriculture et agroalimentaire Canada/Agriculture and Agri-Food Canada
>Téléphone/Telephone: 819-780-7319
>Télécopieur/Facsimile: 819-564-5507
>2000, Rue Collège / 2000, College Street
>Sherbrooke (Québec)  J1M 0C8
>Canada
>roger.martineau using canada.ca