R/nodesplit-class.R

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

##############################################
#### Functions for class("nodesplit") ####
##############################################


#' @describeIn mb.nodesplit Plot outputs from nodesplit models
#'
#' @param x An object of `class("nodesplit")`
#' @param plot.type A character string that can take the value of `"forest"` to plot
#' only forest plots, `"density"` to plot only density plots, or left as `NULL` (the
#' default) to plot both types of plot.
#' @param params A character vector corresponding to a time-course parameter(s) for which to plot results.
#' If left as `NULL` (the default), nodes-split results for all time-course parameters will be plotted.
#' @param ... Arguments to be sent to `ggplot2::ggplot()` or `ggdist::stat_halfeye()`
#'
#' @details The S3 method `plot()` on an `mb.nodesplit` object generates either
#' forest plots of posterior medians and 95\\% credible intervals, or density plots
#' of posterior densities for direct and indirect evidence.
#'
#' @return Plots the desired graph(s) and returns an object (or list of objects if
#' `plot.type=NULL`) of `class(c("gg", "ggplot"))`, which can be edited using `ggplot` commands.
#'
#' @export
plot.nodesplit <- function(x, plot.type=NULL, params=NULL, ...) {

  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(x, "nodesplit", add=argcheck)
  checkmate::assertChoice(plot.type, choices = c("density", "forest"), null.ok=TRUE, add=argcheck)
  checkmate::assertCharacter(params, null.ok=TRUE, add=argcheck)
  checkmate::reportAssertions(argcheck)

  # Declare global variables
  Parameter <- NULL
  Comparison <- NULL
  value <- NULL
  Source <- NULL

  # Set colours
  cols <- RColorBrewer::brewer.pal(3, "Set1")

  plotlist <- list()

  # Check if params are within nodesplit
  if (!is.null(params)) {
    for (i in seq_along(params)) {
      if (!(params[i] %in% names(x[[1]]))) {
        stop(paste0(params[i], " is not a time-course parameter in `x` for which node-split results are available"))
      }
    }
  } else if (is.null(params)) {
    params <- names(x[[1]])
  }


  if (is.null(plot.type)) {
    plot.type <- c("density", "forest")
  }

  plot.df <- data.frame(value=NA, Source=NA, Parameter=NA, Comparison=NA)
  for (i in seq_along(x)) {
    for (k in seq_along(params)) {
      temp <- x[[i]][[k]]$mcmc[,c("value", "Source")]
      temp$Parameter <- rep(params[k], nrow(temp))
      temp$Comparison <- rep(names(x)[i], nrow(temp))
      plot.df <- rbind(plot.df, temp)
    }
  }
  plot.df <- plot.df[-1,]

  # Replace "pred"
  plot.df$Parameter <- gsub("time", "time = ", plot.df$Parameter)

  # Factorise
  plot.df <- plot.df %>% dplyr::mutate(
    Parameter = factor(Parameter),
    Comparison = factor(Comparison)
  )

  # Forest plots
  if ("forest" %in% plot.type) {
    gg <-
      ggplot2::ggplot(data=plot.df, ggplot2::aes(x=value, y=Source,
                                                 linetype=Source, shape=Source, color=Source)) +
      ggdist::stat_halfeye(.width=0.95) +
      ggplot2::ylab("") + ggplot2::xlab("Treatment effect (95% CrI)") +
      ggplot2::scale_color_manual(values=cols) +
      ggplot2::facet_wrap(Parameter~Comparison, scales = "free_x")

    gg <- gg + theme_mbnma()
    #gg <- do.call(ggdist::stat_halfeye, args=list(...))

    graphics::plot(gg)
    plotlist[[length(plotlist)+1]] <- gg
  }


  # Density plots (with shaded area of overlap)
  if ("density" %in% plot.type) {
    dens.df <- subset(plot.df, Source!="NMA")

    dens <- ggplot2::ggplot(dens.df, ggplot2::aes(x=value, linetype=Source, fill=Source), ...) +
      ggplot2::geom_density(alpha=0.2, linewidth=0.8) +
      ggplot2::xlab("Treatment effect") +
      ggplot2::ylab("Posterior density") +
      ggplot2::scale_fill_manual(values=cols) +
      ggplot2::facet_wrap(Parameter~Comparison, scales = "free_x") +
      theme_mbnma()

    graphics::plot(dens)
    plotlist[[length(plotlist)+1]] <- dens
  }

  return(invisible(plotlist))
}





#' Prints basic results from a node-split to the console
#'
#' @param x An object of class `"nodesplit"` generated by `mb.nodeplit()`
#' @param groupby A character object that can take the value `"time.param"` to present
#' results grouped by time-course parameter (the default) or `"comparison"` to present
#' results grouped by treatment comparison.
#' @param ... arguments to be sent to `knitr::kable()`
#'
#' @return Prints summary details of nodesplit results to the console
#'
#' @export
print.nodesplit <- function(x, groupby="time.param", ...) {

  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(x, "nodesplit", add=argcheck)
  checkmate::assertChoice(groupby, choices=c("time.param", "comparison"), add=argcheck)
  checkmate::reportAssertions(argcheck)

  cat(crayon::bold("========================================\nNode-splitting analysis of inconsistency\n========================================\n"))

  sum.df <- summary.nodesplit(x)

  if (groupby=="time.param") {
    params <- names(x[[1]])

    for (i in seq_along(params)) {

      temp <- params[i]
      temp <- gsub("time", "time = ", temp)

      cat(paste0("\n", crayon::bold(crayon::underline(temp)), "\n\n"))

      param.df <- sum.df[sum.df$Time.Param==params[i],]

      comparisons <- unique(param.df$Comparison)

      out.df <- param.df[1,]
      for (k in seq_along(comparisons)) {
        head <- param.df[param.df$Comparison==comparisons[k],][1,]
        head$Median <- NA
        head$`2.5%` <- NA
        head$`97.5%` <- NA

        tail <- param.df[param.df$Comparison==comparisons[k],]
        tail$Comparison <- c("-> direct", "-> indirect")
        tail$p.value <- rep(NA,2)

        tail <- rbind(tail, rep(NA, ncol(tail)))

        out.df <- rbind(out.df, rbind(head,tail))
      }
      out.df <- out.df[-1, c("Comparison", "p.value", "Median", "2.5%", "97.5%")]

      rownames(out.df) <- NULL

      out <- knitr::kable(out.df, col.names = c("Comparison", "p-value", "Median", "2.5%", "97.5%"), ...)
      cat(gsub('\\bNA\\b', '  ', out), sep='\n')
    }

  } else if (groupby=="comparison") {
    params <- names(x)

    for (i in seq_along(params)) {
      cat(paste0("\n", crayon::bold(crayon::underline(params[i])), "\n\n"))

      param.df <- sum.df[sum.df$Comparison==params[i],]
      timeparams <- unique(param.df$Time.Param)

      out.df <- param.df[1,]
      for (k in seq_along(timeparams)) {
        head <- param.df[param.df$Time.Param==timeparams[k],][1,]
        head$Median <- NA
        head$`2.5%` <- NA
        head$`97.5%` <- NA

        tail <- param.df[param.df$Time.Param==timeparams[k],]
        tail$Time.Param <- c("-> direct", "-> indirect")
        tail$p.value <- rep(NA,2)

        tail <- rbind(tail, rep(NA, ncol(tail)))

        out.df <- rbind(out.df, rbind(head,tail))
      }
      out.df <- out.df[-1, c("Time.Param", "p.value", "Median", "2.5%", "97.5%")]

      rownames(out.df) <- NULL

      out <- knitr::kable(out.df, col.names = c("Time-course parameter", "p-value", "Median", "2.5%", "97.5%"), ...)
      cat(gsub('\\bNA\\b', '  ', out), sep='\n')

    }
  }
}




#' Takes node-split results and produces summary data frame
#'
#' @param object An object of class `"nodesplit"` generated by `mb.nodeplit()`
#' @param ... further arguments passed to or from other methods
#'
#' @return
#' A data frame of summary node-split results with the following variables:
#' * `Comparison` The treatment comparison on which a node-split has been performed
#' * `Time.Param` The time-course parameter on which a node-split has been performed
#' * `Evidence` The evidence contribution for the given comparison (either "Direct" or "Indirect")
#' * `Median` The posterior median
#' * `2.5%` The lower 95% credible interval limit
#' * `97.5%` The upper 95% credible interval limit
#' * `p.value` The Bayesian p-value for the overlap between direct and indirect evidence for
#' the given comparison (it will therefore have an identical value for direct and indirect evidence
#' within a particular comparison and time-course parameter)
#'
#' @export
summary.nodesplit <- function(object, ...) {
  checkmate::assertClass(object, "nodesplit")

  sum.mat <- matrix(ncol=3)
  comp <- vector()
  time.param <- vector()
  evidence <- vector()
  pvals <- vector()

  for (i in seq_along(object)) {
    for (k in seq_along(object[[i]])) {
      post <- object[[i]][[k]]$quantiles

      sum.mat <- rbind(sum.mat, post$direct)
      evidence <- c(evidence, "Direct")

      sum.mat <- rbind(sum.mat, post$indirect)
      evidence <- c(evidence, "Indirect")

      pvals <- c(pvals, rep(object[[i]][[k]]$p.values, 2))
      time.param <- c(time.param, rep(names(object[[i]])[k], 2))
      comp <- c(comp, rep(names(object)[i], 2))
    }
  }
  sum.mat <- round(sum.mat[-1,], digits = max(3L, getOption("digits") - 5L))
  pvals <- round(pvals, max(3L, getOption("digits") - 5L))

  sum.df <- data.frame(comp, time.param,
                       evidence, sum.mat[,2],
                       sum.mat[,1], sum.mat[,3],
                       pvals, ...
  )

  names(sum.df) <- c("Comparison", "Time.Param", "Evidence", "Median",
                     "2.5%", "97.5%", "p.value")

  return(sum.df)

}