-
virgile.baudrot authoredvirgile.baudrot authored
plot.survFitTT.R 18.77 KiB
#' Plotting method for \code{survFitTT} objects
#'
#' This is the generic \code{plot} S3 method for the \code{survFitTT} class. It
#' plots concentration-response fit under target time survival analysis.
#'
#' The fitted curve represents the \strong{estimated survival probability} at
#' the target time as a function of the concentration of chemical compound;
#' When \code{adddata = TRUE} the black dots depict the \strong{observed survival
#' probability} at each tested concentration. 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
#' probability (by default the grey area around the fitted curve) and 95\% binomial confidence
#' intervals for the observed survival probability (as black segments 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 intervals.
#' If 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 (10\% of the MCMC chains are randomly
#' taken for this sample).
#'
#' @param x an object of class \code{survFitTT}
#' @param xlab a label for the \eqn{X}-axis, default is \code{Concentration}
#' @param ylab a label for the \eqn{Y}-axis, default is \code{Survival probability}
#' @param main main title for the plot
#' @param fitcol color of the fitted curve
#' @param fitlty line type of the fitted curve
#' @param fitlwd width of the fitted curve
#' @param spaghetti if \code{TRUE}, the credible interval is represented by
#' multiple curves
#' @param cicol color of the 95 \% credible interval limits
#' @param cilty line type for the 95 \% credible interval limits
#' @param cilwd width of the 95 \% credible interval limits
#' @param ribcol color of the ribbon between lower and upper credible limits.
#' Transparent if \code{NULL}
#' @param adddata if \code{TRUE}, adds the observed data with confidence intervals
#' to the plot
#' @param addlegend if \code{TRUE}, adds a default legend to the plot
#' @param log.scale if \code{TRUE}, displays \eqn{X}-axis in log-scale
#' @param style graphical backend, can be \code{'generic'} or \code{'ggplot'}
#' @param \dots Further arguments to be passed to generic methods
#' @note When \code{style = "ggplot"}, the function calls function
#' \code{\link[ggplot2]{ggplot}} and returns an object of class \code{ggplot}.
#'
#' @return a plot of class \code{ggplot}
#'
#' @examples
#'
#' # (1) Load the data
#' data(cadmium1)
#'
#' # (2) Create an object of class "survData"
#' dat <- survData(cadmium1)
#'
#' \donttest{
#' # (3) Run the survFitTT function with the log-logistic
#' # binomial model
#' out <- survFitTT(dat, lcx = c(5, 10, 15, 20, 30, 50, 80),
#' quiet = TRUE)
#'
#' # (4) Plot the fitted curve
#' plot(out, log.scale = TRUE, adddata = TRUE)
#'
#' # (5) Plot the fitted curve with ggplot style
#' plot(out, xlab = expression("Concentration in" ~ mu~g.L^{-1}),
#' fitcol = "blue", adddata = TRUE, cicol = "blue",
#' style = "ggplot")
#' }#'
#' @keywords plot
#'
#' @import grDevices
#' @import ggplot2
#' @importFrom gridExtra grid.arrange arrangeGrob
#' @importFrom grid grid.rect gpar
#' @importFrom graphics plot axis legend lines par points polygon segments
#' @importFrom stats aggregate
#' @importFrom reshape2 melt
#'
#' @export
plot.survFitTT <- function(x,
xlab = "Concentration",
ylab = "Survival probability",
main = NULL,
fitcol = "orange",
fitlty = 1,
fitlwd = 1,
spaghetti = FALSE,
cicol = "orange",
cilty = 2,
cilwd = 1,
ribcol = "grey70",
adddata = FALSE,
addlegend = FALSE,
log.scale = FALSE,
style = "ggplot", ...) {
# plot the fitted curve estimated by survFitTT
# INPUTS
# - x: survFitTt object
# - xlab : label x
# - ylab : label y
# - main : main title
# - fitcol : color fitted curve
# - fitlty : type line fitted curve
# - fitlwd : width line fitted curve
# - cicol : color ci ribbon
# - cilty : type line ci ribbon
# - cilwd : width line ci ribbon
# - addlegend : boolean
# - log.scale : x log option
# - style : generic or ggplot
# OUTPUT:
# - plot of fitted regression
# Selection of datapoints that can be displayed given the type of scale
sel <- if(log.scale) x$dataTT$conc > 0 else TRUE
dataTT <- x$dataTT[sel, ]
dataTT$resp <- dataTT$Nsurv / dataTT$Ninit
# data points are systematically pooled, since our model does not
# take individual variation into account
dataTT <- aggregate(resp ~ conc, dataTT, mean)
transf_data_conc <- optLogTransform(log.scale, dataTT$conc)
# Concentration values used for display in linear scale
display.conc <- (function() {
x <- optLogTransform(log.scale, dataTT$conc)
s <- seq(min(x),max(x), length = 100)
if(log.scale) exp(s) else s
})()
# Possibly log transformed concentration values for display
curv_conc <- optLogTransform(log.scale, display.conc)
conf.int <- survLlbinomConfInt(x, log.scale)
cred.int <- survMeanCredInt(x, display.conc)
spaghetti.CI <- if (spaghetti) { survSpaghetti(x, display.conc) } else NULL
dataCIm <- if (spaghetti) {melt(cbind(curv_conc, spaghetti.CI), id.vars = c("curv_conc", "conc"))} else NULL
curv_resp <- data.frame(conc = curv_conc, resp = cred.int[["q50"]],
Line = "loglogistic")
if (style == "generic") {
survFitPlotGenericCredInt(x,
dataTT$conc, transf_data_conc, dataTT$resp,
curv_conc, curv_resp,
conf.int, cred.int, spaghetti.CI, dataCIm,
xlab, ylab, fitcol, fitlty, fitlwd,
main, addlegend, adddata,
cicol, cilty, cilwd, ribcol, log.scale)
}
else if (style == "ggplot") {
survFitPlotGG(x,
dataTT$conc, transf_data_conc, dataTT$resp,
curv_conc, curv_resp,
conf.int, cred.int, spaghetti.CI, dataCIm,
xlab, ylab, fitcol, fitlty, fitlwd,
main, addlegend, adddata,
cicol, cilty, cilwd / 2, ribcol)
}
else stop("Unknown style")
}
#' @importFrom stats aggregate binom.test
survLlbinomConfInt <- function(x, log.scale) {
# create confidente interval on observed data for the log logistic
# binomial model by a binomial test
# INPUT:
# - x : object of class survFitTT
# - log.scale : boolean
# OUTPUT:
# - ci : confidente interval
x <- cbind(aggregate(Nsurv ~ time + conc, x$dataTT, sum),
Ninit = aggregate(Ninit ~ time + conc, x$dataTT, sum)$Ninit)
ci <- apply(x, 1, function(x) {
binom.test(x["Nsurv"], x["Ninit"])$conf.int
})
rownames(ci) <- c("qinf95", "qsup95")
colnames(ci) <- x$conc
if (log.scale) ci <- ci[ ,colnames(ci) != 0]
return(ci)
}
#' @importFrom stats quantile
survMeanCredInt <- function(fit, x) {
# create the parameters for credible interval for the log logistic binomial
# model
# INPUT:
# - fit : object of class survFitTT
# - x : vector of concentrations values (x axis)
# OUTPUT:
# - ci : credible limit
mctot <- do.call("rbind", fit$mcmc)
k <- nrow(mctot)
# parameters
if (fit$det.part == "loglogisticbinom_3") {
d2 <- mctot[, "d"]
}
log10b2 <- mctot[, "log10b"]
b2 <- 10^log10b2
log10e2 <- mctot[, "log10e"]
e2 <- 10^log10e2
# quantiles
qinf95 = NULL
q50 = NULL
qsup95 = NULL
for (i in 1:length(x)) {
# llbinom 2 parameters
if (fit$det.part == "loglogisticbinom_2") {
theomean <- 1 / (1 + (x[i] / e2)^(b2)) # mean curve
}
# llbinom 3 parameters
else if (fit$det.part == "loglogisticbinom_3") {
theomean <- d2 / (1 + (x[i] / e2)^(b2)) # mean curve
}
# IC 95%
qinf95[i] <- quantile(theomean, probs = 0.025, na.rm = TRUE)
q50[i] <- quantile(theomean, probs = 0.5, na.rm = TRUE)
qsup95[i] <- quantile(theomean, probs = 0.975, na.rm = TRUE)
}
# values for CI
ci <- data.frame(qinf95 = qinf95,
q50 = q50,
qsup95 = qsup95)
return(ci)
}
survSpaghetti <- function(fit, x) {
mctot <- do.call("rbind", fit$mcmc)
sel <- sample(nrow(mctot))[1:ceiling(nrow(mctot) / 10)]
# parameters
if (fit$det.part == "loglogisticbinom_3") {
d2 <- mctot[, "d"][sel]
}
log10b2 <- mctot[, "log10b"][sel]
b2 <- 10^log10b2
log10e2 <- mctot[, "log10e"][sel]
e2 <- 10^log10e2
# all theorical
dtheo <- array(data = NA, dim = c(length(x), length(e2)))
for (i in 1:length(e2)) {
# llbinom 2 parameters
if (fit$det.part == "loglogisticbinom_2") {
dtheo[, i] <- 1 / (1 + (x / e2[i])^(b2[i])) # mean curve
}
# llbinom 3 parameters
else if (fit$det.part == "loglogisticbinom_3") {
dtheo[, i] <- d2[i] / (1 + (x / e2[i])^(b2[i])) # mean curve
}
}
dtheof <- as.data.frame(cbind(x, dtheo))
names(dtheof) <- c("conc", paste0("X", 1:length(sel)))
return(dtheof)
}
survFitPlotGenericCredInt <- function(x,
data_conc, transf_data_conc, data_resp,
curv_conc, curv_resp,
conf.int, cred.int, spaghetti.CI, dataCIm,
xlab, ylab, fitcol, fitlty, fitlwd,
main, addlegend, adddata,
cicol, cilty, cilwd, ribcol, log.scale)
{
# plot the fitted curve estimated by survFitTT # with generic style with credible interval
plot(transf_data_conc, data_resp,
xlab = xlab,
ylab = ylab,
main = main,
xaxt = "n",
yaxt = "n",
ylim = c(0, 1.01),
type = "n")
# axis
axis(side = 2, at = pretty(c(0, max(c(conf.int["qsup95",],
cred.int[["qsup95"]])))))
axis(side = 1,
at = transf_data_conc,
labels = data_conc)
# Plotting the theoretical curve
# CI ribbon + lines
if (!is.null(spaghetti.CI)) {
color <- "gray"
color_transparent <- adjustcolor(color, alpha.f = 0.05)
by(dataCIm, dataCIm$variable, function(x) {
lines(x$curv_conc, x$value, col = color_transparent)
})
} else if(!is.null(ribcol)) {
polygon(c(curv_conc, rev(curv_conc)), c(cred.int[["qinf95"]],
rev(cred.int[["qsup95"]])),
col = ribcol, border = NA)
}
lines(curv_conc, cred.int[["qsup95"]], type = "l", col = cicol, lty = cilty,
lwd = cilwd)
lines(curv_conc, cred.int[["qinf95"]], type = "l", col = cicol, lty = cilty,
lwd = cilwd)
if (adddata) {
# segment CI
segments(transf_data_conc, data_resp,
transf_data_conc, conf.int["qsup95", ])
segments(transf_data_conc, data_resp,
transf_data_conc, conf.int["qinf95", ])
# points
points(transf_data_conc, data_resp, pch = 16)
}
# fitted curve
lines(curv_conc, curv_resp[, "resp"], type = "l", col = fitcol, lty = fitlty,
lwd = fitlwd)
# legend
if (addlegend) {
legend("bottomleft", pch = c(ifelse(adddata, 16, NA), NA, NA, NA),
lty = c(NA, ifelse(adddata, 1, NA), cilty, fitlty),
lwd = c(NA, ifelse(adddata,1, NA), cilwd, fitlwd),
col = c(ifelse(adddata, 1, NA), 1, cicol, fitcol),
legend = c(ifelse(adddata, "Observed values", NA),
ifelse(adddata, "Confidence interval", NA),
"Credible limits", x$det.part),
bty = "n")
}
}
#' @importFrom grid arrow unit
survFitPlotGGCredInt <- function(x, data, curv_resp, conf.int, cred.int,
spaghetti.CI, dataCIm, cilty, cilwd,
valCols, fitlty, fitlwd, ribcol, xlab, ylab, main,
adddata) { # IC
data.three <- data.frame(conc = data$transf_conc,
qinf95 = conf.int["qinf95",],
qsup95 = conf.int["qsup95",],
Conf.Int = "Confidence interval")
data.four <- data.frame(conc = curv_resp$conc,
qinf95 = cred.int[["qinf95"]],
qsup95 = cred.int[["qsup95"]],
Cred.Lim = "Credible limits")
if (adddata) {
plt_3 <- ggplot(data) +
geom_segment(aes(x = conc, xend = conc, y = qinf95, yend = qsup95,
linetype = Conf.Int), data.three,
color = valCols$cols3) +
scale_linetype(name = "") +
theme_minimal()
}
plt_302 <- if (!is.null(spaghetti.CI)) {
ggplot(data) + geom_line(data = dataCIm, aes(x = curv_conc, y = value,
group = variable),
col = "gray", alpha = 0.05)
} else {
ggplot(data) + geom_ribbon(data = data.four, aes(x = conc, ymin = qinf95,
ymax = qsup95),
fill = ribcol, col = NA, alpha = 0.4)
}
plt_32 <- plt_302 +
geom_line(data = data.four, aes(conc, qinf95, color = Cred.Lim),
linetype = cilty, size = cilwd) +
geom_line(data = data.four, aes(conc, qsup95, color = Cred.Lim),
linetype = cilty, size = cilwd) +
scale_color_manual("", values = valCols$cols4) +
theme_minimal()
# plot IC
# final plot
if (!is.null(spaghetti.CI)) {
plt_40 <- ggplot(data) +
geom_line(data = dataCIm, aes(x = curv_conc, y = value, group = variable),
col = "gray", alpha = 0.05)
} else {
plt_40 <- ggplot(data) + geom_ribbon(data = data.four, aes(x = conc,
ymin = qinf95,
ymax = qsup95),
fill = ribcol,
col = NA, alpha = 0.4)
}
plt_401 <- plt_40 +
geom_line(data = data.four, aes(conc, qinf95),
linetype = cilty, size = cilwd, color = valCols$cols4) +
geom_line(data = data.four, aes(conc, qsup95),
linetype = cilty, size = cilwd, color = valCols$cols4) +
geom_line(data = curv_resp, aes(conc, resp),
linetype = fitlty, size = fitlwd, color = valCols$cols2) +
scale_color_discrete(guide = "none") +
ylim(0, 1) +
labs(x = xlab, y = ylab) +
ggtitle(main) + theme_minimal()
if (adddata) {
plt_4 <- plt_401 + geom_point(data = data, aes(transf_conc, resp)) +
geom_segment(aes(x = conc, xend = conc, y = qinf95, yend = qsup95),
data.three, col = valCols$cols3)
} else {
plt_4 <- plt_401
}
return(list(plt_3 = if (adddata) plt_3 else NULL,
plt_32 = plt_32,
plt_4 = plt_4))
}
survFitPlotGG <- function(x,
data_conc, transf_data_conc, data_resp,
curv_conc, curv_resp,
conf.int, cred.int, spaghetti.CI, dataCIm,
xlab, ylab, fitcol, fitlty, fitlwd,
main, addlegend, adddata,
cicol, cilty, cilwd, ribcol) {
if (Sys.getenv("RSTUDIO") == "") {
dev.new() # create a new page plot
# when not use RStudio
}
# dataframes points (data) and curve (curv)
data <- data.frame(conc = data_conc, transf_conc = transf_data_conc,
resp = data_resp, Points = "Observed values")
# colors
valCols <- fCols(data, fitcol, cicol)
if (adddata) {
# points (to create the legend)
plt_1 <- ggplot(data) +
geom_point(data = data, aes(transf_conc, resp, fill = Points),
col = valCols$cols1) + scale_fill_hue("") +
theme_minimal()
}
# curve (to create the legend)
plt_2 <- ggplot(data) +
geom_line(data = curv_resp, aes(conc, resp, colour = Line),
linetype = fitlty, size = fitlwd) +
scale_colour_manual("", values = valCols$cols2) +
theme_minimal()
plt_4 <-
survFitPlotGGCredInt(x, data, curv_resp, conf.int, cred.int, spaghetti.CI,
dataCIm, cilty, cilwd, valCols, fitlty, fitlwd, ribcol,
xlab, ylab, main, adddata)$plt_4
if (addlegend) { # legend yes
# create legends
mylegend_1 <- if (adddata) { legendGgplotFit(plt_1) } else NULL # points legend
mylegend_2 <- legendGgplotFit(plt_2) # mean line legend
plt_5 <- plt_4 + scale_x_continuous(breaks = data$transf_conc,
labels = data$conc)
plt_3 <- survFitPlotGGCredInt(x, data, curv_resp, conf.int, cred.int,
spaghetti.CI, dataCIm, cilty, cilwd,
valCols, fitlty, fitlwd, ribcol, xlab, ylab, main,
adddata)$plt_3
plt_32 <- survFitPlotGGCredInt(x, data, curv_resp, conf.int, cred.int,
spaghetti.CI, dataCIm, cilty, cilwd,
valCols, fitlty, fitlwd, ribcol, xlab, ylab, main,
adddata)$plt_32
mylegend_3 <- if (adddata) { legendGgplotFit(plt_3) } else NULL
mylegend_32 <- legendGgplotFit(plt_32)
if (adddata) {
grid.arrange(plt_5, arrangeGrob(mylegend_1, mylegend_3, mylegend_32,
mylegend_2, nrow = 6), ncol = 2, widths = c(6, 2))
} else {
grid.arrange(plt_5, arrangeGrob(mylegend_32,
mylegend_2, nrow = 6), ncol = 2,
widths = c(6, 2))
}
}
else { # no legend
plt_5 <- plt_4 + scale_x_continuous(breaks = data$transf_conc,
labels = data$conc)
return(plt_5)
}
}