R/compare_models_by_treatment.R

In rumenGP: Rumen Gas Production Modeling, Comparison, and Visualization

#' Compare Models by Treatment
#'
#' Calculates model performance separately
#' for each treatment.
#'
#' Performance metrics are computed using
#' treatment-level predictions and observations,
#' allowing direct comparison of competing models
#' within each treatment.
#'
#' Typical metrics include:
#'
#' \itemize{
#'   \item R-squared (R²)
#'   \item Root Mean Squared Error (RMSE)
#'   \item Residual Sum of Squares (RSS)
#'   \item Akaike Information Criterion (AIC)
#'   \item Bayesian Information Criterion (BIC)
#' }
#'
#' This function is useful for determining whether
#' different treatments are best described by
#' different kinetic models.
#'
#' @param ... Fitted model objects.
#'
#' @examples
#'
#' files <- example_data()
#'
#' raw_data <- read_ankom(
#'   files$ankom
#' )
#'
#' metadata <- read_metadata(
#'   files$metadata
#' )
#'
#' gp <- process_ankom(
#'   raw_data,
#'   metadata,
#'   headspace_ml = 210,
#'   temperature_c = 39
#' )
#'
#' groot_fit <- fit_groot(
#'   gp
#' )
#'
#' gompertz_fit <- fit_gompertz(
#'   gp
#' )
#'
#' compare_models_by_treatment(
#'   Groot = groot_fit,
#'   Gompertz = gompertz_fit
#' )
#'
#' @return A data frame containing treatment-level
#' performance metrics for each fitted model.
#'
#' @seealso
#' \code{\link{compare_models}},
#' \code{\link{rank_models_by_treatment}},
#' \code{\link{best_model_by_treatment}},
#' \code{\link{model_win_frequency}}
#'
#' @export
compare_models_by_treatment <- function(...) {

  fits <- list(...)

  purrr::imap_dfr(
    fits,
    function(fit, model_name) {

      fit$diagnostics |>
        dplyr::group_by(
          Treatment
        ) |>
        dplyr::summarise(

          Model = model_name,

          Mean_R2 =
            mean(
              R2,
              na.rm = TRUE
            ),

          Mean_RMSE =
            mean(
              RMSE,
              na.rm = TRUE
            ),

          Mean_AIC =
            mean(
              AIC,
              na.rm = TRUE
            ),

          Mean_BIC =
            mean(
              BIC,
              na.rm = TRUE
            ),

          .groups = "drop"

        )

    }

  )

}