Nothing
#' Plotting method for \code{survFitTKTD} objects
#'
#' This is the generic \code{plot} S3 method for the
#' \code{survFitTKTD}. It plots the fit obtained for each
#' concentration of chemical compound in the original dataset.
#'
#' The fitted curves represent the \strong{estimated survival probablity} as a function
#' of time for each concentration
#' When \code{adddata = TRUE} the black dots depict the \strong{observed survival
#' probablity} at each time point. Note that since our model does not take
#' inter-replicate variability into consideration, replicates are systematically
#' pooled in this plot.
#' The function plots both 95\% credible intervals for the estimated survival
#' probablity (by default the grey area around the fitted curve) and 95\% binomial confidence
#' intervals for the observed survival probablity (as black error bars if
#' \code{adddata = TRUE}).
#' Both types of intervals are taken at the same level. Typically
#' a good fit is expected to display a large overlap between the two types of intervals.
#' If \code{spaghetti = TRUE}, the credible intervals are represented by two
#' dotted lines limiting the credible band, and a spaghetti plot is added to this band.
#' This spaghetti plot consists of the representation of simulated curves using parameter values
#' sampled in the posterior distribution (2\% of the MCMC chains are randomly
#' taken for this sample).
#'
#' @param x An object of class \code{survFitTKTD}.
#' @param xlab A label for the \eqn{X}-axis, by default \code{Time}.
#' @param ylab A label for the \eqn{Y}-axis, by default \code{Survival probablity}.
#' @param main A main title for the plot.
#' @param concentration A numeric value corresponding to some specific concentration in
#' \code{data}. If \code{concentration = NULL}, draws a plot for each concentration.
#' @param spaghetti if \code{TRUE}, draws a set of survival curves using
#' parameters drawn from the posterior distribution
#' @param one.plot if \code{TRUE}, draws all the estimated curves in
#' one plot instead of one plot per concentration.
#' @param adddata if \code{TRUE}, adds the observed data to the plot
#' with (frequentist binomial) confidence intervals
#' @param addlegend if \code{TRUE}, adds a default legend to the plot.
#' @param style graphical backend, can be \code{'generic'} or \code{'ggplot'}
#' @param \dots Further arguments to be passed to generic methods.
#'
#' @keywords plot
#'
#' @return a plot of class \code{ggplot}
#'
#' @examples
#'
#' # (1) Load the survival data
#' data(propiconazole)
#'
#' # (2) Create an object of class "survData"
#' dataset <- survData(propiconazole)
#'
#' \donttest{
#' # (3) Run the survFitTKTD function ('SD' model only)
#' out <- survFitTKTD(dataset)
#'
#' # (4) Plot the fitted curves in one plot
#' plot(out)
#'
#' # (5) Plot one fitted curve per concentration with credible limits as
#' # spaghetti, data and confidence intervals
#' # and with a ggplot style
#' plot(out, spaghetti = TRUE , adddata = TRUE, one.plot = FALSE,
#' style = "ggplot")
#'
#' # (6) Plot fitted curve for one specific concentration
#' plot(out, concentration = 36, style = "ggplot")
#' }
#'
#' @export
#'
#' @import ggplot2
#' @import grDevices
#' @importFrom dplyr filter
#' @importFrom gridExtra grid.arrange arrangeGrob
#' @importFrom grid grid.rect gpar
#' @importFrom graphics plot
#'
plot.survFitTKTD <- function(x,
xlab = "Time",
ylab = "Survival probablity",
main = NULL,
concentration = NULL,
spaghetti = FALSE,
one.plot = FALSE,
adddata = FALSE,
addlegend = FALSE,
style = "ggplot", ...) {
if (one.plot && !is.null(concentration))
one.plot <- FALSE
if ((addlegend && is.null(concentration)) ||
(addlegend && !one.plot))
warning("The legend is available only if [one.plot] is TRUE or if [concentration] is not NULL !", call. = FALSE)
if (!is.null(concentration) && !any(x$transformed.data$conc == concentration))
stop("The [concentration] argument is not one of the possible concentration !")
if (one.plot)
warning("The credible limits and confidence intervals are not drawn when 'one.plot' = TRUE.", call. = FALSE)
conf.int <- survTKTDConfInt(x)
dobs <- data.frame(conc = x$transformed.data$conc,
time = x$transformed.data$time,
psurv = x$transformed.data$N_alive / x$transformed.data$N_init,
Points = "Observed values",
color = as.numeric(as.factor(x$transformed.data$conc)),
conf.int)
# remove time 0 in dobs
dobs <- filter(dobs, time != 0)
data.credInt <- survFitPlotCITKTD(x)
if (style == "generic") {
survFitPlotTKTDGeneric(data.credInt, dobs, xlab, ylab, main, concentration,
one.plot, spaghetti,
adddata, addlegend)
} else if (style == "ggplot") {
survFitPlotTKTDGG(data.credInt, dobs, xlab, ylab, main, concentration,
one.plot, spaghetti, adddata, addlegend)
} else stop("Unknown style")
}
Surv <- function(Cw, time, ks, kd, NEC, m0)
{
S <- exp(-m0*time)
x <- ifelse(Cw > NEC, 1 - NEC/Cw, NA)
tNEC <- -(1/kd)*log(x)
y <- ifelse(time > tNEC,
exp( ks/kd*Cw*(exp(-kd*tNEC) -exp(-kd*time))
- ks*(Cw-NEC)*(time - tNEC)),
NA)
return(ifelse(!is.na(x) & !is.na(y), S * y, S))
}
#' @importFrom stats aggregate binom.test
survTKTDConfInt <- function(x) {
# create confidente interval on observed data
# binomial model by a binomial test
# INPUT:
# - x : object of class survFitTT
# OUTPUT:
# - ci : confidente interval
ci <- apply(x$transformed.data, 1, function(x) {
binom.test(x["N_alive"], x["N_init"])$conf.int
})
ci <- as.data.frame(t(ci))
colnames(ci) <- c("qinf95", "qsup95")
ci$Conf.Int <- "Confidence interval"
return(ci)
}
survFitPlotCITKTD <- function(x) {
# INPUT
# x : An object of class survFitTKTD
# OUTPUT
# A list of - dobs : observed values
# - dtheo : estimated values
npoints <- 100
concobs <- unique(x$transformed.data$conc)
tfin <- seq(0, max(x$jags.data$t), length.out = npoints)
# prameters
mctot <- do.call("rbind", x$mcmc)
ks <- 10^mctot[, "log10ks"]
kd <- 10^mctot[, "log10kd"]
m0 <- 10^mctot[, "log10m0"]
nec <- 10^mctot[, "log10NEC"]
# all theorical
k <- 1:length(concobs)
j <- 1:npoints
dtheo <- lapply(k, function(y) { # conc
sapply(j, function(z) { # time
Surv(Cw = concobs[y], time = tfin[z],
ks = ks, kd = kd,
NEC = nec,
m0 = m0)
})
})
# transpose dtheo
dtheo <- do.call("rbind", lapply(dtheo, t))
# quantile
qinf95 <- apply(dtheo, 1, quantile, probs = 0.025, na.rm = TRUE)
qsup95 <- apply(dtheo, 1, quantile, probs = 0.975, na.rm = TRUE)
q50 <- apply(dtheo, 1, quantile, probs = 0.5, na.rm = TRUE)
dtheo <- as.data.frame(dtheo)
colnames(dtheo) <- paste0("X", 1:length(kd))
# divide number of mcmc by 50
sel <- sample(ncol(dtheo))[1:ceiling(ncol(dtheo) / 50)]
dtheo <- dtheo[, sel]
dtheo$conc <- rep(concobs, rep(npoints, length(concobs)))
dtheo$time <- rep(tfin, length(concobs))
# add credible limits
dtheo$qinf95 <- qinf95
dtheo$qsup95 <- qsup95
dtheo$q50 <- q50
# names for legend plots
dtheo$Cred.Lim <- "Credible limits"
dtheo$Mean.C <- "Mean curve"
# vector color
dtheo$color <- as.numeric(as.factor(dtheo$conc))
return(dtheo)
}
survFitPlotTKTDGeneric <- function(data, dobs, xlab, ylab, main, concentration,
one.plot, spaghetti, adddata,
addlegend) {
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
if (one.plot) {
survFitPlotTKTDGenericOnePlot(data, dobs, xlab, ylab, main, adddata,
addlegend)
} else if (!one.plot && is.null(concentration)) {
par(mfrow = plotMatrixGeometry(length(unique(dobs$conc))))
survFitPlotTKTDGenericNoOnePlot(data, dobs, xlab, ylab, spaghetti,
adddata, concentration)
par(mfrow = c(1, 1))
} else {
survFitPlotTKTDGenericNoOnePlot(data, dobs, xlab, ylab, spaghetti,
adddata, concentration, addlegend)
}
}
survFitPlotTKTDGenericOnePlot <- function(data, dobs, xlab, ylab, main, adddata,
addlegend) {
plot(c(0, dobs$time),
c(dobs$psurv, 1),
xlab = xlab,
ylab = ylab,
type = "n",
main = main)
by(data, list(data$conc),
function(x) {
lines(x$time, x$q50, # lines
col = x$color)
})
# points
if (adddata) {
points(dobs$time,
dobs$psurv,
pch = 20,
col = dobs$color)
}
if (addlegend) {
legend("bottomleft",
legend = c(ifelse(adddata, "Observed values", NA),
unique(dobs$conc)),
pch = c(ifelse(adddata, 20, NA),
rep(NA, length(unique(dobs$conc)))),
lty = c(NA, rep(1, length(unique(dobs$conc)))),
bty = "n",
cex = 1,
ncol = 2,
col = c(ifelse(adddata, "black", NA),
unique(dobs$color)),
title = "Concentrations")
}
}
#' @importFrom reshape2 melt
survFitPlotTKTDGenericNoOnePlot <- function(data, dobs, xlab, ylab, spaghetti,
adddata, concentration,
addlegend = FALSE) {
delta <- 0.01 * (max(dobs$time) - min(dobs$time))
#remove unused column before melt
dataTmp <- data
dataTmp[, c("Cred.Lim", "Mean.C")] <- list(NULL)
dtheoQm <- melt(dataTmp,
id.vars = c("conc", "time", "color"))
rm(dataTmp)
if (is.null(concentration)) {
dobs <- split(dobs, dobs$conc)
dtheoQ <- split(data, data$conc)
dtheoQm <- split(dtheoQm, dtheoQm$conc)
} else {
dobs <- list(filter(dobs, conc == concentration))
dtheoQ <- list(filter(data, conc == concentration))
dtheoQm <- list(filter(dtheoQm, conc == concentration))
}
mapply(function(x, y, z) {
plot(x[, "time"],
x[, "q50"],
xlab = xlab,
ylab = ylab,
type = "n",
ylim = c(0, 1),
main = paste0("Concentration = ", unique(x[, "conc"])))
if (spaghetti) {
color <- "grey"
color_transparent <- adjustcolor(color, alpha.f = 0.05)
by(z, z$variable, function(x) {
lines(x[, "time"], x[, "value"], col = color_transparent)
})
} else {
polygon(c(x[, "time"], rev(x[, "time"])), c(x[, "qinf95"],
rev(x[, "qsup95"])),
col = "grey70", border = NA)
}
lines(x[, "time"], x[, "q50"], # lines
col = "orange")
lines(x[, "time"], x[, "qinf95"],
col = "gray70")
lines(x[, "time"], x[, "qsup95"],
col = "gray70")
if (adddata) {
points(y[, "time"],
y[, "psurv"],
pch = 20,
col = "black") # points
segments(y[, "time"], y[, "qinf95"],
y[, "time"], y[, "qsup95"],
col = "black")
}
if (addlegend) {
legend("bottomleft", pch = c(ifelse(adddata, 20, NA), NA, NA, NA),
lty = c(NA, ifelse(adddata, 1, NA), 1, 1),
col = c(ifelse(adddata, "black", NA),
ifelse(adddata, "black", NA),
"orange", "grey70"),
legend = c(ifelse(adddata, "Observed values", NA),
ifelse(adddata, "Confidence interval", NA),
"Mean curve", "Credible limits"),
bty = "n")
}
}, x = dtheoQ, y = dobs, z = dtheoQm)
}
survFitPlotTKTDGG <- function(data, dobs, xlab, ylab, main, concentration,
one.plot, spaghetti, adddata, addlegend) {
if (one.plot) {
survFitPlotTKTDGGOnePlot(data, dobs, xlab, ylab, main, adddata, addlegend)
} else if (!one.plot && is.null(concentration)) {
survFitPlotTKTDGGNoOnePlot(data, dobs, xlab, ylab, main, spaghetti,
adddata, concentration)
} else {
survFitPlotTKTDGGNoOnePlot(data, dobs, xlab, ylab, main, spaghetti,
adddata, concentration, addlegend)
}
}
survFitPlotTKTDGGOnePlot <- function(data, dobs, xlab, ylab, main, adddata,
addlegend) {
gf <- ggplot(dobs) +
geom_line(aes(x = time, y = q50, colour = factor(conc)),
data = data) +
labs(x = xlab, y = ylab) + ggtitle(main) +
ylim(c(0, 1)) +
theme_minimal()
if (adddata) {
gf <- gf + geom_point(aes(x = time, y = psurv, colour = factor(conc)),
data = dobs)
}
if (addlegend) {
gf + scale_color_discrete("Concentrations")
} else {
gf + scale_color_discrete(guide = "none")
}
}
#' @importFrom dplyr filter
#' @importFrom reshape2 melt
survFitPlotTKTDGGNoOnePlot <- function(data, dobs, xlab, ylab, main, spaghetti,
adddata, concentration,
addlegend = FALSE) {
# colors
valCols <- fCols(data, "orange", "gray70")
dataTmp <- data
dataTmp[, c("Cred.Lim", "Mean.C")] <- list(NULL)
dtheoQm <- melt(dataTmp,
id.vars = c("conc", "time", "color"))
rm(dataTmp)
if (!is.null(concentration)) {
dobs <- filter(dobs, conc == concentration)
data <- filter(data, conc == concentration)
dtheoQm <- filter(dtheoQm, conc == concentration)
curv_resp <- data.frame(time = data$time,
resp = data$q50,
Line = "Mean curve")
}
if (spaghetti) {
gf <- ggplot(dobs) +
geom_line(data = dtheoQm, aes(x = time, y = value, group = variable),
alpha = 0.05)
} else {
gf <- ggplot(dobs) +
geom_ribbon(data = data, aes(x = time, ymin = qinf95,
ymax = qsup95),
fill = valCols$cols4,
# color = valCols$cols4,
alpha = 0.4)
}
if (!is.null(concentration)) {
gf <- gf + geom_line(data = curv_resp, aes(x = time, y = resp,
color = Line)) +
scale_colour_hue(name = "")
} else {
gf <- gf + geom_line(data = data, aes(x = time, y = q50),
color = valCols$cols5)
}
gf <- gf +
# geom_line(data = data, aes(x = time, y = qinf95), color = "grey70") +
# geom_line(data = data, aes(x = time, y = qsup95), color = "grey70") +
facet_wrap(~conc) +
scale_linetype(name = "") +
labs(x = xlab, y = ylab) + ggtitle(main) +
ylim(c(0, 1)) +
theme_minimal()
if (adddata) {
# dataframes points (data) and curve (curv)
gf <- gf +
geom_point(aes(x = time, y = psurv, fill = Points),
data = dobs, col = valCols$cols1) +
geom_segment(aes(x = time, xend = time, y = qinf95, yend = qsup95,
linetype = Conf.Int),
dobs, col = valCols$cols3,
size = 0.5) +
scale_fill_hue("")
} else {
gf <- gf
}
if (addlegend) {
gf
} else {
gf + theme(legend.position = "none")
}
}
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.