R/lpx.R

In morseTKTD: Bayesian Inference of TKTD Models

#' @name LPxt
#' 
#' @title Lethal Profile calculation
#' 
#' @description
#' Predict the Lethal Profile factor leading to $x$ % of reduction in survival
#' at a specific time $t$.
#' 
#' Generic method for \code{LPxt}, a function denoted \eqn{LP(x,t)} for 
#' \eqn{x}\% Multiplication Factor at time \eqn{t}.
#' 
#' @param fit An object used to select a method
#' @param object An object of class \code{LPx}
#' @param x rate of individuals dying (e.g., \eqn{0.5} for
#'  \eqn{LP_{50}}, \eqn{0.1} for \eqn{LP_{10}}, ...).
#' @param t A number giving the time at which  \eqn{LP_{x}} has to be estimated. 
#' If NULL, the latest time point of the experiment is used.
#' @param display.exposure A vector of the exposure porfile
#' @param interpolate_length of time point in the range of concentration between 0 
#' and the maximal concentration. 100 by default. description.
#' @param max.steps max steps to find the LPxt
#' @param accuracy accuracy of the LPxt algorithm (stop when reaching the accuracy).
#' @param \dots Further arguments to be passed to generic methods
#' 
#' @param \dots Further arguments to be passed to generic methods
#' 
#' @return returns an object of class \code{LPxt}
#' 
#' @export
#' 
lpxt <- function(fit, x, ...){
    UseMethod("lpxt")
}

#' @name LPxt
#' @export
lpxt.SurvFit <- function(fit, x = 0.5, t = NULL,
                         display.exposure = NULL,
                         interpolate_length = NULL,
                         max.steps = 100,
                         accuracy = 0.01,
                         ...){
    
    # EXTRACT FROM FIT
    df.param <- extract_param(fit)
    if (is.null(t)) {
        t = max(display.exposure[["time"]])
    } else{
        if (!(t %in% display.exposure[["time"]])) {
            stop("[t] must be in time of [display.exposure].")
        }
    }
    ##########
    #
    # DO THE BINARY SEARCH ON THE q50 ONLY TO SPEEDUP THE SEARCH!
    #
    df.param_q50 <- apply(df.param, 2, quantile, probs = 0.5)
    df.param_q50 <- as.data.frame(t(df.param_q50))
    if (fit$model_type == "SD") {
        df_q50 <- predict_varSD(display.exposure = display.exposure,
                                 display.parameters = df.param_q50,
                                 interpolate_length = NULL)
    }
    if (fit$model_type == "IT") {
        df_q50 <- predict_varIT(display.exposure = display.exposure,
                                 display.parameters = df.param_q50,
                                 interpolate_length = NULL)
    }
    x_get = df_q50[df_q50[["time"]] == t, ]$Psurv_1
    diff = x_get / x
    MFx_ = c(0, 1)
    x_get_ = c(NA, x_get)
    
    while (abs(diff - 1 ) > accuracy) {
        if (diff > 1) {
            if (max(MFx_) == MFx_[length(MFx_)]) {
                MFx = MFx_[length(MFx_)] * 10
            } else{
                MFx = MFx_[length(MFx_)] + 
                    abs(MFx_[length(MFx_)] - MFx_[length(MFx_) - 1]) / 2
            }
        } else{
            MFx = abs(MFx_[length(MFx_)] - MFx_[length(MFx_) - 1]) / 2
        }
        MFx_ = append(MFx_, MFx)
        
        df <- display.exposure
        df$conc <- df$conc * MFx
        if (fit$model_type == "SD") {
            df_q50 <- predict_varSD(display.exposure = df,
                                     display.parameters = df.param_q50,
                                     interpolate_length = NULL)
        }
        if (fit$model_type == "IT") {
            df_q50 <- predict_varIT(display.exposure = df,
                                     display.parameters = df.param_q50,
                                     interpolate_length = NULL)
        }
        x_get = df_q50[nrow(df_q50), ]$Psurv_1
        x_get_ = append(x_get_, x_get)
        diff = x_get / x
        if (length(MFx_) > max.steps) {
            MFx_test <- NULL
            warning(paste("The number of iterations reached the threshold number of iterations of", max.steps))
            break
        }
    }
    # OUT
    object_MFx <- list(
        fit, x = x,
        display.exposure = display.exposure,
        interpolate_length = interpolate_length,
        x_get_ = x_get_,
        MFx_ = MFx_,
        df_q50 = df_q50
    )
    
    class(object_MFx) <- append("LPxt", class(object_MFx))
    return(object_MFx)
}

#' @name LPxt
#' @export 
update.LPxt <- function(object, accuracy = 0.01, max.steps = 100, ...){
    # RECOVER ATTRIBUTE
    fit = object$fit
    x = object$x
    display.exposure = object$display.exposure
    interpolate_length = object$interpolate_length
    x_get_ = object$x_get_
    x_get = x_get_[length(x_get_)]
    MFx_ = object$MFx_
    MFx = MFx_[length(MFx_)]
    # RUN 
    diff = x_get / x
    while (abs(diff - 1 ) > accuracy) {
        cat("\n")
        if (diff > 1) {
            if (max(MFx_) == MFx_[length(MFx_)]) {
                MFx = MFx_[length(MFx_)] * 10
            } else{
                MFx = MFx_[length(MFx_)] + 
                    abs(MFx_[length(MFx_)] - MFx_[length(MFx_) - 1]) / 2
            }
        } else{
            MFx = abs(MFx_[length(MFx_)] - MFx_[length(MFx_) - 1]) / 2
        }
        MFx_ = append(MFx_, MFx)
        df <- display.exposure
        df$conc <- df$conc * MFx

        if (fit$model_type == "SD") {
            df_q50 <- predict_varSD(display.exposure = df,
                                     display.parameters = df.param_q50,
                                     interpolate_length = NULL)
        }
        if (fit$model_type == "IT") {
            df_q50 <- predict_varIT(display.exposure = df,
                                     display.parameters = df.param_q50,
                                     interpolate_length = NULL)
        }
        x_get = df_q50[nrow(df_q50), ]$Psurv_1
        x_get_ = append(x_get_, x_get)
        diff = x_get / x
        if (length(MFx_) > max.steps) {
            MFx_test <- NULL
            warning(paste("The number of iterations reached the threshold number of iterations of", max.steps))
            break
        }
    }
    
    # OUT
    object_MFx <- list(
        fit, x = 0.5,
        display.exposure = display.exposure,
        interpolate_length = interpolate_length,
        x_get_ = x_get_,
        MFx_ = MFx_,
        df_q50 = df_q50
    )
    
    class(object_MFx) <- append("LPxt", class(object_MFx))
    return(object_MFx)
}

#' @name PlotLPxt
#' 
#' @title Plot of the LPxt object
#' 
#' @description
#' Method for plotting output of loxt function returning object
#'  of class \code{LPxt}.
#' 
#' @param x an object of class \code{LPxt}
#' @param plot style of the plot (default is curve)
#' @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.LPxt <- function(x,
                      plot = "curve",
                      xlab = "Time",
                      ylab = "Survival probability",
                      main = NULL, ...){
    
    plt_curve <- ggplot(data = x$df_q50) +
        theme_minimal() +
        labs(x = xlab, y = ylab, title = main) +
        lims(y = c(0,1)) +
        # geom_ribbon(aes(x = conc, ymin = qinf95, ymax = qsup95), fill = fillci) +
        geom_line(aes(x = time, y = Psurv_1), color = colmed) +
        geom_hline(aes(yintercept = x$x), linetype = 2)
    
    plt_MFx <- ggplot() + 
        theme_minimal() +
        labs(y = "x% of survival", x = "Lethal Profile for x%") +
        geom_point(aes(x = x$MFx_, y = x$x_get_), color = colmed) + 
        geom_line(aes(x = x$MFx_, y = x$x_get_), color = colmed, linetype = 2) + 
        geom_hline(aes(yintercept = x$x), linetype = 2) + 
        geom_vline(aes(xintercept = x$MFx_[length(x$MFx_)]), linetype = 2)
    
    if (plot == "curve") {
        plt = plt_curve
    } else{
        plt = plt_MFx
    }
    return(plt)
}