# =============================================================================
# ANALYSE TRADEOFF : Impact de Jobs, Tracks, Site sur NbBranch et TimeRatio
# Approche en deux étapes :
#   Étape 1 : Modéliser les ÉCHECS (régression logistique)
#   Étape 2 : Sur les cas non-échoués, modéliser NbBranch et TimeRatio
# =============================================================================

# --- Packages -----------------------------------------------------------------
pkgs <- c("dplyr", "ggplot2", "tidyr", "lme4", "MASS", "nnet", "pROC",
          "ggeffects", "patchwork", "broom")
for (p in pkgs) if (!require(p, character.only = TRUE, quietly = TRUE))
  install.packages(p)

library(dplyr); library(ggplot2); library(tidyr); library(MASS)
library(nnet); library(pROC); library(ggeffects); library(patchwork); library(broom)

# --- 1. Lecture et préparation ------------------------------------------------
data_raw <- read.csv2("export_tradeoff.csv", header = TRUE,
                      stringsAsFactors = FALSE)

# Conversion numérique
data_raw <- data_raw %>%
  mutate(across(everything(), ~ suppressWarnings(as.numeric(.))))

# Variables factorielles (Site, Jobs, Tracks sont des niveaux discrets)
data_raw <- data_raw %>%
  mutate(
    Site_f    = factor(Site),
    Jobs_f    = factor(Jobs),
    Tracks_f  = factor(Tracks)
  )

cat("=== APERÇU DES DONNÉES ===\n")
cat("Dimensions :", nrow(data_raw), "lignes x", ncol(data_raw), "colonnes\n")
cat("\nValeurs de Jobs   :", sort(unique(data_raw$Jobs)), "\n")
cat("Valeurs de Tracks :", sort(unique(data_raw$Tracks)), "\n")
cat("Valeurs de Site   :", sort(unique(data_raw$Site)), "\n")
cat("Valeurs de Rames  :", sort(unique(data_raw$Rames)), "\n")

# =============================================================================
# ÉTAPE 0 : Vue d'ensemble des taux d'échec
# =============================================================================
data_raw <- data_raw %>%
  mutate(Echec = ifelse(NbBranch == -1, 1L, 0L))

cat("\n=== TAUX D'ÉCHEC GLOBAL ===\n")
cat(sprintf("Échecs : %d / %d (%.1f%%)\n",
            sum(data_raw$Echec), nrow(data_raw),
            100 * mean(data_raw$Echec)))

# Taux d'échec par facteur
for (var in c("Jobs_f", "Tracks_f", "Site_f")) {
  cat(sprintf("\n--- Taux d'échec par %s ---\n", var))
  tbl <- data_raw %>%
    group_by(.data[[var]]) %>%
    summarise(N = n(), Echecs = sum(Echec),
              Taux = round(100 * mean(Echec), 1), .groups = "drop")
  print(tbl)
}

# =============================================================================
# ÉTAPE 1 : MODÉLISATION DES ÉCHECS (Régression logistique)
# =============================================================================
cat("\n\n========== ÉTAPE 1 : MODÈLE LOGISTIQUE (Échec vs Succès) ==========\n")

# Modèle avec Jobs, Tracks, Site comme facteurs
mod_logit <- glm(Echec ~ Jobs_f + Tracks_f + Site_f,
                 data = data_raw, family = binomial(link = "logit"))

cat("\n--- Résumé du modèle logistique ---\n")
print(summary(mod_logit))

cat("\n--- Odds-Ratios (exp(coef)) avec IC 95% ---\n")
or_table <- exp(cbind(OR = coef(mod_logit), confint(mod_logit)))
print(round(or_table, 3))

# Test de déviance (ANOVA du modèle)
cat("\n--- Test de déviance (importance de chaque variable) ---\n")
print(anova(mod_logit, test = "Chisq"))

# AUC / ROC
roc_obj <- roc(data_raw$Echec, fitted(mod_logit), quiet = TRUE)
cat(sprintf("\nAUC = %.3f (pouvoir discriminant du modèle logistique)\n", auc(roc_obj)))

# Modèle avec effets d'interaction Jobs × Tracks
mod_logit_inter <- glm(Echec ~ Jobs_f * Tracks_f + Site_f,
                       data = data_raw, family = binomial(link = "logit"))
cat("\n--- Comparaison modèle sans vs avec interaction Jobs×Tracks ---\n")
print(anova(mod_logit, mod_logit_inter, test = "Chisq"))

# =============================================================================
# VISUALISATION ÉTAPE 1
# =============================================================================

# Graphe 1 : Taux d'échec moyen par Jobs et Tracks
df_heatmap <- data_raw %>%
  group_by(Jobs_f, Tracks_f) %>%
  summarise(TauxEchec = mean(Echec) * 100, N = n(), .groups = "drop")

p1 <- ggplot(df_heatmap, aes(x = Jobs_f, y = Tracks_f, fill = TauxEchec)) +
  geom_tile(color = "white", linewidth = 0.5) +
  geom_text(aes(label = sprintf("%.0f%%\n(n=%d)", TauxEchec, N)), size = 3) +
  scale_fill_gradient2(low = "steelblue", mid = "lightyellow", high = "firebrick",
                       midpoint = 50, name = "% Échec") +
  labs(title = "Taux d'échec (%) par Jobs et Tracks",
       x = "Jobs", y = "Tracks") +
  theme_minimal(base_size = 12)

# Graphe 2 : Taux d'échec par Site
df_site <- data_raw %>%
  group_by(Site_f) %>%
  summarise(TauxEchec = mean(Echec) * 100, N = n(),
            se = sqrt(mean(Echec) * (1 - mean(Echec)) / n()), .groups = "drop")

p2 <- ggplot(df_site, aes(x = Site_f, y = TauxEchec, fill = Site_f)) +
  geom_col(width = 0.6) +
  geom_errorbar(aes(ymin = TauxEchec - 1.96 * se * 100,
                    ymax = TauxEchec + 1.96 * se * 100), width = 0.2) +
  scale_fill_brewer(palette = "Set2", guide = "none") +
  labs(title = "Taux d'échec par Site",
       x = "Site", y = "% Échec") +
  theme_minimal(base_size = 12)

# Graphe 3 : Probabilités prédites par le modèle logistique
pred_data <- expand.grid(
  Jobs_f   = levels(data_raw$Jobs_f),
  Tracks_f = levels(data_raw$Tracks_f),
  Site_f   = levels(data_raw$Site_f)
) %>%
  mutate(prob_echec = predict(mod_logit, newdata = ., type = "response"))

p3 <- ggplot(pred_data, aes(x = Jobs_f, y = prob_echec, color = Tracks_f, group = Tracks_f)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  facet_wrap(~ Site_f, labeller = label_both) +
  scale_color_brewer(palette = "RdYlGn", name = "Tracks", direction = -1) +
  labs(title = "Probabilité d'échec prédite (modèle logistique)",
       x = "Jobs", y = "P(Échec)") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "right")

# =============================================================================
# ÉTAPE 2a : Sous-ensemble des non-échecs
# =============================================================================
data_ok <- data_raw %>%
  filter(NbBranch != -1, TimeRatio != -1) %>%
  mutate(
    Statut = case_when(
      TimeRatio <= 1 ~ "Succès",
      TimeRatio >  1 ~ "Semi-échec"
    ),
    Statut_f    = factor(Statut, levels = c("Succès", "Semi-échec")),
    logTimeRatio = log(TimeRatio)
  )

cat(sprintf("\n\nSous-ensemble non-échoués : %d lignes\n", nrow(data_ok)))
cat(sprintf("  Succès     (TimeRatio ≤ 1) : %d (%.1f%%)\n",
            sum(data_ok$Statut == "Succès"),
            100 * mean(data_ok$Statut == "Succès")))
cat(sprintf("  Semi-échec (TimeRatio > 1)  : %d (%.1f%%)\n",
            sum(data_ok$Statut == "Semi-échec"),
            100 * mean(data_ok$Statut == "Semi-échec")))

# =============================================================================
# ÉTAPE 2b : MODÈLE POUR NbBranch (régression binomiale négative)
# =============================================================================
cat("\n\n========== ÉTAPE 2b : NbBranch (régression binomiale négative) ==========\n")

mod_nb <- glm.nb(NbBranch ~ Jobs_f + Tracks_f + Site_f, data = data_ok)
print(summary(mod_nb))

cat("\n--- IRR (Incidence Rate Ratio = exp(coef)) avec IC 95% ---\n")
irr_table <- exp(cbind(IRR = coef(mod_nb), confint(mod_nb)))
print(round(irr_table, 3))

cat("\n--- Test de déviance (importance de chaque variable) ---\n")
print(anova(mod_nb))

# =============================================================================
# ÉTAPE 2c : MODÈLE POUR TimeRatio (régression log-linéaire)
# =============================================================================
cat("\n\n========== ÉTAPE 2c : TimeRatio (régression log-linéaire) ==========\n")

mod_tr <- lm(logTimeRatio ~ Jobs_f + Tracks_f + Site_f, data = data_ok)
print(summary(mod_tr))

cat("\n--- ANOVA (importance de chaque variable) ---\n")
print(anova(mod_tr))

cat("\n--- Effets sur l'échelle originale (exp(coef)) ---\n")
coef_table <- exp(cbind(Multiplicateur = coef(mod_tr),
                        confint(mod_tr)))
print(round(coef_table, 3))

# =============================================================================
# ÉTAPE 2d : DISTINCTION SUCCÈS vs SEMI-ÉCHEC (régression logistique)
# =============================================================================
cat("\n\n========== ÉTAPE 2d : Succès vs Semi-échec (logistique) ==========\n")

mod_logit2 <- glm(Statut_f ~ Jobs_f + Tracks_f + Site_f,
                  data = data_ok, family = binomial(link = "logit"))
print(summary(mod_logit2))

cat("\n--- Odds-Ratios ---\n")
or2 <- exp(cbind(OR = coef(mod_logit2), confint(mod_logit2)))
print(round(or2, 3))

# =============================================================================
# VISUALISATION ÉTAPE 2
# =============================================================================

# Graphe 4 : Distribution de TimeRatio par Jobs et Tracks (non-échecs)
p4 <- ggplot(data_ok, aes(x = Jobs_f, y = TimeRatio, fill = Jobs_f)) +
  geom_boxplot(outlier.size = 1, alpha = 0.7) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red", linewidth = 0.8) +
  facet_wrap(~ Tracks_f, labeller = label_both, nrow = 1) +
  scale_fill_brewer(palette = "Blues", guide = "none") +
  coord_cartesian(ylim = c(0, quantile(data_ok$TimeRatio, 0.95))) +
  labs(title = "TimeRatio par Jobs et Tracks (cas non-échoués)\n[ligne rouge = seuil 1 : succès / semi-échec]",
       x = "Jobs", y = "TimeRatio") +
  theme_minimal(base_size = 11)

# Graphe 5 : Distribution de NbBranch par Jobs et Tracks
p5 <- ggplot(data_ok, aes(x = Jobs_f, y = NbBranch, fill = Jobs_f)) +
  geom_boxplot(outlier.size = 1, alpha = 0.7) +
  facet_wrap(~ Tracks_f, labeller = label_both, nrow = 1) +
  scale_fill_brewer(palette = "Greens", guide = "none") +
  labs(title = "NbBranch par Jobs et Tracks (cas non-échoués)",
       x = "Jobs", y = "NbBranch") +
  theme_minimal(base_size = 11)

# Graphe 6 : % semi-échec par Jobs et Tracks
df_semi <- data_ok %>%
  group_by(Jobs_f, Tracks_f) %>%
  summarise(PctSemi = 100 * mean(Statut == "Semi-échec"), N = n(), .groups = "drop")

p6 <- ggplot(df_semi, aes(x = Jobs_f, y = Tracks_f, fill = PctSemi)) +
  geom_tile(color = "white") +
  geom_text(aes(label = sprintf("%.0f%%\n(n=%d)", PctSemi, N)), size = 3) +
  scale_fill_gradient2(low = "steelblue", mid = "lightyellow", high = "orange",
                       midpoint = 50, name = "% Semi-échec") +
  labs(title = "% Semi-échec (TimeRatio > 1) parmi les non-échoués",
       x = "Jobs", y = "Tracks") +
  theme_minimal(base_size = 12)

# =============================================================================
# AFFICHAGE DES GRAPHES
# =============================================================================
if (interactive()) {
  print(p1 + p2)
  print(p3)
  print(p4)
  print(p5)
  print(p6)
}

# Sauvegarde des graphes
ggsave("plot1_taux_echec_heatmap.png",  p1, width = 8, height = 6, dpi = 150)
ggsave("plot2_echec_par_site.png",      p2, width = 5, height = 4, dpi = 150)
ggsave("plot3_proba_echec_predit.png",  p3, width = 10, height = 5, dpi = 150)
ggsave("plot4_timeratio_boxplot.png",   p4, width = 12, height = 5, dpi = 150)
ggsave("plot5_nbbranch_boxplot.png",    p5, width = 12, height = 5, dpi = 150)
ggsave("plot6_semiechec_heatmap.png",   p6, width = 8,  height = 6, dpi = 150)

# =============================================================================
# SYNTHÈSE NUMÉRIQUE FINALE
# =============================================================================
cat("\n\n========== SYNTHÈSE FINALE ==========\n")

cat("\n[ÉTAPE 1 - ÉCHEC]\n")
cat("Variables significatives pour prédire l'échec (p < 0.05) :\n")
coef_logit <- tidy(mod_logit) %>% filter(p.value < 0.05, term != "(Intercept)")
print(coef_logit %>% select(term, estimate, p.value) %>%
        mutate(OR = round(exp(estimate), 2), p = round(p.value, 4)))

cat("\n[ÉTAPE 2b - NbBranch]\n")
cat("Variables significatives pour NbBranch (p < 0.05) :\n")
coef_nb <- tidy(mod_nb) %>% filter(p.value < 0.05, term != "(Intercept)")
print(coef_nb %>% select(term, estimate, p.value) %>%
        mutate(IRR = round(exp(estimate), 2), p = round(p.value, 4)))

cat("\n[ÉTAPE 2c - TimeRatio log]\n")
cat("Variables significatives pour log(TimeRatio) (p < 0.05) :\n")
coef_tr <- tidy(mod_tr) %>% filter(p.value < 0.05, term != "(Intercept)")
print(coef_tr %>% select(term, estimate, p.value) %>%
        mutate(Multiplicateur = round(exp(estimate), 2), p = round(p.value, 4)))

cat("\n[ÉTAPE 2d - Succès vs Semi-échec]\n")
cat("Variables significatives pour Statut (p < 0.05) :\n")
coef_l2 <- tidy(mod_logit2) %>% filter(p.value < 0.05, term != "(Intercept)")
print(coef_l2 %>% select(term, estimate, p.value) %>%
        mutate(OR = round(exp(estimate), 2), p = round(p.value, 4)))

cat("\nScript terminé. Graphes sauvegardés dans le répertoire courant.\n")