R/mbnma-class.R

In MBNMAdose: Dose-Response MBNMA Models

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

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

#' Forest plot for results from dose-response MBNMA models
#'
#' Generates a forest plot for dose-response parameters.
#'
#' @param x An S3 object of class `"mbnma"` generated by running
#'   a dose-response MBNMA model
#' @param params A character vector of dose-response parameters to plot.
#' Parameters must be given the same name as monitored nodes in `mbnma` and must be
#' modelled as relative effects (`"rel"`). Can be set to
#' `NULL` to include all available dose-response parameters estimated by `mbnma`.
#' @param ... Arguments to be passed to methods, such as graphical parameters
#'
#' @return A forest plot of class `c("gg", "ggplot")` that has separate panels for
#' different dose-response parameters. Results are plotted on the link scale.
#'
#' @examples
#' \donttest{
#' # Using the triptans data
#' network <- mbnma.network(triptans)
#'
#' # Run an exponential dose-response MBNMA and generate the forest plot
#' exponential <- mbnma.run(network, fun=dexp())
#' plot(exponential)
#'
#' # Plot only Emax parameters from an Emax dose-response MBNMA
#' emax <- mbnma.run(network, fun=demax(), method="random")
#' plot(emax, params=c("emax"))
#'
#'
#' #### Forest plots including class effects ####
#' # Generate some classes for the data
#' class.df <- triptans
#' class.df$class <- ifelse(class.df$agent=="placebo", "placebo", "active1")
#' class.df$class <- ifelse(class.df$agent=="eletriptan", "active2", class.df$class)
#' netclass <- mbnma.network(class.df)
#' emax <- mbnma.run(netclass, fun=demax(), method="random",
#'             class.effect=list("ed50"="common"))
#'
#' }
#'
#' @export
plot.mbnma <- function(x, params=NULL, ...) {

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

  # Declare global variables
  labs <- NULL

  fun <- x$model.arg$fun

  # Cannot plot ume model results
  if (x$model.arg$UME==TRUE) {
    stop("UME model results cannot be plotted")
  }

  # Check that specified params are monitored in model
  if (!all(params %in% x[["parameters.to.save"]])) {
    setparam <- x$parameters.to.save[x$parameters.to.save %in% c(x$model.arg$fun$params, toupper(x$model.arg$fun$params))]
    stop(paste0("Variable 'params': Must contain elements of set {", paste(setparam, collapse = ", "), "}"))
  }

  # Check that specified params are modelled using relative effects
  for (i in seq_along(params)) {
    if (!length(grep(paste0("^", params[i], "\\["), rownames(x$BUGSoutput$summary))>1)) {
      stop(paste0(params[i], " does not vary by agent or class"))
    }
  }

  # Add all available params if is.null(params)
  if (is.null(params)) {
    params <- vector()

    relparams <- names(x$model.arg$fun$apool)[which(fun$apool=="rel")]

    # Add d
    params <- append(params,
                     x[["parameters.to.save"]][x[["parameters.to.save"]] %in% relparams])

    # Add D
    params <- append(params,
                     x[["parameters.to.save"]][x[["parameters.to.save"]] %in% toupper(relparams)])

    # Add non-parametric
    params <- append(params,
                     x[["parameters.to.save"]][x[["parameters.to.save"]] %in% "d.1"])

    if (length(params)==0) {
      stop("No dose-response relative effects can be identified from the model")
    }
  }


  if ("nonparam" %in% fun$name) {

    np.df <- as.data.frame(x$BUGSoutput$summary[grepl("d\\.1", rownames(x$BUGSoutput$summary)), c(3,5,7)])

    np.df$agent <- as.numeric(gsub("(d\\.1\\[)([0-9]+)\\,([0-9]+)\\]", "\\3", rownames(np.df)))
    np.df$dose <- as.numeric(gsub("(d\\.1\\[)([0-9]+)\\,([0-9]+)\\]", "\\2", rownames(np.df)))

    np.df$agent <- x$network$agents[np.df$agent]
    np.df <- subset(np.df, !(agent!="Placebo" & dose==1))
    np.df$treatment <- x$network$treatments

    np.df$dose <-
      sapply(np.df$treatment, FUN = function(x) {as.numeric(strsplit(x, split="_", fixed = TRUE)[[1]][2])})

    if (np.df$`50%`[1]!=0 & np.df$`2.5%`[1]!=0) {
      row <- np.df[0,]
      row[,1:3] <- 0
      row$treatment <- "Placebo_0"
      row$agent <- "Placebo"
      row$dose <- 1
      np.df <- rbind(row, np.df)
    }

    # Plot faceted by agent as dose-response splitplot
    # Add intercept for all agents
    agents <- unique(np.df$agent)
    agents <- agents[agents!="Placebo"]
    for (i in seq_along(agents)) {
      row <- np.df[np.df$agent=="Placebo",]
      row$agent <- agents[i]
      np.df <- rbind(row, np.df)
    }
    np.df <- np.df[np.df$agent!="Placebo",]


    # Generate forest plot
    g <- ggplot2::ggplot(np.df, ggplot2::aes(x=dose, y=`50%`))+#, ...) +
      ggplot2::geom_point() +
      ggplot2::geom_errorbar(ggplot2::aes(ymin=`2.5%`, ymax=`97.5%`, width=.05)) +
      ggplot2::facet_wrap(~factor(agent), scales = "free_x") +
      ggplot2::xlab("Dose") +
      ggplot2::ylab("Effect size on link scale") +
      ggplot2::labs(caption = paste0(x$model.arg$fun$direction, " monotonic dose-response function")) +
      theme_mbnma()


  } else {

    # Compile parameter data into one data frame
    mcmc <- x$BUGSoutput$sims.list
    plotdata <- data.frame(Var2=NA, value=NA, doseparam=NA)
    for (i in seq_along(params)) {
      paramdata <- reshape2::melt(mcmc[[params[i]]])

      # For agent-specific DR functions
      if (length(fun$name)>1) {

        subcodes <- mult2agent(fun, params[i])
        if ("Placebo" %in% x$network$agents) {
          subcodes <- subcodes[!subcodes==1]
          subcodes <- subcodes-1
        }

        if (ncol(paramdata)>1) {
          paramdata <- paramdata[,2:3]
          paramdata$Var2 <- subcodes[paramdata$Var2]
        } else {
          paramdata <- cbind(rep(subcodes, nrow(paramdata)), paramdata)
          colnames(paramdata) <- c("Var2", "value")
        }
      } else {
        paramdata <- paramdata[,2:3]
      }

      paramdata$doseparam <- rep(params[i], nrow(paramdata))
      plotdata <- rbind(plotdata, paramdata)
    }
    plotdata <- plotdata[-1,]

    if ("Placebo" %in% x$network$agents) {
      plotdata[["param"]] <- plotdata[["Var2"]] + 1
    } else {
      plotdata[["param"]] <- plotdata[["Var2"]]
    }
    # plotdata[["param"]] <- as.numeric(gsub("(.+\\[)([0-9]+)(\\])", "\\2", rownames(plotdata)))

    plotdata$level <- "agent"
    plotdata$level[grepl("[A-Z]", plotdata$doseparam)] <- "class"
    # plotdata$level[grepl("A-Z", rownames(plotdata))] <- "class"

    temp.df <- plotdata %>%
      subset(level=="agent") %>%
      dplyr::mutate(labs=x$network$agents[param])

    if ("class" %in% plotdata$level) {
      class.df <- plotdata %>%
        subset(level=="class") %>%
        dplyr::mutate(labs=x$network$classes[param])

      temp.df <- rbind(temp.df, class.df)
    }
    plotdata <- temp.df

    if (any(is.na(plotdata$labs))) {
      stop("Cannot identify labels for agents/classes in 'x'")
    }

    g <- ggplot2::ggplot(plotdata, ggplot2::aes(y=value, x=labs)) +
      ggdist::stat_halfeye() +
      ggplot2::facet_wrap(~doseparam, scales="free") +
      ggplot2::coord_flip()

    # Axis labels
    g <- g + ggplot2::xlab("Agent / Class") +
      ggplot2::ylab("Effect size") +
      theme_mbnma()
  }

  graphics::plot(g)
  return(invisible(g))
}





#' Rank parameter estimates
#'
#' Only parameters that vary by agent/class can be ranked.
#'
#' @param lower_better Indicates whether negative responses are better (`TRUE`) or positive responses are better (`FALSE`)
#' @param to.rank A numeric vector containing the codes for the agents/classes you wish to rank.
#' If left `NULL` then all agents/classes (depending on the value assigned to `level`) in
#' the model will be ranked. Included codes must be greater than
#' `2` if placebo has been modelled, since placebo cannot be included in the ranking
#' @param level Can be set to `"agent"` to rank across different agents or `"class"` to rank
#' across different classes.
#' @param params A character vector of named parameters in the model that vary by either agent
#' or class (depending on the value assigned to `level`). If left as `NULL` (the default), then
#' ranking will be calculated for all available parameters that vary by agent/class.
#' @param ... Arguments to be passed to methods
#' @inheritParams predict.mbnma
#' @inheritParams rank
#'
#' @details Ranking cannot currently be performed on non-parametric dose-response MBNMA
#'
#' @return An object of `class("mbnma.rank")` which is a list containing a summary data
#' frame, a matrix of rankings for each MCMC iteration, a matrix of probabilities
#' that each agent has a particular rank, and a matrix of cumulative ranking probabilities
#' for each agent, for each parameter that has been ranked.
#'
#' @examples
#' \donttest{
#' # Using the triptans data
#' network <- mbnma.network(triptans)
#'
#' # Rank selected agents from a log-linear dose-response MBNMA
#' loglin <- mbnma.run(network, fun=dloglin())
#' ranks <- rank(loglin, to.rank=c("zolmitriptan", "eletriptan", "sumatriptan"))
#' summary(ranks)
#'
#' # Rank only ED50 parameters from an Emax dose-response MBNMA
#' emax <- mbnma.run(network, fun=demax(), method="random")
#' ranks <- rank(emax, params="ed50")
#' summary(ranks)
#'
#'
#' #### Ranking by class ####
#' # Generate some classes for the data
#' class.df <- triptans
#' class.df$class <- ifelse(class.df$agent=="placebo", "placebo", "active1")
#' class.df$class <- ifelse(class.df$agent=="eletriptan", "active2", class.df$class)
#' netclass <- mbnma.network(class.df)
#' emax <- mbnma.run(netclass, fun=demax(), method="random",
#'             class.effect=list("ed50"="common"))
#'
#' # Rank by class, with negative responses being worse
#' ranks <- rank(emax, level="class", lower_better=FALSE)
#' print(ranks)
#'
#' # Print and generate summary data frame for `mbnma.rank` object
#' summary(ranks)
#' print(ranks)
#'
#' # Plot `mbnma.rank` object
#' plot(ranks)
#' }
#'
#' @export
rank.mbnma <- function(x, params=NULL, lower_better=TRUE, level="agent", to.rank=NULL, ...) {

  # Checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(x, classes="mbnma", add=argcheck)
  checkmate::assertCharacter(params, null.ok=TRUE, add=argcheck)
  checkmate::assertLogical(lower_better, add=argcheck)
  checkmate::assertChoice(level, choices = c("agent","class"), add=argcheck)
  checkmate::reportAssertions(argcheck)

  fun <- x$model.arg$fun
  if ("nonparam" %in% fun$name) {
    stop("Ranking cannot currently be performed for non-parametric models")
  }

  # Cannot rank ume model results
  if (x$model.arg$UME==TRUE) {
    stop("rank() does not work for UME models")
  }

  # Cannot rank dmulti() models
  if (length(x$model.arg$fun$name)>1) {
    stop("Ranking cannot be performed for models with agent-specific dose-response functions\nTry ranking of relative effects instead (generated using get.relative())")
  }

  # Change agent/class to agents/classes
  levels <- ifelse(level=="agent", "agents", "classes")

  if (level=="class") {
    if (is.null(x[["model"]][["data"]]()[["class"]])) {
      stop("`level` has been set to `class` but classes have not been used in the model")
    }
    if (!is.null(to.rank)) {
      warning("Codes in `to.rank` correspond to class codes rather than treatment codes")
    }
  }

  # If treats have not been specified then select all of them - WONT WORK IF PLACEBO NOT INCLUDED
  starttrt <- ifelse(x$network$agents[1]=="Placebo", 2, 1)
  codes.mod <- c(starttrt:max(x[["model"]][["data"]]()[[level]], na.rm=TRUE))
  if (is.null(to.rank)) {
    to.rank <- codes.mod
  } else if (is.numeric(to.rank)) {
    if (!all(to.rank %in% seq(1:max(x[["model"]][["data"]]()[[level]], na.rm=TRUE)))) {
      stop("`to.rank` codes must match those in the dataset for either `agent` or `class`")
    }
  } else if (is.character(to.rank)) {
    if (!all(to.rank %in% x$network[[levels]])) {
      stop("`to.rank` agent/class names must match those in the network for either `agent` or `class`")
    }
    to.rank <- as.numeric(factor(to.rank, levels=x$network[[levels]]))
  }

  # Check if placebo included in rankings
  if (x$network$agents[1]=="Placebo") {
    to.rank <- to.rank-1
    if (any(to.rank==0)) {
      warning("Placebo (d[1] or D[1]) cannot be included in the ranking for relative effects and will therefore be excluded")
      to.rank <- to.rank[to.rank!=0]
    }
    agents <- x$network[[levels]][to.rank+1]
  } else {
    agents <- x$network[[levels]][to.rank]
  }

  # Check parameters to rank
  if (is.null(params)) {
    for (i in seq_along(fun$params)) {
      mcmc <- x$BUGSoutput$sims.list[[names(fun$apool)[i]]]
      pass <- TRUE
      if (is.vector(mcmc)) {
        pass <- FALSE
      }
      if (is.matrix(mcmc)) {
        if (ncol(mcmc)<=1) {
          pass <- FALSE
        }
      }

      if (pass==TRUE) {

        if (level=="agent" & fun$params[i] %in% x$parameters.to.save) {
          params <- append(params, fun$params[i])
        }
        if (level=="class" & toupper(fun$params[i]) %in% x$parameters.to.save) {
          params <- append(params, toupper(fun$params[i]))
        }
      }
    }
  } else {
    for (i in seq_along(params)) {
      if (!(params[i] %in% x[["parameters.to.save"]])) {
        stop(paste0(params[i], " has not been monitored by the model. `params` can only include model parameters that have been monitored and vary by agent/class"))
      }
    }
  }

  rank.result <- list()
  for (i in seq_along(params)) {
    if (params[i] %in% x[["parameters.to.save"]]) {

      if (length(fun$name)==1) {
        # subcodes <- codes.mod
        subrank <- to.rank
        subagents <- agents
      } else {
        subcodes <- mult2agent(fun, params[i])

        if ("Placebo" %in% x$network$agents) {
          subcodes <- subcodes-1
        }
        subrank <- 1:length(subcodes)
        subagents <- agents[subcodes]
      }

      param.mod <- x[["BUGSoutput"]][["sims.list"]][[params[i]]]

      # Check that selected parameter is different over multiple treatments
      #if (!is.matrix(param.mod) | ncol(param.mod)!=length(subrank)) {
      if (!is.matrix(param.mod) | ncol(param.mod)==1) {
        msg <- paste0(params[i], " does not vary by ", level, " and therefore cannot be ranked")
        stop(msg)
      }

      param.mod <- param.mod[,subrank]
      rank.mat <- t(apply(param.mod, MARGIN=1, FUN=function(x) {
        order(order(x, decreasing = lower_better), decreasing=FALSE)
      }))
      colnames(rank.mat) <- subagents

      # Calculate ranking probabilities
      prob.mat <- calcprob(rank.mat, treats=subagents)

      # Calculate cumulative ranking probabilities
      cum.mat <- apply(prob.mat, MARGIN=2,
                       FUN=function(col) {cumsum(col)})

      rank.result[[params[i]]] <-
        list("summary"=sumrank(rank.mat),
             "prob.matrix"=prob.mat,
             "rank.matrix"=rank.mat,
             "cum.matrix"=cum.mat)

    }
  }

  if (!is.null(x$model.arg$regress.vars)) {
    regress.vals <- rep(0, length(x$model.arg$regress.vars))
    names(regress.vals) <- x$model.arg$regress.vars
  } else {
    regress.vals <- NULL
  }

  attributes(rank.result) <- list("class"="mbnma.rank",
                                  "names"=names(rank.result),
                                  "lower_better"=lower_better,
                                  "level"=level,
                                  "regress.vals"=regress.vals
                                  )

  if (length(rank.result)==0) {
    stop(paste0("There are no parameters saved in the model that vary by ", level))
  }

  return(rank.result)
}




#' Predict responses for different doses of agents in a given population based on MBNMA
#' dose-response models
#'
#' Used to predict responses for different doses of agents or to predict
#' the results of a new study. This is calculated by combining
#' relative treatment effects with a given reference treatment response
#' (specific to the population of interest).
#'
#' @param object An S3 object of class `"mbnma"` generated by running
#'   a dose-response MBNMA model
#'
#' @param n.doses A number indicating the number of doses at which to make predictions
#'   within each agent. The default is `30`.
#' @param exact.doses A list of numeric vectors. Each named element in the list corresponds to an
#'   agent (either named similarly to agent names given in the data, or named
#'   correspondingly to the codes for agents given in `mbnma`) and each number within the vector
#'   for that element corresponds to a dose of the agent for which to predict responses.
#'   Doses can only take positive values. For models fitted using `dspline()` making predictions at only a very small
#'   number of doses for each agent may throw an error since it can make the spline difficult to identify.
#' @param E0 An object to indicate the value(s) to use for the response at dose = 0 (i.e.
#'   placebo) in the prediction. This can take a number of different formats depending
#'   on how it will be used/calculated. The default is `0.2` since a default of `0` will typically lead
#'   to non-sensical predictions unless an identify link function has been used for the MBNMA model in `object`.
#'   * `numeric()` A single numeric value representing the deterministic response at dose = 0,
#'   given on the natural scale - so for binomial data, proportions should be given and
#'   for Poisson data, a rate should be given.
#'   * `character()` A single string representing a stochastic distribution for the response
#'   at dose = 0, given on the natural scale - so for binomial data, proportions should be given and
#'   for Poisson data, a rate should be given. This is specified as a random number generator
#'   (RNG) given as a string, and can take any RNG distribution for which a function exists
#'   in R. For example: `"rnorm(n, 7, 0.5)"`.
#'   * `data.frame()` A data frame containing data in the long format (one row per study arm) to be meta-analysed
#'   to estimate the dose = 0 (placebo) response. This could be a set of observational
#'   studies that are specific to the population on which to make
#'   predictions, or it can be a subset of the study arms within the MBNMA dataset
#'   that investigate placebo. See [ref.synth()]
#' @param synth A character object that can take the value `"fixed"` or `"random"` to
#'   specify the the type of pooling to use for synthesis of `E0` if a data frame
#'   has been provided for it. Using `"random"` rather
#'   than `"fixed"` for `synth` will result in wider 95\\% CrI for predictions.
#' @param lim Specifies calculation of either 95% credible intervals (`lim="cred"`) or 95% prediction intervals (`lim="pred"`).
#' @param regress.vals A named numeric vector of effect modifier values at which results should
#'   be predicted. Named elements must match variable names specified in `regress.vars` within
#'   the MBNMA model.
#'
#' @param ... Arguments to be sent to [R2jags::jags()] for synthesis of the network
#'   reference treatment effect (using [ref.synth()])
#'
#'
#' @return An S3 object of class `mbnma.predict` that contains the following
#'   elements:
#'
#' * `predicts` A named list of
#'   matrices. Each matrix contains the MCMC results of predicted responses at
#'   follow-up times specified in `times` for each treatment specified in
#'   `treats`
#'
#' * `likelihood` The likelihood used in the MBNMA model `object`
#'
#' * `link` The link function used in the MBNMA model `object`
#'
#' * `network` The dataset in `mbnma.network` format
#'
#' * `E0` A numeric vector of value(s) used for E0 in the prediction, on the
#'   link scale.
#'
#'
#' @details
#' The range of doses on which to make predictions can be specified in one of two ways:
#'
#' 1. Use `max.dose` and `n.doses` to specify the maximum dose for each agent and the
#' number of doses within that agent for which to predict responses. Doses will be chosen
#' that are equally spaced from zero to the maximum dose for each agent. This is useful
#' for generating plots of predicted responses (using `[plot-mbnma.predict]`) as it will
#' lead to fitting a smooth dose-response curve (provided `n.doses` is sufficiently high).
#'
#' 2. Use `exact.doses` to specify the exact doses for which to predict responses for each
#' agent. This may be more useful when ranking different predicted responses using
#' `[rank-mbnma.predict]`
#'
#' @examples
#' \donttest{
#' # Using the triptans data
#' network <- mbnma.network(triptans)
#'
#' # Run an Emax dose-response MBNMA
#' emax <- mbnma.run(network, fun=demax(), method="random")
#'
#'
#' ###########################
#' ###### Specifying E0 ######
#' ###########################
#'
#' #### Predict responses using deterministic value for E0 ####
#' # Data is binomial so we specify E0 on the natural scale as a probability
#' pred <- predict(emax, E0 = 0.2)
#'
#' # Specifying non-sensical values will return an error
#' #pred <- predict(emax, E0 = -10)
#' ### ERROR ###
#'
#' #### Predict responses using stochastic value for E0 ####
#' # Data is binomial so we might want to draw from a beta distribution
#' pred <- predict(emax, E0 = "rbeta(n, shape1=1, shape2=5)")
#'
#' # Misspecifying the RNG string will return an error
#' #pred <- predict(emax, E0 = "rbeta(shape1=1, shape2=5)")
#' ### ERROR ###
#'
#'
#' #### Predict responses using meta-analysis of dose = 0 studies ####
#'
#' # E0 is assigned a data frame of studies to synthesis
#' # Can be taken from placebo arms in triptans dataset
#' ref.df <- network$data.ab[network$data.ab$agent==1,]
#'
#' # Synthesis can be fixed/random effects
#' pred <- predict(emax, E0 = ref.df, synth="random")
#'
#'
#'
#' ######################################################################
#' #### Specifying which doses/agents for which to predict responses ####
#' ######################################################################
#'
#' # Change the number of predictions for each agent
#' pred <- predict(emax, E0 = 0.2, n.doses=20)
#' pred <- predict(emax, E0 = 0.2, n.doses=3)
#'
#' # Specify several exact combinations of doses and agents to predict
#' pred <- predict(emax, E0 = 0.2,
#'             exact.doses=list("eletriptan"=c(0:5), "sumatriptan"=c(1,3,5)))
#' plot(pred) # Plot predictions
#'
#' # Print and summarise `mbnma.predict` object
#' print(pred)
#' summary(pred)
#'
#' # Plot `mbnma.predict` object
#' plot(pred)
#' }
#'
#' @export
predict.mbnma <- function(object, n.doses=30, exact.doses=NULL,
                          E0=0.2, synth="fixed", lim="cred",
                          regress.vals=NULL,
                          ...) {
  ######## CHECKS ########

  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(object, "mbnma", add=argcheck)
  # checkmate::assertList(max.doses, types="numeric", null.ok=TRUE, add=argcheck)
  checkmate::assertList(exact.doses, types="numeric", null.ok=TRUE, add=argcheck)
  checkmate::assertInt(n.doses, lower=2, add=argcheck)
  checkmate::assertChoice(synth, choices=c("random", "fixed"), add=argcheck)
  checkmate::assertChoice(lim, choices=c("cred", "pred"), add=argcheck)
  checkmate::assertNumeric(regress.vals, names = "named", null.ok=TRUE, add=argcheck)
  checkmate::reportAssertions(argcheck)

  agents <- object$model$data()$agent
  mbnma.agents <- object$network[["agents"]]

  # Check regress.vals
  check.predreg(mbnma=object, regress.vals=regress.vals)

  # Checks for doses
  doses <- NULL
  max.doses <- NULL
  if (!is.null(exact.doses) & !is.null(max.doses)) {
    warning("`exact.dose` and `max.doses` have both been specified in the arguments, yet only one of these can be used. Deferring to using argument given for `exact.doses`")
    max.doses <- NULL
  }
  if (!is.null(exact.doses)) {
    doses <- exact.doses

    # Add placebo to exact doses if missing (required for plotting of predictions etc.)
    if (!"Placebo" %in% names(doses) & length(doses)!=length(object$network$agents)) {
      doses <- c("Placebo"=0, doses)
    }

    # Merge list elements if they have the same agent name
    if (!is.null(names(doses))) {
      for (i in seq_along(doses)) {
        if (names(doses)[i] %in% names(doses)[-i]) {

          drop <- which(names(doses) %in% names(doses)[i])[-1]

          for (k in seq_along(drop)) {
            doses[[i]] <- sort(append(doses[[i]], doses[[drop[k]]]))
          }
          doses <- doses[-drop]
        }
      }
    }

  } else if (!is.null(max.doses)) {
    doses <- max.doses

    for (i in seq_along(doses)) {
      doses[[i]] <- signif(seq(0,
                               max(doses[[i]]),
                               length.out=n.doses), 2)
    }

  }

  if (!is.null(doses)) {
    if (length(doses)>max(agents, na.rm=TRUE)) {
      stop("A greater number of agents have been supplied in either `max.doses` or `exact.doses` than are present in the model")
    }

    # If named elements of list are not numeric, check if they match agents in mbnma
    match.pass <- TRUE
    if (is.null(names(doses))) {
      if (length(doses)!=length(mbnma.agents)) {
        stop("If elements in `max.doses` or `exact.doses` are not named then there must be the same number of elements as there are agents in the model, so that they correspond to agent codes")
      }
      names(doses) <- mbnma.agents
      agent.num <- 1:length(mbnma.agents)
    } else if (!is.null(names(doses))) {
      if (any(is.na(suppressWarnings(as.numeric(names(doses)))))) {
        if (!all(names(doses)[names(doses)!="Placebo"] %in% mbnma.agents)) {
          match.pass <- FALSE
        }
        agent.num <- match(names(doses), mbnma.agents)
      } else {
        # if (!all(names(doses) %in% c(1:max(agents, na.rm=TRUE)))) {
        #   match.pass <- FALSE
        # }
        #agent.num <- as.numeric(names(doses)) # Add an agent numerical identifier for included agents
        agent.num <- 1:length(unique(names(doses)))
      }
      if (match.pass==FALSE) {
        stop("Element names in `doses` must correspond either to agent names in data or agent codes in `object`")
      }
    }



    for (i in seq_along(doses)) {
      if (!all(doses[[i]]>=0)) {
        stop(paste0("Doses given in `doses` must be positive. Agent ", names(doses)[i], " contains negative values."))
      }
    }
  } else {
    # Automatically generate doses list for treatments included in data
    # if (any(c("rcs", "bs", "ns", "ls") %in% object$model.arg$fun)) {
    #   dose <- as.vector(object$model$data()$spline[,,1])
    # } else {
    #
    # }
    dose <- as.vector(object$model.arg$jagsdata$dose)

    agent <- as.vector(agents)

    df <- data.frame("agent"=agent, "dose"=dose)
    df <- unique(df[stats::complete.cases(df),])
    df <- dplyr::arrange(df, agent, dose)
    df$agent <- factor(df$agent, labels=mbnma.agents)

    doses <- list()
    for (i in seq_along(mbnma.agents)) {
      doses[[mbnma.agents[i]]] <- signif(seq(0,
                                             max(df$dose[df$agent==mbnma.agents[i]]),
                                             length.out=n.doses), 2)
    }
    agent.num <- c(1:length(mbnma.agents))
  }



  # Set model arguments
  # if (length(object$model.arg$class.effect)>0) {
  #   stop("predict() currently does not work with models that use class effects")
  # }
  if ("nonparam" %in% object$model.arg$fun$name) {
    stop("predict() does not work with non-parametric dose-response functions")
  }
  if (object$model.arg$UME==TRUE) {
    stop("predict() does not work with UME models")
  }
  if (length(object$model.arg$fun$name)>1) {

    # Add posvec here?

    # funs <- c(NA, 1,1,2,3, rep(object$model.arg$knots-1, 3))
    # names(funs) <- c("user", "linear", "exponential", "emax", "emax.hill", "rcs", "bs", "ns")
    # funs <- funs[names(funs) %in% object$model.arg$fun]
    #
    # # Find the location of the agents within the vector of agent names in network
    # funi <- which(object$network$agents %in% names(doses))
    # X <- sapply(object$model.arg$fun[funi], function(x) which(x==names(funs)))
  }

  link <- object$model.arg$link

  n <- object$BUGSoutput$n.keep * object$BUGSoutput$n.chains

  # Write dose-response function used in the model
  DR <- suppressMessages(
    write.dose.fun(fun=object$model.arg$fun,
                   effect="abs"
    )[[1]])
  DR <- gsub("(^.+<-)(.+)", "\\2", DR)

  if (length(object$model.arg$fun$name)>1) {
    DR <- DR[-1]
  }

  # Get dose-response parameter estimates
  betaparams <- get.model.vals(object)

  # Get regression parameter estimates and multiply by regress.vals
  if (!is.null(regress.vals)) {
    regress <- get.regress.vals(object, regress.vals, sum=TRUE)
  }

  # Identify E0
  if (is.null(E0)) {
    stop("`E0` has not been given a value. Responses cannot be predicted without a value(s) for dose = 0 (placebo)")
  }
  if (!any(class(E0) %in% c("numeric", "character", "data.frame"))) {
    stop("`E0` can only be of type `numeric()`, `character()` or `data.frame()`")
  }

  if ((is.numeric(E0) | is.character(E0)) & length(E0)!=1) {
    stop("`E0` must take a single numeric (deterministic E0) or character (stochastic E0) value, or be provided with a data frame so that it may be estimated from the data")
  }
  if (is.character(E0)) {
    if (!grepl("^r[a-z]+\\(n,.+\\)", E0)) {
      stop("Distribution for `E0` does not match stochastic random number generator format.\nFormat must take any R random number generator function")
    }
    E0 <- eval(parse(text=E0))
  }

  if ((is.numeric(E0) | is.character(E0))) {
    E0 <- rescale.link(E0, direction="link", link=link)
  } else if (is.data.frame(E0)) {
    E0 <- ref.synth(data.ab=E0, mbnma=object, synth=synth, ...)

    if (!("sd.mu" %in% names(E0))) {
      E0 <- E0$m.mu
    } else {
      E0 <- stats::dnorm(E0$m.mu, E0$sd.mu)
    }
  }

  # Ensure prediction intervals are used where appropriate
  if (lim=="pred" & object$model.arg$method=="random") {
    addsd <- TRUE

    if (!"sd" %in% object$parameters.to.save) {
      stop(crayon::red("'sd' not included in parameters.to.save - cannot calculate prediction intervals"))
    }

    message("Bayesian prediction intervals to be calculated")
  } else {
    addsd <- FALSE
  }


  predict.result <- list()

  # Add spline basis matrix
  splineopt <- c("rcs", "bs", "ns", "ls")
  fun <- object$model.arg$fun
  if (any(splineopt %in% fun$name)) {

    splinedoses <- doses

    index <- match(names(doses), object$network$agents)

    # If there are multiple spline functions
    if ("posvec" %in% names(fun)) {
      posvec <- fun$posvec
    } else {
      posvec <- rep(1, length(index))
    }
    for (i in seq_along(index)) {
      if (fun$name[posvec[index[i]]] %in% splineopt) {
        #print(agent.num[i])
        #print((object$network$data.ab$dose[object$network$data.ab$agent==agent.num[i]]))
        splinedoses[[i]] <- t(genspline(splinedoses[[i]],
                                   spline = fun$name[posvec[index[i]]],
                                   knots=fun$knots[[posvec[index[i]]]],
                                   degree = fun$degree[posvec[index[i]]],
                                   max.dose=max(object$network$data.ab$dose[object$network$data.ab$agent==agent.num[i]])
                                   ))

      } else {
        # Generate empty matrix for non-spline DR functions
        splinedoses[[i]] <- matrix(0, ncol=length(splinedoses[[i]]), nrow=2)
      }
    }
  }

  # Replace segments of dose-response function string with values for prediction
  for (i in seq_along(doses)) {
    predict.result[[names(doses)[i]]] <- list()
    for (k in seq_along(doses[[i]])) {
      if (names(doses)[i] %in% c("1", "Placebo") | doses[[i]][k]==0) {
        # Ensures reference agent (placebo) takes E0
        # This should be changed if intercept is relaxed
        pred <- E0

      } else {
        if (length(object$model.arg$fun$name)>1) {
          pos <- object$model.arg$fun$posvec[which(object$network$agents==names(doses)[i])]
          tempDR <- DR[pos]
        } else {
          tempDR <- DR
        }

        tempDR <- gsub("\\[agent\\[i,k\\]\\]", "", tempDR)
        tempDR <- gsub("\\[i,k\\]", "", tempDR)
        tempDR <- gsub("(\\[i,k,)([0-9\\])", "[\\2", tempDR) # For splines

        dose <- doses[[i]][k]
        if (any(c("rcs", "bs", "ns", "ls") %in% object$model.arg$fun$name)) {
          spline <- splinedoses[[i]][,k]
        }

        for (param in seq_along(betaparams)) {
          if (is.vector(betaparams[[param]])) {
            assign(paste0("s.", names(betaparams)[param]),
                   betaparams[[param]])
          } else if (is.matrix(betaparams[[param]])) {
            if (ncol(betaparams[[param]])>1) {
              # Look for correct column index for each beta param
              if (length(object$model.arg$fun$name)>1) {
                vec <- object$model.arg$fun$posvec[1:which(object$network$agents==names(doses)[i])]
                vec <- table(vec)[names(table(vec))==pos]
                colnum <- vec

                if (names(vec)=="1") {
                  # Check if placebo in dataset
                  if (object$network$agents[1]=="Placebo") {
                    colnum <- colnum - 1
                  }
                }

              } else {
                colnum <- which(object$network$agents==names(doses)[i])

                # Check if placebo in dataset
                if (object$network$agents[1]=="Placebo") {
                  colnum <- colnum - 1
                }
              }

            } else {
              colnum <- 1
            }

            # Check for if param assignment is valid
            if (colnum<=ncol(betaparams[[param]])) {
              assign(paste0("s.", names(betaparams)[param]),
                     betaparams[[param]][,colnum]
              )
            } else {
              assign(paste0("s.", names(betaparams)[param]), NULL)
            }

          }
        }

        # Add regression
        if (!is.null(regress.vals)) {
          tempDR <- paste0(tempDR, " + reg")
          reg <- regress[,colnum]
        }

        # Evaluate dose-response string for prediction
        chunk <- eval(parse(text=tempDR))
        if (is.list(chunk)) {
          chunk <- chunk[[1]]
        }

        # Incorporate between-study SD
        if (addsd==TRUE) {
          mat <- matrix(nrow=length(chunk), ncol=2)
          mat[,1] <- chunk
          mat[,2] <- stats::median(object$BUGSoutput$sims.list[["sd"]])
          chunk <- apply(mat, MARGIN=1, FUN=function(x) stats::rnorm(1, x[1], x[2]))
        }


        pred <- E0 + chunk
        if (length(pred)<=1) {stop()}
      }

      # Convert to natural scale using link function
      pred <- rescale.link(pred, direction="natural", link=link)

      predict.result[[names(doses)[i]]][[as.character(doses[[i]][k])]] <-
        pred
    }
  }

  output <- list("predicts"=predict.result,
                 "likelihood"=object$model.arg$likelihood, "link"=object$model.arg$link,
                 "network"=object$network, "lim"=lim, "E0"=E0)

  # Add 0 values for missing regression values
  if (is.null(regress.vals) & !is.null(object$model.arg$regress)) {
    regress.vals <- rep(0, ncol(object$model.arg$regress.mat))
    names(regress.vals) <- colnames(object$model.arg$regress.mat)
  }
  output$regress.vals <- regress.vals

  class(output) <- "mbnma.predict"

  return(output)
}




#' Print summary of MBNMA results to the console
#' @param object An S3 object of class `"mbnma"` generated by running
#'   a dose-response MBNMA model
#' @param digits The maximum number of digits for numeric columns
#' @param ... additional arguments affecting the summary produced
#'
#' @export
summary.mbnma <- function(object, digits=4, ...) {
  checkmate::assertClass(object, "mbnma")

  # State that function does not work if "parameters.to.save" has been specified
  if (!is.null(object$model.arg$parameters.to.save)) {
    stop("Cannot use `summary()` method if `parameters.to.save` have been assigned. Use `print()` instead.")
  }
  if (any(object$model.arg$fun$name %in% c("nonparam"))) {
    stop("Cannot use `summary()` method for non-parametric dose-response functions. Use `print()` instead.")
  }

  # Check for rhat < 1.02
  rhat.warning(object)

  ##### Overall section #####
  overall.str(object)

  # Print method section
  cat(method.str(object))

  # Print treatment-level section
  treat.str(object, digits=digits)

  # Class-effect section
  class.str(object, digits=digits)

  # Print regression section
  regress.str(object, digits=digits)

  # Model fit statistics section
  modfit.sect <- modfit.str(object)

  output <- paste("\n\n", modfit.sect, sep="")
  cat(output, ...)
}