Nothing
##############################################
#### Functions for class("mbnma") ####
##############################################
#' Rank parameters from a time-course MBNMA
#'
#' Ranks desired parameters saved from a time-course MBNMA model from "best" to
#' "worst".
#'
#' @param x An S3 object of `class("mbnma")` generated by running
#' a time-course MBNMA model
#' @param param A character object containing any model parameter monitored
#' in `mbnma` for which ranking is desired (e.g. `"beta.1"`, `"emax"`).
#' Parameters must vary by treatment for ranking to be possible. Can also be specified as
#' `"auc"` (see details).
#' @param treats A character vector of treatment/class names (depending on the parameter to be ranked) or
#' a numeric vector of treatment/class codes (as coded in `mbnma`)
#' that indicate which treatments/classes to calculate rankings for. If left `NULL``
#' then rankings will be calculated for all treatments/classes.
#' @param lower_better Indicates whether negative responses are better (`lower_better=TRUE`) or
#' positive responses are better (`lower_better=FALSE`)
#' @param int.range A numeric vector with two elements that indicates the range
#' over which to calculate AUC. Takes the form c(lower bound, upper bound). If left
#' as `NULL` (the default) then the range will be between zero and the maximum follow-up
#' time in the dataset.
#' @param n.iter The number of iterations for which to calculate AUC (if `"auc"` is included in `params`).
#' Must be a positive integer. Default is the value used in `mbnma`.
#' @param ... Arguments to be sent to `integrate()`
#'
#' @return A named list whose elements include:
#' * `summary.rank` A data frame containing
#' mean, sd, and quantiles for the ranks of each treatment given in `treats`
#' * `prob.matrix` A matrix of the proportions of MCMC results for which each
#' treatment/class in `treats` ranked in which position for the given parameter
#' * `rank.matrix` A matrix of the ranks of MCMC results for each treatment/class in
#' `treats` for the given parameter.
#'
#' @details `"auc"` can be specified in `param` to rank treatments based on
#' Area Under the Curve (AUC). This accounts for the effect of multiple
#' time-course parameters simultaneously on the treatment response, but will
#' be impacted by the range of time over which AUC is calculated (`int.range`).
#' This requires integration over `int.range` and can take some time to run (particularly)
#' for spline functions as this uses the trapezoid method rather than adaptive quadrature).
#' Note that `"auc"` can only be calculated at the treatment-level in class effect models.
#'
#' As with other post-estimation functions, `rank()` should only be performed on
#' models which have successfully converged. Note that rankings can be very sensitive to
#' even small changes in treatment effects and therefore failure to converge in only
#' one parameter may have substantial impact on rankings.
#'
#' @examples
#' \donttest{
#' # Create an mb.network object from a dataset
#' network <- mb.network(alog_pcfb)
#'
#' # Run an MBNMA model with an Emax time-course
#' emax <- mb.run(network,
#' fun=temax(pool.emax="rel", method.emax="common",
#' pool.et50="rel", method.et50="random"),
#' intercept=FALSE)
#'
#' # Rank treatments by time-course parameter from the model with lower scores being better
#' rank(emax, param=c("emax"), lower_better=TRUE)
#'
#' # Rank treatments 1-3 by AUC
#' rank(emax, param="auc", treats=c(1:3), lower_better=TRUE,
#' int.range=c(0,20))
#' }
#'
#' @export
rank.mbnma <- function(x, param="auc", lower_better=FALSE, treats=NULL,
int.range=NULL,
n.iter=x$BUGSoutput$n.sims,
...) {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(x, "mbnma", add=argcheck)
checkmate::assertCharacter(param, any.missing=FALSE, unique=TRUE, len=1, add=argcheck)
checkmate::assertLogical(lower_better, null.ok=FALSE, len=1, add=argcheck)
checkmate::assertNumeric(int.range, lower=0, finite=TRUE, any.missing=FALSE, len=2, null.ok=TRUE,
sorted=TRUE, add=argcheck)
#checkmate::assertChoice(level, choices=c("treatment", "class"), add=argcheck)
checkmate::assertInt(n.iter, lower=1, upper=x$BUGSoutput$n.sims, add=argcheck)
checkmate::reportAssertions(argcheck)
# Stop UME models
if (!FALSE %in% x$model.arg$UME) {
stop("rank.mbnma() cannot be used with UME models")
}
# Check param is monitored
if (!"auc" %in% param) {
if (!param %in% x$parameters.to.save) {
stop(paste0(param, " is not a valid paramter saved from the MBNMA model"))
}
}
# Set level
if (grepl("[[:upper:]]", param)) {
level <- "classes"
} else {
level <- "treatments"
}
# if (level=="class" & length(x$model.arg$class.effect)==0) {
# stop("`level` has been specified as `class` yet `x` is not a class effect model")
# }
# level <- ifelse(level=="treatment", "treatments", "classes")
# # Ensure AUC is the last estimate to be called
# if ("auc" %in% param) {
# params <- c(params[params!="auc"], "auc")
# }
# If treats have not been specified then select all of them
if (is.null(treats)) {
treats <- x$network[[level]]
} else if (!is.null(treats)) {
if (is.character(treats)) {
if (!all(treats %in% x$network[[level]])) {
stop("`treats` includes treatments/classes not included in `x`")
}
} else if (is.numeric(treats)) {
if (any(treats > x[["model"]][["data"]]()[["NT"]] | any(treats<1))) {
stop("If given as numeric treatment/class codes, `treats` must be numbered similarly to treatment/class codes in `x`")
}
treats <- x$network[[level]][treats]
}
}
# Provide int.range default values
if ("auc" %in% param) {
if (is.null(int.range)) {
# treatsnum <- which(x$network[[level]] %in% treats)
# fupdata <- x$model.arg$jagsdata
#
# int.max <- max(fupdata$time,
# na.rm = TRUE)
int.max <- round(max(x$network$data.ab$time, na.rm = TRUE))
int.range <- c(0, int.max)
}
}
# Change beta to d (if present) so that it is identified in mcmc output
for (i in 1:4) {
if (paste0("beta.",i) %in% param) {
if (!paste0("beta.",i) %in% x[["parameters.to.save"]]) {
param[which(param==paste0("beta.",i))] <- paste0("d.",i)
}
}
}
if (level=="class" & !grepl("[[:upper:]]", param)) {
stop(paste0("`level` has been specified as `class` but ", param, " is not modelled at the class level"))
}
if (param %in% x[["parameters.to.save"]]) {
param.mod <- x[["BUGSoutput"]][["sims.list"]][[param]]
# Check that selected parameter is different over multiple treatments
if (!is.matrix(param.mod) | ncol(param.mod)<=1) {
msg <- paste0(param, " does not vary by treatment and therefore cannot be ranked by treatment")
stop(msg)
}
param.mod <- param.mod[,which(x$network[[level]] %in% treats)]
rank.mat <- t(apply(param.mod, MARGIN=1, FUN=function(x) {
order(order(x, decreasing = !lower_better), decreasing=FALSE)
}))
colnames(rank.mat) <- treats
# Ranking probabilityes
prob.mat <- calcprob(rank.mat, treats=treats)
# Calculate cumulative ranking probabilities
cum.mat <- apply(prob.mat, MARGIN=2,
FUN=function(col) {cumsum(col)})
rank.result <-
list("param"=param,
"summary"=sumrank(rank.mat),
"prob.matrix"=prob.mat,
"rank.matrix"=rank.mat,
"cum.matrix"=cum.mat,
"lower_better"=lower_better
)
} else if (param=="auc") {
auc <- rankauc(x, lower_better=lower_better,
treats=treats, level=level,
int.range=int.range, n.iter=n.iter, ...)
# Calculate cumulative ranking probs for aux
auc[["cum.matrix"]] <- apply(auc$prob.matrix, MARGIN=2,
FUN=function(col) {cumsum(col)})
auc[["lower_better"]] <- lower_better
auc[["param"]] <- "auc"
rank.result <- auc
}
# rank.result <- list()
# for (i in seq_along(params)) {
# if (level=="class" & !grepl("[[:upper:]]", params[i])) {
# stop(paste0("`level` has been specified as `class` but ", param[i], " is not modelled at the class level"))
# }
#
# if (params[i] %in% x[["parameters.to.save"]]) {
# 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)<=1) {
# msg <- paste0(params[i], " does not vary by treatment and therefore cannot be ranked by treatment")
# stop(msg)
# }
#
# param.mod <- param.mod[,which(x$network[[level]] %in% treats)]
#
# rank.mat <- t(apply(param.mod, MARGIN=1, FUN=function(x) {
# order(order(x, decreasing = !lower_better), decreasing=FALSE)
# }))
# colnames(rank.mat) <- treats
#
# # Ranking probabilityes
# prob.mat <- calcprob(rank.mat, treats=treats)
#
# # 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,
# "lower_better"=lower_better
# )
#
# } else if (params[i]=="auc") {
#
# auc <- rankauc(x, lower_better=lower_better,
# treats=treats, level=level,
# int.range=int.range, n.iter=n.iter, ...)
#
# # Calculate cumulative ranking probs for aux
# auc[["cum.matrix"]] <- apply(auc$prob.matrix, MARGIN=2,
# FUN=function(col) {cumsum(col)})
# auc[["lower_better"]] <- lower_better
#
# rank.result[["auc"]] <- auc
# } else {
# stop(paste0(params[i],
# " is not a valid paramter saved from the MBNMA model"))
# }
# }
class(rank.result) <- "mb.rank"
return(rank.result)
}
#' Print summary MBNMA results to the console
#'
#' @inheritParams predict.mbnma
#' @param ... further arguments passed to `knitr::kable`
#'
#' @return Prints summary details of the model results to the console
#'
#' @export
summary.mbnma <- function(object, ...) {
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.")
}
# Overall section
overall.str(object)
cat("\n\n")
# Print treatment-level section
treat.str(object, ...)
# Class-effect section
class.str(object, ...)
# Correlation section
cor.str(object, ...)
# Model fit statistics section
cat(modfit.str(object))
# Check for rhat < 1.02
cat("\n")
rhat.warning(object)
}
#' Predict effects over time in a given population based on MBNMA time-course
#' models
#'
#' Used to predict effects over time for different treatments or to predict
#' the results of a new study. For MBNMA models that include consistency
#' relative effects on time-course parameters, this is calculated by combining
#' relative treatment effects with a given reference treatment response
#' (specific to the population of interest).
#'
#' By default the network reference treatment baseline (`E0`) and time-course
#' parameter values are set to zero so that `predict()` estimates mean differences
#' (/relative treatment effects) over time versus the network reference treatment.
#'
#' @param object An S3 object of `class("mbnma")` generated by running
#' a time-course MBNMA model
#' @param times A sequence of positive numbers indicating which time points to
#' predict mean responses for (or at which to conduct a node-split if used with `mb.nodesplit()`)
#' @param E0 An object to indicate the value(s) to use for the response at time = 0
#' in the prediction. This can take a number of different formats depending
#' on how it will be used/calculated. The default is 0 but this may lead
#' to non-sensical predictions if Ratio of Means are modeled.
#' * `numeric()` A single numeric value representing the deterministic response at time = 0
#' * `formula()` A formula representing a stochastic distribution for the response
#' at time = 0. 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)`.
#' @param treats A character vector of treatment/class names or a numeric vector of treatment/class codes (as coded
#' in `mbnma`) that indicates which treatments/classes to calculate predictions for. If left as `NULL` then
#' predictions will be calculated for all treatments/classes. Whether the vector should correspond to treatments or
#' classes depends on the value of `level`.
#' @param level Can take either `"treatment"` to make predictions for treatments, or `"class"` to make predictions for classes (in
#' which case `object` must be a class effect model).
#' @param ref.resp An object to indicate the value(s) to use for the reference treatment response in MBNMA models
#' in which the reference treatment response is not estimated within the model (i.e. those that model any time-
#' course parameters using `pool="rel"`). This can take a number of different formats depending
#' on how it will be used/calculated. There are two approaches for this:
#'
#' 1. The reference response can be estimated from a dataset of studies investigating the reference
#' treatment using meta-analysis. This dataset could be a set of observational
#' studies that are specific to the population on which to make
#' predictions, or it could be a subset of the study arms within the MBNMA dataset
#' that investigate the reference treatment. The data should be provided to `ref.resp` as a
#' `data.frame()` containing the data in long format (one row per observation). See [ref.synth()]
#'
#' 2. Values for the reference treatment response can be assigned to different time-course parameters
#' within the model that have been modelled using consistency relative effects (`pool="rel"`).
#' These are given as a list, in which each named element corresponds to a time-course
#' parameter modelled in `mbnma`, specified on the corresponding scale (i.e. specified on the log scale if modelled
#' on the log scale using Ratios of Means). Their values can be either of the following:
#' * `numeric()` A numeric value representing the deterministic value of the time-course parameter in
#' question in individuals given the reference treatment. `0` is used as the default, which assumes no
#' effect of time on the reference treatment (i.e. mean differences / relative effects versus the
#' reference treatment are modeled).
#' * `formula()` A formula representing a stochastic distribution for the value of the time-course
#' parameter in question. This is specified as a random number generator (RNG) given as a formula,
#' and can take any RNG distribution for which a function exists in R. For example: `~rnorm(n, -3, 0.2)`.
#' @param synth A character object that can take the value `"common"` or `"random"` that
#' specifies the the type of pooling to use for synthesis of `ref.resp`. Using `"random"` rather
#' than `"common"` 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 ... Arguments to be sent to R2jags for synthesis of the network
#' reference treatment effect (using [ref.synth()])
#'
#'
#' @return An S3 object of class `mb.predict` that contains the following
#' elements:
#' * `summary` A named list of data frames. Each data frame contains
#' a summary of predicted responses at follow-up times specified in `times`
#' for each treatment specified in `treats`
#' * `pred.mat` 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`
#'
#' @details `ref.resp` only needs to be specified if `mbnma` has
#' been estimated using consistency relative effects (`pool="rel"`) for
#' any time-course parameters, as these inform the absolute values of the
#' network reference treatment parameters which can then be added to the
#' relative effects to calculate specific predictions.
#'
#' @examples
#' \donttest{
#' # Create an mb.network object from a dataset
#' network <- mb.network(osteopain)
#'
#' # Run an MBNMA model with an Emax time-course
#' emax <- mb.run(network,
#' fun=temax(pool.emax="rel", method.emax="common",
#' pool.et50="abs", method.et50="common"))
#'
#' # Predict responses using a stochastic baseline (E0) and a distribution for the
#' #network reference treatment
#' preds <- predict(emax, times=c(0:10),
#' E0=~rnorm(n, 7, 0.5),
#' ref.resp=list(emax=~rnorm(n, -0.5, 0.05)))
#' summary(preds)
#'
#' # Predict responses using the original dataset to estimate the network reference
#' #treatment response
#' paindata.ref <- osteopain[osteopain$treatname=="Placebo_0",]
#' preds <- predict(emax, times=c(5:15),
#' E0=10,
#' ref.resp=paindata.ref)
#' summary(preds)
#'
#' # Repeat the above prediction but using a random effects meta-analysis of the
#' #network reference treatment response
#' preds <- predict(emax, times=c(5:15),
#' E0=10,
#' ref.resp=paindata.ref,
#' synth="random")
#' summary(preds)
#' }
#'
#' @export
predict.mbnma <- function(object, times=seq(0, max(object$model.arg$jagsdata$time, na.rm=TRUE), length.out=30),
E0=0,
treats = NULL, level="treatment",
ref.resp=NULL, synth="common",
lim="cred",
...) {
######## CHECKS ########
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(object, "mbnma", add=argcheck)
checkmate::assertNumeric(times, lower=0, finite=TRUE, any.missing=FALSE, unique=TRUE,
sorted=TRUE, add=argcheck)
checkmate::assertChoice(level, choices=c("treatment", "class"), add=argcheck)
checkmate::assertChoice(synth, choices=c("random", "common"), add=argcheck)
checkmate::assertChoice(lim, choices=c("cred", "pred"), add=argcheck)
#checkmate::assertClass(treats, classes=c("numeric", "character"), null.ok=TRUE, add=argcheck)
checkmate::reportAssertions(argcheck)
# if (object$model.arg$link=="log") {
# stop("'predict()' cannot currently be used with MBNMAs modelled using ratios of means (link='log')")
# }
# Check if level="class" that class effect model was fitted
if (level=="class") {
if (length(object[["model.arg"]][["class.effect"]])==0) {
stop(crayon::red(crayon::bold("`level` has been set to `class` but no class effect models were not used")))
}
if (!isTRUE(all.equal(
object$model.arg$fun$params[object$model.arg$fun$apool %in% "rel"],
names(object$model.arg$class.effect)
))) {
stop(crayon::red(crayon::bold("To predict level='class' all relative effects must be modelled with class effects")))
}
level <- "classes"
} else if (level=="treatment") {
level <- "treatments"
}
# Check whether UME has been used and stop if criteria for it are not met
if (!(FALSE %in% object[["model.arg"]][["UME"]])) {
if (length(treats)!=2) {
stop(crayon::red(crayon::bold("UME model can only be used for prediction of direct estimates for a single treatment comparison: length(treats) must equal 2")))
}
object <- prepare.ume.predict(mbnma=object, treats=treats, level=level)
}
# Check ref.resp has been specified correctly
if ("rel" %in% object$model.arg$fun$apool) {
if (is.null(ref.resp)) {
# If ref.resp is not given then assign 0 to all rel time-course parameters
ref.resp <- list()
rels <- names(object$model.arg$fun$apool)[object$model.arg$fun$apool %in% "rel"]
for (i in seq_along(rels)) {
ref.resp[[rels[i]]] <- 0
}
# # If ref.resp is not given then assign mbnma MCMC value to all abs time-course parameters
# abs <- names(object$model.arg$fun$apool)[object$model.arg$fun$apool %in% "abs"]
# for (i in seq_along(abs)) {
# ref.resp[[abs[i]]] <- mbnma$BUGSoutput$sims.list[[abs[i]]]
# }
} else {
# If ref.resp is given ensure it is of the correct class
if (!(any(class(ref.resp) %in% c("data.frame", "tibble", "list")))) {
stop(crayon::red(crayon::bold("`object` includes time-course parameters modelled using relative effects (pool=`rel`).
The reference treatment response for them must be provided to `ref.resp` as a list,
or estimated from a dataset of reference treatment studies by providing a data frame.")))
}
}
}
# If treats have not been specified then select all of them
NT <- ifelse(level=="treatments", object$model.arg$jagsdata$NT, object$model.arg$jagsdata$Nclass)
if (is.null(treats)) {
#treats <- c(1:object[["model"]][["data"]]()[["NT"]])
treats <- object$network[[level]]
} else if (!is.null(treats)) {
if (is.numeric(treats)) {
if (any(treats > NT | any(treats<1))) {
stop(crayon::red(crayon::bold("If given as numeric treatment/class codes, `treats` must be numbered similarly to treatment/class codes in `object`")))
}
# Sort treats into order in object$network
treats <- sort(treats)
treats <- object$network[[level]][treats]
}
if (is.character(treats)) {
if (!all(treats %in% object$network[[level]])) {
stop(crayon::red(crayon::bold("`treats` includes treatments/classes not included in `object`")))
}
# Sort treats into order in object$network
treats <- object$network[[level]][sort(match(treats, object$network[[level]]))]
}
}
#### Check E0 ####
if (is.null(E0)) {
stop(crayon::red(crayon::bold("E0 has not been defined")))
}
# Check that distribution for E0 is of the correct format
if (inherits(E0,"formula")) {
E0 <- as.character(E0)[2]
if (grepl("r[A-z]+\\(n,.+\\)", E0)==FALSE) {
stop(crayon::red(crayon::bold("Stochastic distribution for E0 must be expressed as a string in the form of a supported R distribution (e.g. '~rnorm(n, 5,2)')")))
}
} else if (is.numeric(E0)) {
if (length(E0)!=1) {
stop(crayon::red(crayon::bold("`E0` can only take a single numeric value if not expressed as a stochastic distribution")))
}
} else {
stop(crayon::red(crayon::bold("'E0' has been incorrectly specified")))
}
# Ensure prediction intervals are used where appropriate
if (lim=="pred" & "random" %in% object$model.arg$a.method) {
addsd <- TRUE
if (!any(grepl("sd\\.", object$parameters.to.save))) {
stop(crayon::red("'sd' for time-course parameters not included in parameters.to.save...\n...cannot calculate prediction intervals"))
}
message("Bayesian prediction intervals to be calculated")
} else {
addsd <- FALSE
}
###### Extract info from mbnma #######
n <- object$BUGSoutput$n.sims
# Initial predict parameters
# timecourse <- list(init.predict(object)[["timecourse"]])
# beta.incl <- init.predict(object)[["beta.incl"]]
# beta.incl <- object$model.arg$fun$params
# timecourse <- paste0("alpha + ", object$model.arg$fun$jags)
# Extract parameter values from MBNMA result
model.vals <- get.model.vals(mbnma=object, E0=E0, level=level, lim=lim, link=object$model.arg$link)
timecourse <- model.vals[["timecourse"]]
time.params <- model.vals[["time.params"]]
########## Get reference treatment effect ###########
#mu.prior <- model.vals[["mu.prior"]]
mu.params <- time.params[grepl("^mu\\.", time.params)]
if (!is.null(ref.resp)) {
# If ref.resp specified as values for each time-course parameter (in a list)
if (any(class(ref.resp)=="list")) {
msg <- paste0("Priors required for: ", paste(mu.params, collapse=", "))
message(msg)
names(ref.resp) <- paste0("mu.", match(names(ref.resp), object$model.arg$fun$params))
if (identical(sort(mu.params), sort(names(ref.resp)))==FALSE) {
msg <- "Named elements of `ref.resp` do not correspond to consistency time-course parameters monitored within the model."
stop(crayon::bold(crayon::red(msg)))
} else {
message(crayon::bold(crayon::green("Success: Elements in prior match consistency time-course treatment effect parameters")))
}
# Assign ref.resp to mu values in model
for (i in seq_along(ref.resp)) {
if (class(ref.resp[[i]]) %in% c("formula", "character")) {
if (class(ref.resp[[i]]) %in% "formula") {
ref.resp[[i]] <- as.character(ref.resp[[i]])[2]
if (grepl("r[A-z]+\\(n,.+\\)", ref.resp[[i]])==FALSE) {
stop(crayon::red("Stochastic distribution for ref.resp must be expressed as a formula in the form of a supported R distribution (e.g. ~rnorm(n, 5,2))"))
}
}
assign(mu.params[which(names(ref.resp)[i]==mu.params)],
eval(parse(text=ref.resp[[i]])))
} else if (class(ref.resp[[i]]) %in% "numeric") {
assign(mu.params[which(names(ref.resp)[i]==mu.params)],
rep(ref.resp[[i]], n))
}
}
} else if (any(class(ref.resp) %in% c("data.frame", "tibble"))) {
### PLACEBO SYNTHESIS MODEL ###
args <- list(...)
synth.result <- do.call(ref.synth, args=c(list(data.ab=ref.resp, mbnma=object, synth=synth), args))
# synth.result <- ref.synth(data.ab=ref.resp, mbnma=object, synth=synth, ...)
synth.result <- synth.result$BUGSoutput$median
synth.result[["deviance"]] <- NULL
# Assign synth.result to mu values in model
for (i in seq_along(mu.params)) {
if (synth=="random") {
assign(mu.params[i],
stats::rnorm(n,
synth.result[[mu.params[i]]],
synth.result[[paste0("sd.", mu.params[i])]])
)
} else if (synth=="common") {
assign(mu.params[i], rep(synth.result[[mu.params[i]]],
n))
} else (stop(crayon::red("synth must be either `common` or `random`")))
}
}
}
# Convert predicted times to splines
if (any(c("bs", "ns", "ls") %in% object$model.arg$fun)) {
timecourse <- gsub("\\[i\\,", "[", timecourse)
spline <- genspline(times, spline=object$model.arg$fun$name, knots=object$model.arg$fun$knots, degree=object$model.arg$fun$degree,
boundaries = c(0, max(object$model.arg$jagsdata$time, na.rm=TRUE)))
}
if ("itp" %in% object$model.arg$fun$name) {
maxtime <- max(c(max(object$model.arg$jagsdata$time, na.rm=TRUE), times))
}
########## Predict responses ###########
# Assign E0 to alpha in model
#alpha <- eval(parse(text=E0)) # TO BE REMOVED
alpha <- model.vals$alpha
beta.params <- time.params[grepl("^beta.", time.params)]
# Assign single beta results to beta values in model
for (i in seq_along(beta.params)) {
if (!is.matrix(model.vals[[beta.params[i]]])) {
assign(beta.params[i], model.vals[[beta.params[i]]])
} else if (is.matrix(model.vals[[beta.params[i]]])) {
if (ncol(model.vals[[beta.params[i]]])==1) {
assign(beta.params[i], model.vals[[beta.params[i]]])
}
}
}
d.params <- time.params[grepl("^d\\.", time.params)]
predicts <- list()
treatsnum <- which(object$network[[level]] %in% treats)
for (treat in seq_along(treatsnum)) {
# Assign d results to d values in model
for (i in seq_along(d.params)) {
assign(d.params[i], model.vals[[d.params[i]]][,treatsnum[treat]])
}
treatpred <- data.frame("pred"=rep(NA,n))
for (m in seq_along(times)) {
time <- times[m]
# Evaluate function
pred <- eval(parse(text=timecourse))
if (any(is.na(pred))) {
pred[is.na(pred)] <- 0
}
#treatpred <- cbind(treatpred, pred)
# if (is.vector(pred)) {
# pred <- as.matrix(pred, ncol=1)
# }
treatpred[[paste0("time", times[m])]] <- pred
}
# predicts[[paste0(treats[treat])]] <- treatpred[,1]
predicts[[paste0(treats[treat])]] <- treatpred[-1]
}
# Generate summary data frame
sumpred <- list()
for (i in seq_along(treats)) {
summary <- data.frame("time"=times)
summary[["mean"]] <- apply(predicts[[as.character(treats[i])]], MARGIN=2,
FUN=function(x) mean(x))
summary[["sd"]] <- apply(predicts[[as.character(treats[i])]], MARGIN=2,
FUN=function(x) stats::sd(x))
quantiles <- apply(predicts[[as.character(treats[i])]], MARGIN = 2,
function(x) stats::quantile(x,
probs=c(0.025, 0.25, 0.5, 0.75, 0.975)))
summary <- cbind(summary, t(quantiles))
sumpred[[as.character(treats[i])]] <- summary
}
#predict.result <- list("summary"=sumpred, "pred.mat"=predicts, "treatments"=object$network$treatments, "mbnma"=object)
predict.result <- list("summary"=sumpred, "pred.mat"=predicts, "network"=object$network,
"times"=times, "link"=object$model.arg$link, "lim"=lim)
class(predict.result) <- "mb.predict"
return(predict.result)
}
#' Forest plot for results from time-course MBNMA models
#'
#' Generates a forest plot for time-course parameters of interest from results from time-course MBNMA models.
#' Posterior densities are plotted above each result using `ggdist:stat_:halfeye()`
#'
#' @param x An S3 object of class `"mbnma"` generated by running
#' a time-course MBNMA model
#' @param params A character vector of time-course parameters to plot.
#' Parameters must be given the same name as monitored nodes in `mbnma` and must vary by treatment or class. Can be set to
#' `NULL` to include all available time-course parameters estimated by `mbnma`.
#' @param treat.labs A character vector of treatment labels. If left as `NULL` (the default) then
#' labels will be used as defined in the data.
#' @param class.labs A character vector of class labels if `mbnma` was modelled using class effects
#' If left as `NULL` (the default) then labels will be used as defined in the data.
#' @param ... Arguments to be sent to `ggdist::stat_halfeye()`
#'
#' @return A forest plot of class `c("gg", "ggplot")` that has separate panels for different time-course parameters
#'
#' @examples
#'\donttest{
#' # Create an mb.network object from a dataset
#' alognet <- mb.network(alog_pcfb)
#'
#' # Run an MBNMA model with an Emax time-course
#' emax <- mb.run(alognet,
#' fun=temax(pool.emax="rel", method.emax="common",
#' pool.et50="rel", method.et50="common"),
#' intercept=FALSE)
#'
#' # Generate forest plot
#' plot(emax)
#'
#' # Plot results for only one time-course parameter
#' plot(emax, params="emax")
#' }
#' @export
plot.mbnma <- function(x, params=NULL, treat.labs=NULL, class.labs=NULL, ...) {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(x, "mbnma", add=argcheck)
checkmate::assertCharacter(params, null.ok=TRUE, add=argcheck)
checkmate::assertCharacter(treat.labs, null.ok=TRUE, add=argcheck)
checkmate::assertCharacter(class.labs, null.ok=TRUE, add=argcheck)
checkmate::reportAssertions(argcheck)
# Declare global variables
timeparam <- NULL
Var2 <- NULL
value <- NULL
ndistinct <- NULL
if (!FALSE %in% x$model.arg$UME) {
stop("plot.mbnma() cannot be used with UME models")
}
# Change beta to d (if present) so that it is identified in mcmc output
for (i in 1:4) {
if (paste0("beta.",i) %in% params) {
if (!paste0("beta.",i) %in% x[["parameters.to.save"]]) {
params[which(params==paste0("beta.",i))] <- paste0("d.",i)
}
}
}
# Check that specified params are monitored in model
if (!all(params %in% x[["parameters.to.save"]])) {
stop(paste0("Variable 'params': Must contain elements of set {", paste(x[["parameters.to.save"]], collapse = ", "), "}"))
}
# Check that specified params vary over treatment
for (i in seq_along(params)) {
if (length(x[["BUGSoutput"]][["median"]][[params[i]]])<=1) {
stop(paste0(params[i], " does not vary over treatment/class and cannot be plotted"))
}
}
# Add all available params if is.null(params)
if (is.null(params)) {
params <- c(x$model.arg$fun$params,
paste0("sd.", x$model.arg$fun$params),
toupper(x$model.arg$fun$params),
paste0("sd.", toupper(x$model.arg$fun$params)),
paste0("d.", 1:4),
paste0("sd.beta.", 1:4),
paste0("D.", 1:4),
paste0("sd.BETA.", 1:4)
)
params <- params[params %in% x$parameters.to.save]
# Selects only parameters which vary by class or agent (with more than 1 element)
drop <- vector()
for (i in seq_along(params)) {
if (length(x[["BUGSoutput"]][["median"]][[params[i]]])<=1) {
drop <- append(drop, i)
}
}
if (length(drop)>0) {
params <- params[-drop]
}
if (length(params)==0) {
stop("No time-course consistency parameters can be identified from the model")
}
}
# Compile parameter data into one data frame
mcmc <- x$BUGSoutput$sims.list
plotdata <- data.frame(Var2=NA, value=NA, timeparam=NA)
for (i in seq_along(params)) {
paramdata <- reshape2::melt(mcmc[[params[i]]])[,2:3]
paramdata$timeparam <- rep(params[i], nrow(paramdata))
plotdata <- rbind(plotdata, paramdata)
}
plotdata <- plotdata[-1,]
# Change param labels for treatments
treatdat <- plotdata[!grepl("[[:upper:]]", plotdata$timeparam),]
if (!is.null(treat.labs)) {
treatcodes <- treatdat$Var2
if (length(treat.labs)!=max(treatcodes)) {
stop("`treat.labs` length does not equal number of treatments that have been modelled for this time-course parameter")
} else {
t.labs <- treat.labs[sort(unique(treatcodes))]
}
} else if ("treatments" %in% names(x$network)) {
t.labs <- x$network[["treatments"]]
} else {
t.labs <- sort(unique(treatdat$param))
}
# Change param labels for classes
classdat <- plotdata[grepl("[[:upper:]]", plotdata$timeparam),]
c.labs <- vector()
if (nrow(classdat)!=0) {
if (!is.null(class.labs)) {
classcodes <- classdat$Var2
c.labs <- class.labs[classcodes]
} else if ("classes" %in% names(x$network)) {
c.labs <- x$network[["classes"]]
} else {
c.labs <- sort(unique(classdat$param))
}
}
# Increase param number (Var2) for classes
ntreat <- ifelse(nrow(treatdat)>0, max(treatdat$Var2), 0)
plotdata$Var2[grepl("[[:upper:]]", plotdata$timeparam)] <-
plotdata$Var2[grepl("[[:upper:]]", plotdata$timeparam)] + ntreat
# Attach labels
if (nrow(treatdat)>0) {
all.labs <- c(t.labs, c.labs)
} else {all.labs <- c.labs}
plotdata$Var2 <- factor(plotdata$Var2, labels=all.labs)
if (any(is.na(levels(plotdata$param)))) {
stop("`treat.labs` or `class.labs` have not been specified correctly")
}
# Drop paramers fixed to zero (i.e. network reference)
plotdata <- plotdata %>% dplyr::group_by(timeparam, Var2) %>%
dplyr::mutate(ndistinct=dplyr::n_distinct(value)) %>%
subset(ndistinct>1)
# Create forest plot
g <- ggplot2::ggplot(plotdata, ggplot2::aes(x = value, y = Var2)) +
ggdist::stat_halfeye() +
ggplot2::facet_wrap(~timeparam, scales="free")
# Axis labels
g <- g + ggplot2::xlab("Effect size") +
ggplot2::ylab("Treatment / Class") +
theme_mbnma()
g <- g + do.call(ggdist::stat_halfeye, args = list(...))
graphics::plot(g)
return(invisible(g))
}
#' Calculates relative effects/mean differences at a particular time-point
#'
#' Uses mbnma time-course parameter estimates to calculate treatment
#' differences between treatments or classes at a particular time-point.
#' Can be used to compare treatments evaluated in studies at different follow-up times, or even
#' to compare treatments in different MBNMA models via a common comparator.
#'
#' @param time A numeric value for the time at which to estimate relative effects/mean differences.
#' @param treats A character vector of treatment names for which to calculate relative effects/mean
#' differences. Must be a subset of `mbnma$network$treatments`.
#' @param classes A character vector of class names for which to calculate relative effects/mean
#' differences. Must be a subset of `mbnma$network$classes`. Only works for class effect models.
#' @param mbnma.add An S3 object of `class("mbnma")` generated by running
#' a time-course MBNMA model. This should only be specified if results from two different MBNMA models
#' are to be combined to perform a 2-stage MBNMA (see Details).
#' @inheritParams predict.mbnma
#' @inheritParams fitplot
#'
#' @return An object of class `"relative.array"` list containing:
#' * The time-point for which results are estimated
#' * Matrices of posterior means, medians, SDs and upper and lower 95% credible intervals for the
#' differences between each treatment
#' * An array containing MCMC results for the differences between all treatments specified in `treats`
#' or all classes specified in `classes`.
#'
#' Results are reported in tables as the row-defined treatment minus the column-defined treatment.
#'
#' @details
#' `get.relative()` can also be used to perform a 2-stage MBNMA that allows synthesis of results
#' from two different MBNMA models via a single common comparator.
#' In an MBNMA model, all treatments must share the same time-course function. However, a 2-stage
#' approach can enable fitting of different time-course functions to different sets ("subnetworks") of
#' treatments. For example, some treatments may have rich time-course information,
#' allowing for a more complex time-course function to be used, whereas others may be sparse,
#' requiring a simpler time-course function.
#'
#' Relative comparisons between treatments in the two datasets at specific follow-up times
#' can then be estimated from MBNMA predicted effects versus a common comparator
#' using the Bucher method and assuming consistency. See the MBNMAtime vignette for further details.
#'
#' @examples
#' \donttest{
#' # Create an mb.network object from a dataset
#' alognet <- mb.network(alog_pcfb)
#'
#' # Run a quadratic time-course MBNMA using the alogliptin dataset
#' mbnma <- mb.run(alognet,
#' fun=tpoly(degree=2,
#' pool.1="rel", method.1="random",
#' pool.2="rel", method.2="common"
#' )
#' )
#'
#' # Calculate differences between all treatments at 20 weeks follow-up
#' allres <- get.relative(mbnma, time=20)
#'
#' # Calculate difference between a subset of treatments at 10 weeks follow-up
#' subres <- get.relative(mbnma, time=10,
#' treats=c("alog_50", "alog_25", "placebo"))
#'
#'
#'
#' ###########################
#' ##### 2-stage MBNMA #####
#' ###########################
#'
#' # Using the osteoarthritis dataset
#' # With placebo (Pl_0) as common comparator between subnetworks
#'
#' #### Sparse model ####
#'
#' # Treatments on which time-course data is limited
#' sparse.trt <- c("Ce_100", "Ce_400", "Du_90", "Lu_200", "Lu_400",
#' "Lu_NA", "Et_5", "Ox_44")
#'
#' # Create a subnetwork of studies comparing these treatments
#' sparse.df <- osteopain %>% dplyr::group_by(studyID) %>%
#' dplyr::filter(any(treatment %in% sparse.trt)) %>%
#' dplyr::ungroup() %>%
#' subset(treatment %in% c("Pl_0", sparse.trt))
#'
#' sparse.net <- mb.network(sparse.df)
#'
#' # Run a ITP MBNMA with a known rate
#' sparse.mbnma <- mb.run(sparse.net, fun=titp(method.rate=0.8, pool.rate="abs"))
#'
#'
#' #### Complex model ####
#'
#' # Treatments on which time-course data is rich
#' rich.trt <- levels(osteopain$treatment)[!levels(osteopain$treatment) %in%
#' c("Pl_0", "Ce_100", "Ce_400", "Du_90", "Lu_200",
#' "Lu_400", "Lu_NA", "Et_5", "Ox_44")]
#'
#' # Create a subnetwork of studies comparing these treatments
#' rich.df <- osteopain %>% dplyr::group_by(studyID) %>%
#' dplyr::filter(any(treatment %in% rich.trt)) %>%
#' dplyr::ungroup() %>%
#' subset(treatment %in% c("Pl_0", rich.trt))
#'
#' rich.net <- mb.network(rich.df)
#'
#' # Run a Emax MBNMA
#' rich.mbnma <- mb.run(rich.net, temax(p.expon = FALSE))
#'
#'
#' #### Calculate relative effects between models ####
#'
#' # At 10 weeks follow-up
#' rels.sparse <- get.relative(sparse.mbnma, time=10)
#' rels.rich <- get.relative(rich.mbnma, time=10)
#'
#' rels.all <- get.relative(mbnma=rich.mbnma,
#' mbnma.add=sparse.mbnma, time=10)
#'
#' print(rels.all$median)
#'
#'
#' }
#' @export
get.relative <- function(mbnma, mbnma.add=NULL, time=max(mbnma$model.arg$jagsdata$time, na.rm=TRUE),
treats=unique(c(mbnma$network$treatments, mbnma.add$network$treatments)),
classes=NULL, lim="cred") {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertNumeric(time, lower=0, len=1, null.ok=FALSE, add=argcheck)
checkmate::assertSubset(treats, choices=c(mbnma$network$treatments,
mbnma.add$network$treatments),
add=argcheck)
checkmate::assertSubset(classes, choices=c(mbnma$network$classes,
mbnma.add$network$classes),
add=argcheck)
checkmate::assertClass(mbnma, classes="mbnma", add=argcheck)
checkmate::assertClass(mbnma.add, classes="mbnma", null.ok = TRUE, add=argcheck)
checkmate::reportAssertions(argcheck)
if (!is.null(classes)) {
treats <- classes
level <- "class"
levels <- "classes"
} else {
level <- "treatment"
levels <- "treatments"
}
if (!FALSE %in% mbnma$model.arg$UME) {
stop("get.relative() cannot be used with UME models\nEstimation of relative effects uses consistency assumption")
}
if (is.null(mbnma.add)) {
pred <- suppressMessages(predict.mbnma(mbnma, times=time,
treats = treats, level = level, lim=lim))
# Matrix of results
mat <- do.call(cbind, pred$pred.mat)
} else if (!is.null(mbnma.add)) {
if (!FALSE %in% mbnma.add$model.arg$UME) {
stop("get.relative() cannot be used with UME models\nEstimation of relative effects uses consistency assumption")
}
# Identify common treatment to be the reference
match1 <- which(treats %in% mbnma$network[[levels]])
match2 <- which(treats %in% mbnma.add$network[[levels]])
ind1 <- which(match1 %in% match2)
ind2 <- which(match2 %in% match1)
ref <- treats[match1[ind1]]
if (length(ind1)!=1 | length(ind2)!=1) {
stop("mbnma and mbnma.add must have a single treatment/class in common in order to perform 2-stage MBNMA")
}
# Split treatments/classes into subnetworks
treats2 <- c(ref, treats[!treats %in% mbnma$network[[levels]]])
treats1 <- c(ref, treats[!treats %in% mbnma.add$network[[levels]]])
# Predict relative effects at specific time
pred <- suppressMessages(predict.mbnma(mbnma, times=time,
treats = treats1, level = level, lim=lim))
pred.add <- suppressMessages(predict.mbnma(mbnma.add, times=time,
treats = treats2, level = level, lim=lim))
# Matrix of results
mat1 <- do.call(cbind, pred$pred.mat)
mat2 <- do.call(cbind, pred.add$pred.mat)
# Ensure matrices are same size
niter <- min(c(nrow(mat1), c(nrow(mat2))))
if (nrow(mat1)>nrow(mat2)) {
mat1 <- mat1[sample(1:nrow(mat1), nrow(mat2), replace = FALSE),]
} else if (nrow(mat1)<nrow(mat2)) {
mat2 <- mat2[sample(1:nrow(mat2), nrow(mat1), replace = FALSE),]
}
# Rearrange matrix to common reference treatment
if (which(names(pred$summary) %in% ref)!=1) {
mat1 <- mat1 - mat1[,ind1]
# Rearrange columns
refcol <- which(names(pred$summary) %in% ref)
mat1 <- cbind(mat1[,refcol], mat1[,-refcol])
}
if (which(names(pred.add$summary) %in% ref)!=1) {
mat2 <- mat2 - mat2[,ind2]
# Rearrange columns
refcol <- which(names(pred.add$summary) %in% ref)
mat2 <- cbind(mat2[,refcol], mat2[,-refcol])
}
# Drop shared reference from mat2
mat2 <- mat2[,-which(treats2 %in% ref)]
# Combine results
mat <<- cbind(mat1, mat2)
treats <- c(treats1, treats2[-which(treats2 %in% ref)])
}
# For lower triangle
outmat <- array(dim=c(length(treats), length(treats), nrow(mat)))
for (i in 1:(ncol(mat)-1)) {
temp <- mat[,-i, drop=FALSE] - mat[,i]
#temp <- apply(temp, MARGIN=2, FUN=function(x) {neatCrI(quantile(x, probs=c(0.025, 0.5, 0.975)), digits = 2)})
outmat[(1+i):dim(outmat)[1],i,] <- t(temp[,i:ncol(temp)])
}
# For upper triangle
for (i in 1:(ncol(mat)-1)) {
temp <- mat[,i] - mat[,-i, drop=FALSE]
outmat[i,(1+i):dim(outmat)[2],] <- t(temp[,i:ncol(temp)])
}
dimnames(outmat)[[1]] <- treats
dimnames(outmat)[[2]] <- treats
######### Summary matrixes ######
xmat <- outmat
meanmat <- matrix(nrow=nrow(xmat), ncol=ncol(xmat))
semat <- meanmat
medmat <- meanmat
l95mat <- medmat
u95mat <- medmat
sum.df <- data.frame(comparison=NA, mean=NA, sd=NA, l95=NA, med=NA, u95=NA)
for (i in 1:nrow(xmat)) {
for (k in 1:ncol(xmat)) {
if (!is.na(xmat[i,k,1])) {
meanmat[i,k] <- mean(xmat[i,k,])
semat[i,k] <- stats::sd(xmat[i,k,])
medmat[i,k] <- stats::median(xmat[i,k,])
l95mat[i,k] <- stats::quantile(xmat[i,k,], probs = 0.025)
u95mat[i,k] <- stats::quantile(xmat[i,k,], probs = 0.975)
if (k>i) {
sum.df <- dplyr::add_row(sum.df,
comparison = paste0(treats[i], " vs ", treats[k]),
mean=meanmat[i,k],
sd=semat[i,k],
l95=l95mat[i,k],
med=medmat[i,k],
u95=u95mat[i,k]
)
}
}
}
}
sum.df <- sum.df[-1,]
names(sum.df) <- c("comparison", "mean", "sd", "2.5%", "50%", "97.5%")
sumlist <- list("mean"=meanmat, "sd"=semat, "median"=medmat, "lower95"=l95mat, "upper95"=u95mat)
for (i in seq_along(sumlist)) {
dimnames(sumlist[[i]])[[1]] <- dimnames(xmat)[[1]]
dimnames(sumlist[[i]])[[2]] <- dimnames(xmat)[[2]]
}
out <- list("time"=time, "data.frame"=sum.df, "relarray"=outmat)
out <- c(out, sumlist, lim=lim)
class(out) <- "relative.array"
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.