evaluation.R

suppressPackageStartupMessages({
  library(broom)
  library(dplyr)
  library(purrr)
  library(tibble)
})

summarise_backtest_results <- function(predictions, group_cols = c("model")) {
  metric_input <- predictions %>%
    filter(is.finite(actual), is.finite(prediction))

  if ("horizon_day" %in% names(metric_input)) {
    metric_input <- metric_input %>%
      mutate(horizon_bucket = horizon_bucket(horizon_day))
  } else if (!"horizon_bucket" %in% names(metric_input)) {
    metric_input <- metric_input %>%
      mutate(horizon_bucket = "All Horizons")
  }

  metric_input %>%
    mutate(
      error = prediction - actual,
      abs_error = abs(error)
    ) %>%
    group_by(across(all_of(group_cols))) %>%
    summarise(
      mae = mean(abs_error, na.rm = TRUE),
      rmse = sqrt(mean(error ^ 2, na.rm = TRUE)),
      bias = mean(error, na.rm = TRUE),
      mape = yardstick::mape_vec(truth = pmax(actual, 1), estimate = pmax(prediction, 1)),
      smape = yardstick::smape_vec(truth = actual, estimate = prediction),
      rsq = yardstick::rsq_vec(truth = actual, estimate = prediction),
      n = n(),
      .groups = "drop"
    )
}

score_model_stability <- function(predictions) {
  predictions %>%
    group_by(model, origin_date) %>%
    summarise(fold_mae = mean(abs(prediction - actual), na.rm = TRUE), .groups = "drop") %>%
    group_by(model) %>%
    summarise(
      monthly_mae_sd = stats::sd(fold_mae, na.rm = TRUE),
      monthly_mae_range = max(fold_mae, na.rm = TRUE) - min(fold_mae, na.rm = TRUE),
      .groups = "drop"
    )
}

select_final_model <- function(backtest_predictions) {
  overall_summary <- summarise_backtest_results(backtest_predictions, group_cols = "model")
  stability_summary <- score_model_stability(backtest_predictions)

  scorecard <- overall_summary %>%
    left_join(stability_summary, by = "model") %>%
    mutate(
      model_class = case_when(
        model %in% c("annual_seasonal_naive", "lag_7") ~ "baseline",
        model == "linear" ~ "transparent",
        model == "glmnet" ~ "regularized",
        model == "xgboost" ~ "tree",
        TRUE ~ "other"
      )
    )

  baseline_thresholds <- scorecard %>%
    filter(model %in% c("annual_seasonal_naive", "lag_7")) %>%
    summarise(max_baseline_mae = max(mae, na.rm = TRUE))

  eligible_candidates <- scorecard %>%
    filter(
      !model %in% c("annual_seasonal_naive", "lag_7"),
      mae < baseline_thresholds$max_baseline_mae
    ) %>%
    arrange(mae, rmse, abs(bias), monthly_mae_sd)

  if (nrow(eligible_candidates) == 0) {
    selected_model <- scorecard %>%
      filter(model == "lag_7") %>%
      slice(1) %>%
      mutate(selection_reason = "No candidate beat both naive benchmarks; retained the 7-day lag baseline.")

    return(list(selected_model = selected_model, scorecard = scorecard))
  }

  selected_model <- eligible_candidates %>% slice(1)

  glmnet_row <- eligible_candidates %>% filter(model == "glmnet")
  xgboost_row <- eligible_candidates %>% filter(model == "xgboost")

  if (nrow(glmnet_row) == 1 && nrow(xgboost_row) == 1) {
    # Prefer the simpler GLMNET unless XGBoost wins clearly and stays stable.
    xgboost_improvement <- (glmnet_row$mae - xgboost_row$mae) / glmnet_row$mae
    xgboost_rmse_ok <- xgboost_row$rmse <= glmnet_row$rmse * 1.01
    xgboost_bias_ok <- abs(xgboost_row$bias) <= abs(glmnet_row$bias) * 1.05 || abs(xgboost_row$bias) <= abs(glmnet_row$bias)

    monthly_wins <- backtest_predictions %>%
      filter(model %in% c("glmnet", "xgboost")) %>%
      group_by(model, origin_date) %>%
      summarise(mae = mean(abs(prediction - actual), na.rm = TRUE), .groups = "drop") %>%
      group_by(origin_date) %>%
      slice_min(mae, n = 1, with_ties = FALSE) %>%
      ungroup() %>%
      count(model, name = "months_won")

    xgboost_months_won <- monthly_wins$months_won[monthly_wins$model == "xgboost"] %||% 0

    if (xgboost_improvement < 0.03 || !xgboost_rmse_ok || !xgboost_bias_ok || xgboost_months_won < 4) {
      selected_model <- glmnet_row
    } else {
      selected_model <- xgboost_row
    }
  }

  if (selected_model$model != "linear") {
    linear_row <- eligible_candidates %>% filter(model == "linear")

    if (nrow(linear_row) == 1) {
      mae_gap <- (linear_row$mae - selected_model$mae) / linear_row$mae

      if (mae_gap <= 0.01 && linear_row$rmse <= selected_model$rmse * 1.01 && abs(linear_row$bias) <= abs(selected_model$bias) * 1.05) {
        selected_model <- linear_row
      }
    }
  }

  selected_model <- selected_model %>%
    mutate(
      selection_reason = case_when(
        model == "linear" ~ "Selected for transparency because accuracy was effectively tied with more complex models.",
        model == "glmnet" ~ "Selected as the best balance of accuracy, stability, and interpretability.",
        model == "xgboost" ~ "Selected because it cleared the 3% MAE improvement threshold over GLMNET without creating a bias or stability penalty.",
        model == "lag_7" ~ "Retained as the operational baseline because no candidate model beat both naive benchmarks.",
        TRUE ~ "Selected by the scorecard decision rule."
      )
    )

  list(selected_model = selected_model, scorecard = scorecard)
}

extract_feature_importance <- function(fitted_workflow, model_name, top_n = 20) {
  if (model_name %in% c("linear", "glmnet")) {
    tidy_output <- broom::tidy(fitted_workflow) %>%
      filter(term != "(Intercept)") %>%
      mutate(importance = abs(estimate)) %>%
      arrange(desc(importance))

    return(tidy_output %>% slice_head(n = top_n))
  }

  if (model_name == "xgboost") {
    if (requireNamespace("vip", quietly = TRUE)) {
      return(vip::vi(fitted_workflow) %>% arrange(desc(Importance)) %>% slice_head(n = top_n))
    }

    xgb_fit <- workflows::extract_fit_engine(fitted_workflow)
    importance <- xgboost::xgb.importance(model = xgb_fit)
    return(as_tibble(importance) %>% slice_head(n = top_n))
  }

  tibble()
}

add_prediction_intervals <- function(predictions, reference_predictions) {
  interval_80 <- calc_empirical_intervals(reference_predictions, level = 0.8)
  interval_95 <- calc_empirical_intervals(reference_predictions, level = 0.95)

  predictions %>%
    mutate(
      lower_80 = pmax(prediction + interval_80[["lower"]], 0),
      upper_80 = pmax(prediction + interval_80[["upper"]], 0),
      lower_95 = pmax(prediction + interval_95[["lower"]], 0),
      upper_95 = pmax(prediction + interval_95[["upper"]], 0)
    )
}