R/subgroup_analyses_mixed_effects_function.R

Defines functions sgame

Documented in sgame

#' 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
#'
#' @importFrom meta forest metagen update.meta
#'
#' @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") {


    # 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

}
MathiasHarrer/dmetar documentation built on March 29, 2020, 7:46 a.m.