R/mbnma.rank-class.R

In MBNMAdose: Dose-Response MBNMA Models

###########################################
#### Functions for class("mbnma.rank") ####
###########################################

## quiets concerns of R CMD check re: the .'s that appear in pipelines
if(getRversion() >= "2.15.1")  utils::globalVariables(c(".", "studyID", "agent", "dose", "Var1", "value",
                                                        "Parameter", "fupdose", "groupvar", "y",
                                                        "network", "a", "param", "med", "l95", "u95", "value",
                                                        "Estimate", "2.5%", "50%", "97.5%", "treatment"))



#' Plot histograms of rankings from MBNMA models
#'
#' @param x An object of class `"mbnma.rank"` generated by `rank.mbnma()`
#' @param treat.labs A vector of treatment labels in the same order as treatment codes.
#' Easiest to use treatment labels stored by `mbnma.network()`
#' @param ... Arguments to be sent to `ggplot::geom_bar()`
#' @inheritParams rank.mbnma
#'
#' @return A series of histograms that show rankings for each treatment/agent/prediction, with a
#' separate panel for each parameter.
#' The object returned is a list containing a separate element for each parameter in `params`
#' which is an object of `class(c("gg", "ggplot"))`.
#'
#' @examples
#' \donttest{
#' # Using the triptans data
#' network <- mbnma.network(triptans)
#'
#' # Estimate rankings  from an Emax dose-response MBNMA
#' emax <- mbnma.run(network, fun=demax(), method="random")
#' ranks <- rank(emax)
#'
#' # Plot rankings for both dose-response parameters (in two separate plots)
#' plot(ranks)
#'
#' # Plot rankings just for ED50
#' plot(ranks, params="ed50")
#'
#' # Plot rankings from prediction
#' doses <- list("eletriptan"=c(0,1,2,3), "rizatriptan"=c(0.5,1,2))
#' pred <- predict(emax, E0 = "rbeta(n, shape1=1, shape2=5)",
#'             exact.doses=doses)
#' rank <- rank(pred)
#' plot(rank)
#'
#'
#' # Trying to plot a parameter that has not been ranked will return an error
#' #### ERROR ####
#' # plot(ranks, params="not.a.parameter")
#' }
#' @export
plot.mbnma.rank <- function(x, params=NULL, treat.labs=NULL, ...) {
  # ... are commands to be sent to geom_histogram

  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(x, "mbnma.rank", add=argcheck)
  checkmate::assertCharacter(params, null.ok=TRUE, add=argcheck)
  checkmate::assertCharacter(treat.labs, null.ok=TRUE, add=argcheck)
  checkmate::reportAssertions(argcheck)

  output <- list()

  if (is.null(params)) {
    params <- names(x)
  } else {
    # Check params is in x
    if (!all(params %in% names(x))) {
      stop("'params' must be a subset of named list elements in 'x' representing ranked parameters")
    }
  }

  # Plot a separate set of histograms for each parameter in param
  for (param in seq_along(params)) {

    rank.mat <- x[[params[param]]]$rank.matrix
    treats <- colnames(rank.mat)

    ranks <- vector()
    treat <- vector()
    for (i in seq_along(treats)) {
      treat <- append(treat, rep(treats[i], nrow(rank.mat)))
      ranks <- append(ranks, rank.mat[,i])
    }
    data <- data.frame("ranks"=ranks, "treat"=treat)

    if (!is.null(treat.labs)) {
      if (length(treat.labs)!=length(unique(data$treat))) {
        stop("`treat.labs` must be a character vector of the same length as the number of ranked tretments/agents")
      }
      data$treat <- factor(data$treat, labels=treat.labs)
    } else {
      data$treat <- factor(data$treat)
    }

    # Plot histograms
    g <- ggplot2::ggplot(data, ggplot2::aes(x=ranks)) +
      ggplot2::geom_bar(...) +
      ggplot2::xlab("Rank (1 = best)") +
      ggplot2::ylab("MCMC iterations") +
      ggplot2::facet_wrap(~treat) +
      ggplot2::ggtitle(params[param]) +
      theme_mbnma()

    graphics::plot(g)

    output[[params[param]]] <- g
  }

  return(invisible(output))
}




#' Generates summary data frames for an mbnma.rank object
#'
#' @param object An object of `class("mbnma.rank")` generated by `rank.mbnma()`
#' @param ... additional arguments affecting the summary produced
#'
#' @return A list in which each element represents a parameter that has been ranked
#' in `mbnma.rank` and contains a data frame of summary ranking results.
#'
#' @export
summary.mbnma.rank <- function(object, ...) {
  checkmate::assertClass(object, "mbnma.rank")

  output <- list(...)
  for (i in seq_along(object)) {
    output[[names(object)[i]]] <- object[[i]]$summary
  }
  return(output)
}




#' Prints summary information about an mbnma.rank object
#'
#' @param x An object of class `"mbnma.rank"` generated by `rank.mbnma()`
#' @param ... further arguments passed to or from other methods
#'
#' @export
print.mbnma.rank <- function(x, ...) {
  checkmate::assertClass(x, "mbnma.rank")

  attrs <- attributes(x)

  head <- crayon::bold("\n================================\nRanking of dose-response MBNMA\n================================")

  if ("predictions" %in% attrs$level) {
    level.str <- "predictions"
  } else if ("agent" %in% attrs$level) {
    level.str <- "agents"
  } else if ("class" %in% attrs$level) {
    level.str <- "classes"
  } else if ("relefs" %in% attrs$level) {
    level.str <- "relefs"
  }

  intro <- c()

  if ("predictions" %in% attrs$level) {
    intro <- c(intro, paste0("Includes ranking of predictions from dose-response MBNMA"))
  } else if ("relefs" %in% attrs$level) {
    intro <- c(intro, paste0("Includes ranking of relative effects"))
  }

  if (any(c("agent", "class") %in% attrs$level)) {
    relef <- vector()
    classef <- vector()
    for (i in seq_along(names(x))) {
      if (grepl("[A-Z]", names(x)[i])) {
        classef <- append(classef, names(x)[i])
      } else {
        relef <- append(relef, names(x)[i])
      }
    }
    if (length(relef)>1) {
      add <- "Includes ranking of relative treatment effects from dose-response MBNMA:"
      add <- paste(add,
                   crayon::bold(paste(relef, collapse="\t")), "",
                   sep="\n")
      intro <- c(intro, add)
    }
    if (length(classef)>0) {
      add <- "Includes ranking of relative class effects from dose-response MBNMA:"
      add <- paste(add,
                   crayon::bold(paste(classef, collapse="\t")), "",
                   sep="\n")
      intro <- c(intro, add)
    }
  }

  # Regression
  if ("regress.vals" %in% names(attrs)) {
    intro <- c(intro,
               paste0("Rankings generated from a model adjusting for the following effect modififers:\n",
                      paste(crayon::bold(names(attrs$regress.vals)), collapse="\t")))
  }

  rankinfo <- paste(nrow(x[[1]]$summary), level.str, "ranked", sep=" ")
  if (attrs$lower_better==FALSE) {
    rankinfo <- paste(rankinfo, "with positive responses being", crayon::green(crayon::bold("`better`")), sep=" ")
  } else if (attrs$lower_better==TRUE) {
    rankinfo <- paste(rankinfo, "with negative responses being", crayon::red(crayon::bold("`worse`")), sep=" ")
  }

  intro <- paste(intro, collapse="\n")
  out <- paste(head, intro, rankinfo, sep="\n\n")
  cat(paste0(out, "\n\n"))

  # Add summary table
  for (param in seq_along(names(x))) {

    head <- paste0(crayon::bold(crayon::underline(paste0(names(x)[param], " ranking"))), " (from best to worst)")

    sumtab <- x[[names(x)[param]]]$summary
    sumtab <- dplyr::arrange(sumtab, mean)
    sumtab <- sumtab[,c(1,2,6,4,8)]

    cat(head)

    print(knitr::kable(sumtab, col.names = c("Treatment", "Mean", "Median", "2.5%", "97.5%"), digits = 2))
    cat("\n\n")
  }

}