R/extract-SurvFit.R

In morseTKTD: Bayesian Inference of TKTD Models

#' @name Extract
#' 
#' @title Extraction methods to recover output of fit object.
#'
#' @description
#' - **extract_Nsurv_ppc**: extract the \code{Nsurv} generated with the sampler.
#' To be used for the Posterior Predictive Check (PPC).
#' 
#' - **extract_Nsurv_sim**: extract the \code{Nsurv} generated with the sampler.
#' To be used for the Simulation (sim).
#' 
#' - **extract_param**: extract parameters of \code{SD} or \code{IT} models.
#' 
#' - **priors_distribution**: Return a \code{data.frame} with prior density 
#' distributions of parameters used in the model.
#' 
#' @param fit An object of class \code{SurvFit}
#' @param size_sample Size of the random generation of the distribution.
#' @param \dots Further arguments to be passed to generic methods
#'
#' @return a \code{data.frame} with the extracted object from stanfit.
#'
#' @export
#' 
extract_Nsurv_ppc <- function(fit){
    ext = rstan::extract(fit$stanfit, pars = "Nsurv_ppc")
    quant = t(apply(ext$Nsurv_ppc, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE))
    dquant = data.frame(quant)
    colnames(dquant) = c("qinf95", "q50", "qsup95")
    dquant$time = fit$standata$time_N
    dquant$replicate = fit$standata$replicate_N
    dquant$Nsurv = fit$standata$Nsurv
    return(dquant)
}

#' @name Extract
#' @export
extract_Nsurv_sim <- function(fit){
    ext = rstan::extract(fit$stanfit, pars = "Nsurv_sim")
    quant = t(apply(ext$Nsurv_sim, 2, quantile, probs = c(0.025, 0.5, 0.975), na.rm = TRUE))
    dquant = data.frame(quant)
    colnames(dquant) = c("qinf95", "q50", "qsup95")
    dquant$time = fit$standata$time_N
    dquant$replicate = fit$standata$replicate_N
    dquant$Nsurv = fit$standata$Nsurv
    return(dquant)
}

#' @name Extract
#' @export
extract_param <- function(fit){
    if (fit$model_type == "SD") {
        ext = rstan::extract(fit$stanfit, pars = c("kd", "hb", "kk", "z"))
    }
    if (fit$model_type == "IT") {
        ext = rstan::extract(fit$stanfit, pars = c("kd", "hb", "alpha", "beta"))
    }
    df <- as.data.frame(ext)
    return(df)
}


#' @name Extract
#' @export
priors_distribution <- function(fit, ...){
    UseMethod("priors_distribution")
}

#' @name Extract
#' @export
priors_distribution.SurvFit <- function(fit,
                                        size_sample = 1e3,
                                        ...){
    
    param <- fit$standata
    
    kd_log10 <- rnorm(size_sample,
                      mean = param$kd_meanlog10,
                      sd = param$kd_sdlog10)
    hb_log10 <- rnorm(size_sample,
                      mean = param$hb_meanlog10,
                      sd = param$hb_sdlog10)
    
    df = data.frame(kd = 10^kd_log10,
                    hb = 10^hb_log10)
    
    if (fit$model_type == "SD") {
        z_log10 <- rnorm(size_sample,
                         mean = param$z_meanlog10,
                         sd = param$z_sdlog10)
        df$z <- 10^z_log10
        kk_log10 <- rnorm(size_sample,
                          mean = param$kk_meanlog10,
                          sd = param$kk_sdlog10)
        df$kk <- 10^kk_log10
    }
    if (fit$model_type == "IT") {
        alpha_log10 <- rnorm(size_sample,
                             mean = param$alpha_meanlog10,
                             sd = param$alpha_sdlog10)
        df$alpha <- 10^alpha_log10
        beta_log10 <- runif(size_sample,
                            min = param$beta_minlog10,
                            max = param$beta_maxlog10)
        df$beta <- 10^beta_log10
    }
    return(df)
}