#' Subgroup analysis using a mixed-effects model
#'
#' This function performs a mixed-effects (random-effects model within subgroups,
#' fixed-effect model between subgroups) subgroup analysis using \code{meta} objects.
#'
#' @usage subgroup.analysis.mixed.effects(x, subgroups, exclude = "none")
#'
#' @param x An object of class \code{meta}, generated by the \code{metabin}, \code{metagen},
#' \code{metacont}, \code{metacor}, \code{metainc}, or \code{metaprop} function.
#' @param subgroups A character vector of the same length as the number of studies within the
#' meta-analysis, with a unique code for the subgroup each study belongs to. Must have the
#' same order as the studies in the \code{meta} object.
#' @param exclude Single string or concatenated array of strings. The name(s) of the subgroup
#' levels to be excluded from the subgroup analysis. If \code{"none"} (default), all subgroup
#' levels are used for the analysis.
#'
#' @details This function conducts a test for differences in effect sizes between subgroups of a meta-analysis.
#' The function implements a mixed-effect model, in which the overall effect size for each subgroup is
#' calculated using a random-effect model, and the test for subgroup differences is conducted using a
#' fixed-effect model. The implementation follows the fixed-effects (plural) model described in Borenstein
#' and Higgins (2013).
#'
#' This model is appropriate for subgroup tests when the subgroup levels under study
#' are assumed to be exhaustive for the characteristic at hand, and are not randomly chosen instances
#' of a "population" of subgroup levels. For example, the fixed-effects (plural) model used in the function
#' is valid when differences between studies published before and after a certain year are considered as a
#' (binary) subgroup level. When subgroup levels can be assumed to be random samples from a distribution of
#' subgroup levels, a random-effects model is more appropriate, and may be calculated using
#' the \code{\link[meta]{update.meta}} function.
#'
#' The function uses the study effect sizes \code{TE} and their standard error \code{seTE} of the provided
#' \code{meta} object to perform the subgroup analyses. Specifications of the summary measure \code{sm} are
#' inherited and used to backtransform log-transformed effect sizes to their original metrics if necessary.
#'
#' Results can be inspected by plugging the function output into the \code{summary} function. Forest plots
#' can be generated using \code{forest}. Additional arguments of the \code{\link[meta]{forest.meta}} function
#' can be passed to the \code{forest} function for additional styling.
#'
#' @references
#'
#' Harrer, M., Cuijpers, P., Furukawa, T.A, & Ebert, D. D. (2019).
#' \emph{Doing Meta-Analysis in R: A Hands-on Guide}. DOI: 10.5281/zenodo.2551803. \href{https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/subgroup.html}{Chapter 7}.
#'
#' Borenstein, M. & Higgins, J. P. T. (2013). Meta-Analysis and Subgroups. \emph{Prevention Science, 14} (2): 134–43.
#'
#' @author Mathias Harrer & David Daniel Ebert
#'
#' @return Returns a \code{list} with five objects:
#' \itemize{
#' \item \code{within.subgroup.results}: The pooled effect size for each subgroup and corresponding measures of heterogeneity (
#' \code{Q} and \code{I2}). If the summary measure \code{sm} is defined as one of
#' \code{"RR"}, \code{"RD"}, \code{"OR"}, \code{"ASD"}, \code{"HR"} or \code{"IRR"} in the
#' \code{meta} object provided in \code{x}, the backtransformed (exponentiated)
#' pooled effect for each subgroup effect size along with the 95\% confidence interval is also provided.
#' \item \code{subgroup.analysis.results}: The results for the \code{Q}-test for subgroup differences, its degrees of freedom \code{df} and
#' \emph{p}-value.
#' \item \code{m.random}: An object of class \code{meta} containing the results of the random-effects model applied
#' for pooling results in each subgroup in the first step.
#' \item \code{method.tau}: The \eqn{\tau^2} estimator used for within-subgroup pooling
#' (inherited from the \code{meta} object provided in \code{x}).
#' \item \code{k}: The total number of included studies.
#' }
#'
#' @aliases sgame
#'
#' @export subgroup.analysis.mixed.effects
#' @export sgame
#'
#' @import meta
#' @import utils
#'
#' @seealso \code{\link{multimodel.inference}}
#'
#' @examples
#' # Example 1: Hedges' g as effect size, precalculated effect sizes
#' suppressPackageStartupMessages(library(dmetar))
#' suppressPackageStartupMessages(library(meta))
#' data("ThirdWave")
#' ThirdWave = ThirdWave[c(1,2,3,5,9,18),]
#'
#' m1 <- metagen(TE = TE,
#' seTE = seTE,
#' studlab = paste(ThirdWave$Author),
#' data=ThirdWave,
#' comb.fixed = FALSE,
#' method.tau = "PM",
#' sm = "SMD")
#'
#' sgame1 = subgroup.analysis.mixed.effects(x = m1, subgroups = ThirdWave$TypeControlGroup)
#' summary(sgame1)
#'
#' # Example 2: Hedges' g as effect size, raw effect data
#' suppressPackageStartupMessages(library(meta))
#' data(amlodipine)
#'
#' # Create an arbitrary subgroup for illustration purposes
#' amlodipine$subgroup = rep(c("A","B"),4)
#'
#' m2 <- metacont(n.amlo, mean.amlo, sqrt(var.amlo),
#' n.plac, mean.plac, sqrt(var.plac),
#' data=amlodipine, studlab=amlodipine$study,
#' sm = "SMD")
#'
#' sgame2 = subgroup.analysis.mixed.effects(x = m2, subgroups = amlodipine$subgroup)
#' summary(sgame2)
#'
#' # Example 3: Risk ratio as effect size, binary outcome data, exlcude one level
#' suppressPackageStartupMessages(library(meta))
#' data(Olkin95)
#'
#' # Create an arbitrary subgroup for illustration purposes
#' Olkin95$subgroup = c(rep(c("A","B"), 30), rep("C",10))
#'
#' m3 <- metabin(event.e, n.e, event.c, n.c,
#' data = Olkin95, studlab = Olkin95$author,
#' method = "Inverse")
#'
#' # Use shorthand
#' sgame3 = sgame(x = m3, subgroups = Olkin95$subgroup,
#' exclude = "B")
#' summary(sgame3)
#'
#'
#' # Example 4: IRR as effect size, incidence data
#' suppressPackageStartupMessages(library(meta))
#' data(smoking)
#'
#' # Create an arbitrary subgroup for illustration purposes
#' smoking$subgroup = c(rep(c("A"), 4), rep(c("B"), 3))
#'
#' m4 <- metainc(d.smokers, py.smokers,
#' d.nonsmokers, py.nonsmokers,
#' data=smoking, studlab=study, sm="IRR")
#'
#' sgame4 = subgroup.analysis.mixed.effects(x = m4, subgroups = smoking$subgroup)
#' summary(sgame4)
#'
#' \dontrun{
#' # Generate Forest Plot
#' # Additional arguments of the meta::forest.meta can be supplied
#' forest(sgame1, col.diamond = "darkgreen")
#' forest(sgame2)
#' forest(sgame3)
#' forest(sgame4)
#' }
subgroup.analysis.mixed.effects = sgame = function(x, subgroups, exclude = "none") {
update.meta = getFromNamespace("update.meta", "meta")
# Define variables
m = x
subgroups = subgroups
exclude = exclude
# Levels of subgroup
subgroups = as.factor(subgroups)
k = as.vector(summary(subgroups))
levels = levels(subgroups)
k.level.df = data.frame(level = levels, k = k)
# Out Loop for wrong input
if (length(subgroups) != length(m$studlab)) {
stop("Subgroup variable does not contain the same number of cases as the 'meta' object. You need to define a variable which provides a subgroup value for each effect size included in your 'meta' results object.")
}
# get 'Exclude' Subgroup level names
if (exclude[1] != "none") {
levels = levels[!levels %in% exclude]
k = k.level.df[(k.level.df$level %in% levels), ]$k
}
# Create Loop for subgroups
list = list()
for (x in levels) {
list[[x]] = which(subgroups %in% c(paste(x)))
}
# Loop over list to generate subgroup results
sg.results = list()
for (x in 1:length(list)) {
sg.results[[x]] = update.meta(m, subset = list[[x]])
}
# Loop over sg.results to get effect size estimates
ES = vector()
SE = vector()
Qsg = vector()
I2sg = vector()
I2sg.lower = vector()
I2sg.upper = vector()
for (x in 1:length(sg.results)) {
ES[x] = sg.results[[x]]$TE.random
SE[x] = sg.results[[x]]$seTE.random
Qsg[x] = sg.results[[x]]$Q
I2sg[x] = sg.results[[x]]$I2
I2sg.lower[x] = sg.results[[x]]$lower.I2
I2sg.upper[x] = sg.results[[x]]$upper.I2
}
me.data = data.frame(Subgroup = levels, TE = ES, seTE = SE)
# Fixed Meta-Analysis betweens subgroups
meta = metagen(TE, seTE, data = me.data,
comb.fixed = TRUE, comb.random = FALSE,
byvar = Subgroup, hakn = FALSE)
# Create full output dataset
if (m$sm %in% c("RR", "RD", "OR", "ASD", "HR", "IRR")){
subgroup.results = data.frame(Subgroup = me.data$Subgroup,
k = k,
TE = me.data$TE,
seTE = me.data$seTE,
sm = exp(me.data$TE),
LLCI = round(exp(meta$lower), 3),
ULCI = round(exp(meta$upper), 3),
p = meta$pval,
Q = Qsg,
I2 = round(I2sg, 2),
I2.lower = round(I2sg.lower, 2),
I2.upper = round(I2sg.upper, 2))
colnames(subgroup.results)[5] = m$sm
} else {
subgroup.results = data.frame(Subgroup = me.data$Subgroup,
k = k,
TE = me.data$TE,
seTE = me.data$seTE,
LLCI = round(meta$lower, 3),
ULCI = round(meta$upper, 3),
p = meta$pval,
Q = Qsg,
I2 = round(I2sg, 2),
I2.lower = round(I2sg.lower, 2),
I2.upper = round(I2sg.upper, 2))
if (m$sm != ""){
colnames(subgroup.results)[3] = m$sm
colnames(subgroup.results)[4] = "SE"
}
}
mixedeffects.results = data.frame(Q = meta$Q,
df = meta$df.Q,
p = meta$pval.Q,
row.names = "Between groups")
# Create Forest plot data
forest.m = update.meta(m, byvar=subgroups, comb.fixed = FALSE)
if (exclude[1] != "none"){
exclude.level = levels(subgroups)[levels(subgroups) %in% exclude]
exclude.nums = which(subgroups %in% exclude.level)
not.exclude.nums = 1:m$k
not.exclude.nums = not.exclude.nums[!(not.exclude.nums %in% exclude.nums)]
forest.m = update.meta(m, byvar=subgroups, comb.fixed = FALSE, subset = not.exclude.nums)
}
forest.m$TE.random = meta$TE.fixed
forest.m$lower.random = meta$lower.fixed
forest.m$upper.random = meta$upper.fixed
forest.m$Q = meta$Q
forest.m$df.Q = meta$df.Q
forest.m$pval.Q = meta$pval.Q
rownames(subgroup.results) = subgroup.results$Subgroup
subgroup.results$Subgroup = NULL
reslist = list(within.subgroup.results = subgroup.results,
subgroup.analysis.results = mixedeffects.results,
m.random = forest.m,
method.tau = m$method.tau,
k = sum(k))
class(reslist) = "subgroup.analysis.mixed.effects"
invisible(reslist)
reslist
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.