R/AIC_functions.R

In bioinactivation: Mathematical Modelling of (Dynamic) Microbial Inactivation

#' Goodness of fit for Isothermal fits
#' 
#' @param object An object of class IsoFitInactivation.
#' 
#' @importFrom purrr %>%
#' @importFrom stats logLik
#' 
#' 
goodness_iso <- function(object) {
    
    residuals <- residuals(object$nls)
    
    my_dof <- summary(object$nls)$df
    
    n_pars <- my_dof[1]
    
    n_obs <- sum(my_dof)
    
    data.frame(ME = mean(residuals),
               RMSE = sqrt(mean(residuals^2)),
               loglik = logLik(object$nls)
    ) %>%
        mutate(AIC = 2*n_pars - 2*.data$loglik) %>%
        mutate(AICc = .data$AIC + 2*n_pars*(n_pars + 1)/(n_obs - n_pars - 1),
               BIC = log(n_obs)*n_pars - 2*.data$loglik) %>%
        mutate(Af = 10^.data$RMSE, Bf = 10^.data$ME)
    
}

#' Goodness of fit for Dynamic fits
#' 
#' @importFrom FME modCost
#' @importFrom dplyr select
#' 
#' @param object An instance of FitInactivation
#' 
goodness_dyna <- function(object) {
    
    my_simulation <- object$best_prediction$simulation %>%
        select("time", "logN")
    
    my_cost <- modCost(model = my_simulation, obs = object$data)
    
    residuals <- residuals(object$fit_results)
    
    n_pars <- length(object$fit_results$par)
    
    n_obs <- nrow(object$data)
    
    data.frame(ME = mean(residuals),
               RMSE = my_cost$model,
               loglik = my_cost$minlogp
    ) %>%
        mutate(AIC = 2*n_pars - 2*.data$loglik) %>%
        mutate(AICc = .data$AIC + 2*n_pars*(n_pars + 1)/(n_obs - n_pars - 1),
               BIC = log(n_obs)*n_pars - 2*.data$loglik) %>%
        mutate(Af = 10^.data$RMSE, Bf = 10^.data$ME)
    
}

#' Goodness of fit for MCMC fits
#' 
#' @importFrom FME modCost
#' @importFrom dplyr select
#' @importFrom rlang .data
#' 
#' @param object An instance of FitInactivationMCMC
#' 
goodness_MCMC <- function(object) {
    
    my_simulation <- object$best_prediction$simulation %>%
        select("time", "logN")
    
    my_cost <- modCost(model = my_simulation, obs = object$data)
    
    residuals <- my_cost$residuals$res
    
    n_pars <- length(object$modMCMC$bestpar)

    n_obs <- nrow(object$data)

    data.frame(ME = mean(residuals),
                      RMSE = my_cost$model,
                      loglik = my_cost$minlogp
    ) %>%
        mutate(AIC = 2*n_pars - 2*.data$loglik) %>%
        mutate(AICc = .data$AIC + 2*n_pars*(n_pars + 1)/(n_obs - n_pars - 1),
               BIC = log(n_obs)*n_pars - 2*.data$loglik) %>%
        mutate(Af = 10^.data$RMSE, Bf = 10^.data$ME)

    
}

#' Goodness of fit for microbial inactivation models
#' 
#' Generates a table with several statistical indexes describing the
#' quality of the model fit: 
#'         \itemize{
#'           \item ME: Mean Error.
#'           \item RMSE: Root Mean Squared Error.
#'           \item loglik: log-likelihood.
#'           \item AIC: Akaike Information Criterion.
#'           \item AICc: Akaike Information Criterion with correction for 
#'           finite sample size.
#'           \item BIC: Bayesian Infromation Criterion.
#'           \item Af: Accuracy factor.
#'           \item Bf: Bias factor.
#'           }
#' 
#' @param object A model fit generated by bioinactivation
#' 
#' @export
#'
goodness_of_fit <- function(object) {
    
    if (is.IsoFitInactivation(object)) {
        
        goodness_iso(object)
        
    } else if (is.FitInactivation(object)) {
        
        goodness_dyna(object)
        
    } else if (is.FitInactivationMCMC(object)) {
        
        goodness_MCMC(object)

        
    } else {
        
        stop("Not recognized object type")
        
    }
}