Nothing
# Graphical functions
# Author: Hugo Pedder
# Date created: 2018-09-10
#' Get large vector of distinct colours using Rcolorbrewer
genmaxcols <- function() {
cols1 <- RColorBrewer::brewer.pal(9, "Set1")
cols2 <- RColorBrewer::brewer.pal(9, "Pastel1")
cols <- c(rbind(cols1, cols2))
cols1 <- RColorBrewer::brewer.pal(8, "Set2")
cols2 <- RColorBrewer::brewer.pal(8, "Pastel2")
cols <- c(cols, c(rbind(cols1, cols2)))
cols <- c(cols, RColorBrewer::brewer.pal(12, "Set3"))
return(cols)
}
#' Calculate position of label with respect to vertex location within a circle
#'
#' Useful for graphs drawn using `igraph` to reposition labels relative to vertices when vertices
#' are laid out in a circle (as is common in network plots). `igraph` interprets position within
#' `vertex.label.degree` as radians, so it is necessary to convert locations into radian values. This
#' is the main role of this function.
#'
#' @param x A numeric vector of positions around a circle, typically sequentially numbered.
#' @param start A number giving the offset from 12 o'clock in radians for the label locations.
#' @param direction Either `1` for clockwise numbering (based on the order of `x`) or `-1` for
#' anti-clockwise.
#'
#' @return A numeric vector of rescaled values
#'
#' @references
#' https://gist.github.com/kjhealy/834774/a4e677401fd6e4c319135dabeaf9894393f9392c
radian.rescale <- function(x, start=0, direction=1) {
c.rotate <- function(x) (x + start) %% (2 * pi) * direction
c.rotate(scales::rescale(x, c(0, 2 * pi), range(x)))
}
#' Overlays observations as shaded regions on a time-course plot
#'
#' @param g An object of `c("gg", "ggplot")` that is a plot of predicted responses from an MBNMA model.
#' @inheritParams plot.mb.predict
#'
#' @noRd
disp.obs <- function(g, predict, col="blue", max.col.scale=NULL) {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(g, c("gg", "ggplot"), add=argcheck)
#checkmate::assertDataFrame(data.ab, add=argcheck)
checkmate::assertClass(predict, "mb.predict", add=argcheck)
checkmate::assertInt(max.col.scale, lower=1, null.ok = TRUE, add=argcheck)
checkmate::reportAssertions(argcheck)
#treats <- as.numeric(names(predict[["summary"]]))
#treats <- which(predict$mbnma$network$treatments %in% names(predict$summary))
treats <- names(predict[["summary"]])
#data.ab <- jagstonetwork(predict$mbnma)
data.ab <- predict$network$data.ab
#raw.data <- network[["data.ab"]]
raw.data <- dplyr::arrange(data.ab, studyID, fupcount, arm)
raw.data$studyID <- as.character(raw.data$studyID)
predict.data <- predict[["summary"]][[1]]
predict.data[["treat"]] <- NA
predict.data[["count"]] <- NA
predict.data <- predict.data[0,]
for (i in seq_along(predict[["summary"]])) {
temp <- predict[["summary"]][[i]]
# temp[["treat"]] <- rep(as.numeric(names(predict[["summary"]][i])),
# nrow(temp))
# temp[["treat"]] <- rep(which(predict$mbnma$network$treatments %in% names(predict$summary)[i]),
# nrow(temp))
temp[["treat"]] <- rep(names(predict$summary)[i],
nrow(temp))
temp[["count"]] <- rep(NA, nrow(temp))
dplyr::arrange(temp, time)
if (temp[["time"]][1]>0) {
temp[["count"]][1] <- nrow(raw.data[raw.data$treatment==which(predict$network$treatments %in% temp[["treat"]][1]) &
raw.data[["time"]]<=temp[["time"]][1] &
raw.data[["time"]]>0
,])
} else {temp[["count"]][1] <- 0 }
for (k in 2:nrow(temp)) {
temp[["count"]][k] <- nrow(raw.data[raw.data$treatment==which(predict$network$treatments %in% temp[["treat"]][k]) &
raw.data[["time"]]<=temp[["time"]][k] &
raw.data[["time"]]>temp[["time"]][k-1]
,])
}
predict.data <- rbind(predict.data, temp)
}
# Check max.col.scale
n.cut <- max(predict.data$count)
if (!is.null(max.col.scale)) {
if (!is.numeric(max.col.scale)) {
stop("max.col.scale must be a number greater than the maximum number of observed counts in the plotted treatments")
}
if (n.cut > max.col.scale) {
stop("max.col.scale must be a number greater than the maximum number of observed counts in the plotted treatments")
}
n.cut <- max.col.scale
}
# Generate colours
#cols <- col.scale(n.cut=max(predict.data$count), col=col)
cols <- alpha.scale(n.cut=n.cut, col=col)
for (treat in seq_along(treats)) {
subset <- predict.data[predict.data$treat==treats[treat],]
# Start assuming lowest time = 0
for (m in 2:nrow(subset)) {
g <- g + ggplot2::geom_ribbon(data=subset[subset$time<=subset$time[m] &
subset$time>=subset$time[m-1]
,],
ggplot2::aes(x=time, ymin=`2.5%`, ymax=`97.5%`),
fill=cols[subset$count[m]+1])
}
}
return(g)
}
#' Generates a scale of colours
#'
#' FUNCTION IS DEPRECATED
#' @noRd
col.scale <- function(n.cut, col="blue",
rgb.start=NULL, rgb.end=NULL) {
# Set colour intensities
if (col=="blue") {
max.rgb <- c(0,0,200)
min.rgb <- c(200,200,250)
} else if (col=="red") {
max.rgb <- c(120,0,0)
min.rgb <- c(250,200,200)
} else if (col=="green") {
max.rgb <- c(0,100,0)
min.rgb <- c(200,250,200)
}
if (!is.null(rgb.start) & !is.null(rgb.end)) {
max.rgb <- rgb.start
min.rgb <- rgb.end
}
#min.rgb <- rgb.end
cut <- (min.rgb-max.rgb)/(n.cut-1)
#rgb <- c(255,255,255)
col.mat <- c(255,255,255)
rgb <- min.rgb
col.mat <- cbind(col.mat, rgb)
for (i in 1:n.cut) {
rgb <- rgb - cut
col.mat <- cbind(col.mat, abs(rgb))
}
#col.mat <- rbind(col.mat, c(0, rep(255, (n.cut+1))))
col.mat <- rbind(col.mat, c(0, rep(255, (ncol(col.mat)-1))))
hexcol <- vector()
for (i in 1:(n.cut)) {
hexcol <- append(hexcol, rgb(col.mat[1,i], col.mat[2,i],
col.mat[3,i], col.mat[4,i], maxColorValue=255
))
}
return(hexcol)
}
#' Generates colours with varying degrees of transparency
#'
#' Identical to MBNMAdose
#'
#' @param ncut A number indicating the number of different counts in the dataset
#' @param col Colour to use for shading
#'
#' @noRd
alpha.scale <- function(n.cut, col="blue") {
# Run checks
checkmate::assertIntegerish(n.cut, lower=1, len=1)
# Set colour intensities
if (is.character(col)) {
if (col=="blue") {
hue <- c(0,0,200)
} else if (col=="red") {
hue <- c(120,0,0)
} else if (col=="green") {
hue <- c(0,100,0)
} else {
stop("Permitted colours are `blue`, `red`, `green`, or an RGB code")
}
} else if (is.numeric(col)) {
if (length(col)!=3) {
stop("Specified RGB code must have length 3")
}
if (any(col>255) | any(col<0)) {
stop("RGB values must be between 0 and 255")
}
hue <- col
}
#min.rgb <- rgb.end
cut <- (255-0)/(n.cut)
alpha <- 0
alpha.vec <- alpha
for (i in 1:n.cut) {
alpha <- alpha+cut
alpha.vec <- append(alpha.vec, alpha)
}
hexcol <- vector()
for (i in 1:(n.cut+1)) {
hexcol <- append(hexcol, grDevices::rgb(hue[1], hue[2],
hue[3], alpha.vec[i], maxColorValue=255
))
}
return(hexcol)
}
#' Predict "lumped" NMA predictions for plotting
#'
#' @param pred An object of `class("mb.predict")`
#' @param timebins A numeric vector defining the boundaries of the time bins within which to perform
#' a standard NMA. Length must be >=2. See details.
#' @inheritParams plot.mb.predict
#' @inheritParams mb.run
#' @inheritParams mb.predict
#'
#' @details `timebins` indicate regions of the data (defined as "time bins") over which it may be reasonable to "lump" different
#' follow-up times from different studies together and assume a standard NMA model.
#'
#' @noRd
overlay.nma <- function(pred, timebins, method="common", link="identity", lim="cred",
plottype="pred", ...) {
# Run checks
argcheck <- checkmate::makeAssertCollection()
if (plottype=="pred") {
checkmate::assertClass(pred, c("mb.predict"), add=argcheck)
}
checkmate::assertNumeric(timebins, lower=0, null.ok=FALSE, add=argcheck)
checkmate::assertChoice(lim, choices=c("cred", "pred"), add=argcheck)
checkmate::assertChoice(plottype, choices=c("rel", "pred"), add=argcheck)
checkmate::reportAssertions(argcheck)
# Declare global variable
cor <- NULL
dif <- NULL
network <- pred$network
outlist <- list()
# Loop across timebins
for (bin in 2:length(timebins)) {
skip <- FALSE
incl.range <- c(timebins[bin-1], timebins[bin])
str.incl.range <- paste(c("time=", incl.range[1], " and time=", incl.range[2]), collapse="")
if (plottype=="pred") {
# Predict NMA results at specific time point
timedif <- abs(pred$times - mean(incl.range))
if (!any(timedif < (mean(incl.range) - incl.range[1]))) {
warning(paste(c("No time-point in predicted values is between ", str.incl.range,
". Cannot overlay.nma for this time bin."), collapse=""))
skip <- TRUE
}
timeindex <- order(timedif)[1]
predref <- pred$pred.mat[[1]][[timeindex]]
} else if (plottype=="rel") {
timeindex <- 1
predref <- 0
}
if (skip==FALSE) {
nmanet <- network$data.ab[network$data.ab$time>incl.range[1] & network$data.ab$time<=incl.range[2],]
# Take the follow-up closes to mean(incl.range) if multiple fups are within the range
nmanet <- nmanet %>% dplyr::group_by(studyID) %>%
#dplyr::mutate(dif=time-mean(incl.range)) %>%
dplyr::mutate(dif=time-ifelse(plottype=="pred", pred$times[timeindex], mean(incl.range))) %>%
dplyr::arrange(dif) %>%
dplyr::group_by(studyID, arm) %>%
dplyr::slice_head()
if (!1 %in% nmanet$treatment) {
message(paste(c("overlay.nma not possible between ", str.incl.range,
". Data for network reference (treatment=1) not available in this time bin"),
collapse = ""))
skip <- TRUE
}
}
if (skip==FALSE) {
# Store max range of data over which NMA is performed
if (plottype=="pred") {
dat.range <- range(c(nmanet$time, pred$times[timeindex]))
} else if (plottype=="rel") {
dat.range <- range(c(nmanet$time))
}
# Check network connectedness
comparisons <- mb.comparisons(nmanet)
nodes <- unique(sort(nmanet$treatment))
nodes <- network$treatments[nodes]
g <- igraph::graph.empty()
g <- g + igraph::vertex(nodes)
ed <- t(matrix(c(comparisons[["t1"]], comparisons[["t2"]]), ncol = 2))
ed <- factor(as.vector(ed), labels=nodes)
edges <- igraph::edges(ed, weight = comparisons[["nr"]], arrow.mode=0)
g <- g + edges
discon <- suppressWarnings(check.network(g))
# Drop treatments disconnected to network reference
if (length(discon)>0) {
# message(paste("The following treatments are disconnected from the network reference for studies in 'incl.range' and will be excluded from overlay.nma:",
# paste(discon, collapse = "\n"), sep="\n"))
drops <- which(network$treatments %in% discon)
nmanet <- nmanet[!nmanet$treatment %in% drops,]
nodes <- unique(sort(nmanet$treatment))
nodes <- network$treatments[nodes]
}
# Recode treatments from 1-X
nmanet$treatment <- as.numeric(factor(nmanet$treatment))
# Print message
message(paste0("Running overlay.nma for ", str.incl.range))
# Run model (incl write priors)
nma <- nma.run(data.ab=nmanet, method=method, link=link,
...)
if (method=="common" | "cred" %in% lim) {
predtrt <- sample(predref, size=nma$BUGSoutput$n.sims, replace=TRUE) + nma$BUGSoutput$sims.list$d[,-1]
} else if (method=="random" & "pred" %in% lim) {
if (is.vector(nma$BUGSoutput$sims.list$d[,-1])) {
cols <- 1
} else {
cols <- ncol(nma$BUGSoutput$sims.list$d[,-1])
}
mat <- array(dim=c(nma$BUGSoutput$n.sims,
cols,
2)
)
mat[,,1] <- nma$BUGSoutput$sims.list$d[,-1]
mat[,,2] <- nma$BUGSoutput$sims.list$sd
predtrt <- apply(mat, MARGIN=c(1,2), FUN=function(x) stats::rnorm(1, x[1], x[2]))
predtrt <- sample(predref, size=nma$BUGSoutput$n.sims, replace=TRUE) + predtrt
}
if (is.vector(predtrt)) {
predtrt <- matrix(predtrt, ncol=1)
}
colnames(predtrt) <- nodes[-1]
pred.df <- as.data.frame(t(apply(predtrt, MARGIN=2, FUN=function(x){stats::quantile(x, probs = c(0.025,0.5,0.975))})))
pred.df$time <- pred$times[timeindex]
pred.df$treat <- rownames(pred.df)
# Select only treatments included in pred
if (plottype=="pred") {
pred.df <- pred.df[pred.df$treat %in% names(pred$summary),]
}
# Add range of data to pred.df
pred.df$tmin <- dat.range[1]
pred.df$tmax <- dat.range[2]
# Create width tolerance for plotting
if (dat.range[1]==dat.range[2]) {
scale <- range(timebins)
width <- (scale[2] - scale[1])/200
pred.df$tmin <- pred.df$tmin - width
pred.df$tmax <- pred.df$tmax + width
}
if (method=="random") {
sd <- nma$BUGSoutput$summary[rownames(nma$BUGSoutput$summary)=="sd",]
sd <- round(sd, 3)
sd <- paste(sd[5], " (", sd[3], ", ", sd[7], ")", sep="")
} else {
sd <- NULL
}
listnam <- paste0("bin", bin-1)
outlist[[listnam]] <- list(pred.df=pred.df,
nma.summary=nma$BUGSoutput$summary,
totresdev=round(nma$BUGSoutput$median$totresdev,1), ndat=nrow(nmanet),
dic=round(nma$BUGSoutput$median$totresdev+nma$BUGSoutput$pD,1),
sd=sd,
incl.range=dat.range)
}
}
return(outlist)
}
#' Check if all nodes in the network are connected (identical to MBNMAdose)
#' @noRd
check.network <- function(g, reference=1) {
connects <- is.finite(igraph::shortest.paths(igraph::as.undirected(g),
to=reference))
treats <- rownames(connects)[connects==FALSE]
if (length(treats>0)) {
warning(paste0("The following treatments/agents are not connected to the network reference:\n",
paste(treats, collapse = "\n")))
}
return(treats)
}
#' Get data for nodesplit plots from mb.nodesplit
#' @noRd
nodesplit.plotdata <- function(nodesplit, type) {
type <- paste0(type, ".plot")
plotdata <- data.frame()
time.param <- vector()
comparison <- vector()
for (param in seq_along(nodesplit[[1]])) {
for (split in seq_along(nodesplit)) {
temp <- nodesplit[[split]][[param]][[type]]$data
temp$time.param <- rep(names(nodesplit[[1]])[param], nrow(temp))
temp$comparison <- rep(names(nodesplit)[split], nrow(temp))
plotdata <- rbind(plotdata, temp)
}
}
return(plotdata)
}
#' Plot raw responses over time by treatment or class
#'
#' @param plotby A character object that can take either `"arm"` to indicate that raw responses
#' should be plotted separately for each study arm, or `"rel"` to indicate that the within-study
#' relative effects/treatment differences should be plotted. In this way the time-course of both the absolute
#' effects and the relative effects can be examined.
#' @param ... Arguments to be sent to `ggplot()`
#' @inheritParams plot.mb.network
#' @inheritParams mb.run
#'
#' @return The function returns an object of `class(c("gg", "ggplot")`. Characteristics
#' of the object can therefore be amended as with other plots generated by `ggplot()`.
#'
#' @details Plots can be faceted by either treatment (`level="treatment"`) or class
#' (`level="class"`) to investigate similarity of treatment responses within classes/agents.
#' Points represent observed responses and lines connect between observations within the
#' same study and arm.
#'
#' @examples
#' \donttest{
#' # Make network
#' goutnet <- mb.network(goutSUA_CFB)
#'
#' # Use timeplot to plot responses grouped by treatment
#' timeplot(goutnet)
#'
#' # Use timeplot ot plot resposes grouped by class
#' timeplot(goutnet, level="class")
#'
#' # Plot matrix of relative effects
#' timeplot(goutnet, level="class", plotby="rel")
#'
#' # Plot using Standardised Mean Differences
#' copdnet <- mb.network(copd)
#' timeplot(copdnet, plotby="rel", link="smd")
#'
#' }
#'
#' @export
timeplot <- function(network, level="treatment",
plotby="arm", link="identity",
...) {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(network, "mb.network", add=argcheck)
checkmate::assertChoice(level, choices = c("treatment", "class"), add=argcheck)
checkmate::assertChoice(plotby, choices = c("arm", "rel"), add=argcheck)
checkmate::assertChoice(link, choices = c("identity", "smd", "log"), add=argcheck)
checkmate::reportAssertions(argcheck)
# Define global variables
poolvar <- NULL
Rx.Name.x <- NULL
Rx.Name.y <- NULL
if (level=="class" & !("classes" %in% names(network))) {
stop("`level` cannot be set to class if there is no class variable included in the dataset/network")
}
plotdata <- network$data.ab
cfb <- network$cfb
if (plotby=="arm") {
studies <- unique(plotdata$studyID)
for (i in seq_along(studies)) {
base.dat <- subset(plotdata, fupcount==1 & studyID==studies[i])
if (cfb[i]==TRUE & any(base.dat$time!=0)) {
base.dat$y <- rep(0, nrow(base.dat))
base.dat$se <- rep(0, nrow(base.dat))
base.dat$time <- rep(0, nrow(base.dat))
plotdata <- rbind(base.dat, plotdata)
}
}
plotdata <- dplyr::arrange(plotdata, studyID, time, treatment)
# Do not run this function with pylr loaded!!
plotdata <- plotdata %>%
dplyr::group_by(studyID, time) %>%
dplyr::mutate(arm = sequence(dplyr::n()))
g <- ggplot2::ggplot(plotdata,
ggplot2::aes(x=plotdata$time, y=plotdata$y,
group=(paste(plotdata$studyID, plotdata$arm, sep="_")))) +
ggplot2::geom_line() +
ggplot2::geom_point(size=1)
if (level=="treatment") {
g <- g + ggplot2::facet_wrap(~factor(.data$treatment, labels=network$treatments))
#g <- g + ggplot2::facet_wrap(ggplot2::vars(treatment))
} else if (level=="class") {
g <- g + ggplot2::facet_wrap(~factor(.data$class, labels=network$classes))
}
g <- g + ggplot2::ylab("Response")
} else if (plotby=="rel") {
diffs <- plotdata %>%
dplyr::mutate(Rx.Name = factor(network$treatments[.data$treatment], levels=network$treatments))
if (link=="identity") {
diffs <- diffs %>%
dplyr::inner_join(diffs, by=c("studyID", "time")) %>%
dplyr::filter(.data$treatment.x < .data$treatment.y) %>%
dplyr::mutate(pairDiff = .data$y.y - .data$y.x)
} else if (link=="log") {
diffs <- diffs %>%
dplyr::inner_join(diffs, by=c("studyID", "time")) %>%
dplyr::filter(.data$treatment.x < .data$treatment.y) %>%
dplyr::mutate(pairDiff = log(.data$y.y/.data$y.x))
} else if (link=="smd") {
if (!"n" %in% names(diffs)) {
stop("'n' must be included in dataset if using link='smd'")
}
diffs <- diffs %>%
dplyr::inner_join(diffs, by=c("studyID", "time")) %>%
dplyr::filter(.data$treatment.x < .data$treatment.y) %>%
dplyr::mutate(
var.y = (.data$se.y * (.data$n.y)^0.5)^2,
var.x = (.data$se.x * (.data$n.x)^0.5)^2,
poolvar = ((.data$var.y * (.data$n.y-1)) + (.data$var.x * (.data$n.x-1))) / (.data$n.y + .data$n.x-2),
pairDiff = (.data$y.y - .data$y.x) / (poolvar^0.5)
)
}
# Handle cfb data
studies <- unique(diffs$studyID)
diffs$bl <- FALSE
for (i in seq_along(cfb)) {
if (cfb[i]==FALSE) {
diffs$bl[diffs$studyID==studies[i] & diffs$fupcount.x==1 & diffs$fupcount.y==1] <- TRUE
}
}
diffs <- diffs %>% dplyr::bind_rows(diffs %>% dplyr::group_by(.data$studyID,
.data$arm.x, .data$arm.y) %>%
dplyr::slice(1) %>%
dplyr::mutate(time=dplyr::case_when(bl==0 ~ 0),
pairDiff=dplyr::case_when(bl==0 ~ 0),
bl=dplyr::case_when(bl==0 ~ 1)))
diffs <- diffs[!is.na(diffs$time),]
if (level=="class") {
diffs <- diffs %>%
dplyr::mutate(Rx.Name.x = network$classkey$class[match(Rx.Name.x, network$classkey$treatment)],
Rx.Name.y = network$classkey$class[match(Rx.Name.y, network$classkey$treatment)]
)
}
g <- ggplot2::ggplot(data=diffs, ggplot2::aes(x=.data$time, y=.data$pairDiff,
# group=.data$studyID),
group=paste(.data$studyID, .data$arm.x, .data$arm.y, sep="_")),
...) +
ggplot2::geom_line() +
ggplot2::geom_point() +
ggplot2::facet_grid(rows=ggplot2::vars(.data$Rx.Name.y), cols=ggplot2::vars(.data$Rx.Name.x))
if (link=="identity") {
g <- g + ggplot2::ylab("Response")
} else if (link=="log") {
g <- g + ggplot2::ylab("Log-Ratio of Means")
} else if (link=="smd") {
g <- g + ggplot2::ylab("Standardised Mean Difference")
}
}
g <- g + ggplot2::xlab("Time") +
theme_mbnma()
graphics::plot(g)
return(invisible(g))
}
#' Plot relative effects from NMAs performed at multiple time-bins
#'
#' @param legend `TRUE`/`FALSE` to indicate whether a legend should be plotted.
#' @inheritParams mb.run
#' @inheritParams plot.mb.predict
#' @inheritParams predict.mbnma
#'
#' @return Plots treatment effects from NMAs performed within discrete time bins. The
#' object returned is a list containing the plot and a sublist of NMA results and
#' predictions from each time bin specified in `overlay.nma`.
#'
#' @details
#' Performs several standard NMAs at different time "bins", time periods within
#' which treatment effects are assumed to be constant over time. Separate NMAs
#' are then performed within each time bin on data points from studies that fall
#' within the time bin (only a single follow-up time is taken from each study
#' to avoid double counting).
#'
#' Note that the wider the time bin boundaries specified by the user, the
#' larger the potential range of included follow-up times and this can
#' introduce heterogeneity or inconsistency.
#'
#' Results are plotted versus the network reference and are plotted on the
#' specified link scale. Each time bin window is marked on the plot by
#' vertical dashed lines. The NMA estimates within each time bin are plotted
#' as a horizontal solid black line (the posterior median) with a shaded region
#' indicating the 95% credible interval (prediction intervals can instead
#' be plotted). The width of these shaded regions is equal to the range of study
#' time-points included in the NMA performed within that timebin, which
#' may therefore be more narrow than the time bin specified in the `binplot()`
#' command due to the follow-up times at which data is available in included
#' studies.
#'
#' @inheritSection plot.mb.predict Overlaying NMA results
#'
#' @examples
#' \donttest{
#' # Create an mb.network object from a dataset
#' alognet <- mb.network(alog_pcfb)
#'
#' # Plot relative effects from NMAs calculated for a single time-bins
#' # Do not plot time-bin boundaries
#' binplot(alognet, overlay.nma=c(0,5), plot.bins=FALSE)
#'
#' # Plot relative effects from NMAs at multiple time-bins
#' # With random treatment effects
#' binplot(alognet, overlay.nma=c(5,10,15,20),
#' method="random")
#' }
#'
#' @export
binplot <- function(network, overlay.nma=c(0, stats::quantile(network$data.ab$time)),
method="common", link="identity", lim="cred", plot.bins=TRUE,
legend=TRUE, ...) {
checkmate::assertNumeric(overlay.nma, min.len = 2)
# Ensure timebins are unique
overlay.nma <- unique(overlay.nma)
# Create pred to match overlay.nma
pred <- list()
pred$times <- seq(0,max(overlay.nma), length.out=100)
pred$pred.mat <- list(matrix(0, ncol=100, nrow=1))
pred$network <- network
class(pred) <- "mb.predict"
nma <- overlay.nma(pred, timebins=overlay.nma, method=method, link=link, lim=lim,
plottype="rel", ...)
if (length(nma)==0) {
stop("No NMA can be performed between time points specified in overlay.nma")
}
# Bind data frames
plot.df <- nma[[1]]$pred.df[0,]
for (i in seq_along(nma)) {
plot.df <- rbind(plot.df, nma[[i]]$pred.df)
}
g <- ggplot2::ggplot(data=plot.df, ggplot2::aes(x=tmin, xend=tmax, ymin=`2.5%`, ymax=`97.5%`, y=`50%`)) +
ggplot2::geom_rect(ggplot2::aes(ymin=`2.5%`, ymax=`97.5%`, xmin=tmin, xmax=tmax,
fill="95% Interval"),
alpha=1) +
ggplot2::geom_segment(ggplot2::aes(y=`50%`, yend=`50%`, x=tmin, xend=tmax, color="Posterior median"),
linewidth=0.8) +
ggplot2::geom_hline(ggplot2::aes(yintercept=0), linetype="dotted")
colorvals <- c("Posterior median"="gray0")
capt <- paste0("Results versus ", network$treatments[1],
"\nWidth of shaded region denotes range of follow-up times included in data for each NMA",
"\nHeight of shaded region denotes the 95% interval of the posterior distribution")
if (plot.bins==TRUE) {
capt <- paste0(capt, "\nVertical dashed lines indicate time bin boundaries")
g <- g + ggplot2::geom_vline(xintercept=overlay.nma, linetype="dashed", alpha=0.5) +
ggplot2::scale_x_continuous(breaks=unique(c(0, overlay.nma)))
}
g <- g +
ggplot2::labs(caption=capt) +
ggplot2::scale_fill_manual(name="", values=c("95% Interval"="lightblue"))
g <- g + ggplot2::facet_wrap(~factor(treat)) +
ggplot2::labs(y="Treatment effect (on link scale)", x="Time") +
ggplot2::scale_color_manual(name="", values=colorvals) +
theme_mbnma() +
ggplot2::theme(legend.position="none")
if (legend==TRUE) {
# legend.png <- png::readPNG("man/figures/binplot_legend.PNG")
legend.png <- png::readPNG(system.file("figures/binplot_legend.PNG", package = "MBNMAtime"))
g2 <- grid::rasterGrob(legend.png, interpolate = TRUE)
g <- gridExtra::arrangeGrob(g, g2, ncol=2, widths=c(10, 1))
}
graphics::plot(g)
out <- list("graph"=g)
out[["overlay.nma"]] <- nma
return(invisible(out))
}
#' Plot deviance contributions from an MBNMA model
#'
#' @param mbnma An S3 object of class `"mbnma"` generated by running
#' a time-course MBNMA model
#' @param dev.type Deviances to plot - can be either residual deviances (`"resdev"`) or deviances (`"dev"`, the default)
#' @param plot.type Deviances can be plotted either as scatter points (`"scatter"` - using
#' `geom_point()`) or as boxplots (`"box"`, the default)
#' @param xaxis A character object that indicates whether deviance contributions should be plotted
#' by time (`"time"`) or by follow-up count (`"fup"`)
#' @param facet A boolean object that indicates whether or not to facet by treatment
#' @param n.iter The number of iterations to update the model whilst monitoring additional parameters (if necessary).
#' Must be a positive integer. Default is the value used in `mbnma`.
#' @param n.thin The thinning rate. Must be a positive integer. Default is the value used in `mbnma`.
#' @param ... Arguments to be sent to `ggplot2::ggplot()`
#'
#' @return Generates a plot of deviance contributions and returns a list containing the
#' plot (as an object of class `c("gg", "ggplot")`), and a data.frame of posterior mean
#' deviance/residual deviance contributions for each observation.
#'
#' @details
#' Deviances should only be plotted for models that have converged successfully. If deviance
#' contributions have not been monitored in `mbnma$parameters.to.save` then additional
#' iterations will have to be run to get results for these.
#'
#' Deviance contributions cannot be calculated for models with a multivariate likelihood (i.e.
#' those that account for correlation between observations) because the covariance matrix in these
#' models is treated as unknown (if `rho="estimate"`) and deviance contributions will be correlated.
#'
#' @examples
#' \donttest{
#' # Make network
#' alognet <- mb.network(alog_pcfb)
#'
#' # Run MBNMA
#' mbnma <- mb.run(alognet, fun=tpoly(degree=2), intercept=FALSE)
#'
#' # Plot residual deviance contributions in a scatterplot
#' devplot(mbnma)
#'
#' # Plot deviance contributions in boxplots at each follow-up measurement
#' # Monitor for 500 additional iterations
#' devplot(mbnma, dev.type="dev", plot.type="box", xaxis="fup", n.iter=500)
#' }
#' @export
devplot <- function(mbnma, dev.type="dev", plot.type="box",
xaxis="time", facet=TRUE,
n.iter=round(mbnma$BUGSoutput$n.iter/4), n.thin=mbnma$BUGSoutput$n.thin,
...) {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(mbnma, "mbnma", add=argcheck)
checkmate::assertChoice(dev.type, choices = c("dev", "resdev"), add=argcheck)
checkmate::assertChoice(plot.type, choices = c("scatter", "box"), add=argcheck)
checkmate::assertLogical(facet, add=argcheck)
checkmate::assertChoice(xaxis, choices = c("time", "fup"), add=argcheck)
checkmate::assertInt(n.iter, lower=1, null.ok = TRUE, add=argcheck)
checkmate::assertInt(n.thin, lower=1, upper=n.iter, null.ok = TRUE, add=argcheck)
checkmate::reportAssertions(argcheck)
# Check if !is.null(rho) (in which case cannot monitor deviances)
# if (!is.null(mbnma$model.arg$rho)) {
if (mbnma$model.arg$covar %in% c("AR1", "CS")) {
stop("Correlation between time points has been modelled using a multivariate likelihood.
Deviances cannot be calculated for this model.")
}
if (!(dev.type %in% mbnma$parameters.to.save)) {
msg <- paste0("`", dev.type, "` not monitored in mbnma$parameters.to.save.\n",
"additional iterations will be run in order to obtain results for `", dev.type, "`")
message(msg)
dev.df <- mb.update(mbnma, param=dev.type, n.iter=n.iter, n.thin=n.thin)
} else {
dev.df <- get.theta.dev(mbnma, param=dev.type)
}
# Changes fup to time if xaxis=time
if (xaxis=="time") {
dev.df$fup <- apply(dev.df, MARGIN=1, FUN=function(x) mbnma$model.arg$jagsdata$time[x[1], x[3]])
xlab <- "Time"
} else if (xaxis=="fup") {
xlab <- "Follow-up count"
}
if (facet==TRUE) {
dev.df$facet <- apply(dev.df, MARGIN=1, FUN=function(x) mbnma$model.arg$jagsdata$treat[x[1], x[2]])
}
# Plots the residual deviances over time grouped by study and arm
if (plot.type=="scatter") {
g <- ggplot2::ggplot(dev.df, ggplot2::aes(x=fup, y=mean), group=(paste(study, arm, sep="_")), ...) +
ggplot2::geom_point()
} else if (plot.type=="box") {
g <- ggplot2::ggplot(dev.df, ggplot2::aes(x=factor(fup), y=mean), ...) +
ggplot2::geom_boxplot()
}
# Add facets
if (facet==TRUE) {
g <- g + ggplot2::facet_wrap(~factor(facet, labels=mbnma$network$treatments), scales="free_x")
}
# Add axis labels
g <- g + ggplot2::xlab(xlab) +
ggplot2::ylab("Posterior mean") +
ggplot2::geom_hline(yintercept = 0, lty="dashed") +
theme_mbnma()
graphics::plot(g)
return(invisible(list("graph"=g, "dev.data"=dev.df)))
}
#' Plot fitted values from MBNMA model
#'
#' @param treat.labs A character vector of treatment labels with which to name graph panels.
#' Can use `mb.network()[["treatments"]]` with original dataset if in doubt.
#' @param disp.obs A boolean object to indicate whether raw data responses should be
#' plotted as points on the graph
#' @param ... Arguments to be sent to `ggplot2::ggplot()`
#' @inheritParams R2jags::jags
#' @inheritParams devplot
#'
#' @return Generates a plot of fitted values from the MBNMA model and returns a list containing
#' the plot (as an object of class `c("gg", "ggplot")`), and a data.frame of posterior mean
#' fitted values for each observation.
#'
#' @details
#' Fitted values should only be plotted for models that have converged successfully.
#' If fitted values (`theta`) have not been monitored in `mbnma$parameters.to.save`
#' then additional iterations will have to be run to get results for these.
#'
#' @examples
#' \donttest{
#' # Make network
#' painnet <- mb.network(osteopain)
#'
#' # Run MBNMA
#' mbnma <- mb.run(painnet,
#' fun=temax(pool.emax="rel", method.emax="common",
#' pool.et50="abs", method.et50="random"))
#'
#' # Plot fitted values from the model
#' # Monitor fitted values for 500 additional iterations
#' fitplot(mbnma, n.iter=500)
#' }
#'
#' @export
fitplot <- function(mbnma, treat.labs=NULL, disp.obs=TRUE,
n.iter=round(mbnma$BUGSoutput$n.iter/4), n.thin=mbnma$BUGSoutput$n.thin, ...) {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(mbnma, "mbnma", add=argcheck)
checkmate::assertLogical(disp.obs, add=argcheck)
checkmate::assertInt(n.iter, lower=1, null.ok = TRUE, add=argcheck)
checkmate::assertInt(n.thin, lower=1, upper=n.iter, null.ok = TRUE, add=argcheck)
checkmate::reportAssertions(argcheck)
if (!("theta" %in% mbnma$parameters.to.save)) {
msg <- paste0("`theta` not monitored in mbnma$parameters.to.save.\n",
"additional iterations will be run in order to obtain results")
message(msg)
theta.df <- mb.update(mbnma, param="theta", n.iter=n.iter, n.thin=n.thin)
} else {
theta.df <- get.theta.dev(mbnma, param="theta")
}
# Get treatment codes
theta.df <- get.treattimes(mbnma, theta.df, type="treat")
# Get times
theta.df <- get.treattimes(mbnma, theta.df, type="time")
# Get y
theta.df <- get.treattimes(mbnma, theta.df, type="y")
# Check if intercept is modelled
base.dat <- theta.df[theta.df$fup==1,]
if (mbnma$model.arg$intercept==FALSE | all(base.dat$time!=0)) {
if (any(base.dat$time!=0)) {
message(cat("Absence of observations with time=0 in network - data assumed to be change from baseline:",
"plotting response=0 at time=0", sep="\n"))
}
base.dat$y <- rep(0, nrow(base.dat))
base.dat$time <- rep(0, nrow(base.dat))
base.dat$fup <- rep(0, nrow(base.dat))
base.dat$mean <- rep(0, nrow(base.dat))
theta.df <- rbind(theta.df, base.dat)
}
# Generate plot
g <- ggplot2::ggplot(theta.df,
ggplot2::aes(x=time, y=mean, group=paste(study, arm, sep="_")), ...) +
ggplot2::geom_line()
# Overlay observed responses
if (disp.obs==TRUE) {
g <- g + ggplot2::geom_point(ggplot2::aes(y=y), size=1)
}
# Add facet labels
if (!is.null(treat.labs)) {
if (length(treat.labs)!=max(theta.df$treat)) {
stop("treat.labs must be the same length as the number of treatment codes in the model")
}
g <- g + ggplot2::facet_wrap(~factor(treat, labels=treat.labs))
} else {
g <- g + ggplot2::facet_wrap(~factor(treat, labels=mbnma$network$treatments))
}
# Add axis labels
g <- g + ggplot2::xlab("Time") +
ggplot2::ylab("Response") +
theme_mbnma()
graphics::plot(g)
}
# Get matching treatments or times for a data.frame from an MBNMA model
get.treattimes <- function(mbnma, add.df, type="treat") {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(mbnma, "mbnma", add=argcheck)
checkmate::assertChoice(type, choices = c("time", "treat", "y"), add=argcheck)
checkmate::assertDataFrame(add.df, add=argcheck)
checkmate::reportAssertions(argcheck)
# Get codes
#index.df <- reshape2::melt(mbnma$model$data()[[type]])
index.df <- reshape2::melt(mbnma$model.arg$jagsdata[[type]])
index.df <- index.df[stats::complete.cases(index.df),]
# Match codes to data frame
if (type=="treat") {
var <- "arm"
} else if (type=="time") {
var <- "fup"
}
if (type %in% c("treat", "time")) {
studyarm.index <- paste(index.df[,1], index.df[,2], sep="_")
studyarm.df <- paste(add.df$study, add.df[[var]], sep="_")
add.df[[type]] <- index.df[,3][match(studyarm.df, studyarm.index)]
} else if (type=="y") {
studyarm.index <- paste(index.df[,1], index.df[,2], index.df[,3], sep="_")
studyarm.df <- paste(add.df$study, add.df$arm, add.df$fup, sep="_")
add.df[[type]] <- index.df[,4][match(studyarm.df, studyarm.index)]
}
return(add.df)
}
#' Extracts fitted values or deviance contributions into a data.frame with indices
#' @noRd
get.theta.dev <- function(mbnma, param="theta") {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(mbnma, "mbnma", add=argcheck)
checkmate::assertChoice(param, choices = c("dev", "resdev", "theta"), add=argcheck)
checkmate::reportAssertions(argcheck)
dev.df <- as.data.frame(
mbnma$BUGSoutput$summary[grep(paste0("^", param, "\\["),
rownames(mbnma$BUGSoutput$summary)),]
)
if (mbnma$type=="time") {
# Takes the study, arm and follow-up measurement identifier for each residual deviance point
id <- matrix(unlist(strsplit(
gsub(
"([a-z]+\\[)([0-9]+)(,)([0-9]+)(,)([0-9]+)(\\])",
"\\2.\\4.\\6", rownames(dev.df)),
split="\\.")),
byrow=TRUE,
ncol=3
)
dev.df$fup <- as.numeric(id[,3])
} else if (mbnma$type=="dose") {
# Takes the study, arm and follow-up measurement identifier for each residual deviance point
id <- matrix(unlist(strsplit(
gsub(
"([a-z]+\\[)([0-9]+)(,)([0-9]+)(\\])",
"\\2.\\4", rownames(dev.df)),
split="\\.")),
byrow=TRUE,
ncol=2
)
}
dev.df$study <- as.numeric(id[,1])
dev.df$arm <- as.numeric(id[,2])
return(dev.df)
}
#' Plot illustrative time-course functions (CURRENTLY NOT FUNCTIONAL FOR ALL TIME_COURSE FUNCTIONS)
#'
#' Can be used to plot potential time-course functions to identify if they may
#' be a suitable fit for the data.
#'
#' @param x An object of `class("timefun")`
#' @inheritParams tuser
#'
#' @return An illustrative plot of the time-course function with the parameters specified
#' in `beta.1`, `beta.2`, `beta.3` and `beta.4`
#'
#' @noRd
plot.timefun <- function(x=tpoly(degree=1), beta.1=0, beta.2=0,
beta.3=0, beta.4=0) {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(x, "timefun", add=argcheck)
checkmate::assertNumeric(beta.1, add=argcheck)
checkmate::assertNumeric(beta.2, add=argcheck)
checkmate::assertNumeric(beta.3, add=argcheck)
checkmate::assertNumeric(beta.4, add=argcheck)
checkmate::reportAssertions(argcheck)
time <- seq(0,100)
funstr <- x$jags
#funstr <- gsub("beta", "tempbeta", funstr)
funstr <- gsub("\\[i\\,k\\]", "", funstr)
funstr <- gsub("\\[i\\,m\\]", "", funstr)
if (any(c("ns", "bs", "ls") %in% x$name)) {
spline <- genspline(time, spline=x$name, knots=x$knots)
funstr <- gsub("spline\\[i\\,m", "spline[", funstr)
}
if (any(c("user", "fpoly") %in% x$name)) {
stop("Plotting for 'tfpoly()' and 'tuser()' not currently supported")
}
nparam <- x$nparam
df <- data.frame("y"=NA, "time"=NA)
y <- eval(parse(text=funstr))
df <- rbind(df, stats::setNames(cbind(y,time), names(df)))
df <- df[-1,]
vals <- vector()
for (i in 1:nparam) {
vals <- append(vals, get(paste0("beta.", i)))
}
g <- ggplot2::ggplot(df, ggplot2::aes(x=time, y=y)) +
ggplot2::geom_line(linewidth=1) +
ggplot2::labs(subtitle=paste(paste0("beta.", 1:nparam, " = ", vals), collapse=" ")) +
ggplot2::xlab("Time") +
theme_mbnma()
graphics::plot(g)
return(invisible(g))
}
plot.invisible <- function(...){
ff <- tempfile()
grDevices::png(filename=ff)
res <- graphics::plot(...)
grDevices::dev.off()
unlink(ff)
return(res)
}
#' MBNMA ggplot2 theme style
#' @noRd
theme_mbnma <- function(...) {
ggplot2::theme_bw(...) +
ggplot2::theme(
# change stuff here
panel.background = ggplot2::element_blank(),
plot.background = ggplot2::element_rect(fill="white", colour=NA),
legend.background = ggplot2::element_rect(fill="transparent", colour=NA),
legend.key = ggplot2::element_rect(fill="transparent", colour=NA),
# From multinma
#panel.border = ggplot2::element_rect(colour = "grey70", fill = NA),
panel.grid.major = ggplot2::element_line(colour = "grey95"),
panel.grid.minor = ggplot2::element_line(colour = "grey95"),
strip.background = ggplot2::element_rect(colour = "black",
fill = "lightsteelblue1"),
strip.text = ggplot2::element_text(colour = "black")
)
}
#' Plot cumulative ranking curves from MBNMA models
#'
#' @param x An object of class `"mb.rank"` generated by `rank.mbnma()`
#' @param sucra A logical object to indicate whether Surface Under Cumulative Ranking Curve (SUCRA)
#' values should be calculated and returned as a data frame. Areas calculated
#' using trapezoid approach.
#' @param ... Arguments to be sent to `ggplot::geom_line()`
#'
#' @return Line plots showing the cumulative ranking probabilities for each agent/class for
#' the ranked dose response paramtere in `x`. The object returned is a list which contains the plot
#' (an object of `class(c("gg", "ggplot")`) and a data frame of SUCRA values
#' if `sucra = TRUE`.
#'
#' @examples
#' \donttest{
#' # Using the alogliptin data
#' network <- mb.network(alog_pcfb)
#'
#' # Estimate rankings from an Emax dose-response MBNMA
#' emax <- mb.run(network, fun=temax())
#' ranks <- rank(emax, param=c("emax"))
#'
#' # Plot cumulative rankings for both dose-response parameters simultaneously
#' # Note that SUCRA values are also returned
#' cumrank(ranks)
#' }
#' @export
cumrank <- function(x, sucra=TRUE, ...) {
# Run checks
argcheck <- checkmate::makeAssertCollection()
checkmate::assertClass(x, "mb.rank", add=argcheck)
checkmate::assertLogical(sucra, null.ok=FALSE, add=argcheck)
checkmate::reportAssertions(argcheck)
output <- list()
# if (is.null(params)) {
# params <- names(x)
# }
params <- x$param
df <- data.frame()
cum.mat <- x$cum.matrix
treats <- colnames(cum.mat)
melt <- reshape2::melt(cum.mat)
melt$param <- params
df <- rbind(df, melt)
# for (param in seq_along(params)) {
# if (!params[param] %in% x$param) {
# stop(paste0(params[param], " is not a ranked parameter in x"))
# }
#
# cum.mat <- x[[params[param]]]$cum.matrix
# treats <- colnames(cum.mat)
#
# melt <- reshape2::melt(cum.mat)
# melt$param <- params[param]
#
# df <- rbind(df, melt)
#
# }
df$Parameter <- factor(df$param)
g <- ggplot2::ggplot(df, ggplot2::aes(x=Var1, y=value, linetype=Parameter, colour=Parameter), ...) +
ggplot2::geom_line(linewidth=1)
g <- g + ggplot2::facet_wrap(~factor(Var2)) +
ggplot2::xlab("Rank (1 = best)") +
ggplot2::ylab("Cumulative probability") +
ggplot2::labs(linetype="Parameter", colour="Parameter") +
theme_mbnma()
graphics::plot(g)
output <- list("cumplot"=g)
# Calculate AUC
if (sucra==TRUE) {
df.auc <- df %>%
dplyr::group_by(df$Var2, df$param) %>%
dplyr::do(data.frame(sucra=calcauc(.)))
df.auc <- dplyr::ungroup(df.auc)
names(df.auc) <- c("treatment", "parameter", "sucra")
output[["sucra"]] <- df.auc
print(df.auc)
}
return(invisible(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.