R/predict-SurvFitVarExp.R

In morseTKTD: Bayesian Inference of TKTD Models

#' @name PredictSurvFit
#' @export
predict_SurvFitVarExp <- function(fit,
                            display.exposure = NULL,
                            hb_value = NULL,
                            interpolate_length = NULL,
                            interpolate_method = "linear",
                            ...) {
    df <- display.exposure
    # PARAMETERS
    df.param <- extract_param(fit)
    # list over replicate
    ls_df <- lapply(unique(df$replicate), function(r){
      df[df$replicate == r, ]
    })
    # RUN PREDICTION
    if (fit$model_type == "SD") {
      df_mcmc <- predict_varSD(display.exposure = df,
                               display.parameters = df.param,
                               hb_value = hb_value,
                               interpolate_length = interpolate_length,
                               interpolate_method = interpolate_method)
    }
    if (fit$model_type == "IT") {
      df_mcmc <- predict_varIT(display.exposure = df,
                               display.parameters = df.param,
                               hb_value = hb_value,
                               interpolate_length = interpolate_length,
                               interpolate_method = interpolate_method)
    }

    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%"]
    }

    return_object <- list(df_quantile = df_quantile,
                          df_mcmc = df_mcmc)
    class(return_object) <- append("SurvPredict", class(return_object))
    return(return_object)
}

#' @name PredictSurvFit
#' @export
predict_varSD <- function(display.exposure = NULL,
                          display.parameters = NULL,
                          hb_value = NULL,
                          interpolate_length = NULL,
                          interpolate_method = NULL){
  df <- display.exposure
  kd <- display.parameters$kd
  hb <- display.parameters$hb
  if (!is.null(hb_value)) {
    hb <- rep(hb_value, length(display.parameters$hb))
  }
  z <- display.parameters$z
  kk <- display.parameters$kk
  ls_df <- lapply(unique(df$replicate), function(r){
    df[df$replicate == r, ]
  })
  ls = lapply(ls_df, function(d) {
    dout = SurvSD_var(Cw = d$conc,
                      time = d$time,
                      kd = kd,
                      hb = hb,
                      z = z,
                      kk = kk,
                      interpolate_length = interpolate_length,
                      interpolate_method = interpolate_method)
    dout$replicate = unique(d$replicate)
    return(dout)
  })
  df_mcmc <- do.call("rbind", c(ls, make.row.names = FALSE))
  return(df_mcmc)
}

#' @name PredictSurvFit
#' @export
predict_varIT <- function(display.exposure = NULL,
                          display.parameters = NULL,
                          hb_value = NULL,
                          interpolate_length = NULL,
                          interpolate_method = NULL){
  df <- display.exposure
  kd <- display.parameters$kd
  hb <- display.parameters$hb
  if (!is.null(hb_value)) {
    hb <- rep(hb_value, length(display.parameters$hb))
  }
  alpha <- display.parameters$alpha
  beta <- display.parameters$beta

  ls_df <- lapply(unique(df$replicate), function(r){
    df[df$replicate == r, ]
  })

  ls = lapply(ls_df, function(d) {
    dout = SurvIT_var(Cw = d$conc,
                      time = d$time,
                      kd = kd,
                      hb = hb,
                      alpha = alpha,
                      beta = beta,
                      interpolate_length = NULL,
                      interpolate_method = interpolate_method)
    dout$replicate = unique(d$replicate)
    return(dout)
  })
  df_mcmc <- do.call("rbind", c(ls, make.row.names = FALSE))
  return(df_mcmc)
}


#' @title Internal predict function
#'
#' @description
#' Survival function for "SD" model with external concentration changing
#'  with time
#'
#'
#' @param Cw A scalar of external concentration
#' @param time A vector of time
#' @param kd a vector of parameter
#' @param hb a vector of parameter
#' @param z a vector of parameter
#' @param kk a vector of parameter
#' @param interpolate_length if \code{display.time} is \code{NULL}, the argument
#' \code{interpolate_length} can be used to provide a sequence from 0 to maximum of
#' the time of exposure in original dataset (used for fitting).
#' @param interpolate_method The interpolation method for concentration.
#'  See package \code{deSolve} for details.
#'
#' @return A data.frame with exposure columns \code{time} and \code{conc} and
#' the resulting probabilisty of survival in \code{Psurv_XX} column where
#' \code{XX} refer to an MCMC iteration
#'
SurvSD_var <- function(Cw, time, kd, hb, z, kk,
                       interpolate_length = NULL,
                       interpolate_method = c("linear", "constant")) {

    interpolate_method <- match.arg(interpolate_method)

    ## external signal with several rectangle impulses
    signal <- data.frame(times = time, import = Cw)
    if (!is.null(interpolate_length)) {
      time_seq <- seq(min(time), max(time), length.out = interpolate_length)
      time <- sort(unique(c(time, time_seq)))
    }

    mcmc_size <- length(kd)
    xstart <- c(rep(c(D = 0), mcmc_size),
                rep(c(H = 0), mcmc_size))
    # ordering of parameters required by compiled function
    parms <- c(mcmc_size, kd, hb, z, kk)
    # solve model
    on.exit(.C("gutsredsd_free")) # clean up
    out <- deSolve::ode(y = xstart,
                 times = time,
                 parms = parms,
                 method = "lsoda",
                 dllname = "morseTKTD",
                 initfunc = "gutsredsd_init",
                 func = "gutsredsd_func",
                 initforc = "gutsredsd_forc",
                 forcings = signal,
                 fcontrol = list(method = interpolate_method,
                                 rule = 2, ties = "ordered"),
                 nout = 1,
                 outnames = "Exposure"
    )
    S <- exp(-out[,grep("H", colnames(out))])
    if (!is.matrix(S)) {
      S <- as.matrix(S)
    }
    df_Psurv <- as.data.frame(S)
    colnames(df_Psurv) = paste("Psurv", 1:ncol(df_Psurv), sep = "_")
    df_Psurv$conc = out[, "Exposure"]
    df_Psurv$time = out[, "time"]
    return(df_Psurv)
}

#' @title Internal predict function
#'
#' @description
#' Survival function for "IT" model with external concentration changing
#'  with time
#'
#' @param Cw A vector of external concentration
#' @param time A vector of time
#' @param kd a vector of parameter
#' @param hb a vector of parameter
#' @param alpha a vector of parameter
#' @param beta a vector of parameter
#' @param interpolate_length if \code{display.time} is \code{NULL}, the argument
#' \code{interpolate_length} can be used to provide a sequence from 0 to maximum of
#' the time of exposure in original dataset (used for fitting).
#' @param interpolate_method The interpolation method for concentration.
#'  See package \code{deSolve} for details.
#'
#' @return A data.frame with exposure columns \code{time} and \code{conc} and
#' the resulting probabilisty of survival in \code{Psurv_XX} column where
#' \code{XX} refer to an MCMC iteration
#'
SurvIT_var <- function(Cw, time, kd, hb, alpha, beta,
                       interpolate_length = NULL,
                       interpolate_method=c("linear","constant")){
    interpolate_method <- match.arg(interpolate_method)

    ## external signal with several rectangle impulses
    signal <- data.frame(times = time, import = Cw)
    if (!is.null(interpolate_length)) {
        time_seq <- seq(min(time), max(time), length.out = interpolate_length)
        time <- sort(unique(c(time, time_seq)))
    }

    ## The parameters
    mcmc_size <- length(kd)
    parms  <- c(mcmc_size, kd, hb)
    ## Start values for steady state
    xstart <- c(rep(c(D = 0), mcmc_size),
                rep(c(H = 0), mcmc_size))
    ## Solve model
    on.exit(.C("gutsredit_free")) # clean up
    out <- deSolve::ode(y = xstart,
                 times = time,
                 parms = parms,
                 method = "lsoda",
                 dllname = "morseTKTD",
                 initfunc = "gutsredit_init",
                 func = "gutsredit_func",
                 initforc = "gutsredit_forc",
                 forcings = signal,
                 fcontrol = list(method = interpolate_method,
                                 rule = 2, ties = "ordered"),
                 nout = 1,
                 outnames = "Exposure"
    )
    D <- out[,grep("D",colnames(out))]
    cumMax_D <- if (is.null(dim(D))) cummax(D) else apply(D, 2, cummax)
    ############### ATTENTION USE LOG-NORMAL RATHER THAN LOG-LOGISTIC !
    thresholdIT <- stats::plnorm(cumMax_D, meanlog = alpha, sdlog = beta)
    # thresholdIT <- t(1 / (1 + (t(cumMax_D) / alpha)^(-beta)))
    S <- (1 - thresholdIT) * exp(time %*% t(-hb))
    df_Psurv <- as.data.frame(S)
    colnames(df_Psurv) = paste("Psurv", 1:ncol(df_Psurv), sep = "_")
    df_Psurv$conc = out[, "Exposure"]
    df_Psurv$time = out[, "time"]
    return(df_Psurv)
}