#' 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 }