netmeta: Network Meta-Analysis using Frequentist Methods

Documented in netmetabin

#' Network meta-analysis of binary outcome data
#' 
#' @description
#' Provides three models for the network meta-analysis of binary data
#' (Mantel-Haenszel method, based on the non-central hypergeometric
#' distribution, and the inverse variance method).
#' 
#' @param event1 Number of events (first treatment).
#' @param n1 Number of observations (first treatment).
#' @param event2 Number of events (second treatment).
#' @param n2 Number of observations (second treatment)
#' @param treat1 Label/Number for first treatment.
#' @param treat2 Label/Number for second treatment.
#' @param studlab An optional - but important! - vector with study
#'   labels (see Details).
#' @param data An optional data frame containing the study
#'   information.
#' @param subset An optional vector specifying a subset of studies to
#'   be used.
#' @param sm A character string indicating underlying summary measure,
#'   i.e., \code{"RD"}, \code{"RR"}, \code{"OR"}, \code{"ASD"}.
#' @param method A character string indicating which method is to be
#'   used for pooling of studies. One of \code{"Inverse"},
#'   \code{"MH"}, or \code{"NCH"}, can be abbreviated.
#' @param cc.pooled A logical indicating whether \code{incr} should be
#'   used as a continuity correction, when calculating the network
#'   meta-analysis estimates.
#' @param incr A numerical value which is added to each cell count,
#'   i.e., to the numbers of events and non-events, of all treatment
#'   arms in studies with zero events or non-events in any of the
#'   treatment arms ("continuity correction").
#' @param allincr A logical indicating whether \code{incr} should be
#'   added to each cell count of all studies if a continuity
#'   correction was used for at least one study (only considered if
#'   \code{method = "Inverse"}). If FALSE (default), \code{incr} is
#'   used as continuity correction only for studies with zero events
#'   or zero non-events in any of the treatment arms.
#' @param addincr A logical indicating whether \code{incr} should be
#'   added to each cell count of all studies, irrespective of zero
#'   cell counts (only considered if \code{method = "Inverse"}).
#' @param allstudies A logical indicating whether studies with zero
#'   events or non-events in all treatment arms should be included in
#'   an inverse variance meta-analysis (applies only if \code{method =
#'   "Inverse"} and \code{sm} is equal to either \code{"RR"} or
#'   \code{"OR"}).
#' @param level The level used to calculate confidence intervals for
#'   individual studies.
#' @param level.ma The level used to calculate confidence intervals
#'   for network estimates.
#' @param common A logical indicating whether a common effects network
#'   meta-analysis should be conducted.
#' @param random A logical indicating whether a random effects network
#'   meta-analysis should be conducted.
#' @param prediction A logical indicating whether a prediction
#'   interval should be printed (only considered if \code{method =
#'   "Inverse"}).
#' @param level.predict The level used to calculate prediction
#'   interval for a new study (only considered if \code{method =
#'   "Inverse"}).
#' @param reference.group Reference treatment.
#' @param baseline.reference A logical indicating whether results
#'   should be expressed as comparisons of other treatments versus the
#'   reference treatment (default) or vice versa. This argument is
#'   only considered if \code{reference.group} has been specified.
#' @param all.treatments A logical or \code{"NULL"}. If \code{TRUE},
#'   matrices with all treatment effects, and confidence limits will
#'   be printed.
#' @param seq A character or numerical vector specifying the sequence
#'   of treatments in printouts.
#' @param tau.preset An optional value for manually setting the
#'   square-root of the between-study variance \eqn{\tau^2} (only
#'   considered if \code{method = "Inverse"}).
#' @param tol.multiarm A numeric for the tolerance for consistency of
#'   treatment estimates in multi-arm studies which are consistent by
#'   design (only considered if \code{method = "Inverse"}).
#' @param tol.multiarm.se A numeric for the tolerance for consistency
#'   of standard errors in multi-arm studies which are consistent by
#'   design (only considered if the argument is not \code{NULL} and
#'   \code{method = "Inverse"}).
#' @param details.chkmultiarm A logical indicating whether treatment
#'   estimates and / or variances of multi-arm studies with
#'   inconsistent results or negative multi-arm variances should be
#'   printed (only considered if \code{method = "Inverse"}).
#' @param details.chkdata A logical indicating whether number of
#'   events and participants of studies with inconsistent data should
#'   be printed.
#' @param sep.trts A character used in comparison names as separator
#'   between treatment labels.
#' @param nchar.trts A numeric defining the minimum number of
#'   characters used to create unique treatment names (see Details).
#' @param func.inverse R function used to calculate the pseudoinverse
#'   of the Laplacian matrix L (see \code{\link{netmeta}}).
#' @param backtransf A logical indicating whether results should be
#'   back transformed in printouts and forest plots. If
#'   \code{backtransf = TRUE}, results for \code{sm = "OR"} are
#'   presented as odds ratios rather than log odds ratios, for
#'   example.
#' @param title Title of meta-analysis / systematic review.
#' @param keepdata A logical indicating whether original data (set)
#'   should be kept in netmeta object.
#' @param warn A logical indicating whether warnings should be printed
#'   (e.g., if studies are excluded from meta-analysis due to zero
#'   standard errors).
#' @param warn.deprecated A logical indicating whether warnings should
#'   be printed if deprecated arguments are used.
#' @param \dots Additional arguments (to catch deprecated arguments).
#' 
#' @details
#' This function implements three models for the network meta-analysis
#' of binary data:
#' \itemize{
#' \item The Mantel-Haenszel network meta-analysis model, as described
#'   in Efthimiou et al. (2019) (\code{method = "MH"});
#' \item a network meta-analysis model using the non-central
#'   hypergeometric distribution with the Breslow approximation, as
#'   described in Stijnen et al. (2010) (\code{method = "NCH"});
#' \item the inverse variance method for network meta-analysis
#'   (\code{method = "Inverse"}), also provided by
#'   \code{\link{netmeta}}.
#' }
#' 
#' Comparisons belonging to multi-arm studies are identified by
#' identical study labels (argument \code{studlab}). It is therefore
#' important to use identical study labels for all comparisons
#' belonging to the same multi-arm study.
#' 
#' Data entry for this function is in \emph{contrast-based} format,
#' that is, each line of the data corresponds to a single pairwise
#' comparison between two treatments (arguments \code{treat1},
#' \code{treat2}, \code{event1}, \code{n1}, \code{event2}, and
#' \code{n2}). If data are provided in \emph{arm-based} format, that
#' is, number of events and participants are given for each treatment
#' arm separately, function \code{\link{pairwise}} can be used to
#' transform the data to \emph{contrast-based} format (see help page
#' of function \code{\link{pairwise}}).
#' 
#' Note, all pairwise comparisons must be provided for a multi-arm
#' study. Consider a multi-arm study of \emph{p} treatments with known
#' variances. For this study, the number of events and observations
#' must be provided for each treatment, for each of \emph{p}(\emph{p}
#' - 1) / 2 possible comparisons in separate lines in the data. For
#' instance, a three-arm study contributes three pairwise comparisons,
#' a four-arm study even six pairwise comparisons. Function
#' \code{\link{pairwise}} automatically calculates all pairwise
#' comparisons for multi-arm studies.
#' 
#' For \code{method = "Inverse"}, both common and random effects
#' models are calculated regardless of values choosen for arguments
#' \code{common} and \code{random}. Accordingly, the network estimates
#' for the random effects model can be extracted from component
#' \code{TE.random} of an object of class \code{"netmeta"} even if
#' argument \code{random = FALSE}.  However, all functions in R
#' package \bold{netmeta} will adequately consider the values for
#' \code{common} and \code{random}. E.g. function
#' \code{\link{print.summary.netmeta}} will not print results for the
#' random effects model if \code{random = FALSE}.
#' 
#' For the random-effects model, the direct treatment estimates are
#' based on the common between-study variance \eqn{\tau^2} from the
#' network meta-analysis.
#'
#' For \code{method = "MH"} and \code{method = "NCH"}, only a common
#' effects model is available.
#' 
#' By default, treatment names are not abbreviated in
#' printouts. However, in order to get more concise printouts,
#' argument \code{nchar.trts} can be used to define the minimum number
#' of characters for abbreviated treatment names (see
#' \code{\link{abbreviate}}, argument \code{minlength}). R function
#' \code{\link{treats}} is utilised internally to create abbreviated
#' treatment names.
#' 
#' Names of treatment comparisons are created by concatenating
#' treatment labels of pairwise comparisons using \code{sep.trts} as
#' separator (see \code{\link{paste}}). These comparison names are
#' used in the covariance matrices \code{Cov.common} and
#' \code{Cov.random} and in some R functions, e.g,
#' \code{\link{decomp.design}}. By default, a colon is used as the
#' separator. If any treatment label contains a colon the following
#' characters are used as separator (in consecutive order):
#' \code{"-"}, \code{"_"}, \code{"/"}, \code{"+"}, \code{"."},
#' \code{"|"}, and \code{"*"}. If all of these characters are used in
#' treatment labels, a corresponding error message is printed asking
#' the user to specify a different separator.
#'
#' @return
#' An object of class \code{netmetabin} and \code{netmeta} with
#' corresponding \code{print}, \code{summary}, \code{forest}, and
#' \code{netrank} functions. The object is a list containing the
#' following components:
#' \item{studlab, treat1, treat2}{As defined above.}
#' \item{n1, n2, event1, event2}{As defined above.}
#' \item{TE}{Estimate of treatment effect, i.e. difference between
#'   first and second treatment (e.g. log odds ratio).}
#' \item{seTE}{Standard error of treatment estimate.}
#' \item{k}{Total number of studies.}
#' \item{m}{Total number of pairwise comparisons.}
#' \item{n}{Total number of treatments.}
#' \item{d}{Total number of designs (corresponding to the unique set
#'   of treatments compared within studies).}
#' \item{trts}{Treatments included in network meta-analysis.}
#' \item{k.trts}{Number of studies evaluating a treatment.}
#' \item{n.trts}{Number of observations receiving a treatment.}
#' \item{events.trts}{Number of events observed for a treatment.}
#' \item{studies}{Study labels coerced into a factor with its levels
#'   sorted alphabetically.}
#' \item{narms}{Number of arms for each study.}
#' \item{designs}{Unique list of designs present in the network. A
#'   design corresponds to the set of treatments compared within a
#'   study.}
#' \item{TE.common, seTE.common}{\emph{n}x\emph{n} matrix with estimated
#'   overall treatment effects and standard errors for common effects
#'   model.}
#' \item{lower.common, upper.common}{\emph{n}x\emph{n} matrices with
#'   lower and upper confidence interval limits for common effects
#'   model.}
#' \item{statistic.common, pval.common}{\emph{n}x\emph{n} matrices with
#'   z-value and p-value for test of overall treatment effect under
#'   common effects model.}
#' \item{TE.random, seTE.random}{\emph{n}x\emph{n} matrix with
#'   estimated overall treatment effects and standard errors for
#'   random effects model (only available if \code{method =
#'   "Inverse"}).}
#' \item{lower.random, upper.random}{\emph{n}x\emph{n} matrices with
#'   lower and upper confidence interval limits for random effects
#'   model (only available if \code{method = "Inverse"}).}
#' \item{statistic.random, pval.random}{\emph{n}x\emph{n} matrices
#'   with z-value and p-value for test of overall treatment effect
#'   under random effects model (only available if \code{method =
#'   "Inverse"}).}
#' \item{TE.direct.common, seTE.direct.common}{\emph{n}x\emph{n} matrix
#'   with estimated treatment effects and standard errors from direct
#'   evidence under common effects model.}
#' \item{lower.direct.common, upper.direct.common}{\emph{n}x\emph{n}
#'   matrices with lower and upper confidence interval limits from
#'   direct evidence under common effects model.}
#' \item{statistic.direct.common, pval.direct.common}{\emph{n}x\emph{n}
#'   matrices with z-value and p-value for test of overall treatment
#'   effect from direct evidence under common effects model.}
#' \item{TE.direct.random, seTE.direct.random}{\emph{n}x\emph{n}
#'   matrix with estimated treatment effects and standard errors from
#'   direct evidence under random effects model (only available if
#'   \code{method = "Inverse"}).}
#' \item{lower.direct.random, upper.direct.random}{\emph{n}x\emph{n}
#'   matrices with lower and upper confidence interval limits from
#'   direct evidence under random effects model (only available if
#'   \code{method = "Inverse"}).}
#' \item{statistic.direct.random,
#'   pval.direct.random}{\emph{n}x\emph{n} matrices with z-value and
#'   p-value for test of overall treatment effect from direct evidence
#'   under random effects model (only available if \code{method =
#'   "Inverse"}).}
#' \item{Q}{Overall heterogeneity / inconsistency statistic. (only
#'   available if \code{method = "Inverse"})}
#' \item{df.Q}{Degrees of freedom for test of heterogeneity /
#'   inconsistency.}
#' \item{pval.Q}{P-value for test of heterogeneity / inconsistency.}
#' \item{I2, lower.I2, upper.I2}{I-squared, lower and upper confidence
#'   limits (only available if \code{method = "Inverse"}).}
#' \item{tau}{Square-root of between-study variance (only available if
#'   \code{method = "Inverse"}).}
#' \item{Q.heterogeneity}{Overall heterogeneity statistic. (only
#'   available if \code{method = "Inverse"})}
#' \item{df.Q.heterogeneity}{Degrees of freedom for test of overall
#'   heterogeneity.}
#' \item{pval.Q.heterogeneity}{P-value for test of overall
#'   heterogeneity.}
#' \item{Q.inconsistency}{Overall inconsistency statistic.}
#' \item{df.Q.inconsistency}{Degrees of freedom for test of overall
#'   inconsistency.}
#' \item{pval.Q.inconsistency}{P-value for test of overall
#'   inconsistency.}
#' \item{A.matrix}{Adjacency matrix (\emph{n}x\emph{n}).}
#' \item{H.matrix}{Hat matrix (\emph{m}x\emph{m})}
#' \item{n.matrix}{\emph{n}x\emph{n} matrix with number of
#'   observations in direct comparisons.}
#' \item{events.matrix}{\emph{n}x\emph{n} matrix with number of events
#'   in direct comparisons.}
#' \item{sm, method, level, level.ma}{As defined above.}
#' \item{incr, allincr, addincr, allstudies, cc.pooled}{As defined
#'   above.}
#' \item{common, random}{As defined above.}
#' \item{prediction, level.predict}{As defined above.}
#' \item{reference.group, baseline.reference, all.treatments}{As
#'   defined above.}
#' \item{seq, tau.preset, tol.multiarm, tol.multiarm.se}{As defined
#'   above.}
#' \item{details.chkmultiarm, details.chkdata}{As defined above.}
#' \item{sep.trts, nchar.trts}{As defined above.}
#' \item{backtransf, title, warn, warn.deprecated}{As defined above.}
#' \item{data}{Data set (in contrast-based format).}
#' \item{data.design}{List with data in arm-based format (each list
#'   element corresponds to a single design).}
#' \item{call}{Function call.}
#' \item{version}{Version of R package netmeta used to create object.}
#' 
#' @author Orestis Efthimiou \email{oremiou@@gmail.com}, Guido
#'   Schwarzer \email{guido.schwarzer@@uniklinik-freiburg.de}
#' 
#' @seealso \code{\link{pairwise}}, \code{\link{netmeta}}
#' 
#' @references
#' Efthimiou O, Rücker G, Schwarzer G, Higgins J, Egger M, Salanti G
#' (2019):
#' A Mantel-Haenszel model for network meta-analysis of rare events.
#' \emph{Statistics in Medicine},
#' \bold{38}, 2992--3012
#' 
#' Senn S, Gavini F, Magrez D, Scheen A (2013):
#' Issues in performing a network meta-analysis.
#' \emph{Statistical Methods in Medical Research},
#' \bold{22}, 169--89
#' 
#' Stijnen T, Hamza TH, Ozdemir P (2010):
#' Random effects meta-analysis of event outcome in the framework of
#' the generalized linear mixed model with applications in sparse
#' data.
#' \emph{Statistics in Medicine},
#' \bold{29}, 3046--67

#' @examples
#' data(Dong2013)
#' 
#' # Only consider first ten studies (to reduce runtime of example)
#' #
#' first10 <- subset(Dong2013, id <= 10)
#' 
#' # Transform data from long arm-based format to contrast-based
#' # format. Argument 'sm' has to be used for odds ratio as summary
#' # measure; by default the risk ratio is used in the metabin
#' # function called internally.
#' #
#' p1 <- pairwise(treatment, death, randomized, studlab = id,
#'   data = first10, sm = "OR")
#' 
#' # Conduct Mantel-Haenszel network meta-analysis (without continuity
#' # correction)
#' #
#' nb1 <- netmetabin(p1, ref = "plac")
#' nb1
#' 
#' # Obtain the league table
#' #
#' netleague(nb1)
#' 
#' \dontrun{
#' # Conduct Mantel-Haenszel network meta-analysis for the whole
#' # dataset
#' #
#' p2 <- pairwise(treatment, death, randomized, studlab = id,
#'   data = Dong2013, sm = "OR")
#' netmetabin(p2, ref = "plac")
#'   
#' # Conduct network meta-analysis using the non-central
#' # hypergeometric model (without continuity correction)
#' #
#' netmetabin(p2, ref = "plac", method = "NCH")
#' 
#' # Conduct Mantel-Haenszel network meta-analysis (with continuity
#' # correction of 0.5; include all studies)
#' #
#' netmetabin(p2, ref = "plac", cc.pooled = TRUE)
#' 
#' data(Gurusamy2011)
#' 
#' p3 <- pairwise(treatment, death, n, studlab = study,
#'   data = Gurusamy2011, sm = "OR")
#' 
#' # Conduct Mantel-Haenszel network meta-analysis (without continuity
#' # correction)
#' #
#' netmetabin(p3, ref = "cont")
#' }
#' 
#' @export netmetabin


netmetabin <- function(event1, n1, event2, n2,
                       treat1, treat2, studlab,
                       data = NULL, subset = NULL,
                       sm,
                       method = "MH",
                       cc.pooled = FALSE,
                       incr, allincr, addincr, allstudies,
                       level = gs("level"),
                       level.ma = gs("level.ma"),
                       common = gs("common"),
                       random = method == "Inverse" &
                         (gs("random") | !is.null(tau.preset)),
                       ##
                       prediction = FALSE,
                       level.predict = gs("level.predict"),
                       ##
                       reference.group = "",
                       baseline.reference = TRUE,
                       all.treatments = NULL,
                       seq = NULL,
                       ##
                       tau.preset = NULL,
                       ##
                       tol.multiarm = 0.001,
                       tol.multiarm.se = NULL,
                       details.chkmultiarm = FALSE,
                       details.chkdata = TRUE,
                       ##
                       sep.trts = ":",
                       nchar.trts = 666,
                       ##
                       func.inverse = invmat,
                       ##
                       backtransf = gs("backtransf"),
                       ##
                       title = "",
                       keepdata = gs("keepdata"),
                       warn = TRUE, warn.deprecated = gs("warn.deprecated"),
                       ...) {
  
  
  ##
  ##
  ## (1) Check arguments
  ##
  ##
  modtext <-
    paste0("must be equal to 'Inverse' (classic network meta-analysis), ",
           "'MH' (Mantel-Haenszel, the default) or ",
           "'NCH' (common-effects non-central hypergeometric).")
  method <- setchar(method, c("Inverse", "MH", "NCH"), modtext)
  chklogical(cc.pooled)
  ##
  chklevel(level)
  chklevel(level.predict)
  ##
  chklogical(prediction)
  ##
  ## Check for deprecated arguments in '...'
  ##
  args  <- list(...)
  chklogical(warn.deprecated)
  ##
  level.ma <- deprecated(level.ma, missing(level.ma), args, "level.comb",
                         warn.deprecated)
  chklevel(level.ma)
  ##
  missing.common <- missing(common)
  common <- deprecated(common, missing.common, args, "comb.fixed",
                       warn.deprecated)
  common <- deprecated(common, missing.common, args, "fixed",
                       warn.deprecated)
  chklogical(common)
  ##
  random <-
    deprecated(random, missing(random), args, "comb.random", warn.deprecated)
  chklogical(random)
  ##
  if (method != "Inverse" & !common) {
    warning("Argument 'common' set to TRUE for Mantel-Haenszel ",
            "method and non-central hypergeometric distribution.",
            call. = FALSE)
    common <- TRUE
  }
  ##
  if (method != "Inverse" & random) {
    warning("Argument 'random' set to FALSE for Mantel-Haenszel ",
            "method and non-central hypergeometric distribution.",
            call. = FALSE)
    random <- FALSE
  }
  ##
  if (method != "Inverse" & prediction) {
    warning("Argument 'prediction' set to FALSE for Mantel-Haenszel ",
            "method and non-central hypergeometric distribution.",
            call. = FALSE)
    prediction <- FALSE
  }
  ##
  chklogical(baseline.reference)
  ##
  if (!is.null(all.treatments))
    chklogical(all.treatments)
  ##
  if (!is.null(tau.preset))
    chknumeric(tau.preset, min = 0, length = 1)
  ##
  chknumeric(tol.multiarm, min = 0, length = 1)
  if (!is.null(tol.multiarm.se))
    chknumeric(tol.multiarm.se, min = 0, length = 1)
  chklogical(details.chkmultiarm)
  chklogical(details.chkdata)
  ##
  missing.sep.trts <- missing(sep.trts)
  chkchar(sep.trts)
  chknumeric(nchar.trts, min = 1, length = 1)
  ##
  chklogical(backtransf)
  ##
  chkchar(title)
  chklogical(keepdata)
  chklogical(warn)
  ##
  ## Check value for reference group
  ##
  if (is.null(all.treatments))
    if (reference.group == "")
      all.treatments <- TRUE
    else
      all.treatments <- FALSE
  ##
  chklogical(baseline.reference)
  ##
  missing.incr <- missing(incr)
  
  
  ##
  ##
  ## (2) Read data
  ##
  ##
  nulldata <- is.null(data)
  sfsp <- sys.frame(sys.parent())
  mc <- match.call()
  ##
  if (nulldata)
    data <- sfsp
  ##
  ## Catch 'event1', 'event2', 'n1', 'n2', 'treat1', 'treat2', and
  ## 'studlab' from data:
  ##
  event1 <- catch("event1", mc, data, sfsp)
  ##
  if (is.data.frame(event1)  &&
      (!is.null(attr(event1, "pairwise")) ||
       inherits(event1, "pairwise"))) {
    is.pairwise <- TRUE
    ##
    if (missing.incr)
      incr <-
        replaceNULL(attr(event1, "incr"), gs("incr"))
    if (missing(allincr))
      allincr <-
        replaceNULL(attr(event1, "allincr"), gs("allincr"))
    if (missing(addincr))
      addincr <-
        replaceNULL(attr(event1, "addincr"), gs("addincr"))
    if (missing(allstudies))
      allstudies <-
        replaceNULL(attr(event1, "allstudies"), gs("allstudies"))
    ##
    if (missing(sm) & method == "Inverse")
      sm <- replaceNULL(attr(event1, "sm"), gs("smbin"))
    else if (method != "Inverse") {
      if (!missing(sm) && tolower(sm) != "or")
        warning("Argument 'sm' set to 'OR'.", call. = FALSE)
      sm <- "OR"
    }
    ##
    n1 <- event1$n1
    event2 <- event1$event2
    n2 <- event1$n2
    treat1 <- event1$treat1
    treat2 <- event1$treat2
    studlab <- event1$studlab
    ##
    pairdata <- event1
    data <- event1
    ##
    event1 <- event1$event1
    ##
    chknull(event1, text = "Variable")
    chknull(n1, text = "Variable")
    chknull(event2, text = "Variable")
    chknull(n2, text = "Variable")
    chknull(treat1, text = "Variable")
    chknull(treat2, text = "Variable")
    chknull(studlab, text = "Variable")
  }
  else {
    is.pairwise <- FALSE
    ##
    if (missing.incr)
      incr <- gs("incr")
    if (missing(allincr))
      allincr <- gs("allincr")
    if (missing(addincr))
      addincr <- gs("addincr")
    if (missing(allstudies))
      allstudies <- gs("allstudies")
    ##
    if (missing(sm) & method == "Inverse") {
      if (!nulldata && !is.null(attr(data, "sm")))
        sm <- attr(data, "sm")
      else
        sm <- ""
    }
    else if (method != "Inverse") {
      if (!missing(sm) && tolower(sm) != "or")
        warning("Argument 'sm' set to 'OR'.", call. = FALSE)
      sm <- "OR"
    }
    ##
    event2 <- catch("event2", mc, data, sfsp)
    ##
    n1 <- catch("n1", mc, data, sfsp)
    n2 <- catch("n2", mc, data, sfsp)
    ##
    treat1 <- catch("treat1", mc, data, sfsp)
    treat2 <- catch("treat2", mc, data, sfsp)
    ##
    studlab <- catch("studlab", mc, data, sfsp)
    ##
    chknull(event1)
    chknull(n1)
    chknull(event2)
    chknull(n2)
    chknull(treat1)
    chknull(treat2)
    chknull(studlab)
  }
  ##
  chknumeric(event1)
  chknumeric(n1)
  chknumeric(event2)
  chknumeric(n2)
  ##
  k.Comp <- length(event1)
  ##
  if (is.factor(treat1))
    treat1 <- as.character(treat1)
  if (is.factor(treat2))
    treat2 <- as.character(treat2)
  if (is.factor(studlab))
    studlab <- as.character(studlab)
  ##
  ## Remove leading and trailing whitespace
  ##
  treat1 <- rmSpace(rmSpace(treat1, end = TRUE))
  treat2 <- rmSpace(rmSpace(treat2, end = TRUE))
  ##
  ## Keep original order of studies
  ##
  .order <- seq_along(studlab)
  ##
  subset <- catch("subset", mc, data, sfsp)
  missing.subset <- is.null(subset)
  
  
  ##
  ##
  ## (2b) Store complete dataset in list object data
  ##
  ##
  if (nulldata & !is.pairwise)
    data <- data.frame(.event1 = event1)
  else if (nulldata & is.pairwise) {
    data <- pairdata
    data$.order <- .order
    data$.event1 <- event1
  }
  else
    data$.event1 <- event1
  ##
  data$.n1 <- n1
  data$.event2 <- event2
  data$.n2 <- n2
  data$.treat1 <- treat1
  data$.treat2 <- treat2
  data$.studlab <- studlab
  data$.order <- .order
  ##
  if (!missing.subset) {
    if (length(subset) == dim(data)[1])
      data$.subset <- subset
    else {
      data$.subset <- FALSE
      data$.subset[subset] <- TRUE
    }
  }
  ##
  chklogical(incr)
  chklogical(allincr)
  chklogical(addincr)
  chklogical(allstudies)
  ##
  m.data <- metabin(event1, n1, event2, n2,
                    sm = sm, method = "Inverse",
                    incr = incr, allincr = allincr,
                    addincr = addincr, allstudies = allstudies,
                    method.tau = "DL", method.tau.ci = "",
                    warn = FALSE,
                    warn.deprecated = FALSE)
  ##
  data$.TE <- m.data$TE
  data$.seTE <- m.data$seTE
  
  
  ##
  ##
  ## (3) Use subset for analysis
  ##
  ##
  if (!missing.subset) {
    if ((is.logical(subset) & (sum(subset) > k.Comp)) ||
        (length(subset) > k.Comp))
      stop("Length of subset is larger than number of studies.",
           call. = FALSE)
    ##
    studlab <- studlab[subset]
    treat1 <- treat1[subset]
    treat2 <- treat2[subset]
    ##
    event1 <- event1[subset]
    event2 <- event2[subset]
    n1 <- n1[subset]
    n2 <- n2[subset]
    ##
    .order <- .order[subset]
  }
  ##
  ## Check for correct treatment order within comparison
  ##
  wo <- treat1 > treat2
  ##
  if (any(wo)) {
    tevent1 <- event1
    event1[wo] <- event2[wo]
    event2[wo] <- tevent1[wo]
    ##
    tn1 <- n1
    n1[wo] <- n2[wo]
    n2[wo] <- tn1[wo]
    ##
    ttreat1 <- treat1
    treat1[wo] <- treat2[wo]
    treat2[wo] <- ttreat1[wo]
  }
  ##
  labels <- sort(unique(c(treat1, treat2)))
  ##
  if (compmatch(labels, sep.trts)) {
    if (!missing.sep.trts)
      warning("Separator '", sep.trts, "' used in at least ",
              "one treatment label. Try to use predefined separators: ",
              "':', '-', '_', '/', '+', '.', '|', '*'.",
              call. = FALSE)
    ##
    if (!compmatch(labels, ":"))
      sep.trts <- ":"
    else if (!compmatch(labels, "-"))
      sep.trts <- "-"
    else if (!compmatch(labels, "_"))
      sep.trts <- "_"
    else if (!compmatch(labels, "/"))
      sep.trts <- "/"
    else if (!compmatch(labels, "+"))
      sep.trts <- "+"
    else if (!compmatch(labels, "."))
      sep.trts <- "-"
    else if (!compmatch(labels, "|"))
      sep.trts <- "|"
    else if (!compmatch(labels, "*"))
      sep.trts <- "*"
    else
      stop("All predefined separators (':', '-', '_', '/', '+', '.', '|', '*') ",
           "are used in at least one treatment label.\n   ",
           "Please specify a different character that should be used ",
           " as separator (argument 'sep.trts').",
           call. = FALSE)
  }
  ##
  if (reference.group != "")
    reference.group <- setref(reference.group, labels)
  ##
  rm(labels)


  ##
  ##
  ## (4) Additional checks
  ##
  ##
  if (any(treat1 == treat2))
    stop("Treatments must be different (arguments 'treat1' and 'treat2').",
         call. = FALSE)
  ##
  if (length(studlab) != 0)
    studlab <- as.character(studlab)
  else {
    if (warn)
      warning("No information given for argument 'studlab'. ",
              "Assuming that comparisons are from independent studies.",
              call. = FALSE)
    studlab <- seq(along = event1)
  }
  ##
  ## Check for correct number of comparisons
  ##
  tabnarms <- table(studlab)
  sel.narms <- !is.wholenumber((1 + sqrt(8 * tabnarms + 1)) / 2)
  ##
  if (sum(sel.narms) == 1)
    stop(paste0("Study '", names(tabnarms)[sel.narms],
                "' has a wrong number of comparisons.",
                "\n  Please provide data for all treatment comparisons",
                " (two-arm: 1; three-arm: 3; four-arm: 6, ...)."),
         call. = FALSE)
  if (sum(sel.narms) > 1)
    stop(paste0("The following studies have a wrong number of comparisons: ",
                paste(paste0("'", names(tabnarms)[sel.narms], "'"),
                      collapse = " - "),
                "\n  Please provide data for all treatment comparisons",
                " (two-arm: 1; three-arm: 3; four-arm: 6, ...)."),
         call. = FALSE)
  ##
  ## Check number of subgraphs
  ##
  n.subnets <- netconnection(treat1, treat2, studlab)$n.subnets
  ##
  if (n.subnets > 1)
    stop(paste0("Network consists of ", n.subnets, " separate sub-networks.\n",
                "  Use R function 'netconnection' to identify sub-networks."),
         call. = FALSE)
  ##
  ## Catch 'incr' from data:
  ##
  if (!missing.incr)
    incr <- catch("incr", mc, data, sfsp)
  chknumeric(incr, min = 0)
  ##
  if (method != "Inverse" & length(incr) > 1) {
    warning("Argument 'incr' must be a single value for ",
            "Mantel-Haenszel and common-effects non-central ",
            "hypergeometric method. Set to zero (default).",
            call. = FALSE)
    incr <- 0
  }
  ##
  if (all(incr == 0) & method != "Inverse" & cc.pooled == TRUE)
    cc.pooled <- FALSE
  ##
  lengthunique <- function(x) length(unique(x))


  ##
  ##
  ## (5) Create analysis dataset
  ##
  ##
  dat.long <- rbind(data.frame(studlab, treat = treat1, event = event1, n = n1,
                               .order,
                               stringsAsFactors = FALSE),
                    data.frame(studlab, treat = treat2, event = event2, n = n2,
                               .order,
                               stringsAsFactors = FALSE))
  ##
  dat.wide <- data.frame(studlab = studlab, treat1 = treat1, treat2 = treat2,
                         event1 = event1, n1 = n1, event2 = event2, n2 = n2,
                         .order,
                         stringsAsFactors = FALSE)
  dat.wide <- dat.wide[order(dat.wide$studlab,
                             dat.wide$treat1, dat.wide$treat2), ]
  ##
  ## Check for correct number of treatments
  ##
  d.long <- unique(dat.long[, c("studlab", "treat", "event", "n")])
  d.long <- d.long[order(d.long$studlab, d.long$treat), , drop = FALSE]
  ##
  tabnarms <- table(d.long$studlab, d.long$treat)
  larger.one <- function(x) any(x > 1)
  sel.study <- apply(tabnarms, 1, larger.one)
  ##
  if (sum(sel.study) == 1) {
    if (details.chkdata) {
      cat("Inconsistent number of events and participants:\n")
      d.l <- d.long[d.long$studlab == rownames(tabnarms)[sel.study], ]
      prmatrix(d.l,
               quote = FALSE, right = TRUE,
               rowlab = rep("", nrow(d.l)))
    }
    stop(paste0("Study '", rownames(tabnarms)[sel.study],
                "' has inconsistent data. Please check the original data."),
         call. = FALSE)
  }
  if (sum(sel.study) > 1) {
    if (details.chkdata) {
      cat("Inconsistent number of events and participants:\n")
      d.l <- d.long[d.long$studlab %in% rownames(tabnarms)[sel.study], ]
      prmatrix(d.l,
               quote = FALSE, right = TRUE,
               rowlab = rep("", nrow(d.l)))
    }
    ##
    stop(paste0("The following studies have inconsistent data: ",
                paste(paste0("'", rownames(tabnarms)[sel.study], "'"),
                      collapse = " - "),
                "\n  Please check the original data."),
         call. = FALSE)
  }
  ##
  get.designs <- function(x) {
    net1 <-
      netmeta(pairwise(studlab = x$studlab, treat = x$treat,
                       event = x$event, n = x$n,
                       sm = "RD"),
              warn = FALSE)
    ##
    if (net1$n == 2)
      res <- data.frame(studlab = net1$studlab, design = net1$designs)
    else
      res <- nma.krahn(net1)$studies[, c("studlab", "design")]
    ##
    res$design <- as.character(res$design)
    res <- unique(res)
    ##
    res
  }
  ##
  data$.drop <- rep(FALSE, nrow(data))
  ##
  ## Add variable 'non.event'
  ##
  dat.long$non.event <- dat.long$n - dat.long$event
  ##
  dat.wide$non.event1 <- dat.wide$n1 - dat.wide$event1
  dat.wide$non.event2 <- dat.wide$n2 - dat.wide$event2
  ##
  ## Remove duplicated rows from dataset (due to multi-arm studies)
  ##
  dupl <- duplicated(dat.long[, c("studlab", "treat", "event", "n")])
  ##
  if (any(dupl))
    dat.long <- dat.long[!dupl, ]
  ##
  rm(dupl)
  ##
  all.designs <- get.designs(dat.long)
  ##
  ## Add variable 'design' with study design
  ##
  tdat.design <- get.designs(dat.long)
  ##
  dat.long <- merge(dat.long, tdat.design, by = "studlab",
                    all.x = TRUE)
  dat.wide <- merge(dat.wide, tdat.design, by = "studlab",
                    all.x = TRUE)
  ##
  dat.long <- dat.long[order(dat.long$.order), ]
  dat.wide <- dat.wide[order(dat.wide$.order), ]
  ##
  names(tdat.design) <- c("studlab", ".design")
  ##
  data <- merge(data, tdat.design,
                by.x = ".studlab", by.y = "studlab",
                all.x = TRUE)
  data <- data[order(data$.order), ]
  ##
  rm(tdat.design)
  
  
  ##
  ##
  ## (6) Stage I: setting up the data (Efthimiou et al., 2019)
  ##
  ##
  ## Step i. Remove all-zero or all-event studies (only MH and NCH methods)
  ##
  if (method != "Inverse") {
    events.study <- with(dat.long, tapply(event, studlab, sum))
    nonevents.study <- with(dat.long, tapply(n - event, studlab, sum))
    ##
    if (any(events.study == 0) | any(nonevents.study == 0)) {
      zeroevents <- events.study == 0
      allevents <- nonevents.study == 0 & events.study != 0
      keep <- !(zeroevents | allevents)
      ##
      if (warn) {
        if (sum(zeroevents) == 1)
          warning("Study '", names(events.study)[zeroevents],
                  "' without any events excluded from network meta-analysis.",
                  call. = FALSE)
        else if (sum(zeroevents) > 1)
          warning("Studies without any events excluded ",
                  "from network meta-analysis: ",
                  paste(paste0("'", names(events.study)[zeroevents], "'"),
                        collapse = " - "),
                  call. = FALSE)
        ##
        if (sum(allevents) == 1)
          warning("Study '", names(nonevents.study)[allevents],
                  "' with all events excluded from network meta-analysis.",
                  call. = FALSE)
        else if (sum(allevents) > 1)
          warning("Studies with all events excluded ",
                  "from network meta-analysis: ",
                  paste(paste0("'", names(nonevents.study)[allevents], "'"),
                        collapse = " - "),
                  call. = FALSE)
      }
      ##
      dat.long <- dat.long[dat.long$studlab %in% names(events.study)[keep], ,
                           drop = FALSE]
      dat.wide <- dat.wide[dat.wide$studlab %in% names(events.study)[keep], ,
                           drop = FALSE]
      ##
      data$.drop <- data$.drop | data$studlab %in% names(events.study)[!keep]
      ##
      rm(zeroevents, allevents, keep)
    }
    ##
    rm(events.study)
  }
  ##
  ## Add variable 'study' with study numbers
  ##
  dat.long$study <- as.numeric(factor(dat.long$studlab,
                                      levels = unique(dat.long$studlab)))
  dat.wide$study <- as.numeric(factor(dat.wide$studlab,
                                      levels = unique(dat.wide$studlab)))
  ##
  ## Sort by study and treatment
  ##
  dat.long <- dat.long[order(dat.long$studlab, dat.long$treat), ]
  dat.wide <- dat.wide[order(dat.wide$studlab,
                              dat.wide$treat1, dat.wide$treat2), ]
  ##
  dat.long$incr <- 0
  dat.wide$incr <- 0
  data$.incr <- 0
  ##
  ## Step iii. Drop treatment arms without or all events from individual
  ##           designs (argument 'cc.pooled' is FALSE) or add
  ##           increment if argument 'cc.pooled' is TRUE and argument
  ##           'incr' is larger than zero
  ##
  if (method != "Inverse") {
    ##
    events.arm <- with(dat.long, tapply(event, list(design, treat), sum))
    nonevents.arm <- with(dat.long, tapply(n - event, list(design, treat), sum))
    ##
    if (any(events.arm == 0, na.rm = TRUE) |
        any(nonevents.arm == 0, na.rm = TRUE)) {
      ##
      ## Identify entries without events
      ##
      zeroevents <- events.arm == 0
      zerocells <- as.data.frame(which(zeroevents, arr.ind = TRUE),
                                 stringsAsFactors = FALSE)
      ##
      zerocells$design <- rownames(zeroevents)[zerocells$row]
      zerocells$treat <- colnames(zeroevents)[zerocells$col]
      ##
      zero.long <- rep(0, nrow(dat.long))
      zero.wide <- rep(0, nrow(dat.wide))
      zero.data <- rep(0, nrow(data))
      ##
      for (i in seq_along(zerocells$design)) {
        if (!cc.pooled)
          zero.long <- zero.long + (dat.long$design == zerocells$design[i] &
                                    dat.long$treat == zerocells$treat[i])
        else
          zero.long <- zero.long + (dat.long$design == zerocells$design[i])
        ##
        zero.wide <- zero.wide + (dat.wide$design == zerocells$design[i] &
                                  (dat.wide$treat1 == zerocells$treat[i] |
                                   dat.wide$treat2 == zerocells$treat[i]))
        ##
        zero.data <- zero.data + (data$.design == zerocells$design[i] &
                                  (data$.treat1 == zerocells$treat[i] |
                                   data$.treat2 == zerocells$treat[i]))
      }
      ##
      zero.long <- zero.long > 0
      zero.wide <- zero.wide > 0
      zero.data <- zero.data > 0
      ##
      ## Identify entries with all events
      ##
      allevents <- nonevents.arm == 0
      allcells <- as.data.frame(which(allevents, arr.ind = TRUE),
                                 stringsAsFactors = FALSE)
      ##
      allcells$design <- rownames(allevents)[allcells$row]
      allcells$treat <- colnames(allevents)[allcells$col]
      ##
      all.long <- rep(0, nrow(dat.long))
      all.wide <- rep(0, nrow(dat.wide))
      all.data <- rep(0, nrow(data))
      ##
      for (i in seq_along(allcells$design)) {
        if (!cc.pooled)
          all.long <- all.long + (dat.long$design == allcells$design[i] &
                                  dat.long$treat == allcells$treat[i])
        else
          all.long <- all.long + (dat.long$design == allcells$design[i])
        ##
        all.wide <- all.wide + (dat.wide$design == allcells$design[i] &
                                (dat.wide$treat1 == allcells$treat[i] |
                                 dat.wide$treat2 == allcells$treat[i]))
        ##
        all.data <- all.data + (data$.design == allcells$design[i] &
                                (data$.treat1 == allcells$treat[i] |
                                 data$.treat2 == allcells$treat[i]))
      }
      ##
      all.long <- all.long > 0
      all.wide <- all.wide > 0
      all.data <- all.data > 0
      ##
      ## Identify entries with sparse binary data (i.e., with zero or
      ## all events)
      ##
      sparse.long <- zero.long | all.long
      sparse.wide <- zero.wide | all.wide
      sparse.data <- zero.data | all.data
      ##
      if (!cc.pooled) {
        if (warn) {
          if (sum(zeroevents, na.rm = TRUE) == 1)
            warning("Treatment arm '", zerocells$treat,
                    "' without events in design '",
                    zerocells$design, "' excluded from network meta-analysis.",
                    call. = FALSE)
          else if (sum(zeroevents, na.rm = TRUE) > 1)
            warning("Treatment arms without events in a design excluded ",
                    "from network meta-analysis:\n    ",
                    paste0("'",
                           paste0(paste0(zerocells$treat, " in "),
                                  zerocells$design),
                           "'",
                           collapse = " - "),
                    call. = FALSE)
          ##
          if (sum(allevents, na.rm = TRUE) == 1)
            warning("Treatment arm '", allcells$treat,
                    "' with all events in design '",
                    allcells$design, "' excluded from network meta-analysis.",
                    call. = FALSE)
          else if (sum(allevents, na.rm = TRUE) > 1)
            warning("Treatment arms with all events in a design excluded ",
                    "from network meta-analysis:\n    ",
                    paste0("'",
                           paste0(paste0(allcells$treat, " in "),
                                  allcells$design),
                           "'",
                           collapse = " - "),
                    call. = FALSE)
        }
        ##
        dat.long <- dat.long[!sparse.long, , drop = FALSE]
        dat.wide <- dat.wide[!sparse.wide, , drop = FALSE]
        data$.drop <- data$.drop | sparse.data
      }
      else {
        dat.long$incr[sparse.long] <- incr
        ##
        dat.long$event[sparse.long] <- dat.long$event[sparse.long] + incr
        dat.long$non.event[sparse.long] <- dat.long$non.event[sparse.long] + incr
        dat.long$n[sparse.long] <- dat.long$n[sparse.long] + 2 * incr
        ##
        dat.wide$incr[sparse.wide] <- incr
        ##
        data$.incr[sparse.data] <- incr
      }
      ##
      rm(zeroevents, zerocells, zero.long, zero.wide, zero.data,
         allevents, allcells, all.long, all.wide, all.data,
         sparse.long, sparse.wide, sparse.data)
    }
    ##
    rm(events.arm, nonevents.arm)
  }
  ##
  ## Step iv. Remove designs with single treatment arm from dataset
  ##
  if (method != "Inverse") {
    d.single <- with(dat.long, tapply(treat, design, lengthunique))
    ##
    if (any(d.single == 1, na.rm = TRUE)) {
      single <- !is.na(d.single) & d.single == 1
      design.single <- names(d.single)[single]
      ##
      if (warn)
        if (sum(single) == 1)
          warning("Design '", design.single,
                  "' with single treatment arm excluded ",
                  "from network meta-analysis.",
                  call. = FALSE)
        else
          warning("Designs with single treatment arm excluded ",
                  "from network meta-analysis: ",
                  paste(paste0("'", design.single, "'"),
                        collapse = " - "),
                  call. = FALSE)
      ##
      dat.long <- dat.long[!(dat.long$design %in% design.single), , drop = FALSE]
      dat.wide <- dat.wide[!(dat.wide$design %in% design.single), , drop = FALSE]
      ##
      data$.drop <- data$.drop | data$.design %in% design.single
      ##
      rm(single, design.single)
    }
    ##
    rm(d.single)
  }
  ##
  dat.long <- dat.long[, c("studlab", "treat",
                           "event", "non.event", "n", "incr",
                           "design", "study", ".order")]
  dat.wide <- dat.wide[, c("studlab", "treat1", "treat2",
                           "event1", "event2", "non.event1", "non.event2",
                           "n1", "n2", "incr",
                           "design", "study", ".order")]
  ##
  dat.long$studlab <- as.character(dat.long$studlab)
  dat.long$design <- as.character(dat.long$design)
  dat.long$treat <- as.character(dat.long$treat)
  ##
  dat.wide$studlab <- as.character(dat.wide$studlab)
  dat.wide$design <- as.character(dat.wide$design)
  dat.wide$treat1 <- as.character(dat.wide$treat1)
  dat.wide$treat2 <- as.character(dat.wide$treat2)
  ##
  o <- order(dat.wide$.order)
  studlab <- dat.wide$studlab[o]
  treat1 <- dat.wide$treat1[o]
  treat2 <- dat.wide$treat2[o]
  ##
  event1 <- dat.wide$event1[o]
  event2 <- dat.wide$event2[o]
  n1 <- dat.wide$n1[o]
  n2 <- dat.wide$n2[o]
  ##
  trts <- sort(unique(c(treat1, treat2)))
  n.treat <- length(trts)
  n.studies <- length(unique(studlab))
  m <- length(studlab)
  ##
  treat1.pos <- match(treat1, trts)
  treat2.pos <- match(treat2, trts)
  ##
  ## Calculate adjacency matrix
  ##
  B <- createB(treat1.pos, treat2.pos, ncol = n.treat)
  ##
  M <- t(B) %*% B    # unweighted Laplacian matrix
  D <- diag(diag(M)) # diagonal matrix
  A <- D - M         # adjacency matrix (n x n)
  ##
  rownames(A) <- colnames(A) <- trts
  ##
  rm(B, M, D)
  ##
  ## Empty matrices for results
  ##
  NAmatrix <- matrix(NA, nrow = n.treat, ncol = n.treat)
  rownames(NAmatrix) <- colnames(NAmatrix) <- trts
  ##
  TE.common <- seTE.common <- NAmatrix
  TE.direct.common <- seTE.direct.common <- NAmatrix
  Q.direct <- tau2.direct <- I2.direct <- NAmatrix
  
  
  ##
  ##
  ## (7) Conduct classic network meta-analysis using inverse variance
  ##     method
  ##
  ##
  dat.iv <- dat.long[order(dat.long$.order), ]
  ##
  if (method == "Inverse") {
    if (missing(warn))
      warn.iv <- gs("warn")
    else
      warn.iv <- warn
    incr.iv <- incr
    allstudies.iv <- allstudies
  }
  else {
    warn.iv <- FALSE
    incr.iv <- incr * cc.pooled
    allstudies.iv <- cc.pooled
  }
  ##
  p.iv <- pairwise(studlab = dat.iv$studlab,
                   treat = dat.iv$treat,
                   event = dat.iv$event,
                   n = dat.iv$n,
                   data = dat.iv,
                   sm = sm,
                   incr = incr.iv,
                   allincr = allincr, addincr = addincr,
                   allstudies = allstudies.iv)
  ##
  net.iv <- netmeta(p.iv,
                    level = level, level.ma = level.ma,
                    common = common, random = random,
                    prediction = prediction, level.predict = level.predict,
                    reference.group = reference.group,
                    baseline.reference = baseline.reference,
                    all.treatments = all.treatments,
                    seq = seq,
                    tau.preset = tau.preset,
                    tol.multiarm = tol.multiarm,
                    tol.multiarm.se = tol.multiarm.se,
                    details.chkmultiarm = details.chkmultiarm,
                    sep.trts = sep.trts,
                    nchar.trts = nchar.trts,
                    backtransf = backtransf,
                    title = title,
                    keepdata = keepdata,
                    warn = warn.iv)
  ##
  if (method == "Inverse")
    return(net.iv)
  else {
    net.iv$Cov.common[!is.na(net.iv$Cov.common)] <- NA
    net.iv$Cov.random[!is.na(net.iv$Cov.random)] <- NA
  }
  
  
  ##
  ##
  ## (8) Stage 2: Direct meta-analyses per design (MH and NCH methods)
  ##
  ##
  n.d <- tapply(dat.long$treat,   dat.long$design, lengthunique)
  k.d <- tapply(dat.long$studlab, dat.long$design, lengthunique)
  ##
  designs <- names(k.d)
  d <- length(designs)
  seq.d <- seq_len(d)
  ##
  dat.design <- vector("list", d)
  names(dat.design) <- designs
  ##
  for (i in seq.d)
    dat.design[[i]] <- dat.long[dat.long$design == designs[i], , drop = FALSE]
  ##
  if (method == "MH") {
    ##
    ## MH method
    ##
    treat.per.design <- vector("list", d)
    studies.in.design <- vector("list", d)
    ##
    c.xy <- vector("list", d)
    d.xy <- vector("list", d)
    C.xy <- vector("list", d)
    L.xy <- vector("list", d)
    ##
    U.xyy <- vector("list", d)
    U.xyz <- vector("list", d)
    ##
    t.pl <- vector("list", d)
    ##
    ps1 <- ps2 <- ps3 <- ps4 <- vector("list", d)
    ##
    L.bar.xy <- vector("list", d)
    ##
    U.plus.xx <- vector("list", d)
    ##
    U.new <- vector("list", d)
    ##
    U.plus.xy <- vector("list", d)
    ##
    Var.Lbar <- vector("list", d)
    ##
    CoVar.Lbar <- vector("list", d)
    ##
    length.y <- sum(n.d - 1)
    ##
    y <- c(rep(0, length.y))
    ##
    V1 <- matrix(rep(0, length.y * length.y), nrow = length.y)
    V2 <- matrix(rep(0, length.y * length.y), nrow = length.y)
    ##
    X <- matrix(0, nrow = length.y, ncol = n.treat - 1)
    ##
    counter1 <- counter2 <- counter3 <- 0
    ##
    list1 <- matrix(0, nrow = length.y, ncol = 2)
    ##
    N.j <- rep(0, d)
    ##
    if (d > 1)
      for (i in 2:d)
        N.j[i] <- N.j[i - 1] + n.d[i - 1] - 1
    ##
    ## Loop for designs
    ##
    for (i in seq.d) {
      treat.per.design[[i]] <- unique(dat.design[[i]]$treat)
      studies.in.design[[i]] <- unique(dat.design[[i]]$study)
      ##
      ## List of studies in each design
      ##
      ## Recode the studies in each design
      ##
      for (j in seq_along(dat.design[[i]]$study))
        for (k in seq_len(k.d[i]))
          if (dat.design[[i]]$study[j] == studies.in.design[[i]][k])
            dat.design[[i]]$id.s[j] <- k
      ##
      ## Recode the treatments in each design
      ##
      for (j in seq_along(dat.design[[i]]$study))
        for (k in seq_len(n.d[i]))
          if (dat.design[[i]]$treat[j] == treat.per.design[[i]][k])
            dat.design[[i]]$id.t[j] <- k
      ##
      ## Calculate total number of patients per studies
      ##
      for (j in seq_along(dat.design[[i]]$studlab))
        dat.design[[i]]$n.study[j] <-
          with(dat.design[[i]], sum(n[which(study == study[j])]))
      ##
      ## Define c.xy and d.xy
      ##
      c.xy[[i]] <- array(rep(0, k.d[i] * n.d[i] * n.d[i]),
                         dim = c(n.d[i], n.d[i], k.d[i]))
      ##
      d.xy[[i]] <- array(rep(0, k.d[i] * n.d[i] * n.d[i]),
                         dim = c(n.d[i], n.d[i], k.d[i]))
      ##
      for (st in seq_len(k.d[i]))
        for (t1 in seq_len(n.d[i]))
          for (t2 in seq_len(n.d[i])) {
            c.xy[[i]][t1, t2, st] <-
              with(dat.design[[i]],
                   sum(    event[which(id.s == st & id.t == t1)]) *
                   sum(non.event[which(id.s == st & id.t == t2)]) /
                   n.study[id.s == st][1])
            ##
            d.xy[[i]][t1, t2, st] <-
              with(dat.design[[i]],
              (sum(    event[which(id.s == st & id.t == t1)]) +
               sum(non.event[which(id.s == st & id.t == t2)])) /
              n.study[id.s == st][1])
          }
      ##
      ## Define C.xy
      ##
      C.xy[[i]] <- matrix(rep(0, n.d[i] * n.d[i]), nrow = n.d[i])
      ##
      for (j in seq_len(n.d[i]))
        for (k in seq_len(n.d[i]))
          C.xy[[i]][j, k] <- sum(c.xy[[i]][j, k, ])
      ##
      ## Define L.xy
      ##
      L.xy[[i]] <- matrix(rep(0, n.d[i] * n.d[i]), nrow = n.d[i])
      for (j in seq_len(n.d[i]))
        for (k in seq_len(n.d[i]))
          L.xy[[i]][j, k] <- log(C.xy[[i]][j, k] / C.xy[[i]][k, j])
      ##
      ## Calculate the variance of L.xy (dimension n.d[i] x n.d[i])
      ##
      U.xyy[[i]] <- matrix(rep(0, n.d[i] * n.d[i]), nrow = n.d[i])
      for (j in seq_len(n.d[i]))
        for (k in seq_len(n.d[i]))
          U.xyy[[i]][j, k] <-
            sum(c.xy[[i]][j, k, ] * d.xy[[i]][j, k, ]) /
            (2 * C.xy[[i]][j, k]^2) +
            sum(c.xy[[i]][j, k, ] * d.xy[[i]][k, j, ] + c.xy[[i]][k, j, ] *
                d.xy[[i]][j, k, ]) / (2 * C.xy[[i]][j, k] * C.xy[[i]][k, j]) +
            sum(c.xy[[i]][k, j, ] * d.xy[[i]][k, j, ]) / (2 * C.xy[[i]][k, j]^2)
      ##
      ## Calculate the covariance matrix U.xyz
      ##
      t.pl[[i]] = array(rep(0, k.d[i] * n.d[i] * n.d[i]),
                        dim = c(n.d[i], n.d[i], k.d[i]))
      ##
      for (j in seq_len(k.d[i]))
        t.pl[[i]][ , , j] <- with(dat.design[[i]], n.study[id.s == j][1])
      ##
      U.xyz[[i]] <- array(rep(0, n.d[i] * n.d[i] * n.d[i]),
                          dim = c(n.d[i], n.d[i], n.d[i]))
      ##
      ## Per study ...
      ##
      for (t1 in seq_len(n.d[i])) {
        for (t2 in seq_len(n.d[i])) {
          for (t3 in seq_len(n.d[i])) {
            for (st in seq_len(k.d[i])) {
              sel1 <- with(dat.design[[i]], which(id.s == st & id.t == t1))
              sel2 <- with(dat.design[[i]], which(id.s == st & id.t == t2))
              sel3 <- with(dat.design[[i]], which(id.s == st & id.t == t3))
              ##
              ps1[[i]][st] <-
                with(dat.design[[i]],
                     event[sel1] /
                     (t.pl[[i]][1, 1, st])^2 *
                     non.event[sel2] * non.event[sel3] *
                     (t1 != t2) * (t1 != t3) * (t2 != t3))
              ##
              ps2[[i]][st] <-
                with(dat.design[[i]],
                     n[sel1] / (t.pl[[i]][1, 1, st])^2 *
                     non.event[sel2] * event[sel3] *
                     (t1 != t2) * (t1 != t3) * (t2 != t3))
              ##
              ps3[[i]][st] <-
                with(dat.design[[i]],
                     n[sel1] / (t.pl[[i]][1, 1, st])^2 *
                     event[sel2] * non.event[sel3] *
                     (t1 != t2) * (t1 != t3) * (t2 != t3))
              ##
              ps4[[i]][st] <-
                with(dat.design[[i]],
                     non.event[sel1] / (t.pl[[i]][1, 1, st])^2 *
                     event[sel2] * event[sel3] *
                     (t1 != t2) * (t1 != t3) * (t2 != t3))
            }
            ##
            U.xyz[[i]][t1, t2, t3] <-
              sum(ps1[[i]][]) / (3 * C.xy[[i]][t1, t2] * C.xy[[i]][t1, t3]) +
              sum(ps2[[i]][]) / (3 * C.xy[[i]][t1, t2] * C.xy[[i]][t3, t1]) +
              sum(ps3[[i]][]) / (3 * C.xy[[i]][t2, t1] * C.xy[[i]][t1, t3]) +
              sum(ps4[[i]][]) / (3 * C.xy[[i]][t2, t1] * C.xy[[i]][t3, t1])
          }
        }
      }
      ##
      ## Calculate L.bar.xy
      ##
      L.bar.xy[[i]] <- matrix(rep(0, n.d[i] * n.d[i]), nrow = n.d[i])
      for (t1 in seq_len(n.d[i]))
        for (t2 in seq_len(n.d[i]))
          L.bar.xy[[i]][t1, t2] <-
            (sum(L.xy[[i]][t1, ]) - sum(L.xy[[i]][t2, ])) / n.d[i]
      ##
      ## Calculate U.plus.xx
      ##
      for (t1 in seq_len(n.d[i]))
        U.xyy[[i]][t1, t1] <- 0
      ##
      U.plus.xx[[i]] <- rep(0, n.d[i])
      for (t1 in seq_len(n.d[i]))
        U.plus.xx[[i]][t1] <-
          sum(U.xyy[[i]][t1, seq_len(n.d[i])]) + sum(U.xyz[[i]][t1, , ])
      ##
      ## Calculate U.plus.xy
      ##
      U.new[[i]] <- array(rep(0, n.d[i] * n.d[i] * n.d[i]),
                          dim = c(n.d[i], n.d[i], n.d[i]))
      ##
      for (t1 in seq_len(n.d[i]))
        for (t2 in seq_len(n.d[i]))
          for (t3 in seq_len(n.d[i]))
            U.new[[i]][t1, t2, t3] <-
              (t1 != t2) * (t1 != t3) * (t2 != t3) *
              U.xyz[[i]][t1, t2, t3] +
              (t1 != t2) * (t2 == t3) * U.xyy[[i]][t1, t2]
      ##
      U.plus.xy[[i]] = matrix(rep(0, n.d[i] * n.d[i]), nrow = n.d[i])
      ##
      for (t1 in seq_len(n.d[i]))
        for (t2 in seq_len(n.d[i]))
          U.plus.xy[[i]][t1, t2] <-
            sum(U.new[[i]][seq_len(n.d[i]), t1, t2]) -
            sum(U.new[[i]][t1, t2, seq_len(n.d[i])]) -
            sum(U.new[[i]][t2, t1, seq_len(n.d[i])]) +
            U.new[[i]][t1, t2, t2]
      ##
      ## Variance of L.bar.xy
      ##
      Var.Lbar[[i]] <- matrix(rep(0, n.d[i] * n.d[i]), nrow = n.d[i])
      for (t1 in seq_len(n.d[i]))
        for (t2 in seq_len(n.d[i]))
          Var.Lbar[[i]][t1, t2] <-
            (U.plus.xx[[i]][t1] - 2 * U.plus.xy[[i]][t1, t2] +
             U.plus.xx[[i]][t2]) / n.d[i]^2
      ##
      ## Covariance of L.bar.xy. Only a subset of covariances are
      ## calculate here, i.e. cov(L.bar_(1, t1), L.bar_(1, t2))
      ##
      CoVar.Lbar[[i]] <- matrix(rep(0, n.d[i]^2), nrow = n.d[i])
      for (t1 in seq_len(n.d[i]))
        for (t2 in seq_len(n.d[i]))
          CoVar.Lbar[[i]][t1, t2] <-
            (U.plus.xx[[i]][1] - U.plus.xy[[i]][1, t2] -
             U.plus.xy[[i]][t1, 1] +
             U.plus.xy[[i]][t1, t2]) / n.d[i]^2
      ##
      ## Define matrix X
      ##
      for (j in seq_along(dat.design[[i]]$studlab))
        for (k in seq_along(trts))
          if (dat.design[[i]]$treat[j] == trts[k])
            dat.design[[i]]$id.t2[j] <- k
      ##
      ##
      ## Create y, the vector of treatment effects from each study. Only
      ## effects vs the first treatment are needed. y is coded as OR
      ## 1vsX
      ##
      for (j in seq_len(n.d[i] - 1))
        y[j + counter1] <- L.bar.xy[[i]][1, j + 1]
      ##
      counter1 <- counter1 + n.d[i] - 1
      ##
      ## Create V, the matrix of covariances of treatment effects from
      ## each study. Only covariances of effects vs the first treatment
      ## are needed.
      ##
      for (j in seq_len(n.d[i] - 1))
        for (k in seq_len(n.d[i] - 1))
          V1[j + counter2, k + counter2] <- (k == j) * Var.Lbar[[i]][1, j + 1]
      ##
      counter2 <- counter2 + n.d[i] - 1
      ##
      for (j in 2:(n.d[i]))
        for (k in 2:(n.d[i]))
          V2[j + counter3 - 1, k + counter3 - 1] <-
            (k != j) * CoVar.Lbar[[i]][j, k]
      ##
      counter3 <- counter3 + n.d[i] - 1
      ##
      for (j in seq_len(n.d[i] - 1)) {
        list1[N.j[i] + j, 1] <- dat.design[[i]]$id.t2[[1]]
        list1[N.j[i] + j, 2] <- dat.design[[i]]$id.t2[[j + 1]]
      }
      ##
      ## Drop unnecessary variables
      ##
      dat.design[[i]]$id.s <- NULL
      dat.design[[i]]$id.t <- NULL
      dat.design[[i]]$id.t2 <- NULL
      dat.design[[i]]$n.study <- NULL
    }
    ##
    V <- V1 + V2
    ##
    basic.contrasts <- 2:n.treat
    ##
    for (i in 1:length.y)
      for (k in 1:(n.treat - 1)) {
        if (list1[i, 1] == basic.contrasts[k])
          X[i, k] = -1
        if (list1[i, 2] == basic.contrasts[k])
          X[i, k] = 1
      }
    ##
    W <- solve(V)
    ##
    TE.basic <- solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% y
  }
  else if (method == "NCH") {
    ##
    ## NCH method
    ##
    dat.long$ttt <- 0
    ##
    for (k in 1:n.treat)
      for (i in seq_along(dat.long$studlab))
        if (dat.long$treat[i] == trts[k])
          dat.long$ttt[i] <- k
    ##
    dat.long$treat <- as.numeric(dat.long$ttt)
    dat.long$ttt <- NULL
    ##
    dat.long$treat <- as.numeric(dat.long$treat)
    ##
    dat.long <- dat.long[order(dat.long$studlab, dat.long$treat), ] # necessary ???
    ##
    for (i in unique(dat.long$studlab)) {
      counter <- 1
      for (j in seq_along(dat.long$studlab))
        if (dat.long$studlab[j] == i) {
          dat.long$count[j] <- counter
          counter <- counter + 1
        }
    }
    ##
    max.arms <- max(dat.long$count)
    ##
    d1 <- reshape(dat.long, idvar = "studlab",
                  timevar = "count", direction = "wide")
    ##
    d1 <- d1[order(d1$studlab), ]
    ##
    d1$N1 <- 0 # d1$n.all = 0
    d1$narms <- 0
    ##
    for (i in seq_len(max.arms)) {
      ev <- paste0("event.", i)
      tot <- paste0("n.", i)
      d1$N1 <- rowSums(cbind(d1$N1, d1[, colnames(d1) == ev]),
                       na.rm = TRUE)
      ## d1$n.all <- rowSums(cbind(d1$n.all, d1[, colnames(d1) == tot]), na.rm = TRUE)
    }
    ##
    for (i in seq_along(d1$studlab))
      for (k in 1:max.arms) {
        ev <- paste0("event.", k)
        if (!is.na(d1[, colnames(d1) == ev][i]))
          d1$narms[i] <- k
      }
    ##
    for (i in seq_along(colnames(d1)))
      for (k in seq_along(colnames(d1))) {
        if (colnames(d1)[k] == paste0("treat.", i))
          colnames(d1)[k] = paste0("t", i)
        ##
        if (colnames(d1)[k] == paste0("event.", i))
          colnames(d1)[k] = paste0("r", i)
        ##
        if (colnames(d1)[k] == paste0("n.", i))
          colnames(d1)[k] = paste0("n", i)
      }
    ##
    ## Likelihood function
    ##
    myLik1 <- function(mypar) {
      x <- mypar
      ##
      myLogLik1 <- myLogLik2 <- myLogLik3 <- rep(0, length(d1$studlab))
      ##
      for (i in seq_along(d1$studlab)) {
        for (k in 2:d1$narms[i]) {
          myLogLik1[i] <-
            myLogLik1[i] +
            d1[, colnames(d1) == paste0("r", k)][i] *
            (x[d1[, colnames(d1) == paste0("t", k)][i]] -
             x[d1$t1[i]] * (d1$t1[i] != 1))
          ##
          myLogLik2[i] <-
            myLogLik2[i] +
            d1[, colnames(d1) == paste0("n", k)][i] *
            exp(x[d1[, colnames(d1) == paste0("t", k)][i]] -
                x[d1$t1[i]] * (d1$t1[i] != 1))
          ##
          myLogLik3[i] <- -d1$N1[i] * log(d1$n1[i] + myLogLik2[i])
        }
      }
      ##
      myLogLik <- sum(myLogLik1 + myLogLik3)
      ##
      myLogLik
    }
    ##
    opt <- optim(rep(0, n.treat), myLik1, method = "L-BFGS-B",
                 lower = -Inf, upper = Inf,
                 control = list(fnscale = -1, maxit = 10000),
                 hessian = TRUE)
    ##
    W <- solve(-opt$hessian[2:n.treat, 2:n.treat])
    ##
    TE.basic <- -(opt$par)[2:n.treat]
  }
  ##
  ## Define H matrix
  ##
  H <- matrix(0,
              nrow = n.treat * ((n.treat - 1) / 2),
              ncol = n.treat - 1)
  ##
  diag(H) <- 1
  ##
  if (n.treat > 2) {
    t1 <- c()
    t2 <- c()
    for (i in 2:(n.treat - 1))
      for (j in (i + 1):(n.treat)) {
        t1 <- rbind(t1, i)
        t2 <- rbind(t2, j)
      }
    ##
    h1 <- matrix(c(t1, t2), nrow = nrow(t1))
    ##
    for (i in 1:((n.treat - 1) * (n.treat - 2) / 2))
      for (j in 1:(n.treat - 1))
        H[i + n.treat - 1, j] <- -(h1[i, 1] == j + 1) + (h1[i, 2] == j + 1)
  }
  ##
  ## Common effects matrix
  ##
  d.hat <- H %*% TE.basic
  ##
  TE.common[lower.tri(TE.common, diag = FALSE)] <-  d.hat
  TE.common <- t(TE.common)
  TE.common[lower.tri(TE.common, diag = FALSE)] <- -d.hat
  diag(TE.common) <- 0
  ##
  ## Matrix with standard errors
  ##
  if (method == "MH")
    cov.d.hat <- H %*% solve(t(X) %*% W %*% X) %*% t(H)
  else
    cov.d.hat <- H %*% W %*% t(H)
  ##
  seTE.common[lower.tri(seTE.common, diag = FALSE)] <- sqrt(diag(cov.d.hat))
  seTE.common <- t(seTE.common)
  seTE.common[lower.tri(seTE.common, diag = FALSE)] <- sqrt(diag(cov.d.hat))
  diag(seTE.common) <- 0
  ##
  ## Confidence intervals
  ##
  ci.f <- ci(TE.common, seTE.common, level = level.ma)
  ##
  ## Inconsistency global
  ##
  if (method == "MH") {
    Q <- as.vector(t(y - X %*% TE.basic) %*% solve(V) %*%
                    (y - X %*% TE.basic))
    df.Q <- sum(n.d - 1) - (n.treat - 1)
    pval.Q <- pvalQ(Q, df.Q)
  }
  else {
    Q <- NA
    df.Q <- NA
    pval.Q <- NA
  }
  
  
  ##
  ##
  ## (9) Inconsistency evaluation: direct MH estimates
  ##
  ##
  B.matrix <- createB(treat1.pos, treat2.pos)
  B.matrix <- unique(B.matrix)
  colnames(B.matrix) <- trts
  ##
  for (i in seq_len(nrow(B.matrix))) {
    sel.treat1 <- colnames(B.matrix)[B.matrix[i, ] ==  1]
    sel.treat2 <- colnames(B.matrix)[B.matrix[i, ] == -1]
    selstud <- treat1 == sel.treat1 & treat2 == sel.treat2
    ##
    m.i <- metabin(event1, n1, event2, n2, subset = selstud,
                   method = "MH", sm = "OR",
                   incr = incr, allincr = allincr, addincr = addincr,
                   allstudies = allstudies, MH.exact = !cc.pooled,
                   method.tau = "DL", method.tau.ci = "",
                   Q.Cochrane = FALSE,
                   warn.deprecated = FALSE)
    ##
    TE.i   <- m.i$TE.common
    seTE.i <- m.i$seTE.common
    Q.i <- m.i$Q
    tau2.i <- m.i$tau2
    I2.i <- m.i$I2
    ##
    TE.direct.common[sel.treat1, sel.treat2]   <- TE.i
    seTE.direct.common[sel.treat1, sel.treat2] <- seTE.i
    Q.direct[sel.treat1, sel.treat2] <- Q.i
    tau2.direct[sel.treat1, sel.treat2] <- tau2.i
    I2.direct[sel.treat1, sel.treat2] <- I2.i
    ##
    TE.direct.common[sel.treat2, sel.treat1]   <- -TE.i
    seTE.direct.common[sel.treat2, sel.treat1] <- seTE.i
    Q.direct[sel.treat2, sel.treat1] <- Q.i
    tau2.direct[sel.treat2, sel.treat1] <- tau2.i
    I2.direct[sel.treat2, sel.treat1] <- I2.i
  }
  ##
  rm(sel.treat1, sel.treat2, selstud, m.i, TE.i, seTE.i)
  ##
  ci.d <- ci(TE.direct.common, seTE.direct.common, level = level.ma)
  
  
  labels <- sort(unique(c(treat1, treat2)))
  ##
  if (!is.null(seq))
    seq <- setseq(seq, labels)
  else {
    seq <- labels
    if (is.numeric(seq))
      seq <- as.character(seq)
  }
  
  
  res <- list(studlab = studlab,
              treat1 = treat1,
              treat2 = treat2,
              ##
              TE = data$.TE[!data$.drop],
              seTE = data$.seTE[!data$.drop],
              seTE.adj = rep(NA, sum(!data$.drop)),
              seTE.adj.common = rep(NA, sum(!data$.drop)),
              seTE.adj.random = rep(NA, sum(!data$.drop)),
              ##
              design = designs(treat1, treat2, studlab)$design,
              ##
              event1 = event1,
              event2 = event2,
              n1 = n1,
              n2 = n2,
              ##
              k = n.studies,
              m = m,
              n = length(trts),
              d = d,
              ##
              trts = trts,
              k.trts = rowSums(A),
              n.trts = NA,
              events.trts = NA,
              ##
              studies = NA,
              narms = NA,
              ##
              designs = designs,
              ##
              TE.common = TE.common,
              seTE.common = seTE.common,
              lower.common = ci.f$lower,
              upper.common = ci.f$upper,
              statistic.common = ci.f$statistic,
              pval.common = ci.f$p,
              ##
              TE.random = NAmatrix,
              seTE.random = NAmatrix,
              lower.random = NAmatrix,
              upper.random = NAmatrix,
              statistic.random = NAmatrix,
              pval.random = NAmatrix,
              ##
              seTE.predict = NAmatrix,
              lower.predict = NAmatrix,
              upper.predict = NAmatrix,
              ##
              prop.direct.common = NA,
              prop.direct.random = NA,
              ##
              TE.direct.common = TE.direct.common,
              seTE.direct.common = seTE.direct.common,
              lower.direct.common = ci.d$lower,
              upper.direct.common = ci.d$upper,
              statistic.direct.common = ci.d$statistic,
              pval.direct.common = ci.d$p,
              ##
              TE.direct.random = NAmatrix,
              seTE.direct.random = NAmatrix,
              lower.direct.random = NAmatrix,
              upper.direct.random = NAmatrix,
              statistic.direct.random = NAmatrix,
              pval.direct.random = NAmatrix,
              ##
              Q.direct = Q.direct,
              tau.direct = sqrt(tau2.direct),
              tau2.direct = tau2.direct,
              I2.direct = I2.direct,
              ##
              TE.indirect.common = NA,
              seTE.indirect.common = NA,
              lower.indirect.common = NA,
              upper.indirect.common = NA,
              statistic.indirect.common = NA,
              pval.indirect.common = NA,
              ##
              TE.indirect.random = NA,
              seTE.indirect.random = NA,
              lower.indirect.random = NA,
              upper.indirect.random = NA,
              statistic.indirect.random = NA,
              pval.indirect.random = NA,
              ##
              Q = NA,
              df.Q = NA,
              pval.Q = NA,
              I2 = NA, lower.I2 = NA, upper.I2 = NA,
              tau = NA,
              ##
              Q.heterogeneity = NA,
              df.Q.heterogeneity = NA,
              pval.Q.heterogeneity = NA,
              Q.inconsistency = Q,
              df.Q.inconsistency = df.Q,
              pval.Q.inconsistency = pval.Q,
              ##
              Q.decomp = NA,
              ##
              A.matrix = A,
              H.matrix = H,
              ##
              n.matrix = NA,
              events.matrix = NA,
              ##
              P.common = NA,
              P.random = NA,
              ##
              Cov.common = net.iv$Cov.common,
              Cov.random = net.iv$Cov.random,
              ##
              sm = sm,
              method = method,
              ##
              incr = incr,
              allincr = allincr,
              addincr = addincr,
              allstudies = allstudies,
              cc.pooled = cc.pooled,
              ##
              level = level,
              level.ma = level.ma,
              common = common,
              random = random,
              ##
              prediction = prediction,
              level.predict = level.predict,
              ##
              reference.group = reference.group,
              baseline.reference = baseline.reference,
              all.treatments = all.treatments,
              seq = seq,
              ##
              tau.preset = tau.preset,
              ##
              tol.multiarm = tol.multiarm,
              tol.multiarm.se = tol.multiarm.se,
              details.chkmultiarm = details.chkmultiarm,
              details.chkdata = details.chkdata,
              ##
              sep.trts = sep.trts,
              nchar.trts = nchar.trts,
              ##
              func.inverse = deparse(substitute(func.inverse)),
              ##
              backtransf = backtransf,
              ##
              title = title,
              ##
              data = if (keepdata) data else NULL,
              data.design = dat.design,
              ##
              warn = warn,
              call = match.call(),
              version = packageDescription("netmeta")$Version
              )
  ##
  class(res) <- c("netmetabin", "netmeta")
  ##
  ## Add results for indirect treatment estimates
  ##
  n <- res$n
  ##
  res$prop.direct.common <-
    netmeasures(res, random = FALSE, warn = warn)$proportion
  if (is.logical(res$prop.direct.common))
    res$prop.direct.common <- as.numeric(res$prop.direct.common)
  ##
  P.common <- P.random <- matrix(NA, n, n)
  colnames(P.common) <- rownames(P.common) <-
    colnames(P.random) <- rownames(P.random) <- trts
  ##
  if (n == 2) {
    ##
    ## For two treatments only direct evidence is available
    ##
    res$prop.direct.common <- 1
    names(res$prop.direct.common) <- paste(labels, collapse = sep.trts)
    ##
    sel <- row(P.common) != col(P.common)
    P.common[sel] <- 1
  }
  else {
    k <- 0
    for (i in 1:(n - 1)) {
      for (j in (i + 1):n) {
        k <- k + 1
        P.common[i, j] <- P.common[j, i] <- res$prop.direct.common[k]
      }
    }
  }
  ##
  ## Set direct evidence estimates to 0 if only indirect evidence is available
  ## (otherwise indirect estimates would be NA as direct estimates are NA)
  ##
  TE.direct.common <- res$TE.direct.common
  ##
  TE.direct.common[abs(P.common) < .Machine$double.eps^0.5] <- 0
  ##
  ## Indirect estimate is NA if only direct evidence is available
  ##
  res$P.common <- P.common
  ##
  P.common[abs(P.common - 1) < .Machine$double.eps^0.5] <- NA
  P.common[P.common > 1] <- NA
  ##
  ## Common effects model
  ##
  ci.if <- ci((res$TE.common - P.common * TE.direct.common) / (1 - P.common),
              sqrt(res$seTE.common^2 / (1 - P.common)),
              level = level)
  ##
  res$TE.indirect.common   <- ci.if$TE
  res$seTE.indirect.common <- ci.if$seTE
  ##
  res$lower.indirect.common <- ci.if$lower
  res$upper.indirect.common <- ci.if$upper
  ##
  res$statistic.indirect.common <- ci.if$statistic
  res$pval.indirect.common <- ci.if$p
  ##
  ## No results for random effects model
  ##
  res$prop.direct.random <- res$prop.direct.common
  res$prop.direct.random[!is.na(res$prop.direct.random)] <- NA
  res$P.random <- P.random
  ##
  res$TE.indirect.random <- res$seTE.indirect.random <-
    res$lower.indirect.random <- res$upper.indirect.random <-
      res$statistic.indirect.random <- res$pval.indirect.random <-
        NAmatrix
  
  
  ##
  ## Study overview
  ##
  p0 <- prepare(rep(1, nrow(dat.wide)),
                rep(1, nrow(dat.wide)),
                dat.wide$treat1,
                dat.wide$treat2,
                dat.wide$studlab,
                func.inverse = func.inverse)
  ##
  tdata <- data.frame(studies = p0$studlab, narms = p0$narms,
                      stringsAsFactors = FALSE)
  ##
  tdata <- tdata[!duplicated(tdata[, c("studies", "narms")]), , drop = FALSE]
  res$studies <- tdata$studies
  res$narms <- tdata$narms
  ##
  if (all(suppressWarnings(!is.na(as.numeric(res$studies))))) {
    o <- order(as.numeric(res$studies))
    res$studies <- res$studies[o]
    res$narms <- res$narms[o]
  }
  ##
  res$data <- merge(res$data,
                    data.frame(.studlab = res$studies,
                               .narms = res$narms),
                    by = ".studlab", all.x = TRUE)
  
  
  res$data <- res$data[order(res$data$.order), ]
  res$data$.order <- NULL
  rownames(res$data) <- seq_len(nrow(res$data))
  
  
  res$n.matrix <- netmatrix(res, n1 + n2, func = "sum")
  ##
  dat.n <- data.frame(studlab = c(studlab, studlab),
                      treat = c(treat1, treat2),
                      n = c(n1, n2))
  dat.n <- dat.n[!duplicated(dat.n[, c("studlab", "treat")]), ]
  dat.n <- by(dat.n$n, dat.n$treat, sum, na.rm = TRUE)
  res$n.trts <- as.vector(dat.n[trts])
  names(res$n.trts) <- trts
  ##
  res$events.matrix <- netmatrix(res, event1 + event2, func = "sum")
  ##
  dat.e <- data.frame(studlab = c(studlab, studlab),
                      treat = c(treat1, treat2),
                      n = c(event1, event2))
  dat.e <- dat.e[!duplicated(dat.e[, c("studlab", "treat")]), ]
  dat.e <- by(dat.e$n, dat.e$treat, sum, na.rm = TRUE)
  res$events.trts <- as.vector(dat.e[trts])
  names(res$events.trts) <- trts
  ##
  ## Backward compatibility
  ##
  res$fixed <- res$common
  res$comb.fixed <- res$common
  res$comb.random <- res$random
  ##
  res$seTE.adj.fixed <- res$seTE.adj.common
  ##
  res$TE.fixed <- res$TE.common
  res$seTE.fixed <- res$seTE.common
  res$lower.fixed <- res$lower.common
  res$upper.fixed <- res$upper.common
  res$statistic.fixed <- res$statistic.common
  res$pval.fixed <- res$pval.common
  ##
  res$prop.direct.fixed <- res$prop.direct.common
  ##
  res$TE.direct.fixed <- res$TE.direct.common
  res$seTE.direct.fixed <- res$seTE.direct.common
  res$lower.direct.fixed <- res$lower.direct.common
  res$upper.direct.fixed <- res$upper.direct.common
  res$statistic.direct.fixed <- res$statistic.direct.common
  res$pval.direct.fixed <- res$pval.direct.common
  ##
  res$TE.indirect.fixed <- res$TE.indirect.common
  res$seTE.indirect.fixed <- res$seTE.indirect.common
  res$lower.indirect.fixed <- res$lower.indirect.common
  res$upper.indirect.fixed <- res$upper.indirect.common
  res$statistic.indirect.fixed <- res$statistic.indirect.common
  res$pval.indirect.fixed <- res$pval.indirect.common
  ##
  res$P.fixed <- res$P.common
  res$Cov.fixed <- res$Cov.common               
  
  
  res
}