R/predict.functions.R

In MBNMAdose: Dose-Response MBNMA Models

# Functions for predicting responses in MBNMAdose
# Author: Hugo Pedder
# Date created: 2019-04-25

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



#' Get MBNMA model values for dose-response parameters
#'
#' @inheritParams predict.mbnma
#' @noRd
get.model.vals <- function(mbnma) {

  fun <- mbnma$model.arg$fun

  betaparams <- list()
  for (i in seq_along(fun$apool)) {
    temp <- vector()
    res.mat <- mbnma$BUGSoutput$sims.list
    if (fun$apool[i] %in% c("rel")) {
      temp <- as.matrix(res.mat[[names(fun$apool)[i]]], ncol=1)
    } else if (fun$apool[i] %in% c("common")) {
      temp <- as.vector(res.mat[[names(fun$apool)[i]]])
    } else if (fun$apool[i] %in% "random") {

      # Incorporates SD from between-study SD for ABSOLUTE pooling
      mat <- matrix(nrow=mbnma$BUGSoutput$n.sims, ncol=2)
      mat[,1] <- res.mat[[names(fun$apool)[i]]]
      mat[,2] <- stats::median(res.mat[[paste0("sd.", names(fun$apool)[i])]])
      mat <- apply(mat, MARGIN=1, FUN=function(x) stats::rnorm(1, x[1], x[2]))

      temp <- as.vector(mat)
    } else if (suppressWarnings(!is.na(as.numeric(fun$apool[i])))) {
      temp <- rep(as.numeric(fun$apool[i]), mbnma$BUGSoutput$n.sims)
    }

    betaparams[[fun$bname[i]]] <- temp
  }

  return(betaparams)

}





#' Synthesise single arm dose = 0 / placebo studies to estimate E0
#'
#' Synthesises single arm studies to estimate E0. Used in predicting responses from a
#' dose-response MBNMA.
#'
#' @inheritParams predict.mbnma
#' @inheritParams R2jags::jags
#' @param data.ab A data frame of arm-level data in "long" format containing the
#'   columns:
#'   * `studyID` Study identifiers
#'   * `y` Numeric data indicating the aggregate response for a continuous outcome. Required for
#'   continuous data.
#'   * `se` Numeric data indicating the standard error for a given observation. Required for
#'   continuous data.
#'   * `r` Numeric data indicating the number of responders within a study arm. Required for
#'   binomial or poisson data.
#'   * `n` Numeric data indicating the total number of participants within a study arm. Required for
#'   binomial data
#'   * `E` Numeric data indicating the total exposure time for participants within a study arm. Required
#'   for poisson data.
#' @param mbnma An S3 object of class `"mbnma"` generated by running
#'   a dose-response MBNMA model
#'
#' @details `data.ab` can be a collection of studies that closely resemble the
#'   population of interest intended for the prediction, which could be
#'   different to those used to estimate the MBNMA model, and could include
#'   single arms of RCTs or observational studies. If other data is not
#'   available, the data used to estimate the MBNMA model can be used by
#'   selecting only the studies and arms that investigate dose = 0 (placebo).
#'
#'   Defaults for `n.iter`, `n.burnin`, `n.thin` and `n.chains` are those used to estimate
#'   `mbnma`.
#'
#' @return A list of named elements corresponding to E0 and the between-study standard deviation for
#'   E0 if `synth="random"`. Each element contains the full MCMC results from the synthesis.
#'
#' @examples
#' \donttest{
#' # Using the triptans data
#' network <- mbnma.network(triptans)
#'
#' # Run an Emax dose-response MBNMA
#' emax <- mbnma.run(network, fun=demax(), method="random")
#'
#' # Data frame for synthesis can be taken from placebo arms
#' ref.df <- triptans[triptans$agent=="placebo",]
#'
#' # Meta-analyse placebo studies using fixed treatment effects
#' E0 <- ref.synth(ref.df, emax, synth="fixed")
#' names(E0)
#'
#' # Meta-analyse placebo studies using random treatment effects
#' E0 <- ref.synth(ref.df, emax, synth="random")
#' names(E0)
#' }
#'
#' @export
ref.synth <- function(data.ab, mbnma, synth="fixed",
                      n.iter=mbnma$BUGSoutput$n.iter,
                      n.burnin=mbnma$BUGSoutput$n.burnin,
                      n.thin=mbnma$BUGSoutput$n.thin,
                      n.chains=mbnma$BUGSoutput$n.chains,
                      ...) {

  # First need to validate data.frame to check dataset is in correct format...maybe another function for this
  # Change it to correct format if it is not already
  data.ab <- E0.validate(data.ab, likelihood=mbnma$model.arg$likelihood)[["data.ab"]]

  # Run checks
  argcheck <- checkmate::makeAssertCollection()

  checkmate::assertClass(mbnma, "mbnma", add=argcheck)
  checkmate::assertChoice(synth, choices=c("random", "fixed"), add=argcheck)
  checkmate::assertInt(n.iter, lower=1, add=argcheck)
  checkmate::assertInt(n.burnin, lower=1, add=argcheck)
  checkmate::assertInt(n.thin, lower=1, add=argcheck)
  checkmate::assertInt(n.chains, lower=1, add=argcheck)
  checkmate::reportAssertions(argcheck)

  # To get model for meta-analysis of placebo must create v similar model
  #to study model
  # Do all the mbnma.write bits but without the consistency bits

  # Calculate outcome measure scale
  om <- calcom(data.ab=data.ab, likelihood=mbnma$model.arg$likelihood, link=mbnma$model.arg$link)

  jagsmodel <- write.E0.synth(synth=synth,
                              likelihood=mbnma$model.arg$likelihood,
                              link=mbnma$model.arg$link,
                              om=om
  )

  parameters.to.save <- c("m.mu", "resdev", "totresdev")
  if (synth=="random") {
    parameters.to.save <- append(parameters.to.save, "sd.mu")
  }

  # Run synthesis
  jags.result <- suppressWarnings(
    mbnma.jags(data.ab, model=jagsmodel,
               parameters.to.save=parameters.to.save,
               likelihood=mbnma$model.arg$likelihood,
               link=mbnma$model.arg$link,
               n.iter=n.iter, n.burnin=n.burnin,
               n.thin=n.thin, n.chains=n.chains,
               ...)[["jagsoutput"]]
  )

  result <- list(jagsmod=jags.result,
                 m.mu=jags.result$BUGSoutput$sims.list[["m.mu"]])
  if (synth=="random") {
    result[["sd.mu"]] <- jags.result$BUGSoutput$sims.list[["sd.mu"]]
  }

  if (any(jags.result$BUGSoutput$summary[,
                                         colnames(jags.result$BUGSoutput$summary)=="Rhat"
                                         ]>1.02)) {
    warning("Rhat values for parameter(s) in reference treatment synthesis model are >1.02. Suggest running for more iterations.")
  }

  return(result)

}





#' Checks the validity of ref.estimate
#'
#' Ensures `E0` takes the correct form to allow for synthesis of network
#' reference treatment response
#'
#' @param likelihood A character object that can take any of `"normal"`, `"binomial"` or `"poisson"`
#'
#' @inheritParams ref.synth
#' @noRd
E0.validate <- function(data.ab, likelihood=NULL) {

  data.ab <- data.ab[,names(data.ab) %in% c("studyID", "dose", "agent", "treatment", "r", "n", "E", "y", "se")]

  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertDataFrame(data.ab, any.missing=FALSE, add=argcheck)
  checkmate::assertNames(names(data.ab), must.include = c("studyID"), add=argcheck)
  checkmate::reportAssertions(argcheck)

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


  # Sort data.ab
  data.ab <- dplyr::arrange(data.ab, studyID)

  data.ab <- data.ab %>%
    dplyr::group_by(studyID) %>%
    dplyr::mutate(narm=dplyr::n())

  if (any(data.ab$narm>1)) {
    stop("Studies in `data.ab` contain >1 arm. ref.synth() only pools single arms.")
  }

  print("Data frame must contain only data from reference treatment")

  #### Prepare data frame ####
  # Add arm index (=1 since only one arm in each study)
  data.ab[["arm"]] <- 1
  data.ab[["narm"]] <- 1

  # Ensuring studies are numbered sequentially
  if (!is.numeric(data.ab[["studyID"]])) {
    print("Studies being recoded to allow sequential numbering")
    data.ab <- transform(data.ab,studyID=as.numeric(factor(studyID, levels=as.character(unique(data.ab$studyID)))))
  } else if (all(abs(diff(data.ab[["studyID"]])) != TRUE)) {
    print("Studies being recoded to allow sequential numbering")
    data.ab <- transform(data.ab,studyID=as.numeric(factor(studyID, levels=as.character(unique(data.ab$studyID)))))
  }
  data.ab <- dplyr::arrange(data.ab, studyID)

  data.ab <- add_index(data.ab)

  return(data.ab)

}




#' Rescale data depending on the link function provided
#'
#' @inheritParams mbnma.run
#' @param x A numeric vector of data to be rescaled
#' @param direction Can take either `"link"` to convert data to a particular scale
#'   as defined by the `link` function, or `"natural"` to return it to the natural scale.
#'
#' @return A rescaled numeric vector
#'
#' @export
rescale.link <- function(x, direction="link", link="logit") {

  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertNumeric(x, add=argcheck)
  checkmate::assertChoice(direction, choices=c("natural", "link"), null.ok = FALSE, add=argcheck)
  checkmate::assertChoice(link, choices=c("logit", "identity", "log", "probit", "cloglog", "smd"), null.ok = FALSE, add=argcheck)
  checkmate::reportAssertions(argcheck)

  if (direction=="link") {
    if (link=="logit") {
      x <- stats::qlogis(x)
    } else if (link=="log") {
      x <- log(x)
    } else if (link=="probit") {
      x <- stats::qnorm(x)
    } else if (link=="cloglog") {
      x <- log(-log(1-x))
    }

  } else if (direction=="natural") {
    if (link=="logit") {
      x <- exp(x) / (1+exp(x))
    } else if (link=="log") {
      x <- exp(x)
    } else if (link=="cloglog") {
      x <- 1-exp(-exp(x))
    } else if (link=="probit") {
      x <- stats::pnorm(x)
    }
  }
  return(x)
}





#' Calculates values for EDx from an Emax model, the dose at which x% of the maximal response (Emax)
#' is reached
#'
#' @inheritParams devplot
#' @param x A numeric value between 0 and 100 for the dose at which x% of the maximal response (Emax)
#' should be calculated
#'
#' @return A data frame of posterior EDx summary values for each agent
#'
#' @export
calc.edx <- function(mbnma, x=50) {

  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(mbnma, "mbnma", add=argcheck)
  checkmate::assertNumeric(x, len = 1, lower = 0, upper = 100, add=argcheck)
  checkmate::reportAssertions(argcheck)

  if (!"emax" %in% mbnma$model.arg$fun$name) {
    stop("mbnma must be estimated using fun=demax()")
  }

  if (length(mbnma$model.arg$fun$name)>1) {
    stop("calc.edx cannot be used for agent-specific MBNMA models")
  }

  if ("hill" %in% mbnma$model.arg$fun$params) {
    edx <- mbnma$BUGSoutput$sims.list$ed50 * ((x/(100-x)) ^ (1/mbnma$BUGSoutput$sims.list$hill))
  } else {
    edx <- mbnma$BUGSoutput$sims.list$ed50 * (x/(100-x))
  }

  agents <- mbnma$network$agents[mbnma$network$agents!="Placebo"]
  output <- data.frame()
  for (i in seq_along(agents)) {
    quant <- stats::quantile(edx[,i], probs=c(0.025,0.25,0.5,0.75,0.975), na.rm=TRUE)
    df <- data.frame("agent"=agents[i],
                     "mean"=mean(edx[,i]),
                     "sd"=stats::sd(edx[,i]), stringsAsFactors = TRUE
    )
    output <- rbind(output, cbind(df, t(quant)))
  }

  return(output)
}





#' Get MBNMA model values for regression parameters
#'
#' @param sum Logical object to indicate whether resulting covariate and regressor
#'  values should be summed so that they can be easily added to predictions.
#' @inheritParams predict.mbnma
#' @noRd
get.regress.vals <- function(mbnma, regress.vals, sum=TRUE) {

  # Run checks
  argcheck <- checkmate::makeAssertCollection()
  checkmate::assertClass(mbnma, "mbnma", add=argcheck)
  checkmate::assertNumeric(regress.vals, names="named", add=argcheck)
  checkmate::reportAssertions(argcheck)

  # Define level
  if ("agent" %in% mbnma$model.arg$regress.effect) {
    labs <- mbnma$network$agents[mbnma$network$agents!="Placebo"]
  } else if ("class" %in% mbnma$model.arg$regress.effect) {
    labs <- mbnma$network$classes[mbnma$network$classes!="Placebo"]
  } else {
    labs <- NULL
  }

  outlist <- list()
  outval <- 0
  for (i in seq_along(regress.vals)) {
    res.mat <- as.matrix(mbnma$BUGSoutput$sims.list[[paste0("B.", names(regress.vals)[i])]])
    colnames(res.mat) <- labs

    if ("random" %in% mbnma$model.arg$regress.effect) {
      sd.reg <- stats::median(mbnma$BUGSoutput$sims.list[[paste0("sd.B.", names(regress.vals)[i])]])

      # Incorporates SD from random covariate-by-treatment interaction
      mat <- matrix(nrow=mbnma$BUGSoutput$n.sims, ncol=2)
      mat[,1] <- res.mat
      mat[,2] <- sd.reg
      res.mat <- as.matrix(apply(mat, MARGIN=1, FUN=function(x) stats::rnorm(1, x[1], x[2])))
    }

    # Multiply matrix cols out so that there is one for each agent
    if ("class" %in% mbnma$model.arg$regress.effect) {
      classvec <- mbnma$network$classkey$class[mbnma$network$classkey$class!="Placebo"]
      res.mat <- res.mat[,as.character(classvec)]

    } else if (any(c("common", "random") %in% mbnma$model.arg$regress.effect)) {
      res.mat <- res.mat[,rep(1, length(mbnma$network$agents[mbnma$network$agents!="Placebo"]))]
    }

    # Multiply covariate x variable
    regef <- regress.vals[i] * res.mat
    outlist[[names(regress.vals)[i]]] <- regef
    outval <- outval + (regef)
  }

  if (sum==TRUE) {
    return(outval)
  } else if (sum==FALSE) {
    return(outlist)
  }
}






#' Check regression parameters are specified correctly
#'
#' @noRd
check.predreg <- function(mbnma, regress.vals) {

  # Check regress.vals
  if (!is.null(regress.vals)) {
    if (is.null(mbnma$model.arg$regress.mat)) {
      stop("'regress.vals' has been specified but MBNMA is not a meta-regression model")
    }
    if (!setequal(colnames(mbnma$model.arg$regress.mat), names(regress.vals))) {
      stop(paste0("'regress.vals' must contain a single named regressor value for each covariate specified in the MBNMA model:\n", paste(colnames(mbnma$model.arg$regress.mat), collapse="\n")))
    }
  }
}