Nothing
#' Plotting method for \code{survFit} objects
#'
#' This is the generic \code{plot} S3 method for the
#' \code{survFit}. It plots the fit obtained for each
#' concentration of chemical compound in the original dataset.
#'
#' The fitted curves represent the \strong{estimated survival probability} as a function
#' of time for each concentration.
#' The black dots depict the \strong{observed survival
#' probability} at each time point. Note that since our model does not take
#' inter-replicate variability into consideration, replicates are systematically
#' pooled in this plot.
#' The function plots both 95\% credible intervals for the estimated survival
#' probability (by default the grey area around the fitted curve) and 95\% binomial confidence
#' intervals for the observed survival probability (as black error bars if
#' \code{adddata = TRUE}).
#' Both types of intervals are taken at the same level. Typically
#' a good fit is expected to display a large overlap between the two types of intervals.
#' If \code{spaghetti = TRUE}, the credible intervals are represented by two
#' dotted lines limiting the credible band, and a spaghetti plot is added to this band.
#' This spaghetti plot consists of the representation of simulated curves using parameter values
#' sampled in the posterior distribution (2\% of the MCMC chains are randomly
#' taken for this sample).
#'
#' @param x An object of class \code{survFit}.
#' @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 concentration A numeric value corresponding to some specific concentrations in
#' \code{data}. If \code{concentration = NULL}, draws a plot for each concentration.
#' @param spaghetti if \code{TRUE}, draws a set of survival curves using
#' parameters drawn from the posterior distribution
#' @param one.plot if \code{TRUE}, draws all the estimated curves in
#' one plot instead of one plot per concentration.
#' @param adddata if \code{TRUE}, adds the observed data to the plot
#' with (frequentist binomial) confidence intervals
#' @param addlegend if \code{TRUE}, adds a default legend to the plot.
#' @param style graphical backend, can be \code{'generic'} or \code{'ggplot'}
#' @param \dots Further arguments to be passed to generic methods.
#'
#' @keywords plot
#'
#' @return a plot of class \code{ggplot}
#'
#' @export
#'
#' @import ggplot2
#' @import grDevices
#' @importFrom dplyr filter
#' @importFrom dplyr select
#' @importFrom dplyr contains
#' @importFrom gridExtra grid.arrange arrangeGrob
#' @importFrom grid grid.rect gpar
#' @importFrom graphics plot
#' @importFrom tidyr tibble
#' @importFrom tidyr as_tibble
#'
plot.survFitCstExp <- function(x,
xlab = "Time",
ylab = "Survival probability",
main = NULL,
concentration = NULL,
spaghetti = FALSE,
one.plot = FALSE,
adddata = TRUE,
addlegend = FALSE,
style = "ggplot", ...) {
if (one.plot && !is.null(concentration)) one.plot <- FALSE
if ((addlegend && is.null(concentration)) ||
(addlegend && !one.plot))
warning("The legend is available only if [one.plot] is TRUE or if [concentration] is not NULL !", call. = FALSE)
if (!is.null(concentration) && !any(x$transformed.data$conc == concentration))
stop("The [concentration] argument is not one of the possible concentration !")
if (one.plot)
warning("The credible limits and confidence intervals are not drawn when 'one.plot' = TRUE.", call. = FALSE)
conf.int <- survTKTDConfInt_CstExp(x)
dobs <- data.frame(conc = x$transformed.data$conc,
time = x$transformed.data$time,
psurv = x$transformed.data$Nsurv / x$transformed.data$Ninit,
Points = "Observed values",
#color = as.numeric(as.factor(x$transformed.data$conc)),
color = x$transformed.data$replicate,
conf.int)
# remove time 0 in dobs
dobs <- dplyr::filter(dobs, time != 0)
data.credInt <- survFitPlotCITKTD_CstExp(x)
if (style == "generic") {
survFitPlotTKTDGeneric(data.credInt, dobs, xlab, ylab, main, concentration,
one.plot, spaghetti,
adddata, addlegend)
} else if (style == "ggplot") {
survFitPlotTKTDGG(data.credInt, dobs, xlab, ylab, main, concentration,
one.plot, spaghetti, adddata, addlegend)
} else stop("Unknown style")
}
#' @importFrom stats aggregate binom.test
survTKTDConfInt_CstExp <- function(x) {
# create confidente interval on observed data
# binomial model by a binomial test
# INPUT:
# - x : object of class survFitTT
# OUTPUT:
# - ci : confidente interval
ci <- apply(x$transformed.data, 1, function(x) {
binom.test(as.numeric(x["Nsurv"]), as.numeric(x["Ninit"]))$conf.int
})
ci <- as.data.frame(t(ci))
colnames(ci) <- c("qinf95", "qsup95")
ci$Conf.Int <- "Confidence interval"
return(ci)
}
Surv_SD <- function(Cw, time, kk, kd, z, hb)
{
S <- exp(-hb*time)
x <- ifelse(Cw > z, 1 - z/Cw, NA)
tz <- -(1/kd)*log(x)
y <- ifelse(time > tz,
exp( kk/kd*Cw*(exp(-kd*tz) -exp(-kd*time))
- kk*(Cw-z)*(time - tz)),
NA)
return(ifelse(!is.na(x) & !is.na(y), S * y, S))
}
Surv_IT <- function(Cw, time, kd, hb, alpha, beta)
{
D <- Cw*(1-exp(-kd %*% t(time)))
D.max <- t(apply(D, 1, cummax))
S <- exp(-hb %*% t(time)) * (1 - plogis(log(D.max), location = log(alpha), scale = 1/beta))
return(S)
}
survFitPlotCITKTD_CstExp <- function(x) {
# INPUT
# x : An object of class survFitTKTD
# OUTPUT
# A list of - dobs : observed values
# - dtheo : estimated values
npoints <- 100
concobs <- unique(x$transformed.data$conc)
## tfin <- seq(0, max(x$jags.data$t), length.out = npoints)
tfin <- seq(0, max(x$jags.data$time), length.out = npoints)
# 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(concobs)
j <- 1:npoints
model_type = x$model_type
if(model_type == "SD"){
z <- 10^mctot[, "z_log10"]
kk <- 10^mctot[, "kk_log10"]
dtheo <- lapply(k, function(kit) { # conc
sapply(j, function(jit) { # time
Surv_SD(Cw = concobs[kit],
time = tfin[jit],
kk = kk,
kd = kd,
z = z,
hb = hb)
})
})
}
if(model_type == "IT"){
alpha <- 10^mctot[, "alpha_log10"]
beta <- 10^mctot[, "beta_log10"]
dtheo <- lapply(k, function(kit) { # concentration pour chaque concentration
Surv_IT(Cw = concobs[kit],
time = tfin,
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)
dtheo <- as.data.frame(dtheo)
colnames(dtheo) <- paste0("X", 1:length(kd))
# divide number of mcmc by 50
sel <- sample(ncol(dtheo))[1:ceiling(ncol(dtheo) / 50)]
dtheo <- dtheo[, sel]
dtheo$conc <- rep(concobs, rep(npoints, length(concobs)))
dtheo$time <- rep(tfin, length(concobs))
# add credible limits
dtheo$qinf95 <- qinf95
dtheo$qsup95 <- qsup95
dtheo$q50 <- q50
# names for legend plots
dtheo$Cred.Lim <- "Credible limits"
dtheo$Mean.C <- "Mean curve"
# vector color
dtheo$color <- as.numeric(as.factor(dtheo$conc))
return(dtheo)
}
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.