Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • mosaic-software/morse
1 result
Show changes
Commits on Source (3)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 3379 additions and 3357 deletions
.Rbuildignore
View file @ 8c5b023a
......@@ -18,6 +18,7 @@
^morse_*.tar.gz
^Dockerfile
^LICENSE.md
^README.md
.gitignore
^/src/morse.so$
^/src/gutsRedIT.o$
......@@ -25,3 +26,5 @@
.devcontainer
^doc$
^Meta$
^paper$
^notebook/*
DESCRIPTION
View file @ 8c5b023a
Package: morse
Type: Package
Title: Modelling Tools for Reproduction and Survival Data in Ecotoxicology
Version: 3.3.2
Title: Modelling Reproduction and Survival Data in Ecotoxicology
Version: 3.3.3
Encoding: UTF-8
Author:
Virgile Baudrot [cre],
Sandrine Charles [aut],
Marie Laure Delignette-Muller [aut],
Wandrille Duchemin [ctb],
Benoit Goussen [ctb],
Nils Kehrein [ctb],
Guillaume Kon-Kam-King [ctb],
Christelle Lopes [ctb],
Philippe Ruiz [aut],
Alexander Singer [ctb],
Philippe Veber [aut]
Authors@R: c(
person("Virgile", "Baudrot", role = c("aut","cre"), email = "virgile.baudrot@qonfluens.com"),
person("Sandrine", "Charles", role = "aut"),
person("Marie Laure", "Delignette-Muller", role = "aut"),
person("Wandrille", "Duchemin", role = "ctb"),
person("Benoit", "Goussen", role = "ctb"),
person("Nils", "Kehrein", role = "ctb"),
person("Guillaume", "Kon-Kam-King", role = "ctb"),
person("Christelle", "Lopes", role = "ctb"),
person("Philippe", "Ruiz", role = "ctb"),
person("Alexander", "Singer", role = "ctb"),
person("Philippe", "Veber", role = "aut")
)
Maintainer: Virgile Baudrot <virgile.baudrot@qonfluens.com>
URL: https://cran.r-project.org/package=morse
BugReports: https://gitlab.com/qonfluens/model/morse
Description: Tools for ecotoxicologists and regulators dedicated to the
mathematical and statistical modelling of toxicity test data. They use advanced and
innovative methods for a valuable quantitative environmental risk assessment. See our companion paper Baudrot and Charles (2021) <doi:10.21105/joss.03200>,
as well as complementary details in Baudrot et al. (2018) <doi:10.1021/acs.est.7b05464> and Delignette-Muller et al. (2017) <doi:10.1021/acs.est.6b05326>.
URL: https://gitlab.in2p3.fr/mosaic-software/morse
BugReports: https://gitlab.in2p3.fr/mosaic-software/morse/-/issues
Description: Advanced methods for a valuable quantitative environmental risk
assessment using Bayesian inference of survival and reproduction Data. Among
others, it facilitates Bayesian inference of the general unified
threshold model of survival (GUTS). See our companion paper
Baudrot and Charles (2021) <doi:10.21105/joss.03200>,
as well as complementary details in Baudrot et al. (2018)
<doi:10.1021/acs.est.7b05464> and Delignette-Muller et al.
(2017) <doi:10.1021/acs.est.6b05326>.
Depends:
R (>= 3.5.0)
SystemRequirements: JAGS (>= 4.0.0) (see http://mcmc-jags.sourceforge.net)
SystemRequirements: JAGS (>= 4.0.0) (see https://mcmc-jags.sourceforge.io)
Imports:
coda,
deSolve,
......@@ -49,5 +54,5 @@ Suggests:
rmarkdown,
testthat
VignetteBuilder: knitr
RoxygenNote: 7.2.1
RoxygenNote: 7.3.1
LazyData: true
Dockerfile
View file @ 8c5b023a
FROM rocker/verse:4.2.0
# Install JAGS
RUN apt-get update
RUN apt-get install -y jags
# Install R packages
RUN R -e "install.packages(c('roxygen2', 'devtools'), dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('coda',dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('deSolve',dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('epitools',dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('reshape2',dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('rjags',dependencies=TRUE, repos='http://cran.rstudio.com/')"
FROM rocker/verse:4.2.0
# Install JAGS
RUN apt-get update
RUN apt-get install -y jags
# Install R packages
RUN R -e "install.packages(c('roxygen2', 'devtools'), dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('coda',dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('deSolve',dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('epitools',dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('reshape2',dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('rjags',dependencies=TRUE, repos='http://cran.rstudio.com/')"
RUN R -e "install.packages('zoo',dependencies=TRUE, repos='http://cran.rstudio.com/')"
\ No newline at end of file
R/LCx.R
View file @ 8c5b023a
#' Predict \eqn{X}\% Lethal Concentration at the maximum time point (default).
#'
#' Predict median and 95\% credible interval of the x\% Lethal Concentration.
#'
#' When class of \code{object} is \code{survFit}, see \link[=LCx.survFit]{LCx.survFit}.
#'
#' @rdname LCX
#'
#' @param object An object used to select a method
#' @param \dots Further arguments to be passed to generic methods
#'
#' @export
#'
LCx <- function(object, ...){
UseMethod("LCx")
}
#' Predict \eqn{X}\% Lethal Concentration at the maximum time point (default).
#'
#' Predict median and 95\% credible interval of the x\% Lethal Concentration.
#'
#' When class of \code{object} is \code{survFit}, see \link[=LCx.survFit]{LCx.survFit}.
#'
#' @rdname LCX
#'
#' @param object An object used to select a method
#' @param \dots Further arguments to be passed to generic methods
#'
#' @return returns an object of class \code{LCx}.
#'
#' @export
#'
LCx <- function(object, ...){
UseMethod("LCx")
}
R/LCx.survFit.R
View file @ 8c5b023a
#' Predict \eqn{x}\% Lethal Concentration at any specified time point for
#' a \code{survFit} object.
#'
#' The function \code{LCx}, \eqn{x}\% Lethal Concentration (\eqn{LC_x}), is use to compute
#' the dose required to kill \eqn{x}\% of the members of a tested population
#' after a specified test duration (\code{time_LCx}) (default is the maximum
#' time point of the experiment).
#'
#' Mathematical definition of \eqn{x}\% Lethal Concentration at time \eqn{t},
#' denoted \eqn{LC(x,t)}, is:
#'
#' \eqn{S(LC(x,t), t) = S(0, t)*(1- x/100)},
#'
#' where \eqn{S(LC(x,t), t)} is the survival probability at concentration
#' \eqn{LC(x,t)} at time \eqn{t}, and \eqn{S(0,t)} is the survival probability at
#' no concentration (i.e. concentration is \eqn{0}) at time \eqn{t} which
#' reflect the background mortality \eqn{h_b}:
#'
#' \eqn{S(0, t) = exp(-hb* t)}.
#'
#' In the function \code{LCx}, we use the median of \eqn{S(0,t)} to rescale the
#' \eqn{x}\% Lethal Concentration at time \eqn{t}.
#'
#' @rdname LCX
#'
#' @param object An object of class \code{survFit}
#' @param X Percentage of individuals dying (e.g., \eqn{50} for \eqn{LC_{50}}, \eqn{10} for \eqn{LC_{10}}, ...)
#' @param time_LCx A number giving the time at which \eqn{LC_{x}} has to be estimated.
#' If NULL, the latest time point of the experiment is used.
#' @param conc_range A vector of length 2 with minimal and maximal value of the
#' range of concentration. If NULL, the range is
#' define between 0 and the highest tested concentration of the experiment.
#' @param npoints Number of time point in \code{conc_range} between 0 and the maximal concentration. 100 by default.
#' @param \dots Further arguments to be passed to generic methods
#'
#' @return The function returns an object of class \code{LCx}, which is a list
#' with the following information:
#' \item{X_prop}{Survival probability of individuals surviving considering the median
#' of the background mortality (i.e. \eqn{S(0, t)*(1- x/100)})}
#' \item{X_prop_provided}{Survival probability of individuals surviving as provided in arguments (i.e. \eqn{(100-X)/100)}}
#' \item{time_LCx}{A number giving the time at which \eqn{LC_{x}} has to be
#' estimated as provided in arguments or if NULL, the latest time point of the
#' experiment is used.}
#' \item{df_LCx}{A \code{data.frame} with quantiles (median, 2.5\% and 97.5\%)
#' of \eqn{LC_{X}} at time \code{time_LCx} for \eqn{X}\% of individuals}
#' \item{df_dose}{A \code{data.frame} with four columns: \code{concentration}, and median \code{q50} and 95\% credible interval
#' (\code{qinf95} and \code{qsup95}) of the survival probability at time \code{time_LCx}}
#'
#'
#' @examples
#'
#' # (1) Load the data
#' data("propiconazole")
#'
#' # (2) Create an object of class 'survData'
#' dataset <- survData(propiconazole)
#'
#' \dontrun{
#' # (3) Run the survFit function with model_type SD (or IT)
#' out_SD <- survFit(dataset, model_type = "SD")
#'
#' # (4) estimate LC50 at time 4
#' LCx(out_SD, X = 50, time_LCx = 4)
#' }
#'
#' @import zoo
#' @importFrom stats approx
#'
#' @export
#'
LCx.survFit <- function(object,
X,
time_LCx = NULL,
conc_range = NULL,
npoints = 100,
...){
if(is.null(conc_range)){
conc_range = seq(0, max(object$jags.data$conc), length.out = npoints)
} else{
if(length(conc_range) != 2){
stop('conc_range must a vector of length 2 with minimal and maximal value of the range of concentration')
}
conc_range = seq(conc_range[1], conc_range[2], length.out = npoints)
}
if(min(conc_range) != 0){
stop("Minimal value of 'conc_range' must be 0.")
}
if(is.null(time_LCx)){
time_LCx = max(object$jags.data$time)
}
df_dose <- doseResponse_survFitCstExp(x = object, time_LCx = time_LCx, conc_range, npoints)
median_backgroundMortality_Conc0 = dplyr::filter(df_dose, concentration == 0)$q50
X_prop_provided <- (100-X)/100
X_prop <- (100-X)/100*median_backgroundMortality_Conc0
df_LCx <- pointsLCx(df_dose, X_prop)
object_LCx <- list(X_prop = X_prop,
X_prop_provided = X_prop_provided,
time_LCx = time_LCx,
df_LCx = df_LCx,
df_dose = df_dose)
class(object_LCx) <- c("LCx", "list")
return(object_LCx)
}
# dose response curve
#
doseResponse_survFitCstExp <- function(x, time_LCx,
conc_range, npoints){
model_type = x$model_type
# prameters
mctot <- do.call("rbind", x$mcmc)
kd <- 10^mctot[, "kd_log10"]
# "hb" is not in survFit object of morse <v3.2.0
if("hb" %in% colnames(mctot)){
hb <- mctot[, "hb"]
} else{ hb <- 10^mctot[, "hb_log10"] }
# all theorical
k <- 1:length(conc_range)
j <- 1:npoints
if(model_type == "SD"){
if(is.null(time_LCx)){
time_LCx = max(x$jags.data$time)
}
z <- 10^mctot[, "z_log10"]
kk <- 10^mctot[, "kk_log10"]
dtheo <- lapply(k, function(kit) { # conc
Surv_SD(Cw = conc_range[kit],
time = time_LCx,
kk = kk,
kd = kd,
z = z,
hb = hb)
})
}
if(model_type == "IT"){
alpha <- 10^mctot[, "alpha_log10"]
beta <- 10^mctot[, "beta_log10"]
if(is.null(time_LCx)){
time_LCx = max(x$jags.data$time)
}
dtheo <- lapply(k, function(kit) { # concentration pour chaque concentration
Surv_IT_LCx(Cw = conc_range[kit],
time = time_LCx,
kd = kd,
hb = hb,
alpha = alpha,
beta = beta)
})
}
# 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)
df_dose_Resp = data.frame(concentration = conc_range,
q50 = q50,
qinf95 = qinf95,
qsup95 = qsup95)
}
Surv_IT_LCx <- function(Cw, time, kd, hb, alpha, beta)
{
D <- Cw*(1-exp(-kd * time))
S <- exp(-hb * time) * ( 1- 1/(1 + (D/alpha)^(- beta)))
return(S)
}
# points for LCx
#
pointsLCx <- function(df_dose, X_prop){
if(min(df_dose$q50) < X_prop & X_prop < max(df_dose$q50)){
LCX_q50 = approx( df_dose$q50, df_dose$concentration, xout = X_prop, ties = mean)$y
} else {
LCX_q50 = NA
warning(paste("No median for survival probability of", X_prop,
" in the range of concentrations under consideration: [",
min(df_dose$concentration), ";", max(df_dose$concentration), "]"))
}
if(min(df_dose$qinf95) < X_prop & X_prop < max(df_dose$qinf95)){
LCX_qinf95 = approx( df_dose$qinf95, df_dose$concentration, xout = X_prop, ties = mean)$y
} else{
LCX_qinf95 = NA
warning(paste("No 95%inf for survival probability of", X_prop ,
" in the range of concentrations under consideration: [",
min(df_dose$concentration), ";", max(df_dose$concentration), "]"))
}
if(min(df_dose$qsup95) < X_prop & X_prop < max(df_dose$qsup95)){
LCX_qsup95 = approx( df_dose$qsup95, df_dose$concentration, xout = X_prop, ties = mean)$y
} else{
LCX_qsup95 = NA
warning(paste("No 95%sup for survival probability of", X_prop,
" in the range of concentrations under consideration: [",
min(df_dose$concentration), ";", max(df_dose$concentration), "]"))
}
df_LCx <- data.frame(quantile = c("median", "quantile 2.5%", "quantile 97.5%"),
LCx = as.numeric(c(LCX_q50, LCX_qinf95, LCX_qsup95)))
# as.numeric is needed here because if all values are NA, LCx has type logical
return(df_LCx)
}
#' Predict \eqn{x}\% Lethal Concentration at any specified time point for
#' a \code{survFit} object.
#'
#' The function \code{LCx}, \eqn{x}\% Lethal Concentration (\eqn{LC_x}), is use to compute
#' the dose required to kill \eqn{x}\% of the members of a tested population
#' after a specified test duration (\code{time_LCx}) (default is the maximum
#' time point of the experiment).
#'
#' Mathematical definition of \eqn{x}\% Lethal Concentration at time \eqn{t},
#' denoted \eqn{LC(x,t)}, is:
#'
#' \eqn{S(LC(x,t), t) = S(0, t)*(1- x/100)},
#'
#' where \eqn{S(LC(x,t), t)} is the survival probability at concentration
#' \eqn{LC(x,t)} at time \eqn{t}, and \eqn{S(0,t)} is the survival probability at
#' no concentration (i.e. concentration is \eqn{0}) at time \eqn{t} which
#' reflect the background mortality \eqn{h_b}:
#'
#' \eqn{S(0, t) = exp(-hb* t)}.
#'
#' In the function \code{LCx}, we use the median of \eqn{S(0,t)} to rescale the
#' \eqn{x}\% Lethal Concentration at time \eqn{t}.
#'
#' @rdname LCX
#'
#' @param object An object of class \code{survFit}
#' @param X Percentage of individuals dying (e.g., \eqn{50} for \eqn{LC_{50}}, \eqn{10} for \eqn{LC_{10}}, ...)
#' @param time_LCx A number giving the time at which \eqn{LC_{x}} has to be estimated.
#' If NULL, the latest time point of the experiment is used.
#' @param conc_range A vector of length 2 with minimal and maximal value of the
#' range of concentration. If NULL, the range is
#' define between 0 and the highest tested concentration of the experiment.
#' @param npoints Number of time point in \code{conc_range} between 0 and the maximal concentration. 100 by default.
#' @param \dots Further arguments to be passed to generic methods
#'
#' @return The function returns an object of class \code{LCx}, which is a list
#' with the following information:
#' \item{X_prop}{Survival probability of individuals surviving considering the median
#' of the background mortality (i.e. \eqn{S(0, t)*(1- x/100)})}
#' \item{X_prop_provided}{Survival probability of individuals surviving as provided in arguments (i.e. \eqn{(100-X)/100)}}
#' \item{time_LCx}{A number giving the time at which \eqn{LC_{x}} has to be
#' estimated as provided in arguments or if NULL, the latest time point of the
#' experiment is used.}
#' \item{df_LCx}{A \code{data.frame} with quantiles (median, 2.5\% and 97.5\%)
#' of \eqn{LC_{X}} at time \code{time_LCx} for \eqn{X}\% of individuals}
#' \item{df_dose}{A \code{data.frame} with four columns: \code{concentration}, and median \code{q50} and 95\% credible interval
#' (\code{qinf95} and \code{qsup95}) of the survival probability at time \code{time_LCx}}
#'
#'
#' @examples
#'
#' # (1) Load the data
#' data("propiconazole")
#'
#' # (2) Create an object of class 'survData'
#' dataset <- survData(propiconazole)
#'
#' \donttest{
#' # (3) Run the survFit function with model_type SD (or IT)
#' out_SD <- survFit(dataset, model_type = "SD")
#'
#' # (4) estimate LC50 at time 4
#' LCx(out_SD, X = 50, time_LCx = 4)
#' }
#'
#' @import zoo
#' @importFrom stats approx
#'
#' @export
#'
LCx.survFit <- function(object,
X,
time_LCx = NULL,
conc_range = NULL,
npoints = 100,
...){
if(is.null(conc_range)){
conc_range = seq(0, max(object$jags.data$conc), length.out = npoints)
} else{
if(length(conc_range) != 2){
stop('conc_range must a vector of length 2 with minimal and maximal value of the range of concentration')
}
conc_range = seq(conc_range[1], conc_range[2], length.out = npoints)
}
if(min(conc_range) != 0){
stop("Minimal value of 'conc_range' must be 0.")
}
if(is.null(time_LCx)){
time_LCx = max(object$jags.data$time)
}
df_dose <- doseResponse_survFitCstExp(x = object, time_LCx = time_LCx, conc_range, npoints)
median_backgroundMortality_Conc0 = dplyr::filter(df_dose, concentration == 0)$q50
X_prop_provided <- (100-X)/100
X_prop <- (100-X)/100*median_backgroundMortality_Conc0
df_LCx <- pointsLCx(df_dose, X_prop)
object_LCx <- list(X_prop = X_prop,
X_prop_provided = X_prop_provided,
time_LCx = time_LCx,
df_LCx = df_LCx,
df_dose = df_dose)
class(object_LCx) <- c("LCx", "list")
return(object_LCx)
}
# dose response curve
#
doseResponse_survFitCstExp <- function(x, time_LCx,
conc_range, npoints){
model_type = x$model_type
# parameters
mctot <- do.call("rbind", x$mcmc)
kd <- 10^mctot[, "kd_log10"]
# "hb" is not in survFit object of morse <v3.2.0
if("hb" %in% colnames(mctot)){
hb <- mctot[, "hb"]
} else{ hb <- 10^mctot[, "hb_log10"] }
# all theorical
k <- 1:length(conc_range)
j <- 1:npoints
if(model_type == "SD"){
if(is.null(time_LCx)){
time_LCx = max(x$jags.data$time)
}
z <- 10^mctot[, "z_log10"]
kk <- 10^mctot[, "kk_log10"]
dtheo <- lapply(k, function(kit) { # conc
Surv_SD(Cw = conc_range[kit],
time = time_LCx,
kk = kk,
kd = kd,
z = z,
hb = hb)
})
}
if(model_type == "IT"){
alpha <- 10^mctot[, "alpha_log10"]
beta <- 10^mctot[, "beta_log10"]
if(is.null(time_LCx)){
time_LCx = max(x$jags.data$time)
}
dtheo <- lapply(k, function(kit) { # concentration pour chaque concentration
Surv_IT_LCx(Cw = conc_range[kit],
time = time_LCx,
kd = kd,
hb = hb,
alpha = alpha,
beta = beta)
})
}
# 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)
df_dose_Resp = data.frame(concentration = conc_range,
q50 = q50,
qinf95 = qinf95,
qsup95 = qsup95)
}
Surv_IT_LCx <- function(Cw, time, kd, hb, alpha, beta)
{
D <- Cw*(1-exp(-kd * time))
S <- exp(-hb * time) * ( 1- 1/(1 + (D/alpha)^(- beta)))
return(S)
}
# points for LCx
#
pointsLCx <- function(df_dose, X_prop){
if(min(df_dose$q50) < X_prop & X_prop < max(df_dose$q50)){
LCX_q50 = approx( df_dose$q50, df_dose$concentration, xout = X_prop, ties = mean)$y
} else {
LCX_q50 = NA
warning(paste("No median for survival probability of", X_prop,
" in the range of concentrations under consideration: [",
min(df_dose$concentration), ";", max(df_dose$concentration), "]"))
}
if(min(df_dose$qinf95) < X_prop & X_prop < max(df_dose$qinf95)){
LCX_qinf95 = approx( df_dose$qinf95, df_dose$concentration, xout = X_prop, ties = mean)$y
} else{
LCX_qinf95 = NA
warning(paste("No 95%inf for survival probability of", X_prop ,
" in the range of concentrations under consideration: [",
min(df_dose$concentration), ";", max(df_dose$concentration), "]"))
}
if(min(df_dose$qsup95) < X_prop & X_prop < max(df_dose$qsup95)){
LCX_qsup95 = approx( df_dose$qsup95, df_dose$concentration, xout = X_prop, ties = mean)$y
} else{
LCX_qsup95 = NA
warning(paste("No 95%sup for survival probability of", X_prop,
" in the range of concentrations under consideration: [",
min(df_dose$concentration), ";", max(df_dose$concentration), "]"))
}
df_LCx <- data.frame(quantile = c("median", "quantile 2.5%", "quantile 97.5%"),
LCx = as.numeric(c(LCX_q50, LCX_qinf95, LCX_qsup95)))
# as.numeric is needed here because if all values are NA, LCx has type logical
return(df_LCx)
}
R/MFx.R
View file @ 8c5b023a
#' Predict the Multiplication Factor leading to x\% of reduction in survival
#' at a specific time.
#'
#' Generic method for \code{MFx}, a function denoted \eqn{MF(x,t)} for
#' \eqn{x}\% Multiplication Factor at time \eqn{t}.
#'
#' When class of \code{object} is \code{survFit}, see \link[=MFx.survFit]{MFx.survFit}.
#'
#' @param object An object used to select a method
#' @param \dots Further arguments to be passed to generic methods
#'
#' @export
#'
MFx <- function(object, ...){
UseMethod("MFx")
}
#' Predict the Multiplication Factor leading to x\% of reduction in survival
#' at a specific time.
#'
#' Generic method for \code{MFx}, a function denoted \eqn{MF(x,t)} for
#' \eqn{x}\% Multiplication Factor at time \eqn{t}.
#'
#' When class of \code{object} is \code{survFit}, see \link[=MFx.survFit]{MFx.survFit}.
#'
#' @param object An object used to select a method
#' @param \dots Further arguments to be passed to generic methods
#'
#' @return returns an object of class \code{MFx}
#'
#' @export
#'
MFx <- function(object, ...){
UseMethod("MFx")
}
R/MFx.survFit.R
View file @ 8c5b023a
This diff is collapsed.
R/MFx_ode.survFit.R
View file @ 8c5b023a
This diff is collapsed.
R/modelData.survDataCstExp.R
View file @ 8c5b023a
#' Create a data set to analyse a \code{survDataCstExp} object.
#'
#' @param x An object of class \code{survData}
#' @param model_type TKTD GUTS model type ('SD' or 'IT')
#'
#'
modelData.survDataCstExp <- function(x, model_type = NULL){
## 1. Gather replicate when there is the same constante concentration
x_gather <- gather_survDataCstExp(x)
## 2. Creation of additional variable require in the Bayesian model
x_dev <- addVariable_survDataCstExp(x_gather)
##
## ====== Construction of modelData
##
### return priors for model
priorsData <- priors_survData(x = x, model_type = model_type)
### observations
dataList <- list(
replicate = x_dev$replicate,
time = x_dev$time,
conc = x_dev$conc,
Nsurv = x_dev$Nsurv,
Nprec = x_dev$Nprec
)
### parameters
dataList$n_data <- nrow(x_dev)
### other parameters specific to model IT
if(model_type == "IT"){
dataList$i_prec <- x_dev$i_prec
dataList$replicate_ID <- x_dev$replicate_ID
dataList$time_ID <- x_dev$time_ID
}
if(model_type == "SD"){
dataList$tprec <- x_dev$tprec
dataList$i_prec <- x_dev$i_prec
}
##
## =========== Object to return
##
OUT_modelDATA <- list(dataList = dataList,
priorsList = priorsData$priorsList,
priorsMinMax = priorsData$priorsMinMax)
return(OUT_modelDATA)
}
# Gather replicates with the same concentration
#
# @param x An object of class \code{survData}
#
# @return A dataframe
#
gather_survDataCstExp <- function(x){
bool_checkTimeReplicate <- checkTimeReplicate(x)
if( bool_checkTimeReplicate ){
### Sum Nsurv data for each (conc, time) couple
x_dev <- x %>%
dplyr::group_by(conc, time) %>%
dplyr::summarise(Nsurv = sum(Nsurv)) %>%
# concate replicate in the same replicate using factor (conc)
dplyr::mutate(replicate = as.character(conc)) %>%
as.data.frame()
} else{
x_dev <- x
}
return(x_dev)
}
# Check the same number of (time, replicate) for a common concentration
#
# @param x
checkTimeReplicate <- function(x){
df_checkTimeReplicate <- x %>%
dplyr::group_by(conc, time) %>%
dplyr::summarise(nbrReplicate_ConcTime = n_distinct(replicate)) %>%
dplyr::group_by(conc) %>%
dplyr::summarise(nbrReplicate_uniqueConc = length(unique(nbrReplicate_ConcTime)))
return(all(df_checkTimeReplicate$nbrReplicate_uniqueConc == 1))
}
# Add variables for Bayesian fitting
#
# @param x An object of class \code{data.frame}
#
# @return A dataframe
#
addVariable_survDataCstExp <- function(x){
## Creation of additional variable:
## - tprec: previous time
## - Nprec: previous number of survivors
## - time_ID: identification of row number inside a group
## - i_row: identification of row number (for every group)
## - i_prec: identification of previous row number (for every group) except when time_ID (in group) is 1
x_dev <- x %>%
# Add an indice of replicate:
# dplyr::mutate(replicate_ID = group_indices(., .dots = "replicate")) %>%
dplyr::group_by(replicate) %>%
dplyr::mutate(replicate_ID = cur_group_id()) %>%
dplyr::arrange(replicate, time) %>%
dplyr::mutate(tprec = ifelse(time == 0, time, dplyr::lag(time)),
Nprec = ifelse( time == 0, Nsurv, dplyr::lag(Nsurv) )) %>%
# remove time = 0 for the analysis
# dplyr::filter(time != 0) %>%
dplyr::mutate(time_ID = row_number()) %>%
dplyr::ungroup() %>%
dplyr::arrange(replicate, time) %>%
dplyr::mutate(i_row = row_number(),
i_prec = ifelse(time_ID == 1, i_row, dplyr::lag(i_row)))
return(x_dev)
}
#' Create a data set to analyse a \code{survDataCstExp} object.
#'
#' @param x An object of class \code{survData}
#' @param model_type TKTD GUTS model type ('SD' or 'IT')
#'
#' @return A list for parameterization of priors for Bayesian inference.
#'
#'
modelData.survDataCstExp <- function(x, model_type = NULL){
## 1. Gather replicate when there is the same constante concentration
x_gather <- gather_survDataCstExp(x)
## 2. Creation of additional variable require in the Bayesian model
x_dev <- addVariable_survDataCstExp(x_gather)
##
## ====== Construction of modelData
##
### return priors for model
priorsData <- priors_survData(x = x, model_type = model_type)
### observations
dataList <- list(
replicate = x_dev$replicate,
time = x_dev$time,
conc = x_dev$conc,
Nsurv = x_dev$Nsurv,
Nprec = x_dev$Nprec
)
### parameters
dataList$n_data <- nrow(x_dev)
### other parameters specific to model IT
if(model_type == "IT"){
dataList$i_prec <- x_dev$i_prec
dataList$replicate_ID <- x_dev$replicate_ID
dataList$time_ID <- x_dev$time_ID
}
if(model_type == "SD"){
dataList$tprec <- x_dev$tprec
dataList$i_prec <- x_dev$i_prec
}
##
## =========== Object to return
##
OUT_modelDATA <- list(dataList = dataList,
priorsList = priorsData$priorsList,
priorsMinMax = priorsData$priorsMinMax)
return(OUT_modelDATA)
}
# Gather replicates with the same concentration
#
# @param x An object of class \code{survData}
#
# @return A dataframe
#
gather_survDataCstExp <- function(x){
bool_checkTimeReplicate <- checkTimeReplicate(x)
if( bool_checkTimeReplicate ){
### Sum Nsurv data for each (conc, time) couple
x_dev <- x %>%
dplyr::group_by(conc, time) %>%
dplyr::summarise(Nsurv = sum(Nsurv)) %>%
# concate replicate in the same replicate using factor (conc)
dplyr::mutate(replicate = as.character(conc)) %>%
as.data.frame()
} else{
x_dev <- x
}
return(x_dev)
}
# Check the same number of (time, replicate) for a common concentration
#
# @param x
checkTimeReplicate <- function(x){
df_checkTimeReplicate <- x %>%
dplyr::group_by(conc, time) %>%
dplyr::summarise(nbrReplicate_ConcTime = n_distinct(replicate)) %>%
dplyr::group_by(conc) %>%
dplyr::summarise(nbrReplicate_uniqueConc = length(unique(nbrReplicate_ConcTime)))
return(all(df_checkTimeReplicate$nbrReplicate_uniqueConc == 1))
}
# Add variables for Bayesian fitting
#
# @param x An object of class \code{data.frame}
#
# @return A dataframe
#
addVariable_survDataCstExp <- function(x){
## Creation of additional variable:
## - tprec: previous time
## - Nprec: previous number of survivors
## - time_ID: identification of row number inside a group
## - i_row: identification of row number (for every group)
## - i_prec: identification of previous row number (for every group) except when time_ID (in group) is 1
x_dev <- x %>%
# Add an indice of replicate:
# dplyr::mutate(replicate_ID = group_indices(., .dots = "replicate")) %>%
dplyr::group_by(replicate) %>%
dplyr::mutate(replicate_ID = cur_group_id()) %>%
dplyr::arrange(replicate, time) %>%
dplyr::mutate(tprec = ifelse(time == 0, time, dplyr::lag(time)),
Nprec = ifelse( time == 0, Nsurv, dplyr::lag(Nsurv) )) %>%
# remove time = 0 for the analysis
# dplyr::filter(time != 0) %>%
dplyr::mutate(time_ID = row_number()) %>%
dplyr::ungroup() %>%
dplyr::arrange(replicate, time) %>%
dplyr::mutate(i_row = row_number(),
i_prec = ifelse(time_ID == 1, i_row, dplyr::lag(i_row)))
return(x_dev)
}
R/modelData.survDataVarExp.R
View file @ 8c5b023a
......@@ -6,6 +6,8 @@
#' interpolation (comprise between time to compute and fitting accuracy)
#' @param \dots Further arguments to be passed to generic methods
#'
#' @return A list for parameterization of priors for Bayesian inference.
#'
modelData.survDataVarExp <- function(x,
model_type = NULL,
extend_time = 100, ...){
......@@ -16,7 +18,8 @@ modelData.survDataVarExp <- function(x,
## - Nprec: previous number of survivors
## - time_ID_red: identification of row number inside a group
## - i_row: identification of row number (for every group)
## - i_prec: identification of previous row number (for every group) exept when time_ID_red (in group) is 1
## - i_prec: identification of previous row number (for every group) except
## when time_ID_red (in group) is 1
x_interpolate <- survData_interpolate(x, extend_time = extend_time) %>%
dplyr::arrange(replicate, time)
......
R/morse.R
View file @ 8c5b023a
This diff is collapsed.
R/msgTable.R
View file @ 8c5b023a
msgTableCreate <- function() {
t <- data.frame(id = character(0), msg = character(0), stringsAsFactors = FALSE)
class(t) <- c("msgTable", "data.frame")
t
}
msgTableAppend <- function(...) {
u <- rbind(...)
class(u) <- c("msgTable", "data.frame")
u
}
msgTableAdd <- function(t, id, msg) {
newlines <- data.frame(id = id, msg = msg, stringsAsFactors = FALSE)
msgTableAppend(t, newlines)
}
msgTableSingleton <- function(id, msg) {
msgTableAdd(msgTableCreate(), id, msg)
}
msgTableIsEmpty <- function(x)
dim(x)[1] == 0
#' @export
print.msgTable <- function(x, ...) {
if (msgTableIsEmpty(x)) {
cat("Correct format\n")
}
else {
cat("Message(s):\n")
for (m in x$msg) {
cat(paste("\t",m,"\n",sep=""))
}
}
}
msgTableCreate <- function() {
t <- data.frame(id = character(0), msg = character(0), stringsAsFactors = FALSE)
class(t) <- c("msgTable", "data.frame")
t
}
msgTableAppend <- function(...) {
u <- rbind(...)
class(u) <- c("msgTable", "data.frame")
u
}
msgTableAdd <- function(t, id, msg) {
newlines <- data.frame(id = id, msg = msg, stringsAsFactors = FALSE)
msgTableAppend(t, newlines)
}
msgTableSingleton <- function(id, msg) {
msgTableAdd(msgTableCreate(), id, msg)
}
msgTableIsEmpty <- function(x)
dim(x)[1] == 0
#' Print \code{msgTables} objects
#'
#' Print in the REPL the \code{msgTables}
#'
#' @param x an object of class \code{msgTables}
#' @param \dots Further arguments to be passed to generic methods
#'
#' @return Print in the REPL the \code{msgTables}
#'
#' @export
print.msgTable <- function(x, ...) {
if (msgTableIsEmpty(x)) {
cat("Correct format\n")
}
else {
cat("Message(s):\n")
for (m in x$msg) {
cat(paste("\t",m,"\n",sep=""))
}
}
}
R/plot.LCx.R
View file @ 8c5b023a
#' Plotting method for \code{LCx} objects
#'
#' This is the generic \code{plot} S3 method for the
#' \\code{LCx} class. It plots the survival probability as a function of concentration.
#'
#'
#' @param x An object of class \code{LCx}.
#' @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{Survival probability median and 95 CI}.
#' @param main A main title for the plot.
#' @param subtitle A subtitle for the plot
#' @param \dots Further arguments to be passed to generic methods.
#'
#' @keywords plot
#'
#' @examples
#'
#' # (1) Load the data
#' data("propiconazole")
#'
#' # (2) Create an object of class 'survData'
#' dataset <- survData(propiconazole)
#'
#' \dontrun{
#' # (3) Run the survFit function with model_type SD (or IT)
#' out_SD <- survFit(dataset, model_type = "SD")
#'
#' # (4) estimate LC50 at time 4
#' LCx_SD <- LCx(out_SD, X = 50, time_LCx = 4)
#'
#' # (5) plot the object of class 'LCx'
#' plot(LCx_SD)
#' }
#'
#' @export
#'
#'
#'
plot.LCx <- function(x,
xlab = "Concentration",
ylab = "Survival probability \n median and 95 CI",
main = NULL,
subtitle = NULL,
...){
df_dose <- x$df_dose
df_LCx <- x$df_LCx
X_prop <- x$X_prop
X_prop_provided <- x$X_prop_provided
time_LCx <- x$time_LCx
# Check if df_LCx are not all NA:
if(all(is.na(df_LCx$LCx))){
warning(paste0("No LCx can be computed at X=", 100-X_prop_provided*100, " and time_LCx=", time_LCx,
"\nSee the grey dotted line on the graph does not cross zone of credible interval.",
"\nUse LCx function 'LCx' with other values for arguments 'time_LCx' (default is the maximum time of the experimental data),
and/or other 'X'."))
} else{
legend.point = data.frame(
x.pts = df_LCx$LCx,
y.pts = rep(X_prop,3),
pts.leg = c(paste("median: ", round(df_LCx$LCx[1],digits = 2)),
paste("quantile 2.5%: ", round(df_LCx$LCx[2],digits = 2)),
paste("quantile 97.5%: ", round(df_LCx$LCx[3],digits = 2)))
)
}
if(is.null(main)){
main <- paste("Concentration-response curve: LC", 100 - X_prop_provided*100, " at time", time_LCx)
}
LCx_plt <- ggplot() + theme_minimal() +
theme(legend.position = "top",
legend.title = element_blank())+
scale_y_continuous(limits = c(0,1)) +
labs(title = main,
subtitle = subtitle,
x = xlab,
y = ylab) +
geom_ribbon(data = df_dose,
aes(x = concentration, ymin = qinf95, ymax = qsup95), fill = "lightgrey") +
geom_line(data = df_dose,
aes( x = concentration, y = q50), color ="orange") +
geom_hline(yintercept = X_prop, col="grey70", linetype=2)
# LCx points
if(!all(is.na(df_LCx$LCx))){
LCx_plt <- LCx_plt +
geom_point(data = legend.point,
aes(x = x.pts, y=y.pts,
color= pts.leg))+
scale_color_manual(values=c("orange", "black", "black"))
}
return(LCx_plt)
}
#' Plotting method for \code{LCx} objects
#'
#' This is the generic \code{plot} S3 method for the
#' \\code{LCx} class. It plots the survival probability as a function of concentration.
#'
#'
#' @param x An object of class \code{LCx}.
#' @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{Survival probability median and 95 CI}.
#' @param main A main title for the plot.
#' @param subtitle A subtitle for the plot
#' @param \dots Further arguments to be passed to generic methods.
#'
#' @keywords plot
#'
#' @return a plot of class \code{ggplot}
#'
#' @examples
#'
#' # (1) Load the data
#' data("propiconazole")
#'
#' # (2) Create an object of class 'survData'
#' dataset <- survData(propiconazole)
#'
#' \donttest{
#' # (3) Run the survFit function with model_type SD (or IT)
#' out_SD <- survFit(dataset, model_type = "SD")
#'
#' # (4) estimate LC50 at time 4
#' LCx_SD <- LCx(out_SD, X = 50, time_LCx = 4)
#'
#' # (5) plot the object of class 'LCx'
#' plot(LCx_SD)
#' }
#'
#' @export
#'
#'
#'
plot.LCx <- function(x,
xlab = "Concentration",
ylab = "Survival probability \n median and 95 CI",
main = NULL,
subtitle = NULL,
...){
df_dose <- x$df_dose
df_LCx <- x$df_LCx
X_prop <- x$X_prop
X_prop_provided <- x$X_prop_provided
time_LCx <- x$time_LCx
# Check if df_LCx are not all NA:
if(all(is.na(df_LCx$LCx))){
warning(paste0("No LCx can be computed at X=", 100-X_prop_provided*100, " and time_LCx=", time_LCx,
"\nSee the grey dotted line on the graph does not cross zone of credible interval.",
"\nUse LCx function 'LCx' with other values for arguments 'time_LCx' (default is the maximum time of the experimental data),
and/or other 'X'."))
} else{
legend.point = data.frame(
x.pts = df_LCx$LCx,
y.pts = rep(X_prop,3),
pts.leg = c(paste("median: ", round(df_LCx$LCx[1],digits = 2)),
paste("quantile 2.5%: ", round(df_LCx$LCx[2],digits = 2)),
paste("quantile 97.5%: ", round(df_LCx$LCx[3],digits = 2)))
)
}
if(is.null(main)){
main <- paste("Concentration-response curve: LC", 100 - X_prop_provided*100, " at time", time_LCx)
}
LCx_plt <- ggplot() + theme_minimal() +
theme(legend.position = "top",
legend.title = element_blank())+
scale_y_continuous(limits = c(0,1)) +
labs(title = main,
subtitle = subtitle,
x = xlab,
y = ylab) +
geom_ribbon(data = df_dose,
aes(x = concentration, ymin = qinf95, ymax = qsup95), fill = "lightgrey") +
geom_line(data = df_dose,
aes( x = concentration, y = q50), color ="orange") +
geom_hline(yintercept = X_prop, col="grey70", linetype=2)
# LCx points
if(!all(is.na(df_LCx$LCx))){
LCx_plt <- LCx_plt +
geom_point(data = legend.point,
aes(x = x.pts, y=y.pts,
color= pts.leg))+
scale_color_manual(values=c("orange", "black", "black"))
}
return(LCx_plt)
}
R/plot.MFx.R
View file @ 8c5b023a
#' Plotting method for \code{MFx} objects
#'
#' This is the generic \code{plot} S3 method for the
#' \code{MFx} class. It plots the survival probability as a function of
#' the multiplication factor applied or as a function of time.
#'
#'
#' @param x An object of class \code{MFx}.
#' @param x_variable A character to define the variable for the \eqn{X}-axis,
#' either \code{"MFx"} or \code{"Time"}. The default is \code{"MFx"}.
#' @param xlab A label for the \eqn{X}-axis, by default \code{NULL} and depend on the
#' argument \code{x_variable}.
#' @param ylab A label for the \eqn{Y}-axis, by default \code{Survival probability median and 95 CI}.
#' @param main A main title for the plot.
#' @param log_scale If \code{TRUE}, the x-axis is log-scaled. Default is \code{FALSE}.
#' @param ncol An interger for the number of columns when several panels are plotted.
#' @param \dots Further arguments to be passed to generic methods.
#'
#' @keywords plot
#'
#' @examples
#'
#' # (1) Load the data
#' data("propiconazole")
#'
#' # (2) Create an object of class 'survData'
#' dataset <- survData(propiconazole)
#'
#' \dontrun{
#' # (3) Run the survFit function with model_type SD (or IT)
#' out_SD <- survFit(dataset, model_type = "SD")
#'
#'# (4) data to predict
#' data_4prediction <- data.frame(time = 1:10, conc = c(0,0.5,3,3,0,0,0.5,3,1.5,0))
#'
#' # (5) estimate MF for 30% reduction of survival at time 4
#' MFx_SD_30.4 <- MFx(out_SD, data_predict = data_4prediction , X = 30, time_MFx = 4)
#'
#' # (6) plot the object of class 'MFx'
#' plot(MFx_SD_30.4)
#'
#' # (6bis) plot with log-scale of x-axis
#' plot(MFx_SD_30.4, log_scale = TRUE)
#'
#' # (6ter) plot with "Time" as the x-axis
#' plot(MFx_SD_30.4, x_variable = "Time")
#'
#' # (7) plot when X = NULL and along a MFx_range from 5 to 10:
#' MFx_SD_range <- MFx(out_SD, data_predict = data_4prediction ,
#' X = NULL, time_MFx = 4, MFx_range = seq(5, 10, length.out = 50))
#' plot(MFx_SD_range)
#' plot(MFx_SD_range, x_variable = "Time", ncol = 10)
#' }
#'
#' @export
#'
#'
#'
plot.MFx <- function(x,
x_variable = "MFx", # other option is "Time"
xlab = NULL,
ylab = "Survival probability \n median and 95 CI",
main = NULL,
log_scale = FALSE,
ncol = 3,
...){
# definition of xlab and check x_variable
if(is.null(xlab)){
if(x_variable == "MFx"){
xlab = "Multiplication Factor"
} else if(x_variable == "Time"){
xlab = "Time"
} else stop("the argument 'x_variable' must be 'MFx' or 'Time'. The default is 'MFx'.")
}
MFx_plt <- ggplot() + theme_minimal() +
theme(legend.position = "top",
legend.title = element_blank())+
scale_y_continuous(limits = c(0,1))
if(x_variable == "MFx"){
if(is.null(main)){
main <- paste("Multiplication Factor for MF", x$X_prop_provided*100, "% at time", x$time_MFx)
}
if(is.null(x$X_prop)) main <- paste("Survival over [", min(x$MFx_tested), ",", max(x$MFx_tested), "] MF range at time", x$time_MFx)
MFx_plt <- MFx_plt +
scale_color_manual(values=c("orange", "black", "black")) +
labs(title = main,
x = xlab,
y = ylab) +
geom_ribbon(data = x$df_dose,
aes(x = MFx, ymin = qinf95, ymax = qsup95),
fill = "lightgrey") +
geom_point(data = x$df_dose,
aes(x = MFx, y = q50), color = "orange", shape = 4) +
geom_line(data = x$df_dose,
aes(x = MFx, y = q50), color = "orange")
if(!is.null(x$X_prop)){
MFx_plt <- MFx_plt +
geom_point(data = dplyr::filter(x$df_dose, id == "q50"),
aes(x = MFx, y = q50), color = "orange", shape = 4) +
geom_point(data = dplyr::filter(x$df_dose, id == "qinf95"),
aes(x = MFx, y = qinf95), color = "grey", shape = 4) +
geom_point(data = dplyr::filter(x$df_dose, id == "qsup95"),
aes(x = MFx, y = qsup95), color = "grey", shape = 4)
legend.point = data.frame(
x.pts = x$df_MFx$MFx,
y.pts = rep(x$X_prop, 3),
pts.leg = c(paste("median: ", round(x$df_MFx$MFx[1],digits = 2)),
paste("quantile 2.5%: ", round(x$df_MFx$MFx[2],digits = 2)),
paste("quantile 97.5%: ", round(x$df_MFx$MFx[3],digits = 2)))
)
MFx_plt <- MFx_plt +
geom_hline(yintercept = x$X_prop, col="grey70", linetype=2) +
geom_point(data = legend.point,
aes(x = x.pts, y = y.pts, color = pts.leg))
warning("This is not an error message:
Just take into account that MFx has been estimated with a binary
search using the 'accuracy' argument. Cross points indicate the
position of evaluated time series. To improve the shape of the curve, you
can use X = NULL, and compute time series around the median MFx, with the
vector `MFx_range`.")
}
if(log_scale == TRUE){
MFx_plt <- MFx_plt +
scale_x_log10()
}
}
if(x_variable == "Time"){
# Plot
if(is.null(main)) main <- paste("Survival over time. Multiplication Factor of", x$X_prop_provided*100, "percent")
if(is.null(x$X_prop)) main <- paste("Survival over time")
MFx = x$MFx_tested
k <- 1:length(x$ls_predict)
#
# Make a dataframe with quantile of all generated time series
#
ls_predict_quantile <- lapply(k, function(kit){
df_quantile <- x$ls_predict[[kit]]$df_quantile
df_quantile$MFx <- rep(MFx[[kit]], nrow(x$ls_predict[[kit]]$df_quantile))
return(df_quantile)
})
predict_MFx_quantile <- do.call("rbind", ls_predict_quantile)
if(!is.null(x$X_prop)){
initial_predict <- dplyr::filter(predict_MFx_quantile, MFx == 1)
final_predict <- dplyr::filter(predict_MFx_quantile, MFx == x$df_MFx$MFx[1])
y_arrow <- as.numeric(dplyr::filter(initial_predict, time == x$time_MFx)$q50)
yend_arrow <- as.numeric(dplyr::filter(final_predict, time == x$time_MFx)$q50)
MFx_plt <- MFx_plt +
labs(title = main,
x = xlab,
y = ylab) +
# final predict
geom_ribbon(data = final_predict,
aes(x = time, ymin = qinf95, ymax = qsup95),
fill = "grey30", alpha = 0.4) +
geom_line(data = final_predict,
aes(x = time, y = q50),
col="orange", size = 1) +
# initial predict
geom_ribbon(data = initial_predict,
aes(x = time, ymin = qinf95, ymax = qsup95),
fill = "grey60", alpha = 0.4) +
geom_line(data = initial_predict,
aes(x = time, y = q50),
col="orange", size = 0.3) +
# arrow
geom_segment(aes(x = x$time_MFx, y = y_arrow,
xend = x$time_MFx, yend = yend_arrow),
arrow = arrow(length = unit(0.2,"cm")))
}
if(is.null(x$X_prop)){
MFx_plt <- MFx_plt +
labs(title = main,
x = xlab,
y = ylab) +
geom_ribbon(data = predict_MFx_quantile,
aes(x = time, ymin = qinf95, ymax = qsup95),
fill = "grey30", alpha = 0.4) +
geom_line(data = predict_MFx_quantile,
aes(x = time, y = q50),
col="orange", size = 1) +
facet_wrap( ~ round(MFx, digits = 2), ncol = ncol)
}
}
return(MFx_plt)
}
#' Plotting method for \code{MFx} objects
#'
#' This is the generic \code{plot} S3 method for the
#' \code{MFx} class. It plots the survival probability as a function of
#' the multiplication factor applied or as a function of time.
#'
#'
#' @param x An object of class \code{MFx}.
#' @param x_variable A character to define the variable for the \eqn{X}-axis,
#' either \code{"MFx"} or \code{"Time"}. The default is \code{"MFx"}.
#' @param xlab A label for the \eqn{X}-axis, by default \code{NULL} and depend on the
#' argument \code{x_variable}.
#' @param ylab A label for the \eqn{Y}-axis, by default \code{Survival probability median and 95 CI}.
#' @param main A main title for the plot.
#' @param log_scale If \code{TRUE}, the x-axis is log-scaled. Default is \code{FALSE}.
#' @param ncol An interger for the number of columns when several panels are plotted.
#' @param \dots Further arguments to be passed to generic methods.
#'
#' @keywords plot
#'
#' @return a plot of class \code{ggplot}
#'
#' @examples
#'
#' # (1) Load the data
#' data("propiconazole")
#'
#' # (2) Create an object of class 'survData'
#' dataset <- survData(propiconazole)
#'
#' \donttest{
#' # (3) Run the survFit function with model_type SD (or IT)
#' out_SD <- survFit(dataset, model_type = "SD")
#'
#'# (4) data to predict
#' data_4prediction <- data.frame(time = 1:10, conc = c(0,0.5,3,3,0,0,0.5,3,1.5,0))
#'
#' # (5) estimate MF for 30% reduction of survival at time 4
#' MFx_SD_30.4 <- MFx(out_SD, data_predict = data_4prediction , X = 30, time_MFx = 4)
#'
#' # (6) plot the object of class 'MFx'
#' plot(MFx_SD_30.4)
#'
#' # (6bis) plot with log-scale of x-axis
#' plot(MFx_SD_30.4, log_scale = TRUE)
#'
#' # (6ter) plot with "Time" as the x-axis
#' plot(MFx_SD_30.4, x_variable = "Time")
#'
#' # (7) plot when X = NULL and along a MFx_range from 5 to 10:
#' MFx_SD_range <- MFx(out_SD, data_predict = data_4prediction ,
#' X = NULL, time_MFx = 4, MFx_range = seq(5, 10, length.out = 50))
#' plot(MFx_SD_range)
#' plot(MFx_SD_range, x_variable = "Time", ncol = 10)
#' }
#'
#' @export
#'
#'
#'
plot.MFx <- function(x,
x_variable = "MFx", # other option is "Time"
xlab = NULL,
ylab = "Survival probability \n median and 95 CI",
main = NULL,
log_scale = FALSE,
ncol = 3,
...){
# definition of xlab and check x_variable
if(is.null(xlab)){
if(x_variable == "MFx"){
xlab = "Multiplication Factor"
} else if(x_variable == "Time"){
xlab = "Time"
} else stop("the argument 'x_variable' must be 'MFx' or 'Time'. The default is 'MFx'.")
}
MFx_plt <- ggplot() + theme_minimal() +
theme(legend.position = "top",
legend.title = element_blank())+
scale_y_continuous(limits = c(0,1))
if(x_variable == "MFx"){
if(is.null(main)){
main <- paste("Multiplication Factor for MF", x$X_prop_provided*100, "% at time", x$time_MFx)
}
if(is.null(x$X_prop)) main <- paste("Survival over [", min(x$MFx_tested), ",", max(x$MFx_tested), "] MF range at time", x$time_MFx)
MFx_plt <- MFx_plt +
scale_color_manual(values=c("orange", "black", "black")) +
labs(title = main,
x = xlab,
y = ylab) +
geom_ribbon(data = x$df_dose,
aes(x = MFx, ymin = qinf95, ymax = qsup95),
fill = "lightgrey") +
geom_point(data = x$df_dose,
aes(x = MFx, y = q50), color = "orange", shape = 4) +
geom_line(data = x$df_dose,
aes(x = MFx, y = q50), color = "orange")
if(!is.null(x$X_prop)){
MFx_plt <- MFx_plt +
geom_point(data = dplyr::filter(x$df_dose, id == "q50"),
aes(x = MFx, y = q50), color = "orange", shape = 4) +
geom_point(data = dplyr::filter(x$df_dose, id == "qinf95"),
aes(x = MFx, y = qinf95), color = "grey", shape = 4) +
geom_point(data = dplyr::filter(x$df_dose, id == "qsup95"),
aes(x = MFx, y = qsup95), color = "grey", shape = 4)
legend.point = data.frame(
x.pts = x$df_MFx$MFx,
y.pts = rep(x$X_prop, 3),
pts.leg = c(paste("median: ", round(x$df_MFx$MFx[1],digits = 2)),
paste("quantile 2.5%: ", round(x$df_MFx$MFx[2],digits = 2)),
paste("quantile 97.5%: ", round(x$df_MFx$MFx[3],digits = 2)))
)
MFx_plt <- MFx_plt +
geom_hline(yintercept = x$X_prop, col="grey70", linetype=2) +
geom_point(data = legend.point,
aes(x = x.pts, y = y.pts, color = pts.leg))
warning("This is not an error message:
Just take into account that MFx has been estimated with a binary
search using the 'accuracy' argument. Cross points indicate the
position of evaluated time series. To improve the shape of the curve, you
can use X = NULL, and compute time series around the median MFx, with the
vector `MFx_range`.")
}
if(log_scale == TRUE){
MFx_plt <- MFx_plt +
scale_x_log10()
}
}
if(x_variable == "Time"){
# Plot
if(is.null(main)) main <- paste("Survival over time. Multiplication Factor of", x$X_prop_provided*100, "percent")
if(is.null(x$X_prop)) main <- paste("Survival over time")
MFx = x$MFx_tested
k <- 1:length(x$ls_predict)
#
# Make a dataframe with quantile of all generated time series
#
ls_predict_quantile <- lapply(k, function(kit){
df_quantile <- x$ls_predict[[kit]]$df_quantile
df_quantile$MFx <- rep(MFx[[kit]], nrow(x$ls_predict[[kit]]$df_quantile))
return(df_quantile)
})
predict_MFx_quantile <- do.call("rbind", ls_predict_quantile)
if(!is.null(x$X_prop)){
initial_predict <- dplyr::filter(predict_MFx_quantile, MFx == 1)
final_predict <- dplyr::filter(predict_MFx_quantile, MFx == x$df_MFx$MFx[1])
y_arrow <- as.numeric(dplyr::filter(initial_predict, time == x$time_MFx)$q50)
yend_arrow <- as.numeric(dplyr::filter(final_predict, time == x$time_MFx)$q50)
MFx_plt <- MFx_plt +
labs(title = main,
x = xlab,
y = ylab) +
# final predict
geom_ribbon(data = final_predict,
aes(x = time, ymin = qinf95, ymax = qsup95),
fill = "grey30", alpha = 0.4) +
geom_line(data = final_predict,
aes(x = time, y = q50),
col="orange", size = 1) +
# initial predict
geom_ribbon(data = initial_predict,
aes(x = time, ymin = qinf95, ymax = qsup95),
fill = "grey60", alpha = 0.4) +
geom_line(data = initial_predict,
aes(x = time, y = q50),
col="orange", size = 0.3) +
# arrow
geom_segment(aes(x = x$time_MFx, y = y_arrow,
xend = x$time_MFx, yend = yend_arrow),
arrow = arrow(length = unit(0.2,"cm")))
}
if(is.null(x$X_prop)){
MFx_plt <- MFx_plt +
labs(title = main,
x = xlab,
y = ylab) +
geom_ribbon(data = predict_MFx_quantile,
aes(x = time, ymin = qinf95, ymax = qsup95),
fill = "grey30", alpha = 0.4) +
geom_line(data = predict_MFx_quantile,
aes(x = time, y = q50),
col="orange", size = 1) +
facet_wrap( ~ round(MFx, digits = 2), ncol = ncol)
}
}
return(MFx_plt)
}
R/plot.reproData.R
View file @ 8c5b023a
#' Plotting method for \code{reproData} objects
#'
#' This is the generic \code{plot} S3 method for the \code{reproData} class.
#' It plots the cumulated number of offspring as a function of time.
#'
#' @param x an object of class \code{reproData}
#' @param xlab label of the \eqn{X}-axis
#' @param ylab label of the \eqn{Y}-axis, by default \code{Cumulated Number of offspring}
#' @param main main title for the plot
#' @param concentration a numeric value corresponding to some concentration in
#' \code{data}. If \code{concentration = NULL}, draws a plot for each concentration
#' @param style graphical backend, can be \code{'ggplot'} or \code{'generic'}
#' @param pool.replicate if \code{TRUE}, the datapoints of each replicate are
#' summed for a same concentration
#' @param addlegend if \code{TRUE}, adds a default legend to the plot
#' @param remove.someLabels if \code{TRUE}, removes 3/4 of X-axis labels in
#' \code{'ggplot'} style to avoid the label overlap
#' @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{gg} and \code{ggplot}, see function \code{\link[ggplot2]{ggplot}}
#'
#' @keywords plot
#'
#' @examples
#' # (1) Load the data
#' data(cadmium1)
#'
#' # (2) Create an object of class 'reproData'
#' cadmium1 <- reproData(cadmium1)
#'
#' # (3) Plot the reproduction data
#' plot(cadmium1)
#'
#' # (4) Plot the reproduction data for a fixed concentration
#' plot(cadmium1, concentration = 4.36, style = "generic")
#'
#' @import ggplot2
#' @import grDevices
#' @importFrom methods is
#' @importFrom stats aggregate
#'
#' @export
plot.reproData <- function(x,
xlab,
ylab = "Cumulated Number of offspring",
main = NULL,
concentration = NULL,
style = "ggplot",
pool.replicate = FALSE,
addlegend = FALSE,
remove.someLabels = FALSE, ...) {
if (!is(x, "reproData"))
stop("plot.reproData: object of class reproData expected")
if (style == "generic" && remove.someLabels)
warning("'remove.someLabels' argument is valid only in 'ggplot' style.",
call. = FALSE)
if (is.null(concentration) && addlegend)
warning("'addlegend' argument is valid only when 'concentration' is not null.",
call. = FALSE)
if (pool.replicate) {
# agregate by sum of replicate
x <- cbind(aggregate(cbind(Nreprocumul, Nsurv, Ninit) ~ time + conc, x, sum),
replicate = 1)
}
if (is.null(concentration)) {
reproDataPlotFull(x, xlab, ylab, style, remove.someLabels)
}
else {
reproDataPlotFixedConc(x, xlab, ylab, main, concentration, style, addlegend,
remove.someLabels)
}
}
reproDataPlotFull <- function(data, xlab, ylab, style = "generic",
remove.someLabels) {
dataPlotFull(data, xlab, ylab, "Nreprocumul", style,
remove.someLabels)
}
reproDataPlotFixedConc <- function(x,
xlab,
ylab,
main,
concentration,
style = "generic",
addlegend = FALSE,
remove.someLabels = FALSE) {
dataPlotFixedConc(x, xlab, ylab, main, "Nreprocumul",
concentration, style, addlegend, remove.someLabels)
}
#' Plotting method for \code{reproData} objects
#'
#' This is the generic \code{plot} S3 method for the \code{reproData} class.
#' It plots the cumulated number of offspring as a function of time.
#'
#' @param x an object of class \code{reproData}
#' @param xlab label of the \eqn{X}-axis
#' @param ylab label of the \eqn{Y}-axis, by default \code{Cumulated Number of offspring}
#' @param main main title for the plot
#' @param concentration a numeric value corresponding to some concentration in
#' \code{data}. If \code{concentration = NULL}, draws a plot for each concentration
#' @param style graphical backend, can be \code{'ggplot'} or \code{'generic'}
#' @param pool.replicate if \code{TRUE}, the datapoints of each replicate are
#' summed for a same concentration
#' @param addlegend if \code{TRUE}, adds a default legend to the plot
#' @param remove.someLabels if \code{TRUE}, removes 3/4 of X-axis labels in
#' \code{'ggplot'} style to avoid the label overlap
#' @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{gg} and \code{ggplot}, see function \code{\link[ggplot2]{ggplot}}
#'
#' @keywords plot
#'
#' @return a plot of class \code{ggplot}
#'
#' @examples
#' # (1) Load the data
#' data(cadmium1)
#'
#' # (2) Create an object of class 'reproData'
#' cadmium1 <- reproData(cadmium1)
#'
#' # (3) Plot the reproduction data
#' plot(cadmium1)
#'
#' # (4) Plot the reproduction data for a fixed concentration
#' plot(cadmium1, concentration = 4.36, style = "generic")
#'
#' @import ggplot2
#' @import grDevices
#' @importFrom methods is
#' @importFrom stats aggregate
#'
#' @export
plot.reproData <- function(x,
xlab,
ylab = "Cumulated Number of offspring",
main = NULL,
concentration = NULL,
style = "ggplot",
pool.replicate = FALSE,
addlegend = FALSE,
remove.someLabels = FALSE, ...) {
if (!is(x, "reproData"))
stop("plot.reproData: object of class reproData expected")
if (style == "generic" && remove.someLabels)
warning("'remove.someLabels' argument is valid only in 'ggplot' style.",
call. = FALSE)
if (is.null(concentration) && addlegend)
warning("'addlegend' argument is valid only when 'concentration' is not null.",
call. = FALSE)
if (pool.replicate) {
# agregate by sum of replicate
x <- cbind(aggregate(cbind(Nreprocumul, Nsurv, Ninit) ~ time + conc, x, sum),
replicate = 1)
}
if (is.null(concentration)) {
reproDataPlotFull(x, xlab, ylab, style, remove.someLabels)
}
else {
reproDataPlotFixedConc(x, xlab, ylab, main, concentration, style, addlegend,
remove.someLabels)
}
}
reproDataPlotFull <- function(data, xlab, ylab, style = "generic",
remove.someLabels) {
dataPlotFull(data, xlab, ylab, "Nreprocumul", style,
remove.someLabels)
}
reproDataPlotFixedConc <- function(x,
xlab,
ylab,
main,
concentration,
style = "generic",
addlegend = FALSE,
remove.someLabels = FALSE) {
dataPlotFixedConc(x, xlab, ylab, main, "Nreprocumul",
concentration, style, addlegend, remove.someLabels)
}
R/plot.reproFitTT.R
View file @ 8c5b023a
This diff is collapsed.
R/plot.survDataCstExp.R
View file @ 8c5b023a
This diff is collapsed.
R/plot.survDataVarExp.R
View file @ 8c5b023a
#' Plotting method for \code{survDataVarExp} objects
#'
#' This is the generic \code{plot} S3 method for the \code{survDataVarC} class.
#' It plots the number of survivors as a function of time.
#'
#' @param x an object of class \code{survDataVarExp}
#' @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{Number of survivors}
#' @param main main title for the plot
#' @param one.plot if \code{TRUE}, draws all the points in one plot instead of one per \code{replicate}
#' @param facetting_level a vector of \code{characters} to rank \code{replicates} in the multi plot (i.e. \code{one.plot == FALSE})
#' @param \dots Further arguments to be passed to generic methods
#'
#' @return an object of class \code{ggplot}, see function \code{\link[ggplot2]{ggplot}}
#'
#' @keywords plot
#'
#' @import ggplot2
#' @importFrom graphics plot axis legend lines par points title
#' @importFrom methods is
#' @importFrom stats aggregate
#'
#' @export
plot.survDataVarExp <- function(x,
xlab = "Time",
ylab = "Number of survivors",
main = NULL,
one.plot = FALSE,
facetting_level = NULL,
...){
data_plt = filter( x , !is.na(Nsurv))
if(!is.null(facetting_level)){
data_plt$replicate = factor(data_plt$replicate, levels = facetting_level)
}
plt = ggplot(data_plt, aes(x = time, y = Nsurv, group = replicate, color = conc)) +
theme_minimal() +
theme(
## facet elements
strip.background = element_rect(fill="white", colour = "white"),
strip.text = element_text(colour = "black"),
## legend
legend.title = element_text(size = 9),
legend.text=element_text(size = 7),
legend.key.width = unit(0.3, "in"),
legend.key.height = unit(0.15, "in"),
legend.position = "top"
) +
labs(title = main,
x = xlab,
y = ylab,
colour = "Concentration" # legend title
) +
scale_colour_gradient(
name="Concentration",
low="grey20", high="orange"
) +
scale_fill_gradient(
name="Concentration",
low="grey20", high="orange"
) +
expand_limits(x = 0, y = 0) +
geom_point() +
geom_line()
if(one.plot == FALSE){
plt <- plt + facet_wrap(~ replicate, ncol = 2)
} else if(one.plot == TRUE){
plt
} else stop("'one.plot' is a logical 'TRUE' or 'FALSE'.")
return(plt)
}
#' Plotting method for \code{survDataVarExp} objects
#'
#' This is the generic \code{plot} S3 method for the \code{survDataVarC} class.
#' It plots the number of survivors as a function of time.
#'
#' @param x an object of class \code{survDataVarExp}
#' @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{Number of survivors}
#' @param main main title for the plot
#' @param one.plot if \code{TRUE}, draws all the points in one plot instead of one per \code{replicate}
#' @param facetting_level a vector of \code{characters} to rank \code{replicates} in the multi plot (i.e. \code{one.plot == FALSE})
#' @param \dots Further arguments to be passed to generic methods
#'
#' @return an object of class \code{ggplot}, see function \code{\link[ggplot2]{ggplot}}
#'
#' @keywords plot
#'
#' @import ggplot2
#' @importFrom graphics plot axis legend lines par points title
#' @importFrom methods is
#' @importFrom stats aggregate
#'
#' @export
plot.survDataVarExp <- function(x,
xlab = "Time",
ylab = "Number of survivors",
main = NULL,
one.plot = FALSE,
facetting_level = NULL,
...){
data_plt = filter( x , !is.na(Nsurv))
if(!is.null(facetting_level)){
data_plt$replicate = factor(data_plt$replicate, levels = facetting_level)
}
plt = ggplot(data_plt, aes(x = time, y = Nsurv, group = replicate, color = conc)) +
theme_minimal() +
theme(
## facet elements
strip.background = element_rect(fill="white", colour = "white"),
strip.text = element_text(colour = "black"),
## legend
legend.title = element_text(size = 9),
legend.text=element_text(size = 7),
legend.key.width = unit(0.3, "in"),
legend.key.height = unit(0.15, "in"),
legend.position = "top"
) +
labs(title = main,
x = xlab,
y = ylab,
colour = "Concentration" # legend title
) +
scale_colour_gradient(
name="Concentration",
low="grey20", high="orange"
) +
scale_fill_gradient(
name="Concentration",
low="grey20", high="orange"
) +
expand_limits(x = 0, y = 0) +
geom_point() +
geom_line()
if(one.plot == FALSE){
plt <- plt + facet_wrap(~ replicate, ncol = 2)
} else if(one.plot == TRUE){
plt
} else stop("'one.plot' is a logical 'TRUE' or 'FALSE'.")
return(plt)
}
R/plot.survFitCstExp.R
View file @ 8c5b023a
This diff is collapsed.
R/plot.survFitPredict.R
View file @ 8c5b023a
#' Plotting method for \code{survFitPredict} objects
#'
#' This is the generic \code{plot} S3 method for the
#' \code{survFitPredict}. It plots the predicted survival probability for each
#' concentration of the chemical compound in the provided dataset.
#'
#' The fitted curves represent the \strong{predicted survival probability} as a function
#' of time for each concentration.
#' The function plots both the 95\% credible band and the predicted survival
#' probability over time.
#' 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 (10\% of the MCMC chains are randomly
#' taken for this sample).
#'
#' @param x An object of class \code{survFitPredict}.
#' @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 probability}.
#' @param main A main title for the plot.
#' @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 mcmc_size A numerical value refering by default to the size of the mcmc in object \code{survFitPredict}.
#' This option is specific to \code{survFitPredict} objects for which computing time may be long.
#' \code{mcmc_size} can be used to reduce the number of mcmc samples in order to speed up
#' the computation.
#'
#' @param \dots Further arguments to be passed to generic methods.
#'
#' @keywords plot
#'
#' @examples
#'
#' # (1) Load the survival data
#' data("propiconazole_pulse_exposure")
#'
#' # (2) Create an object of class "survData"
#' dataset <- survData(propiconazole_pulse_exposure)
#'
#' \dontrun{
#' # (3) Run the survFit function
#' out <- survFit(dataset , model_type = "SD")
#'
#' # (4) Create a new data table for prediction
#' data_4prediction <- data.frame(time = 1:10, conc = c(0,5,5,5,0,0,5,5,5,5),
#' replicate= rep("predict", 10))
#'
#' # (5) Predict on a new dataset
#' predict_out <- predict(out, data_predict = data_4prediction, spaghetti = TRUE)
#'
#' # (6) Plot the predicted curve
#' plot(predict_out)
#' plot(predict_out, spaghetti = TRUE)
#' }
#'
#' @export
#'
#' @importFrom tidyr gather
#'
plot.survFitPredict <- function(x,
xlab = "Time",
ylab = "Survival probability",
main = NULL,
spaghetti = FALSE,
one.plot = FALSE,
mcmc_size = NULL,
...) {
df_prediction <- x$df_quantile
df_spaghetti <- x$df_spaghetti
# Plot
plt <- ggplot() +
theme_minimal() +
scale_x_continuous(name = xlab) +
scale_y_continuous(name = ylab,
limits = c(0,1)) +
theme(legend.position = "top")
# spaghetti
if(spaghetti == TRUE){
df_spaghetti_gather <- df_spaghetti %>%
tidyr::gather(survRate_key, survRate_value, -c(time,conc,replicate))
plt <- plt +
geom_line(data = df_spaghetti_gather,
aes(x = time, y = survRate_value, group = interaction(survRate_key, replicate)),
alpha = 0.2) +
geom_line(data = df_prediction,
aes(x = time, y= qinf95, group = replicate),
color = "orange", linetype = 2) +
geom_line(data = df_prediction,
aes(x = time, y = qsup95, group = replicate),
color = "orange", linetype = 2)
}
if(spaghetti != TRUE){
plt <- plt +
geom_ribbon(data = df_prediction,
aes(x = time, ymin = qinf95,ymax = qsup95, group = replicate),
fill = "grey70", alpha = 0.4)
}
# Prediction
plt <- plt +
geom_line(data = df_prediction,
aes(x = time, y = q50, group = replicate),
col="orange", size = 1)
# facetting
if(one.plot == FALSE){
plt <- plt + facet_wrap(~ replicate)
}
return(plt)
}
#' Plotting method for \code{survFitPredict} objects
#'
#' This is the generic \code{plot} S3 method for the
#' \code{survFitPredict}. It plots the predicted survival probability for each
#' concentration of the chemical compound in the provided dataset.
#'
#' The fitted curves represent the \strong{predicted survival probability} as a function
#' of time for each concentration.
#' The function plots both the 95\% credible band and the predicted survival
#' probability over time.
#' 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 (10\% of the MCMC chains are randomly
#' taken for this sample).
#'
#' @param x An object of class \code{survFitPredict}.
#' @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 probability}.
#' @param main A main title for the plot.
#' @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 mcmc_size A numerical value refering by default to the size of the mcmc in object \code{survFitPredict}.
#' This option is specific to \code{survFitPredict} objects for which computing time may be long.
#' \code{mcmc_size} can be used to reduce the number of mcmc samples in order to speed up
#' the computation.
#'
#' @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_pulse_exposure")
#'
#' # (2) Create an object of class "survData"
#' dataset <- survData(propiconazole_pulse_exposure)
#'
#' \donttest{
#' # (3) Run the survFit function
#' out <- survFit(dataset , model_type = "SD")
#'
#' # (4) Create a new data table for prediction
#' data_4prediction <- data.frame(time = 1:10, conc = c(0,5,5,5,0,0,5,5,5,5),
#' replicate= rep("predict", 10))
#'
#' # (5) Predict on a new dataset
#' predict_out <- predict(out, data_predict = data_4prediction, spaghetti = TRUE)
#'
#' # (6) Plot the predicted curve
#' plot(predict_out)
#' plot(predict_out, spaghetti = TRUE)
#' }
#'
#' @export
#'
#' @importFrom tidyr gather
#'
plot.survFitPredict <- function(x,
xlab = "Time",
ylab = "Survival probability",
main = NULL,
spaghetti = FALSE,
one.plot = FALSE,
mcmc_size = NULL,
...) {
df_prediction <- x$df_quantile
df_spaghetti <- x$df_spaghetti
# Plot
plt <- ggplot() +
theme_minimal() +
scale_x_continuous(name = xlab) +
scale_y_continuous(name = ylab,
limits = c(0,1)) +
theme(legend.position = "top")
# spaghetti
if(spaghetti == TRUE){
df_spaghetti_gather <- df_spaghetti %>%
tidyr::gather(survRate_key, survRate_value, -c(time,conc,replicate))
plt <- plt +
geom_line(data = df_spaghetti_gather,
aes(x = time, y = survRate_value, group = interaction(survRate_key, replicate)),
alpha = 0.2) +
geom_line(data = df_prediction,
aes(x = time, y= qinf95, group = replicate),
color = "orange", linetype = 2) +
geom_line(data = df_prediction,
aes(x = time, y = qsup95, group = replicate),
color = "orange", linetype = 2)
}
if(spaghetti != TRUE){
plt <- plt +
geom_ribbon(data = df_prediction,
aes(x = time, ymin = qinf95,ymax = qsup95, group = replicate),
fill = "grey70", alpha = 0.4)
}
# Prediction
plt <- plt +
geom_line(data = df_prediction,
aes(x = time, y = q50, group = replicate),
col="orange", size = 1)
# facetting
if(one.plot == FALSE){
plt <- plt + facet_wrap(~ replicate)
}
return(plt)
}