R/plot.functions.R

In MBNMAdose: Dose-Response MBNMA Models

# Functions for plotting in MBNMAdose
# Author: Hugo Pedder
# Date created: 2019-04-18

## 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"))




#' Calculate position of label with respect to vertex location within a circle
#'
#' Useful for graphs drawn using `igraph` to reposition labels relative to vertices when vertices
#' are laid out in a circle (as is common in network plots). `igraph` interprets position within
#' `vertex.label.degree` as radians, so it is necessary to convert locations into radian values. This
#' is the main role of this function.
#'
#' @param x A numeric vector of positions around a circle, typically sequentially numbered.
#' @param start A number giving the offset from 12 o'clock in radians for the label locations.
#' @param direction Either `1` for clockwise numbering (based on the order of `x`) or `-1` for
#' anti-clockwise.
#'
#' @examples
#' radian.rescale(c(1:10), start=0, direction=1)
#'
#' @noRd
#'
#' @references
#' https://gist.github.com/kjhealy/834774/a4e677401fd6e4c319135dabeaf9894393f9392c
radian.rescale <- function(x, start=0, direction=1) {
  c.rotate <- function(x) (x + start) %% (2 * pi) * direction
  c.rotate(scales::rescale(x, c(0, 2 * pi), range(x)))
}




#' Get large vector of distinct colours using Rcolorbrewer
#' @noRd
genmaxcols <- function() {

  cols1 <- RColorBrewer::brewer.pal(9, "Set1")
  cols2 <- RColorBrewer::brewer.pal(9, "Pastel1")
  cols <- c(rbind(cols1, cols2))

  cols1 <- RColorBrewer::brewer.pal(8, "Set2")
  cols2 <- RColorBrewer::brewer.pal(8, "Pastel2")
  cols <- c(cols, c(rbind(cols1, cols2)))

  cols <- c(cols, RColorBrewer::brewer.pal(12, "Set3"))

  return(cols)
}





#' Overlays observations as shaded regions on a time-course
#'
#' @inheritParams mbnma.run
#' @inheritParams plot.mbnma.predict
#' @param g An object of `class("ggplot")`
#' @param max.col.scale The maximum rgb numeric value to use for the colour scale
#' @noRd
disp.obs <- function(g, network, predict, col="red", max.col.scale=NULL) {
  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(g, c("gg", "ggplot"), add=argcheck)
  checkmate::assertClass(network, "mbnma.network", add=argcheck)
  checkmate::assertClass(predict, "mbnma.predict", add=argcheck)
  checkmate::assertInt(max.col.scale, lower=1, null.ok = TRUE, add=argcheck)
  checkmate::reportAssertions(argcheck)


  raw.data <- network[["data.ab"]]
  raw.data$agent <- factor(raw.data$agent, labels=network[["agents"]])
  predict.data <- summary(predict)
  predict.data[["count"]] <- NA
  predict.data[["cum.count"]] <- NA

  # Identify counts of doses in raw data
  for (i in 1:nrow(predict.data)) {
    if (i==1) {
      dose.low <- -0.01
    } else if (predict.data$agent[i-1]!=predict.data$agent[i]) {
      dose.low <- -0.01
    } else {
      dose.low <- predict.data$dose[i-1]
    }

    dose.high <- predict.data$dose[i]
    p.agent <- predict.data$agent[i]

    count <- length(raw.data$studyID[as.character(raw.data$agent)==p.agent &
                                     raw.data$dose>dose.low & raw.data$dose<=dose.high])
    cum.count <- length(raw.data$studyID[as.character(raw.data$agent)==p.agent &
                                       raw.data$dose<=dose.high])

    predict.data$count[i] <- count
    predict.data$cum.count[i] <- cum.count
  }

  # Identify if placebo is in data
  if (any(unique(predict.data$agent) %in% c("1", "Placebo"))) {
    plac.count <- predict.data$count[predict.data$agent %in% c("1", "Placebo")]
    message(paste0(plac.count, " placebo arms in the dataset are not shown within the plots"))
    agents <- unique(predict.data$agent)[!(unique(predict.data$agent) %in% c("1", "Placebo"))]
  } else {
    agents <- unique(predict.data$agent)
  }


  # Check max.col.scale
  n.cut <- max(predict.data$count[predict.data$dose!=0])
  if (!is.null(max.col.scale)) {
    if (!is.numeric(max.col.scale)) {
      stop("max.col.scale must be a number greater than the maximum number of observed counts in the plotted treatments")
    }
    if (n.cut > max.col.scale) {
      stop("max.col.scale must be a number greater than the maximum number of observed counts in the plotted treatments")
    }
    n.cut <- max.col.scale
  }

  # Generate colours
  cols <- alpha.scale(n.cut=n.cut, col=col)

  for (agent in seq_along(agents)) {
    subset <- predict.data[predict.data$agent==agents[agent],]
    subset$agent <- factor(subset$agent, labels=levels(g$data$agent)[agents[agent]])

    # Start assuming lowest time = 0
    for (m in 2:nrow(subset)) {
      g <- g + ggplot2::geom_ribbon(data=subset[subset$dose<=subset$dose[m] &
                                                  subset$dose>=subset$dose[m-1]
                                                ,],
                                    ggplot2::aes(x=dose,
                                                 ymin=`2.5%`,
                                                 ymax=`97.5%`),
                                    fill=cols[subset$count[m]+1])
    }

  }

  return(list(g, cols))
}




#' Generates colours with varying degrees of transparency
#'
#' Identical to function in `MBNMAtime` package
#'
#' @param ncut A number indicating the number of different counts in the dataset
#' @param col Colour to use for shading
#' @noRd
alpha.scale <- function(n.cut, col="blue") {
  # Run checks
  checkmate::assertIntegerish(n.cut, lower=0, len=1)

  # Set colour intensities
  if (is.character(col)) {
    if (col=="blue") {
      hue <- c(0,0,200)
    } else if (col=="red") {
      hue <- c(120,0,0)
    } else if (col=="green") {
      hue <- c(0,100,0)
    } else {
      stop("Permitted colours are `blue`, `red`, `green`, or an RGB code")
    }
  } else if (is.numeric(col)) {
    if (length(col)!=3) {
      stop("Specified RGB code must have length 3")
    }
    if (any(col>255) | any(col<0)) {
      stop("RGB values must be between 0 and 255")
    }
    hue <- col
  }

  cut <- (255-0)/(n.cut)
  alpha <- 0
  alpha.vec <- alpha

  for (i in 1:n.cut) {
    alpha <- alpha+cut
    alpha.vec <- append(alpha.vec, alpha)
  }

  hexcol <- vector()
  for (i in 1:(n.cut+1)) {
    hexcol <- append(hexcol, grDevices::rgb(hue[1], hue[2],
                                 hue[3], alpha.vec[i], maxColorValue=255
    ))
  }

  return(hexcol)
}




#' Overlays results of split (treatment-level) NMA on a plot
#'
#' Requires automatically running an NMA
#'
#' @inheritParams plot.mbnma.predict
#' @inheritParams predict.mbnma
#' @inheritParams disp.obs
#' @noRd
overlay.split <- function(g, network, E0=unique(g$data$`50%`[g$data$dose==0]), method="common",
                          likelihood="binomial", link="logit", lim="cred", ...) {

  # Check/assign link and likelihood
  likelink <- check.likelink(network$data.ab, likelihood=likelihood, link=link)
  likelihood <- likelink[["likelihood"]]
  link <- likelink[["link"]]

  splitNMA <- nma.run(network=network, method=method,
                      likelihood=likelihood, link=link, ...)

  splitmat <- splitNMA[["jagsresult"]]$BUGSoutput$sims.list$d

  if (splitNMA[["trt.labs"]][1]!="Placebo_0") {
    warning("Split NMA network does not include placebo:\nresults cannot be displayed over MBNMA predictions")
    return(g)
  }

  trtvec <- splitNMA[["trt.labs"]]
  agentvec <- as.vector(sapply(splitNMA[["trt.labs"]],
                           function(x) strsplit(x, split="_", fixed=TRUE)[[1]][1]))
  dosevec <- as.vector(as.numeric(sapply(splitNMA[["trt.labs"]],
                                     function(x) strsplit(x, split="_", fixed=TRUE)[[1]][2])))
  agentvec <- factor(agentvec, levels=network$agents)

  # Ensure only agents whose predictions are plotted are included
  splitmat <- splitmat[, agentvec %in% g$data$agent]
  dosevec <- dosevec[agentvec %in% g$data$agent]
  trtvec <- trtvec[agentvec %in% g$data$agent]
  agentvec <- agentvec[agentvec %in% g$data$agent]

  # Add between-study SD
  if (lim=="pred" & method=="random") {
    mat <- array(dim = c(nrow(splitmat), ncol(splitmat), 2))
    mat[,,1] <- splitmat
    mat[,,2] <- stats::median(splitNMA$jagsresult$BUGSoutput$sims.list$sd)
    splitmat <- apply(mat, MARGIN=c(1,2), FUN=function(x) stats::rnorm(1, x[1], x[2]))
  }

  # Ensure E0 is correct length to add to splitmat
  if (length(E0)<nrow(splitmat)) {
    E0 <- rep(E0, length.out=nrow(splitmat))
  } else if (length(E0)>nrow(splitmat)) {
    E0 <- replicate(nrow(splitmat), sample(E0))
  }

  splitmat <- splitmat + E0

  quants <- apply(splitmat, MARGIN=2, FUN=function(x) {
    stats::quantile(x, probs=c(0.025,0.5,0.975))
  })

  # Return to natural scale
  quants <- rescale.link(quants, direction="natural", link=link)

  split.df <- data.frame("agent"=agentvec, "dose"=dosevec, "treatment"=trtvec,
                         "2.5%"=quants[1,], "50%"=quants[2,], "97.5%"=quants[3,])

  names(split.df)[4:6] <- c("2.5%", "50%", "97.5%")


  # Add split NMAs to plot
  g <- g + ggplot2::geom_point(data=split.df, ggplot2::aes(x=dose, y=`50%`)) +
    ggplot2::geom_errorbar(data=split.df, ggplot2::aes(x=dose, ymin=`2.5%`, ymax=`97.5%`), width=.05)


  #### Report split NMA model fit statistics ####
  nma.msg <- vector()
  nma.msg <- c(nma.msg, paste0("Split NMA residual Deviance: ",
                               round(splitNMA$jagsresult$BUGSoutput$median$totresdev, 1)))

  if (method=="random") {
    sd <- splitNMA$jagsresult$BUGSoutput$summary[rownames(splitNMA$jagsresult$BUGSoutput$summary)=="sd",]
    sd <- paste0(signif(sd[5], 3), " (", signif(sd[3], 3), ", ", signif(sd[7], 3), ")")

    nma.msg <- c(nma.msg, paste0("Split NMA between-study SD: ", sd
                                 ))
  }

  message(cat(paste(nma.msg, collapse="\n")))

  return(g)

}







#' Plot deviance contributions from an MBNMA model
#'
#' @param mbnma An S3 object of class `"mbnma"` generated by running
#'   a dose-response MBNMA model
#' @param dev.type *STILL IN DEVELOPMENT FOR MBNMAdose!* Deviances to plot - can be either residual
#' deviances (`"resdev"`, the default) or deviances (`"dev"`)
#' @param plot.type Deviances can be plotted either as scatter points (`"scatter"`)
#' or as boxplots (`"box"`)
#' @param facet A boolean object that indicates whether or not to facet (by agent for `MBNMAdose`
#' and by treatment for `MBNMAtime`)
#' @param ... Arguments to be sent to `ggplot2::geom_point()` or `ggplot2::geom_boxplot`
#' @inheritParams predict.mbnma
#' @inheritParams R2jags::jags
#'
#' @return Generates a plot of deviance contributions and returns a list containing the
#' plot (as an object of `class(c("gg", "ggplot"))`), and a data.frame of posterior mean
#' deviance/residual deviance contributions for each observation.
#'
#' @details
#' Deviances should only be plotted for models that have converged successfully. If deviance
#' contributions have not been monitored in `mbnma$parameters.to.save` then additional
#' iterations will have to be run to get results for these.
#'
#' For `MBNMAtime`, deviance contributions cannot be calculated for models with a multivariate likelihood (i.e.
#' those that account for correlation between observations) because the covariance matrix in these
#' models is treated as unknown (if `rho = "estimate"`) and deviance contributions will be correlated.
#'
#' @examples
#' \donttest{
#' # Using the triptans data
#' network <- mbnma.network(triptans)
#'
#' # Run an Emax dose-response MBNMA and predict responses
#' emax <- mbnma.run(network, fun=demax(), method="random")
#'
#' # Plot deviances
#' devplot(emax)
#'
#' # Plot deviances using boxplots
#' devplot(emax, plot.type="box")
#'
#' # Plot deviances on a single scatter plot (not facetted by agent)
#' devplot(emax, facet=FALSE, plot.type="scatter")
#'
#' # A data frame of deviance contributions can be obtained from the object
#' #returned by `devplot`
#' devs <- devplot(emax)
#' head(devs$dev.data)
#'
#' # Other deviance contributions not currently implemented but in future
#' #it will be possible to plot them like so
#' #devplot(emax, dev.type="dev")
#' }
#'
#' @export
devplot <- function(mbnma, plot.type="box", facet=TRUE, dev.type="resdev",
                    n.iter=mbnma$BUGSoutput$n.iter/2, n.thin=mbnma$BUGSoutput$n.thin,
                    ...) {
  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(mbnma, "mbnma", add=argcheck)
  checkmate::assertChoice(dev.type, choices = c("dev", "resdev"), add=argcheck)
  checkmate::assertChoice(plot.type, choices = c("scatter", "box"), add=argcheck)
  checkmate::assertInt(n.iter, lower=1, null.ok = TRUE, add=argcheck)
  checkmate::assertInt(n.thin, lower=1, upper=n.iter, null.ok = TRUE, add=argcheck)
  checkmate::reportAssertions(argcheck)

  if (!is.null(mbnma$model.arg$rho)) {
    stop("Correlation between time points has been modelled using a multivariate likelihood.
         Deviances cannot be calculated for this model.")
  }

  if (!(dev.type %in% mbnma$parameters.to.save)) {
    msg <- paste0("`", dev.type, "` not monitored in mbnma$parameters.to.save.\n",
                  "additional iterations will be run in order to obtain results for `", dev.type, "`")
    message(msg)

    dev.df <- mbnma.update(mbnma, param=dev.type, n.iter=n.iter, n.thin=n.thin)

  } else {

    dev.df <- get.theta.dev(mbnma, param=dev.type)
  }

  # Add studyID in addition to study index from model and sort data frame
  studvec <- vector()
  for (i in seq_along(mbnma$model.arg$jagsdata$studyID)) {
    studvec <- append(studvec, rep(mbnma$model.arg$jagsdata$studyID[i], mbnma$model.arg$jagsdata$narm[i]))
  }

  dev.df$studyID <- studvec
  # dev.df$studyID <- mbnma$model.arg$jagsdata$studyID
  dev.df <- dplyr::arrange(dev.df, dev.df$study, dev.df$arm)

  # Plots the residual deviances
  agents <- mbnma$network$agents

  # Remove placebo results if they are present
  if (mbnma$network$agents[1]=="Placebo") {
    dev.df <- dev.df[dev.df$facet!=1,]
    agents <- agents[-1]
  }

  xlab <- "Dose"
  facetscale <- "free_x"

  if ("agents" %in% names(mbnma$network)) {
    dev.df$facet <- factor(dev.df$facet, labels=agents)
  } else if ("treatments" %in% names(mbnma$network)) {
    dev.df$facet <- factor(dev.df$facet, labels=mbnma$treatments)
  }

  if (plot.type=="scatter") {
    g <- ggplot2::ggplot(dev.df, ggplot2::aes(x=fupdose, y=mean), group=groupvar) +
      ggplot2::geom_point(...)
  } else if (plot.type=="box") {
    g <- ggplot2::ggplot(dev.df, ggplot2::aes(x=factor(fupdose), y=mean)) +
      ggplot2::geom_boxplot(...)
  }

  # Add axis labels
  g <- g + ggplot2::xlab(xlab) +
    ggplot2::ylab("Posterior mean") +
    theme_mbnma()

  if (facet==TRUE) {
    g <- g + ggplot2::facet_wrap(~facet, scales = facetscale) +
      ggplot2::expand_limits(x=0)
  }

  graphics::plot(g)
  return(invisible(list("graph"=g, "dev.data"=dev.df)))
}





#' Extracts fitted values or deviance contributions into a data.frame with indices
#'
#' @inheritParams predict.mbnma
#' @inheritParams mbnma.update
#' @param The parameter for which to extract values - can take either `"theta"`, `"dev"` or `"resdev"`
#' @noRd
get.theta.dev <- function(mbnma, param="theta", armdat=TRUE) {
  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(mbnma, "rjags", add=argcheck)
  checkmate::assertChoice(param, choices = c("dev", "resdev", "theta"), add=argcheck)
  checkmate::assertLogical(armdat, add=argcheck)
  checkmate::reportAssertions(argcheck)

  dev.df <- as.data.frame(
    mbnma$BUGSoutput$summary[grep(paste0("^", param, "\\["),
                                  rownames(mbnma$BUGSoutput$summary)),]
  )

  if (mbnma$type=="time") {
    # Takes the study, arm and follow-up measurement identifier for each residual deviance point
    id <- matrix(unlist(strsplit(
      gsub(
        "([a-z]+\\[)([0-9]+)(,)([0-9]+)(,)([0-9]+)(\\])",
        "\\2.\\4.\\6", rownames(dev.df)),
      split="\\.")),
      byrow=TRUE,
      ncol=3
    )

    dev.df$fupdose <- as.numeric(id[,3])
    dev.df$groupvar <- paste(as.numeric(id[,1]), as.numeric(id[,2]), sep="_")

    # Treatment as facet
    temp <- replicate(max(dev.df$fupdose), mbnma$model$data()$treatment)
    dev.df$facet <- as.vector(temp)[
      stats::complete.cases(as.vector(temp))
      ]

  } else if (mbnma$type=="dose") {
    # Takes the study, arm and follow-up measurement identifier for each residual deviance point
    id <- matrix(unlist(strsplit(
      gsub(
        "([a-z]+\\[)([0-9]+)(,)([0-9]+)(\\])",
        "\\2.\\4", rownames(dev.df)),
      split="\\.")),
      byrow=TRUE,
      ncol=2
    )

    if (armdat==TRUE) {
      # Agent as facet
      dev.df$facet <- as.vector(mbnma$model.arg$jagsdata$agent)[
        stats::complete.cases(as.vector(mbnma$model.arg$jagsdata$agent))
      ]

      dev.df$fupdose <- as.vector(mbnma$model.arg$jagsdata$dose)[
        stats::complete.cases(as.vector(mbnma$model.arg$jagsdata$dose))
      ]
      dev.df$groupvar <- as.numeric(id[,1])
    }

  }

  dev.df$study <- as.numeric(id[,1])
  dev.df$arm <- as.numeric(id[,2])


  return(dev.df)
}





#' Plot fitted values from MBNMA model
#'
#' @param mbnma An S3 object of class `"mbnma"` generated by running
#'   a dose-response MBNMA model
#' @param disp.obs A boolean object to indicate whether raw data responses should be
#' plotted as points on the graph
#' @param ... Arguments to be sent to `ggplot2::geom_point()` or `ggplot2::geom_line()`
#' @inheritParams predict.mbnma
#' @inheritParams R2jags::jags
#'
#' @return Generates a plot of fitted values from the MBNMA model and returns a list containing
#' the plot (as an object of `class(c("gg", "ggplot"))`), and a data.frame of posterior mean
#' fitted values for each observation.
#'
#' @details
#' Fitted values should only be plotted for models that have converged successfully.
#' If fitted values (`theta`) have not been monitored in `mbnma$parameters.to.save`
#' then additional iterations will have to be run to get results for these.
#'
#' @examples
#' \donttest{
#' # Using the triptans data
#' network <- mbnma.network(triptans)
#'
#' # Run an Emax dose-response MBNMA and predict responses
#' emax <- mbnma.run(network, fun=demax(), method="random")
#'
#' # Plot fitted values and observed values
#' fitplot(emax)
#'
#' # Plot fitted values only
#' fitplot(emax, disp.obs=FALSE)
#'
#' # A data frame of fitted values can be obtained from the object
#' #returned by `fitplot`
#' fits <- fitplot(emax)
#' head(fits$fv)
#' }
#'
#' @export
fitplot <- function(mbnma, disp.obs=TRUE,
                    n.iter=mbnma$BUGSoutput$n.iter, n.thin=mbnma$BUGSoutput$n.thin, ...) {
  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(mbnma, "mbnma", add=argcheck)
  checkmate::assertLogical(disp.obs, add=argcheck)
  checkmate::assertInt(n.iter, lower=1, null.ok = TRUE, add=argcheck)
  checkmate::assertInt(n.thin, lower=1, upper=n.iter, null.ok = TRUE, add=argcheck)
  checkmate::reportAssertions(argcheck)

  if (!("theta" %in% mbnma$parameters.to.save)) {
    msg <- paste0("`theta` not monitored in mbnma$parameters.to.save.\n",
                  "additional iterations will be run in order to obtain results")
    message(msg)

    theta.df <- mbnma.update(mbnma, param="theta", n.iter=n.iter, n.thin=n.thin)
  } else {
    theta.df <- get.theta.dev(mbnma, param="theta")
  }

  # Add studyID in addition to study index from model
  studvec <- vector()
  for (i in seq_along(mbnma$model.arg$jagsdata$studyID)) {
    studvec <- append(studvec, rep(mbnma$model.arg$jagsdata$studyID[i], mbnma$model.arg$jagsdata$narm[i]))
  }

  theta.df$studyID <- studvec
  # theta.df$studyID <- mbnma$model.arg$jagsdata$studyID


  # Obtain raw responses to plot over fitted
  if (mbnma$model.arg$likelihood=="binomial") {
    theta <- mbnma$model$data()$r / mbnma$model$data()$n
  } else if (mbnma$model.arg$likelihood=="poisson") {
    theta <- mbnma$model$data()$r / mbnma$model$data()$E
  } else if (mbnma$model.arg$likelihood=="normal") {
    theta <- mbnma$model$data()$y
  }

  theta <- rescale.link(theta, direction="link", link=mbnma$model.arg$link)
  theta.df$y <- as.vector(theta)[stats::complete.cases(as.vector(theta))]

  # Add E0 for each agent and drop dose==0
  E0 <- mean(theta.df$mean[theta.df$fupdose==0])
  theta.df <- theta.df[theta.df$fupdose!=0,]

  for (i in seq_along(unique(theta.df$study))) {
    subset <- theta.df[theta.df$study==unique(theta.df$study)[i],]
    for (k in seq_along(unique(subset$facet))) {
      row <- subset[subset$facet==unique(subset$facet)[k],][1,]
      row$mean <- E0
      row$y <- NA
      row$fupdose <- 0
      theta.df <- rbind(theta.df, row)
    }
  }

  theta.df <- dplyr::arrange(theta.df, theta.df$study, theta.df$facet, theta.df$fupdose)


  # Axis labels
  xlab <- "Dose"
  facetscale <- "free_x"
  ylab <- "Response on link scale"

  # Add facet labels
  if ("agents" %in% names(mbnma$network)) {
    if (mbnma$network$agents[1]=="Placebo" & mbnma$network$treatments[1]=="Placebo_0") {
      labs <- mbnma$network$agents[-1]
    } else {labs <- mbnma$network$agents}
    theta.df$facet <- factor(theta.df$facet, labels=labs)
  } else if ("treatments" %in% names(mbnma$network)) {
    labs <- mbnma$network$treatments[-1]
    theta.df$facet <- factor(theta.df$facet, labels=labs)
  }

  # Generate plot
  g <- ggplot2::ggplot(theta.df,
                       ggplot2::aes(x=fupdose, y=mean, group=groupvar)) +
    ggplot2::geom_line(...)

  # Overlay observed responses
  if (disp.obs==TRUE) {
    g <- g + ggplot2::geom_point(ggplot2::aes(y=y), size=1, ...)
  }

  # Add facets
  g <- g + ggplot2::facet_wrap(~facet, scales = facetscale)

  # Add axis labels
  g <- g + ggplot2::xlab(xlab) +
    ggplot2::ylab(ylab) +
    theme_mbnma()

  suppressWarnings(graphics::plot(g))

  return(invisible(list("graph"=g, "fv"=theta.df)))
}







plot.invisible <- function(...){
  ff <- tempfile()
  grDevices::png(filename=ff)
  res <- graphics::plot(...)
  grDevices::dev.off()
  unlink(ff)
  return(res)
}






#' Plot cumulative ranking curves from MBNMA models
#'
#' @param x An object of class `"mbnma.rank"` generated by `rank.mbnma()`
#' @param sucra A logical object to indicate whether Surface Under Cumulative Ranking Curve (SUCRA)
#' values should be calculated and returned as a data frame. Areas calculated
#' using trapezoid approach.
#' @param ... Arguments to be sent to `ggplot::geom_line()`
#' @inheritParams rank.mbnma
#'
#' @return Line plots showing the cumulative ranking probabilities for each agent/class and
#' dose-response parameter in `x`. The object returned is a list which contains the plot
#' (an object of `class(c("gg", "ggplot")`) and a data frame of SUCRA values
#' if `sucra = TRUE`.
#'
#' @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 cumulative rankings for both dose-response parameters simultaneously
#' # Note that SUCRA values are also returned
#' cumrank(ranks)
#' }
#' @export
cumrank <- function(x, params=NULL, sucra=TRUE, ...) {

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

  output <- list()

  if (is.null(params)) {
    params <- names(x)
  }

  df <- data.frame()
  for (param in seq_along(params)) {
    if (!params[param] %in% names(x)) {
      stop(paste0(params[param], " is not a ranked parameter in x"))
    }

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

    melt <- reshape2::melt(cum.mat)
    melt$param <- params[param]

    df <- rbind(df, melt)

  }

  df$Parameter <- factor(df$param)

  g <- ggplot2::ggplot(df, ggplot2::aes(x=Var1, y=value, linetype=Parameter, colour=Parameter), ...) +
    ggplot2::geom_line(size=1)

  g <- g + ggplot2::facet_wrap(~factor(Var2)) +
    ggplot2::xlab("Rank (1 = best)") +
    ggplot2::ylab("Cumulative probability") +
    ggplot2::labs(linetype="Parameter", colour="Parameter") +
    theme_mbnma()

  graphics::plot(g)

  output <- list("cumplot"=g)

  # Calculate AUC
  if (sucra==TRUE) {
    df.auc <- df %>%
      dplyr::group_by(df$Var2, df$param) %>%
      dplyr::do(data.frame(sucra=calcauc(.)))

    df.auc <- dplyr::ungroup(df.auc)

    names(df.auc) <- c("agent", "parameter", "sucra")

    output[["sucra"]] <- df.auc

    print(df.auc)
  }

  return(invisible(output))
}




#' Returns a forest plot of nodesplit results
#'
#' @param ... Optional arguments to be passed to `forestplot::forestplot()`
#' @inheritParams plot.nodesplit
#' @noRd
forest.splits <- function(x, ...) {

  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(x, "nodesplit", add=argcheck)
  checkmate::reportAssertions(argcheck)


  df <- summary.nodesplit(x)

  df$Evidence[df$Evidence %in% c("MBNMA", "NMA")] <- "Overall"

  comps <- lapply(x, FUN=function(y) {paste(y$comparison, collapse=" vs \n")})

  match <- lapply(x, FUN=function(y) {paste(y$comparison, collapse=" vs ")})

  df$Comparison <- unlist(comps)[match(df$Comparison, match)]

  df <- dplyr::arrange(df, Comparison, Evidence)

  names(df)[4:5] <- c("l95", "u95")


  df$p.value <- format(df$p.value)
  df$p.value[df$p.value=="0.000"] <- "<0.001"

  # Remove surplus labels
  df$Comparison[df$Evidence!="Direct"] <- NA
  df$p.value[df$Evidence!="Indirect"] <- NA


  # Add blank between groups
  temp <- df[0,]
  for (i in 1:(nrow(df)/3)) {
    seg <- rbind(df[((i-1)*3+1):(i*3),], rep(NA,ncol(df)), rep(NA,ncol(df)))
    temp <- rbind(temp, seg)
  }
  df <- temp

  # Plot forest plot
  forestplot::forestplot(labeltext=cbind(df$Comparison, df$Evidence, df$p.value),
                         mean=df$Median, lower=df$l95, upper=df$u95,
                         boxsize=0.2, graph.pos=3,
                         xlab="Effect size (95% CrI)", hrzl_lines = TRUE,
                         ...)

}





#' Dev-dev plot for comparing deviance contributions from two models
#'
#' Plots the deviances of two model types for comparison. Often used to assess
#' consistency by comparing consistency (NMA or MBNMA) and unrelated mean effects (UME)
#' models (see \insertCite{pedder2021cons;textual}{MBNMAdose}). Models must be run
#' on the *same set of data* or the deviance comparisons will not be valid.
#'
#' @param mod1 First model for which to plot deviance contributions
#' @param mod2 Second model for which to plot deviance contributions
#' @inheritParams devplot
#'
#' @examples
#' \donttest{
#' # Using the triptans data
#' network <- mbnma.network(triptans)
#'
#' # Run an poorly fitting linear dose-response
#' lin <- mbnma.run(network, fun=dpoly(degree=1))
#'
#' # Run a better fitting Emax dose-response
#' emax <- mbnma.run(network, fun=demax())
#'
#' # Run a standard NMA with unrelated mean effects (UME)
#' ume <- nma.run(network, UME=TRUE)
#'
#' # Compare residual deviance contributions from linear and Emax
#' devdev(lin, emax) # Suggests model fit is very different
#'
#' # Compare deviance contributions from Emax and UME
#' devdev(emax, ume) # Suggests model fit is similar
#'
#' }
#'
#' @export
devdev <- function(mod1, mod2, dev.type="resdev",
                   n.iter=2000, n.thin=1,
                   ...) {

  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertChoice(dev.type, choices = c("dev", "resdev"), add=argcheck)
  checkmate::reportAssertions(argcheck)

  if (!any(c("nma", "mbnma") %in% class(mod1))) {
    stop("mod1 must be of class 'mbnma' or 'nma'")
  }
  if (!any(c("nma", "mbnma") %in% class(mod2))) {
    stop("mod2 must be of class 'mbnma' or 'nma'")
  }

  modlist <- list(mod1, mod2)
  devlist <- list()

  for (mod in seq_along(modlist)) {
    # Reformat NMA data if necessary
    if ("nma" %in% class(modlist[[mod]])) {
      modlist[[mod]] <- modlist[[mod]]$jagsresult
    }

    if (!(dev.type %in% modlist[[mod]]$parameters.to.save)) {
      msg <- paste0("`", dev.type, "` not monitored in mod", mod, "$parameters.to.save.\n",
                    "additional iterations will be run in order to obtain results for `", dev.type, "`")
      message(msg)

      devlist[[length(devlist)+1]] <- mbnma.update(modlist[[mod]], param=dev.type, n.iter=n.iter, n.thin=n.thin, armdat = FALSE)

    } else {

      devlist[[length(devlist)+1]] <- get.theta.dev(modlist[[mod]], param=dev.type)
    }
  }

  if (nrow(devlist[[1]])!=nrow(devlist[[2]])) {
    stop("Number of data points differ between mod1 and mod2\nsuggests models have not been run on the same dataset")
  }

  dev.df <- cbind(devlist[[1]], devlist[[2]]$mean)

  dev.df <- dev.df %>% dplyr::rename(mod1.mean="mean", mod2.mean="devlist[[2]]$mean") %>%
    dplyr::select(study, arm, mod1.mean, mod2.mean)

  g <- ggplot2::ggplot(dev.df, ggplot2::aes(x=mod1.mean, y=mod2.mean)) +
    ggplot2::geom_point(...)

  # Add axis labels
  g <- g + ggplot2::xlab(paste0("Posterior mean of ", dev.type, " for mod1")) +
    ggplot2::ylab(paste0("Posterior mean of ", dev.type, " for mod2")) +
    theme_mbnma()

  graphics::plot(g)
  return(invisible(list("graph"=g, "dev.data"=dev.df)))

}




#' MBNMA ggplot2 theme style
#' @noRd
theme_mbnma <- function(...) {
  ggplot2::theme_bw(...) +
    ggplot2::theme(
      # change stuff here
      panel.background  = ggplot2::element_blank(),
      plot.background = ggplot2::element_rect(fill="white", colour=NA),
      legend.background = ggplot2::element_rect(fill="transparent", colour=NA),
      legend.key = ggplot2::element_rect(fill="transparent", colour=NA),

      # From multinma
      #panel.border = ggplot2::element_rect(colour = "grey70", fill = NA),
      panel.grid.major = ggplot2::element_line(colour = "grey95"),
      panel.grid.minor = ggplot2::element_line(colour = "grey95"),
      strip.background = ggplot2::element_rect(colour = "black",
                                               fill = "lightsteelblue1"),
      strip.text = ggplot2::element_text(colour = "black")
    )
}