Nothing
#' Plotting method for \code{reproFitTT} objects
#'
#' This is the generic \code{plot} S3 method for the \code{reproFitTT} class.
#' It plots the concentration-effect fit under target time reproduction
#' analysis.
#'
#' The fitted curve represents the \strong{estimated reproduction rate} at the target time
#' as a function of the chemical compound concentration.
#' The function plots 95\% credible intervals for the estimated reproduction
#' rate (by default the grey area around the fitted curve). Typically
#' a good fit is expected to display a large overlap between the two types of 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.
#' It 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{reproFitTT}
#' @param xlab a label for the \eqn{X}-axis, by default \code{Concentration}
#' @param ylab a label for the \eqn{Y}-axis, by default \code{Nb of offspring per ind/day}
#' @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 limits
#' @param cilty line type of the 95 \% credible limits
#' @param cilwd width of the 95 \% credible limits
#' @param ribcol color of the ribbon between lower and upper credible limits.
#' Transparent if \code{NULL}
#' @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{'ggplot'} or \code{'generic'}
#' @param \dots Further arguments to be passed to generic methods
#'
#' @note When \code{style = "generic"}, the function calls the generic function
#' \code{\link[graphics]{plot}}
#' @note When \code{style = "ggplot"}, the function return an object of class
#' \code{ggplot}, see function \code{\link[ggplot2]{ggplot}}
#'
#'
#' @import ggplot2
#' @import grDevices
#' @importFrom gridExtra grid.arrange arrangeGrob
#' @importFrom grid grid.rect gpar
#' @importFrom graphics plot axis legend lines par points polygon
#' segments title
#' @importFrom reshape2 melt
#'
#' @keywords plot
#'
#' @return a plot of class \code{ggplot}
#'
#' @examples
#'
#' # (1) Load the data
#' data(cadmium1)
#'
#' # (2) Create an object of class "reproData"
#' dataset <- reproData(cadmium1)
#'
#' \donttest{
#' # (3) Run the reproFitTT function with the log-logistic gamma-Poisson model
#' out <- reproFitTT(dataset, stoc.part = "gammapoisson",
#' ecx = c(5, 10, 15, 20, 30, 50, 80), quiet = TRUE)
#'
#'
#' # (4) Plot the fitted curve with generic style
#' plot(out, xlab = expression("Concentration in" ~ mu~g.L^{-1}),
#' fitcol = "blue", cicol = "lightblue",
#' main = "Log-logistic response to concentration")
#' }
#'
#' @export
plot.reproFitTT <- function(x,
xlab = "Concentration",
ylab = "Nb of offspring per ind/day",
main = NULL,
fitcol = "orange",
fitlty = 1,
fitlwd = 1,
spaghetti = FALSE,
cicol = "orange",
cilty = 2,
cilwd = 1,
ribcol = "grey70",
addlegend = FALSE,
log.scale = FALSE,
style = "ggplot", ...) {
# plot the fitted curve estimated by reproFitTT
# INPUTS
# - x: reproFitTT 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
# - cilty : type line ci
# - cilwd : width line ci
# - addlegend : boolean
# - log.scale : x log option
# - style : generic ou 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$Nreprocumul / dataTT$Nindtime
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)
cred.int <- reproMeanCredInt(x, display.conc)
spaghetti.CI <- if (spaghetti) { reproSpaghetti(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")
# ylim
ylim_CI <- if (spaghetti) { max(dataCIm$value, cred.int$qsup95)
} else {
max(cred.int$qsup95)
}
if (style == "generic") {
reproFitPlotGenericCredInt(x, dataTT$conc, transf_data_conc, dataTT$resp,
curv_conc, curv_resp,
cred.int, spaghetti.CI, dataCIm,
xlab, ylab, fitcol, fitlty, fitlwd,
main, addlegend,
cicol, cilty, cilwd, ribcol, log.scale, ylim_CI)
}
else if (style == "ggplot") {
reproFitPlotGG(x, dataTT$conc, transf_data_conc, dataTT$resp,
curv_conc, curv_resp,
cred.int, spaghetti.CI, dataCIm,
xlab, ylab, fitcol, fitlty, fitlwd,
main, addlegend,
cicol, cilty, cilwd, ribcol, log.scale, ylim_CI)
}
else stop("Unknown style")
}
#' @importFrom stats quantile rgamma
reproMeanCredInt <- function(fit, x) {
# create the parameters for credible interval for the log logistic model
# INPUT:
# - fit : object of class reproFitTT
# - x : vector of concentrations values (x axis)
# OUTPUT:
# - ci : credible limit
mctot <- do.call("rbind", fit$mcmc)
k <- nrow(mctot)
# parameters
d2 <- mctot[, "d"]
log10b2 <- mctot[, "log10b"]
b2 <- 10^log10b2
log10e2 <- mctot[, "log10e"]
e2 <- 10^log10e2
# quantiles
qinf95 = NULL
q50 = NULL
qsup95 = NULL
# poisson
if (fit$model.label == "P") {
for (i in 1:length(x)) {
theomean <- d2 / (1 + (x[i] / e2)^(b2)) # mean curve
# IC 95%
qinf95[i] <- quantile(theomean, probs = 0.025, na.rm = TRUE)
qsup95[i] <- quantile(theomean, probs = 0.975, na.rm = TRUE)
q50[i] <- quantile(theomean, probs = 0.5, na.rm = TRUE)
}
}
# gamma poisson
else if (fit$model.label == "GP") {
# parameters
log10omega2 <- mctot[, "log10omega"]
omega2 <- 10^(log10omega2)
for (i in 1:length(x)) {
theomean <- d2 / (1 + (x[i] / e2)^(b2)) # mean curve
theo <- rgamma(n = k, shape = theomean / omega2, rate = 1 / omega2)
# IC 95%
qinf95[i] <- quantile(theo, probs = 0.025, na.rm = TRUE)
qsup95[i] <- quantile(theo, probs = 0.975, na.rm = TRUE)
q50[i] <- quantile(theo, probs = 0.5, na.rm = TRUE)
}
}
# values for cred.int
ci <- data.frame(qinf95 = qinf95,
q50 = q50,
qsup95 = qsup95)
return(ci)
}
reproSpaghetti <- function(fit, x) {
mctot <- do.call("rbind", fit$mcmc)
sel <- sample(nrow(mctot))[1:ceiling(nrow(mctot) / 10)]
k <- nrow(mctot[sel,])
# parameters
d2 <- mctot[, "d"][sel]
log10b2 <- mctot[, "log10b"][sel]
b2 <- 10^log10b2
log10e2 <- mctot[, "log10e"][sel]
e2 <- 10^log10e2
if (fit$model.label == "GP") {
log10omega2 <- mctot[, "log10omega"][sel]
omega2 <- 10^(log10omega2)
}
# all theorical
dtheo <- array(data = NA, dim = c(length(x), length(e2)))
if (fit$model.label == "GP") dtheotemp <- dtheo
for (i in 1:length(e2)) {
if (fit$model.label == "P") {
dtheo[, i] <- d2[i] / (1 + (x / e2[i])^(b2[i])) # mean curve
}
else if (fit$model.label == "GP") {
dtheotemp[, i] <- d2[i] / (1 + (x / e2[i])^(b2[i])) # mean curve
dtheo[, i] <- rgamma(n = length(x), shape = dtheotemp[, i] / omega2[i],
rate = 1 / omega2[i])
}
}
dtheof <- as.data.frame(cbind(x, dtheo))
names(dtheof) <- c("conc", paste0("X", 1:length(sel)))
return(dtheof)
}
#' @importFrom epitools pois.exact
reproFitPlotGenericCredInt <- function(x, data_conc, transf_data_conc, data_resp,
curv_conc, curv_resp,
cred.int, spaghetti.CI, dataCIm,
xlab, ylab, fitcol, fitlty, fitlwd,
main, addlegend,
cicol, cilty, cilwd, ribcol, log.scale, ylim_CI) {
# plot the fitted curve estimated by reproFitTT
# 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, ylim_CI + 0.01),
type = "n")
# axis
axis(side = 2, at = pretty(c(0, ylim_CI)))
axis(side = 1,
at = transf_data_conc,
labels = data_conc)
# Plotting the theoretical curve
# cred.int 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 {
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)
# fitted curve
lines(curv_conc, curv_resp[, "resp"], col = fitcol,
lty = fitlty, lwd = fitlwd, type = "l")
# legend
if(addlegend) {
legend("bottomleft",
lty = c(cilty, fitlty),
lwd = c(cilwd, fitlwd),
col = c(cicol, fitcol),
legend = c("Credible limits", "loglogistic"),
bty = "n")
}
}
reproFitPlotGGCredInt <- function(curv_resp, cred.int, spaghetti.CI, dataCIm,
cicol, cilty, cilwd, valCols, fitlty, fitlwd, ribcol,
xlab, ylab, main, ylim_CI) {
# IC
data.three <- data.frame(conc = curv_resp$conc,
qinf95 = cred.int[["qinf95"]],
qsup95 = cred.int[["qsup95"]],
Cred.Lim = "Credible limits")
plt_31 <- if (!is.null(spaghetti.CI)) {
ggplot(data.three) + geom_line(data = dataCIm, aes(x = curv_conc, y = value,
group = variable),
col = "gray", alpha = 0.05)
} else {
ggplot(data.three) + geom_ribbon(data = data.three, aes(x = conc,
ymin = qinf95,
ymax = qsup95),
fill = ribcol, col = NA,
alpha = 0.4)
}
plt_3 <- plt_31 +
geom_line(data = data.three, aes(conc, qinf95, color = Cred.Lim),
linetype = cilty, size = cilwd) +
geom_line(data = data.three, 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.three) +
geom_line(data = dataCIm, aes(x = curv_conc, y = value, group = variable),
col = "gray", alpha = 0.05)
} else {
plt_40 <- ggplot(data.three) + geom_ribbon(data = data.three,
aes(x = conc,
ymin = qinf95,
ymax = qsup95),
fill = ribcol,
col = NA, alpha = 0.4)
}
plt_4 <- plt_40 +
geom_line(data = data.three, aes(conc, qinf95),
linetype = cilty, size = cilwd, color = valCols$cols4) +
geom_line(data = data.three, aes(conc, qsup95),
linetype = cilty, size = cilwd, color = valCols$cols4) +
geom_line(aes(conc, resp), curv_resp,
linetype = fitlty, size = fitlwd, color = valCols$cols2) +
ylim(0, ylim_CI + 0.2) +
labs(x = xlab, y = ylab) +
ggtitle(main) + theme_minimal()
return(list(plt_3 = plt_3,
plt_4 = plt_4))
}
reproFitPlotGG <- function(x, data_conc, transf_data_conc, data_resp,
curv_conc, curv_resp,
cred.int, spaghetti.CI, dataCIm,
xlab, ylab, fitcol, fitlty, fitlwd,
main, addlegend,
cicol, cilty, cilwd, ribcol, log.scale, ylim_CI) {
if (Sys.getenv("RSTUDIO") == "") {
dev.new() # create a new page plot
# when not use RStudio
}
# dataframes points (data) and curve (curv)
# colors
valCols <- fCols(curv_resp, fitcol, cicol)
plt_4 <-
reproFitPlotGGCredInt(curv_resp, cred.int, spaghetti.CI, dataCIm,
cicol, cilty, cilwd, valCols, fitlty, fitlwd, ribcol, xlab,
ylab, main, ylim_CI)$plt_4
if (addlegend) {
# create legends
# curve (to create the legend)
plt_2 <- ggplot(curv_resp) +
geom_line(data = curv_resp, aes(conc, resp, colour = Line),
linetype = fitlty, size = fitlwd) +
scale_color_manual("", values = valCols$cols2) +
theme_minimal()
mylegend_2 <- legendGgplotFit(plt_2) # mean line legend
plt_5 <- plt_4 + scale_x_continuous(breaks = transf_data_conc,
labels = data_conc)
plt_3 <- reproFitPlotGGCredInt(curv_resp, cred.int, spaghetti.CI, dataCIm,
cicol, cilty, cilwd, valCols, fitlty,
fitlwd, ribcol, xlab, ylab, main, ylim_CI)$plt_3
mylegend_3 <- legendGgplotFit(plt_3)
grid.arrange(plt_5, arrangeGrob(mylegend_2, mylegend_3,
nrow = 6), ncol = 2,
widths = c(6, 2))
}
else { # no legend
plt_5 <- plt_4 + scale_x_continuous(breaks = transf_data_conc,
labels = data_conc)
return(plt_5)
}
}
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.