R/sucra.R

In MathiasHarrer/dmetar: Companion R Package for the Guide 'Doing Meta-Analysis in R'

#' Calculate the Surface Under the Cumulative Ranking score of from a network meta-analysis
#'
#' This function calculates the SUCRA (Surface Under the Cumulative Ranking) score from a rank
#' probability matrix or an object of class \code{mtc.rank.probability} generated by the \code{\link[gemtc]{rank.probability}}
#' function.
#'
#' @usage sucra(x, lower.is.better = FALSE)
#'
#' @param x An object of class \code{mtc.rank.probability} generated by the \code{\link[gemtc]{rank.probability}} function
#' or a matrix/data.frame in which the rows correspond to the treatment, and columns to the
#' probability of a specific treatment having this rank (see Details). Rownames of the matrix should
#' contain the name of the specific treatment.
#' @param lower.is.better Logical. Do lower (i.e., more negative) effect sizes mean that effects are
#' higher? \code{FALSE} by default. Use the default when the provided matrix already contains the
#' correct rank probability for each treatment, and values ought not to be inverted.
#'
#' @details The SUCRA score is a metric to evaluate which treatment in a network is likely
#' to be the most efficacious in the context of network meta-analyses. The SUCRA score is calculated
#' in the function using the formula described in Salanti, Ades and Ioannidis (2011):
#' \deqn{SUCRA_j = \frac{\sum_{b=1}^{a-1}cum_{jb}}{a-1}}
#' Where \eqn{j} is some treatment, \eqn{a} are all competing treatments, \eqn{b} are the
#' \eqn{b = 1, 2, ..., a-1} best treatments, and \eqn{cum} represents the cumulative probability
#' of a treatment being among the \eqn{b} best treatments.
#'
#' Other than an object of class \code{mtc.rank.probability} for argument \code{x}, the function can also be provided
#' with a \eqn{m \times n} matrix where \eqn{m} are rows corresponding to each treatment in the
#' network meta-analysis, and the \eqn{n} columns correspond to each rank (1st, 2nd, etc.). Rank probabilities
#' should be provided as a value from 0 to 1. Rownames of the matrix should correspond to the treatment names.
#' Here is an example rank probability matrix for eight treatments:
#'
#' \tabular{lrrrrrrrr}{
#' . \tab [,1] \tab [,2] \tab [,3] \tab [,4] \tab [,5] \tab [,6] \tab [,7] \tab [,8]\cr
#' CBT \tab 0.000000 \tab 0.000000 \tab 0.000000 \tab 0.000000 \tab 0.000000 \tab 0.001275 \tab 0.087400 \tab 0.911325\cr
#' IPT \tab 0.000000 \tab 0.000000 \tab 0.000000 \tab 0.000000 \tab 0.000000 \tab 0.179400 \tab 0.745875 \tab 0.074725\cr
#' PDT \tab 0.000000 \tab 0.000000 \tab 0.000225 \tab 0.020300 \tab 0.978025 \tab 0.001450 \tab 0.000000 \tab 0.000000\cr
#' PLA \tab 0.002825 \tab 0.551175 \tab 0.262525 \tab 0.181550 \tab 0.001925 \tab 0.000000 \tab 0.000000 \tab 0.000000\cr
#' PST \tab 0.000000 \tab 0.000000 \tab 0.000000 \tab 0.000025 \tab 0.001450 \tab 0.817850 \tab 0.166725 \tab 0.013950\cr
#' SUP \tab 0.000000 \tab 0.216450 \tab 0.398700 \tab 0.383950 \tab 0.000900 \tab 0.000000 \tab 0.000000 \tab 0.000000\cr
#' TAU \tab 0.000375 \tab 0.229200 \tab 0.338525 \tab 0.414175 \tab 0.017700 \tab 0.000025 \tab 0.000000 \tab 0.000000\cr
#' WLC \tab 0.996800 \tab 0.003175 \tab 0.000025 \tab 0.000000 \tab 0.000000 \tab 0.000000 \tab 0.000000 \tab 0.000000
#' }
#'
#' @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/bayesian-network-meta-analysis.html}{Chapter 11.2}.
#'
#'Salanti, G., Ades, A. E. & Ioannidis, J.P.A. (2011). Graphical Methods and Numerical Summaries for
#'Presenting Results from Multiple-Treatment Meta-Analysis: An Overview and Tutorial.
#'\emph{Journal of Clinical Epidemiology, 64} (2): 163–71.
#'
#' @author Mathias Harrer & David Daniel Ebert
#'
#' @import ggplot2
#'
#' @export sucra
#'
#' @seealso
#' \code{\link{direct.evidence.plot}}
#'
#' @examples
#' \dontrun{
#' # Example1 : conduct NMA using gemtc, calculate SUCRAs
#' suppressPackageStartupMessages(library(gemtc))
#' suppressPackageStartupMessages(library(igraph))
#' data("NetDataGemtc")
#'
#' network = suppressWarnings(mtc.network(data.re = NetDataGemtc))
#'
#' plot(network, layout = layout.fruchterman.reingold)
#'
#' model = mtc.model(network, linearModel = "fixed",
#'                    n.chain = 4,
#'                    likelihood = "normal",
#'                    link = "identity")
#'
#' mcmc = mtc.run(model, n.adapt = 5000, n.iter = 100000, thin = 10)
#'
#' rp = rank.probability(mcmc)
#'
#' sucra = sucra(rp, lower.is.better = TRUE)
#' sucra
#' plot(sucra)}
#'
#'
#'
#' # Example 2: construct rank proabability matrix, then use sucra function
#' rp = rbind(CBT = c(0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.001500, 0.088025, 0.910475),
#'            IPT = c(0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.176975, 0.748300, 0.074725),
#'            PDT = c(0.000000, 0.000000, 0.000250, 0.021725, 0.976525, 0.001500, 0.000000, 0.000000),
#'            PLA = c(0.003350, 0.546075, 0.266125, 0.182125, 0.002325, 0.000000, 0.000000, 0.000000),
#'            PST = c(0.000000, 0.000000, 0.000000, 0.000000, 0.001500, 0.820025, 0.163675, 0.014800),
#'            SUP = c(0.000000, 0.217450, 0.403950, 0.378000, 0.000600, 0.000000, 0.000000, 0.000000),
#'            TAU = c(0.000225, 0.232900, 0.329675, 0.418150, 0.019050, 0.000000, 0.000000, 0.000000),
#'            WLC = c(0.996425, 0.003575, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000))
#'
#' sucra(rp, lower.is.better = TRUE)
#' plot(sucra(rp, lower.is.better = TRUE))


sucra = function(x, lower.is.better = FALSE) {

    rank.probability = x


    # Convert rank.probability to matrix
    mat = as.matrix(rank.probability)

    # Loop over treatments, for each treatment: calculate SUCRA
    a = ncol(mat)
    j = nrow(mat)
    names = rownames(mat)

    sucra = numeric()
    for (x in 1:j) {
        sucra[x] = sum(cumsum(mat[x, 1:(a - 1)]))/(a - 1)
    }

    # If condition for lower.is.better
    if (lower.is.better == TRUE) {
        sucra = numeric()
        for (x in 1:j) {
            sucra[x] = 1 - sum(cumsum(mat[x, 1:(a - 1)]))/(a - 1)
        }
    }

    # Make data.frame
    res = data.frame(Treatment = names, SUCRA = sucra)

    # Order
    res = res[order(-res$SUCRA), ]
    rownames(res) = 1:j

    rownames(res) = res$Treatment
    res$Treatment = NULL

    class(res) = "sucra"

    invisible(res)

    res

}