R/lcx.R

In morseTKTD: Bayesian Inference of TKTD Models

#' @name LCxt
#' 
#' @title 
#' Predict Lethal Concentration at any specified time point.
#' 
#' @description
#' 
#' Predict the 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}.
#'  
#' @param fit An object used to select a method
#' @param x rate of individuals dying (e.g., \eqn{0.5} for
#'  \eqn{LC_{50}}, \eqn{0.1} for \eqn{LC_{10}}, ...).
#' @param t 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 exposure_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 interpolate_length of time point in the range of concentration between 0 
#' and the maximal concentration. 100 by default. description.
#' @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: 
#'  \itemize{
#' \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}.
#'  }
#' @export
#' 
lcxt <- function(fit, x, t, ...){
    UseMethod("lcxt")
}

#' @name LCxt
#' @export
lcxt.SurvFit <- function(fit, x = 0.5, t = NULL,
                         exposure_range = NULL,
                         interpolate_length = 50,
                        ...){
    
    if (is.null(t)) {
        t = max(fit$standata$time_N)
    }
    if (!is.null(exposure_range)) {
        exposure = seq(min(exposure_range), max(exposure_range),
                       length.out =  interpolate_length)
    } else{
        exposure = seq(0, max(fit$standata$conc), length.out =  interpolate_length)
    }
    df <- data.frame(
        time = t,
        conc = exposure,
        replicate = 1:interpolate_length
    )
    # EXTRACT FROM FIT
    df.param <- extract_param(fit)
    # RUN PREDICTION
    if (fit$model_type == "SD") {
        df_mcmc <- predict_cstSD(display.exposure = df,
                                 display.parameters = df.param,
                                 interpolate_length = NULL)
    }
    if (fit$model_type == "IT") {
        df_mcmc <- predict_cstIT(display.exposure = df,
                                 display.parameters = df.param,
                                 interpolate_length = NULL)
    }
    Psurv.col.position = grep("Psurv", colnames(df_mcmc))
    if (length(Psurv.col.position) == 1) {
        df_quantile = df_mcmc[, c("conc", "time", "replicate")]
        df_quantile$q50 = df_mcmc$Psurv_1
        df_quantile$qinf95 = df_quantile$q50 
        df_quantile$qsup95 = df_quantile$q50 
    } else{
        df_Psurv = df_mcmc[, Psurv.col.position]
        dfquant <- t(apply(df_Psurv, 1, quantile, probs = c(0.5,0.025,0.975), names = TRUE, na.rm = TRUE))
        
        df_quantile = df_mcmc[,grep("Psurv",colnames(df_mcmc), invert = TRUE)]
        df_quantile$qinf95 = dfquant[, "2.5%"]
        df_quantile$q50 = dfquant[, "50%"]
        df_quantile$qsup95 = dfquant[, "97.5%"]
    }
    
    dfQ <- df_quantile
    x_hb_rescale <- x * dfQ[dfQ$conc == 0, "q50"]
    
    LCX_q50 = linear_interpolate(x_hb_rescale, dfQ$q50, dfQ$conc)
    LCX_qinf95 = linear_interpolate(x_hb_rescale, dfQ$qinf95, dfQ$conc)
    LCX_qsup95 = linear_interpolate(x_hb_rescale, dfQ$qsup95, dfQ$conc)
    
    LCx <- list(
        q50 = LCX_q50,
        qinf95 = LCX_qinf95,
        qsup95 = LCX_qsup95)

    object_LCx <- list(x = x, t = t,
                       x_hb_rescale = x_hb_rescale,
                       df_quantile = df_quantile,
                       LCx = LCx)
    
    class(object_LCx) <- append("LCxt", class(object_LCx))
    
    return(object_LCx)
}

linear_interpolate <- function(target, x, conc){
    if (min(x) < target & target < max(x)) {
        interpol = approx( x, conc, xout = target, ties = mean)$y
    } else {
        interpol = NA
    }
    return(interpol)
}

#' @name PlotLCxt
#' 
#' @title Plot of the LCxt object
#' 
#' @description
#' Method for plotting output of lcxt function returning object
#'  of class \code{LCxt}.
#'  
#' @param x an object of class \code{LCxt}
#' @param xlab argument for the label of the x-axis
#' @param ylab argument for the label of the y-axis
#' @param main argument for the title of the graphic
#' @param \dots Further arguments to be passed to generic methods
#' 
#' @return an object of class \code{ggplot}, see function 
#' \code{\link[ggplot2]{ggplot}}
#' 
#' @export
plot.LCxt <- function(x,
                      xlab = "Concentration",
                      ylab = "Survival probability",
                      main = NULL, ...){
    
    plt <- ggplot(data = x$df_quantile) +
        theme_minimal() +
        labs(x = xlab, y = ylab, title = main) +
        geom_ribbon(aes(x = conc, ymin = qinf95, ymax = qsup95), fill = fillci) +
        geom_line(aes(x = conc, y = q50), color = colmed) +
        geom_hline(aes(yintercept = x$x_hb_rescale), linetype = 2)
    
        if (!is.na(x$LCx[["q50"]])) {
            plt <- plt + geom_point(color = cinter,
                aes(y = x$x_hb_rescale, x = x$LCx[["q50"]]))
        }
        if (!is.na(x$LCx[["qinf95"]])) {
            plt <- plt + geom_point(color = cinter,
                aes(y = x$x_hb_rescale, x = x$LCx[["qinf95"]]))
        }
        if (!is.na(x$LCx[["qsup95"]])) {
            plt <- plt + geom_point(color = cinter,
                aes(y = x$x_hb_rescale, x = x$LCx[["qsup95"]]))
        }
            
    return(plt)
}