R/rank.functions.R

In MBNMAtime: Run Time-Course Model-Based Network Meta-Analysis (MBNMA) Models

# Functions for ranking treatment responses
# Author: Hugo Pedder
# Date created: 2019-07-24


#' Set rank as a method
#'
#' @param x An object on which to apply the rank method
#' @param ... Arguments to be passed to methods
#'
#' @return Uses the rank method
#'
#' @export
rank <- function (x, ...) {
  UseMethod("rank", x)
}





#' Calculates ranking probabilities for AUC from a time-course MBNMA
#'
#' @inheritParams predict.mbnma
#' @inheritParams rank.mbnma
#' @param mbnma An S3 object of class `"mbnma"` generated by running
#' a time-course MBNMA model
#' @param treats A character vector of treatment/class names (depending on the value of `level`). If left `NULL``
#'   then rankings will be calculated for all treatments/classes. Note that unlike `rank.mbnma()` this argument
#'   cannot take a numeric vector.
#' @param subdivisions The number of subdivisions over which to integrate (see \code{\link[stats]{integrate}})
#'
#' @inherit rank.mbnma return
#' @inherit rank.mbnma details
#'
rankauc <- function(mbnma, lower_better=FALSE, treats=NULL, level="treatments",
                    int.range=c(0,max(mbnma$network$data.ab$time)),
                    n.iter=mbnma$BUGSoutput$n.sims, subdivisions=100, ...) {

  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(mbnma, "mbnma", add=argcheck)
  checkmate::assertCharacter(treats, null.ok = FALSE)
  checkmate::assertChoice(level, choices = c("treatments", "classes"), add=argcheck)
  checkmate::assertLogical(lower_better, any.missing=FALSE, len=1, add=argcheck)
  checkmate::assertIntegerish(subdivisions, lower=2, add=argcheck)
  checkmate::assertIntegerish(int.range, lower=0, any.missing=FALSE, len=2, sorted=TRUE,
                              add=argcheck)
  #checkmate::assertInt(subdivisions, lower=1, add=argcheck)
  checkmate::reportAssertions(argcheck)

  # Initial predict parameters
  #nsims <- mbnma$BUGSoutput$n.sims
  # timecourse <- init.predict(mbnma)[["timecourse"]]
  # beta.incl <- init.predict(mbnma)[["beta.incl"]]

  # Extract parameter values from MBNMA result
  model.vals <- get.model.vals(mbnma, E0=0,
                               level=level)

  # Create vector of parameters in expanded time-course function
  timecourse <- model.vals[["timecourse"]]
  timecourse <- gsub("i\\.", "", timecourse)
  time.params <- model.vals[["time.params"]]

  # Switch spline in timecourse and generate spline matrix
  if (any(c("ns", "bs", "ls") %in% mbnma$model.arg$fun$name)) {
    timecourse <- gsub("\\[i,m\\,", "[,", timecourse)

    seg <- seq(from=int.range[1], to=int.range[2], length.out=subdivisions)
    spline <- genspline(seg,
                        spline=mbnma$model.arg$fun$name, knots=mbnma$model.arg$fun$knots, degree=mbnma$model.arg$fun$degree)
  }

  # Replace mu with 0
  for (i in seq_along(time.params)) {
    if (grepl("^mu", time.params[i])) {
      timecourse <- gsub(time.params[i], 0, timecourse)
      time.params <- time.params[-i]
    }
  }


  # NEW SECTION
  auc.result <- matrix(ncol=length(treats))
  rank.mat <- matrix(ncol=length(treats))
  pb <- utils::txtProgressBar(0, n.iter, style = 3)
  for (mcmc in 1:n.iter) {
    utils::setTxtProgressBar(pb, mcmc)

    # Initialise variables
    auc <- vector()
    rank <- matrix(nrow=length(treats), ncol=length(treats))

    treatsnum <- which(mbnma$network[[level]] %in% treats)
    for (treat in seq_along(treatsnum)) {
      time.mcmc <- timecourse

      #for (i in seq_along(params)) {
      for (i in seq_along(time.params)) {
        # Replace parameter in time-course with value for given treat and MCMC
        #timecourse <- gsub(names(params)[i], params[[i]][mcmc], timecourse)
        temp <- model.vals[[time.params[i]]]
        if (is.vector(temp)) {temp <- matrix(temp, ncol=1)}
        time.mcmc <- gsub(time.params[i],
                          ifelse(is.matrix(temp) & ncol(temp)>1, temp[mcmc,treatsnum[treat]], temp[mcmc]),
                          time.mcmc)
      }

      if (any(c("ns", "bs", "ls") %in% mbnma$model.arg$fun$name)) {

        # Using trapezoid method for spline function
        y <- eval(parse(text=time.mcmc))
        auc <- append(auc, sum(diff(seg)*zoo::rollmean(y,2)))

      } else {
        temp <- paste("int.fun <- function(time) {",
                      time.mcmc,
                      "}",
                      sep=" ")

        eval(parse(text=temp))

        integral <- stats::integrate(int.fun,
                                     lower=int.range[1], upper=int.range[2],
                                     subdivisions=subdivisions)

        auc <- append(auc, integral$value)
      }

    }

    # Need to name columns with treats
    auc.result <- rbind(auc.result, auc)

    # Ranking
    #rank.mat <- rbind(rank.mat, order(auc, decreasing = decreasing))
  }
  auc.result <- auc.result[-1,]
  #rank.mat <- rank.mat[-1,]

  # Ranking
  rank.mat <- t(apply(auc.result, MARGIN=1, FUN=function(x) {
    order(order(x, decreasing = !lower_better), decreasing=FALSE)
  }))

  colnames(auc.result) <- treats
  colnames(rank.mat) <- treats

  summary.rank <- sumrank(rank.mat)

  return(list("summary"=summary.rank,
              "prob.matrix"=calcprob(rank.mat, treats=treats),
              "rank.matrix"=rank.mat,
              "auc.int"=auc.result
  ))
}




sumrank <- function(rank.mat) {
  if (is.null(colnames(rank.mat))) {
    colnames(rank.mat) <- c(1:ncol(rank.mat))
  }

  quantiles.rank <- apply(X=rank.mat, MARGIN = 2,
                          function(x) stats::quantile(x, probs=c(0.025, 0.25, 0.5, 0.75, 0.975)))
  summary.rank <- data.frame(
    "treatment"=colnames(rank.mat),
    "mean"= apply(X=rank.mat, MARGIN = 2, function(x) {base::mean(x)}),
    "sd"= apply(X=rank.mat, MARGIN = 2, function(x) {stats::sd(x)})
  )
  summary.rank <- cbind(summary.rank, t(quantiles.rank))
  rownames(summary.rank) <- NULL

  return(summary.rank)
}


#' Calculates a matrix of ranking probabilities from a matrix of treatment/agent/class
#' rankings
#'
#' @noRd
calcprob <- function(rank.mat, treats=NULL) {
  NT <- ncol(rank.mat)
  rank.prob <- vector(length=NT)

  for (c in 1:NT) {
    pos.vec <- vector()
    for (r in 1:NT) {
      pos.vec <- append(pos.vec,
                        length(rank.mat[rank.mat[,c]==r,c])/nrow(rank.mat))
    }
    rank.prob <- cbind(rank.prob, pos.vec)
  }
  rank.prob <- rank.prob[,-1]

  if (!is.null(treats)) {
    colnames(rank.prob) <- treats
  }

  return("rank.prob"=rank.prob)
}





calcauc <- function(df) {

  id <- order(df$Var1)
  auc <- sum(diff(df$Var1[id])*zoo::rollmean(df$value[id],2))

  return(auc)
}