Nothing
######################################
#### 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, ...)
}
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.