PRÉPARATION DES DONNÉES

iMPORTATION DES DONNÉES

#Set working directory or use Session menu
setwd("C:/Users/theow/OneDrive - Universite de Montreal/Documents/PhD/Data/Data_theo/Dendro_Theo/Analyses_preliminaires") # VERRIFIER
           #
load("Etabl_corrige_py.RData")   

# 🔶 Événements (orange)
feu <- data.frame(
                   Site = c("T87","T52","T94"),
                   Start = c(1900, 1920, 1944),
                   End   = c(1900, 1920, 1944))

load("tbe.RData")  # (BD tbe
load("coupes.RData")  # (BD coupes


# Téléchargez le code glmm_funs.R pour vérifier la surdispersion
source(file = "glmm_funs.R")

# On charge le package coefplot2 qui n'est pas dans CRAN

if (!require("coefplot2"))
  remotes::install_github("palday/coefplot2", 
                          subdir = "pkg", 
                          upgrade = "always", 
                          quiet = TRUE)
library(coefplot2)

FONCTION QUI TRANSFORME LES BD DES PERTURBATIONS EN DONNÉES BINAIRES

TRANSFORMATION DE LA BD POUR VOIR LES ABONDANCES SUR TOUTE LA PLAGE D’INTERVAL DE VIE DES SITES

Etabl2 <- Etabl_corrige_py %>%
              # Filtrer les enregistrements où l'année d'établissement est valide (non nulle)
              filter(Year_etab_corr != 0) %>%
            
              # Regrouper les données par ces 4 variables pour étendre chaque combinaison à toutes les années
              group_by(Site, Placette, Esp, Id) %>%
            
              # Créer un nouveau dataframe avec toutes les années de 1770 à 2010 pour chaque groupe
              reframe(
                Year = 1770:2010,  # Génère toutes les années de 1770 à 2010
                Year_etab_corr = first(Year_etab_corr),  # Conserve la valeur originale de l'année d'établissement pour chaque groupe
                Abond = as.integer(1770:2010 == first(Year_etab_corr))  # Initialise Abond à 1 pour l'année d'établissement, 0 sinon
              ) %>% 
             # Corriger la colonne Abond pour s'assurer qu'elle vaut 1 uniquement pour les années valides
              mutate(Abond = ifelse(Year_etab_corr != 0 & Year == Year_etab_corr, 1, 0)) %>%
              # Supprimer les anciennes colonnes binaires
                dplyr::select(-Year_etab_corr)

RAJOUTER LES PERTURBATIONS BINAIRES AUX BASES DE DONNÉES

Etabl3 <- Etabl2 %>%
              
              # Ajouter des colonnes pour différents types de perturbations environnementales
              # Chaque fonction add_perturbation_flag() ajoute une colonne binaire indiquant la présence de la perturbation
              add_perturbation_flag(feu,    "feu") %>%    # Ajoute une colonne 'feu' (1 si perturbation par feu, 0 sinon)
              add_perturbation_flag(tbe,    "tbe") %>%    # Ajoute une colonne 'tbe' (1 si perturbation par tbe, 0 sinon)
              add_perturbation_flag(coupes, "coupes")     # Ajoute une colonne 'coupes' (1 si perturbation par coupes, 0 sinon)

GROUPER LES DONNÉES PAR ANNÉE ET

Etabl4 <- Etabl3 %>%

              group_by(Site,Esp,Year) %>% 
              summarise(
                        Abond_sum = sum(Abond, na.rm = TRUE),         # Ajoute la somme d'Abond
                        feu=first(feu),
                        coupes=first(coupes),
                        tbe=first(tbe),
                        .groups = 'drop'
                      ) %>%
                      ungroup() 

On tranforme les perturbation de binaire en variable catégorielle.

Etabl5 <- Etabl4 %>% 
                mutate(
                        Pert_type = case_when(
                          # Aucune perturbation (RAS)
                          feu == 0 & coupes == 0 & tbe == 0 ~ "RAS",
                    
                          # Perturbations simples
                          feu == 1 & coupes == 0 & tbe == 0 ~ "feu",
                          feu == 0 & coupes == 1 & tbe == 0 ~ "coupes",
                          feu == 0 & coupes == 0 & tbe == 1 ~ "tbe",
                    
                          # Chevauchement spécifique mentionné (tbe_coupe)
                          feu == 0 & coupes == 1 & tbe == 1 ~ "tbe_coupe",
                    
                          # Autres combinaisons de deux perturbations
                          feu == 1 & coupes == 1 & tbe == 0 ~ "feu_coupes",
                          feu == 1 & coupes == 0 & tbe == 1 ~ "feu_tbe",
                    
                          # Toutes les perturbations présentes
                          feu == 1 & coupes == 1 & tbe == 1 ~ "feu_coupes_tbe",
                    
                          # Par défaut (au cas où)
                          TRUE ~ "inconnu"
                        )
                ) #%>%
                # Supprimer les anciennes colonnes binaires
#                dplyr::select(-feu, -coupes, -tbe)

TRansformation de la BD

# Assurez-vous que votre base de données s'appelle bien Etabl5
# Etabl5 <- ...

# 1. Identifier les années de perturbations par site
evenements <- Etabl5 %>%
  group_by(Site) %>%
  summarise(
    feu_y = min(Year[feu == 1], na.rm = TRUE),
    coupes_y = min(Year[coupes == 1], na.rm = TRUE),
    tbe_start = min(Year[tbe == 1], na.rm = TRUE),
    tbe_end = max(Year[tbe == 1], na.rm = TRUE)
  ) %>%
  # Remplacer les Inf (qui apparaissent si un site n'a pas eu la perturbation) par NA
  mutate(across(everything(), ~ifelse(is.infinite(.), NA, .)))

# 2. Générer les fenêtres temporelles de 10 ans avant/après par perturbation
fenetres <- evenements %>%
  pmap_dfr(function(Site, feu_y, coupes_y, tbe_start, tbe_end) {
    
    liste <- list()
    
    # Pour le feu (1 an) : 10 ans avant l'année, 10 ans après l'année
    if(!is.na(feu_y)) {
      liste$feu <- data.frame(
        Site = Site,
        Year = c((feu_y-10):(feu_y-1), (feu_y+1):(feu_y+10)),
        Perturbation = "feu",
        Position = c(rep("Before", 10), rep("After", 10))
      )
    }
    
    # Pour les coupes (1 an) : 10 ans avant, 10 ans après
    if(!is.na(coupes_y)) {
      liste$coupes <- data.frame(
        Site = Site,
        Year = c((coupes_y-10):(coupes_y-1), (coupes_y+1):(coupes_y+10)),
        Perturbation = "coupes",
        Position = c(rep("Before", 10), rep("After", 10))
      )
    }
    
    # Pour le tbe (intervalle d'années) : 10 ans avant le début, 10 ans après la fin
    if(!is.na(tbe_start)) {
      liste$tbe <- data.frame(
        Site = Site,
        Year = c((tbe_start-10):(tbe_start-1), (tbe_end+1):(tbe_end+10)),
        Perturbation = "tbe",
        Position = c(rep("Before", 10), rep("After", 10))
      )
    }
    
    # Combiner les perturbations de ce site
    bind_rows(liste)
  })

# 3. Joindre les fenêtres à la base d'abondance et calculer les sommes
BD_finale_10ans <- Etabl5 %>%
  dplyr::select(Site, Esp, Year, Abond_sum) %>%
  # On ne garde que les années qui tombent dans une fenêtre de 10 ans
  inner_join(fenetres, by = c("Site", "Year")) %>%
  # Groupage final : Site, Espèce, Perturbation, Position (Before/After)
  group_by(Site, Esp, Perturbation, Position) %>%
  summarise(
    Abond_totale = sum(Abond_sum, na.rm = TRUE),
    Nb_annees_dispo = n(), # Colonnes pour vérifier si vous avez bien 10 ans de données
    .groups = "drop"
  )

# Visualiser le résultat
head(BD_finale_10ans)
# A tibble: 6 × 6
  Site  Esp   Perturbation Position Abond_totale
  <chr> <fct> <chr>        <chr>           <dbl>
1 T52   BOP   coupes       After               8
2 T52   BOP   coupes       Before              7
3 T52   BOP   feu          After               1
4 T52   BOP   feu          Before              1
5 T52   BOP   tbe          After               0
6 T52   BOP   tbe          Before              6
# ℹ 1 more variable: Nb_annees_dispo <int>

ERR

PRÉPARATION DES DONNÉES

BD_finale_10ans_ERR <- BD_finale_10ans %>% 
                       filter(Esp=="ERR")


# Imaginons que vous vouliez forcer l'ordre "before" puis "after" pour la variable Position
BD_finale_10ans_ERR$Site <- factor(BD_finale_10ans_ERR$Site, 
                                        levels =  c("T87","T77","T83", "T70", "T58","T52",  "T93",  "T94"))

BD_finale_10ans_ERR$Position <- factor(BD_finale_10ans_ERR$Position, 
                                        levels =  c("Before","After"))

ggplot(BD_finale_10ans_ERR, aes(x = Position, y = Abond_totale)) +
            geom_boxplot(aes(color = Position)) +
            theme(legend.position = "top")

ggplot(BD_finale_10ans_ERR, aes(x = Site, y = Abond_totale)) +
            geom_boxplot(aes(color = Position)) +
            theme(legend.position = "top")

TEST DE CHI2

CALCUL DES FRÉQUENCES PONDERÉE

# 1. Calculer la somme des abondances par Perturbation et Position
freq_table <- BD_finale_10ans_ERR %>%
  group_by(Perturbation, Position) %>%
  summarise(Somme_Abond = sum(Abond_totale, na.rm = TRUE)) %>%
  pivot_wider(
    names_from = Position,
    values_from = Somme_Abond,
    values_fill = 0  # Remplit les valeurs manquantes par 0
  ) %>% 
   column_to_rownames("Perturbation")  # Convertit "Perturbation" en noms de lignes


# Convertir en matrice (sans noms)
matrice <- as.matrix(freq_table)  # On exclut la colonne Perturbation
print(matrice)
       Before After
coupes     78    92
feu         8    12
tbe        36     2

test Khi2 d’indépendance

resultat_chi2 <- chisq.test(matrice)
print(resultat_chi2)

    Pearson's Chi-squared test

data:  matrice
X-squared = 31.406, df = 2, p-value =
1.515e-07

Conclusion : L’association entre la perturbation et la position n’est pas due au hasard. Il existe une différence systématique dans la répartition des abondances selon les positions After et Before.

Verifier la contribution des cellules du tableau

resultat_chi2 <- chisq.test(matrice)
print(resultat_chi2$stdres)  # Résidus standardisés
          Before     After
coupes -3.952761  3.952761
feu    -1.268144  1.268144
tbe     5.581841 -5.581841

Interpretation : Résidus > |2| = cellules contribuant fortement au χ². tbe (After) : -5.58 → Abondance beaucoup plus faible que prévu en position After. tbe (Before) : 5.2 → Abondance beaucoup plus élevée que prévu en position Before. feu (After) : 2.3 → Abondance légèrement plus élevée que prévu en position After.

Attendus vs observés

#résultats attendus/théoriques

resultat_chi2$expected
         Before     After
coupes 90.96491 79.035088
feu    10.70175  9.298246
tbe    20.33333 17.666667
# résultats observés

resultat_chi2$observed
       Before After
coupes     78    92
feu         8    12
tbe        36     2
# différence entre résultats attendus et observés

resultat_chi2$residual
           Before      After
coupes -1.3593542  1.4583428
feu    -0.8258827  0.8860237
tbe     3.4743400 -3.7273425

Graphique

# Mosaic plot (montre les écarts par rapport à l'indépendance)
Assosplot_chi2 <- mosaicplot(resultat_chi2$observed, main = "Association Perturbation × Position")

Assosplot_chi2
NULL
# Barplot des abondances (par perturbation)
barplot(t(matrice), beside = TRUE,
               legend.text = colnames(matrice),
               ylab = "Fréquence", xlab = "Perturbation",
               main = "Fréquence par Perturbation")

Determiner le pourcentage de contribution des résidus

# Observed counts
observed_counts <- resultat_chi2$observed

# Expected counts
expected_counts <- resultat_chi2$expected


print(observed_counts)
       Before After
coupes     78    92
feu         8    12
tbe        36     2
# Calculate contribution to chi-square statistic
contributions <- (observed_counts - expected_counts)^2 / expected_counts

# Calculate percentage contributions
total_chi_square <- resultat_chi2$statistic
percentage_contributions <- 100 * contributions / total_chi_square

# Print percentage contributions
print("Percentage Contributions:")
[1] "Percentage Contributions:"
print(round(percentage_contributions, 2))
       Before After
coupes   5.88  6.77
feu      2.17  2.50
tbe     38.44 44.24

graphique de la contribution des résidus

# Préparation des données
df_heatmap <- melt(percentage_contributions)
colnames(df_heatmap) <- c("Row", "Column", "Value")

# Création du heatmap avec valeurs affichées
ggplot(df_heatmap, aes(x = Column, y = Row, fill = Value)) +
  geom_tile() +
  geom_text(aes(label = round(Value, 1)),  # Affiche les valeurs arrondies à 1 décimale
            color = "black",  # Couleur du texte (noir pour contraste)
            size = 3.5,       # Taille du texte
            vjust = -0.5) +   # Positionne le texte au-dessus des tuiles
  scale_fill_gradientn(colors = c("gray80", "gray60", "gray40"),
                       name = "Contribution (%)") +
  theme_minimal() +
  ggtitle("Percentage Contribution to Chi-Square Statistic") +
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5),  # Étiquettes horizontales
    plot.title = element_text(hjust = 0.5, face = "bold")  # Titre centré
  )

# Sauvegarde du graphique
ggsave("heatmap_contributions.png", width = 18, height = 10, dpi = 300)

BOP

PRÉPARATION DES DONNÉES

BD_finale_10ans_BOP <- BD_finale_10ans %>% 
                       filter(Esp=="BOP")


# Imaginons que vous vouliez forcer l'ordre "before" puis "after" pour la variable Position
BD_finale_10ans_BOP$Site <- factor(BD_finale_10ans_BOP$Site, 
                                        levels =  c("T87","T77","T83", "T70", "T58","T52",  "T93",  "T94"))

BD_finale_10ans_BOP$Position <- factor(BD_finale_10ans_BOP$Position, 
                                        levels =  c("Before","After"))

ggplot(BD_finale_10ans_BOP, aes(x = Position, y = Abond_totale)) +
            geom_boxplot(aes(color = Position)) +
            theme(legend.position = "top")

ggplot(BD_finale_10ans_BOP, aes(x = Site, y = Abond_totale)) +
            geom_boxplot(aes(color = Position)) +
            theme(legend.position = "top")

TEST DE CHI2

CALCUL DES FRÉQUENCES PONDERÉE

# 1. Calculer la somme des abondances par Perturbation et Position
freq_table_BOP <- BD_finale_10ans_BOP %>%
  group_by(Perturbation, Position) %>%
  summarise(Somme_Abond = sum(Abond_totale, na.rm = TRUE)) %>%
  pivot_wider(
    names_from = Position,
    values_from = Somme_Abond,
    values_fill = 0  # Remplit les valeurs manquantes par 0
  ) %>% 
   column_to_rownames("Perturbation")  # Convertit "Perturbation" en noms de lignes


# Convertir en matrice (sans noms)
matrice_BOP <- as.matrix(freq_table_BOP)  # On exclut la colonne Perturbation
print(matrice_BOP)
       Before After
coupes     49    74
feu        14    30
tbe        47     4

test Khi2 d’indépendance

resultat_chi2_BOP <- chisq.test(matrice_BOP)
print(resultat_chi2_BOP)

    Pearson's Chi-squared test

data:  matrice_BOP
X-squared = 47.14, df = 2, p-value =
5.803e-11

Conclusion : L’association entre la perturbation et la position n’est pas due au hasard. Il existe une différence systématique dans la répartition des abondances selon les positions After et Before.

Verifier la contribution des cellules du tableau

resultat_chi2_BOP <- chisq.test(matrice_BOP)
print(resultat_chi2_BOP$stdres)  # Résidus standardisés
          Before     After
coupes -3.568993  3.568993
feu    -2.768131  2.768131
tbe     6.804875 -6.804875

Interpretation : Résidus > |2| = cellules contribuant fortement au χ². tbe (After) : -6.8 → Abondance beaucoup plus faible que prévu en position After. tbe (Before) : 6.8 → Abondance beaucoup plus élevée que prévu en position Before. feu (After) : 2.7 → Abondance légèrement plus élevée que prévu en position After.

Attendus vs observés

#résultats attendus/théoriques

resultat_chi2_BOP$expected
         Before    After
coupes 62.06422 60.93578
feu    22.20183 21.79817
tbe    25.73394 25.26606
# résultats observés

resultat_chi2_BOP$observed
       Before After
coupes     49    74
feu        14    30
tbe        47     4
# différence entre résultats attendus et observés

resultat_chi2_BOP$residual
          Before     After
coupes -1.658299  1.673583
feu    -1.740671  1.756714
tbe     4.192120 -4.230758

Graphique

# Mosaic plot (montre les écarts par rapport à l'indépendance)
Assosplot_chi2_BOP <- mosaicplot(resultat_chi2_BOP$observed, main = "Association Perturbation × Position")

Assosplot_chi2_BOP
NULL
# Barplot des abondances (par perturbation)
barplot(t(matrice_BOP), beside = TRUE,
               legend.text = colnames(matrice),
               ylab = "Fréquence", xlab = "Perturbation")

Determiner le pourcentage de contribution des résidus

# Observed counts
observed_counts_BOP <- resultat_chi2_BOP$observed

# Expected counts
expected_counts_BOP <- resultat_chi2_BOP$expected


print(observed_counts_BOP)
       Before After
coupes     49    74
feu        14    30
tbe        47     4
# Calculate contribution to chi-square statistic
contributions_BOP <- (observed_counts_BOP - expected_counts_BOP)^2 / expected_counts_BOP

# Calculate percentage contributions
total_chi_square_BOP <- resultat_chi2_BOP$statistic
percentage_contributions_BOP <- 100 * contributions_BOP / total_chi_square_BOP

# Print percentage contributions
print("Percentage Contributions:")
[1] "Percentage Contributions:"
print(round(percentage_contributions_BOP, 2))
       Before After
coupes   5.83  5.94
feu      6.43  6.55
tbe     37.28 37.97

graphique de la contribution des résidus

# Préparation des données
df_heatmap_BOP <- melt(percentage_contributions_BOP)
colnames(df_heatmap_BOP) <- c("Row", "Column", "Value")

# Création du heatmap avec valeurs affichées
ggplot(df_heatmap_BOP, aes(x = Column, y = Row, fill = Value)) +
  geom_tile() +
  geom_text(aes(label = round(Value, 1)),  # Affiche les valeurs arrondies à 1 décimale
            color = "black",  # Couleur du texte (noir pour contraste)
            size = 3.5,       # Taille du texte
            vjust = -0.5) +   # Positionne le texte au-dessus des tuiles
  scale_fill_gradientn(colors = c("gray80", "gray60", "gray40"),
                       name = "Contribution (%)") +
  theme_minimal() +
  ggtitle("Percentage Contribution to Chi-Square Statistic") +
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5),  # Étiquettes horizontales
    plot.title = element_text(hjust = 0.5, face = "bold")  # Titre centré
  )

# Sauvegarde du graphique
ggsave("heatmap_contributions_BOP.png", width = 18, height = 10, dpi = 300)

SAB

PRÉPARATION DES DONNÉES

BD_finale_10ans_SAB <- BD_finale_10ans %>% 
                       filter(Esp=="SAB")


# Imaginons que vous vouliez forcer l'ordre "before" puis "after" pour la variable Position
BD_finale_10ans_SAB$Site <- factor(BD_finale_10ans_SAB$Site, 
                                        levels =  c("T87","T77","T83", "T70", "T58","T52",  "T93",  "T94"))

BD_finale_10ans_SAB$Position <- factor(BD_finale_10ans_SAB$Position, 
                                        levels =  c("Before","After"))

ggplot(BD_finale_10ans_SAB, aes(x = Position, y = Abond_totale)) +
            geom_boxplot(aes(color = Position)) +
            theme(legend.position = "top")

ggplot(BD_finale_10ans_SAB, aes(x = Site, y = Abond_totale)) +
            geom_boxplot(aes(color = Position)) +
            theme(legend.position = "top")

TEST DE CHI2

CALCUL DES FRÉQUENCES PONDERÉE

# 1. Calculer la somme des abondances par Perturbation et Position
freq_table_SAB <- BD_finale_10ans_SAB %>%
  group_by(Perturbation, Position) %>%
  summarise(Somme_Abond = sum(Abond_totale, na.rm = TRUE)) %>%
  pivot_wider(
    names_from = Position,
    values_from = Somme_Abond,
    values_fill = 0  # Remplit les valeurs manquantes par 0
  ) %>% 
   column_to_rownames("Perturbation")  # Convertit "Perturbation" en noms de lignes


# Convertir en matrice (sans noms)
matrice_SAB <- as.matrix(freq_table_SAB)  # On exclut la colonne Perturbation
print(matrice_SAB)
       Before After
coupes     30    53
feu         3    20
tbe        82     6

test Khi2 d’indépendance

resultat_chi2_SAB <- chisq.test(matrice_SAB)
print(resultat_chi2_SAB)

    Pearson's Chi-squared test

data:  matrice_SAB
X-squared = 80.673, df = 2, p-value <
2.2e-16

Conclusion : L’association entre la perturbation et la position n’est pas due au hasard. Il existe une différence systématique dans la répartition des abondances selon les positions After et Before.

Verifier la contribution des cellules du tableau

resultat_chi2_SAB <- chisq.test(matrice_SAB)
print(resultat_chi2_SAB$stdres)  # Résidus standardisés
          Before     After
coupes -5.671059  5.671059
feu    -4.807017  4.807017
tbe     8.757349 -8.757349

Interpretation : Résidus > |2| = cellules contribuant fortement au χ². tbe (After) : -8.7 → Abondance beaucoup plus faible que prévu en position After. tbe (Before) : 8.7 → Abondance beaucoup plus élevée que prévu en position Before. feu (After) : 4.8 → Abondance légèrement plus élevée que prévu en position After.

Attendus vs observés

#résultats attendus/théoriques

resultat_chi2_SAB$expected
         Before     After
coupes 49.20103 33.798969
feu    13.63402  9.365979
tbe    52.16495 35.835052
# résultats observés

resultat_chi2_SAB$observed
       Before After
coupes     30    53
feu         3    20
tbe        82     6
# différence entre résultats attendus et observés

resultat_chi2_SAB$residual
          Before     After
coupes -2.737395  3.302728
feu    -2.879954  3.474729
tbe     4.130831 -4.983940

Graphique

# Mosaic plot (montre les écarts par rapport à l'indépendance)
Assosplot_chi2_SAB <- mosaicplot(resultat_chi2_SAB$observed, main = "Association Perturbation × Position")

Assosplot_chi2_SAB
NULL
# Barplot des abondances (par perturbation)
barplot(t(matrice_SAB), beside = TRUE,
               legend.text = colnames(matrice),
               ylab = "Fréquence", xlab = "Perturbation")

Determiner le pourcentage de contribution des résidus

# Observed counts
observed_counts_SAB <- resultat_chi2_SAB$observed

# Expected counts
expected_counts_SAB <- resultat_chi2_SAB$expected


print(observed_counts_SAB)
       Before After
coupes     30    53
feu         3    20
tbe        82     6
# Calculate contribution to chi-square statistic
contributions_SAB <- (observed_counts_SAB - expected_counts_SAB)^2 / expected_counts_SAB

# Calculate percentage contributions
total_chi_square_SAB <- resultat_chi2_SAB$statistic
percentage_contributions_SAB <- 100 * contributions_SAB / total_chi_square_SAB

# Print percentage contributions
print("Percentage Contributions:")
[1] "Percentage Contributions:"
print(round(percentage_contributions_SAB, 2))
       Before After
coupes   9.29 13.52
feu     10.28 14.97
tbe     21.15 30.79

graphique de la contribution des résidus

# Préparation des données
df_heatmap_SAB <- melt(percentage_contributions_SAB)
colnames(df_heatmap_SAB) <- c("Row", "Column", "Value")

# Création du heatmap avec valeurs affichées
ggplot(df_heatmap_SAB, aes(x = Column, y = Row, fill = Value)) +
  geom_tile() +
  geom_text(aes(label = round(Value, 1)),  # Affiche les valeurs arrondies à 1 décimale
            color = "black",  # Couleur du texte (noir pour contraste)
            size = 3.5,       # Taille du texte
            vjust = -0.5) +   # Positionne le texte au-dessus des tuiles
  scale_fill_gradientn(colors = c("gray80", "gray60", "gray40"),
                       name = "Contribution (%)") +
  theme_minimal() +
  ggtitle("Percentage Contribution to Chi-Square Statistic") +
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5),  # Étiquettes horizontales
    plot.title = element_text(hjust = 0.5, face = "bold")  # Titre centré
  )

# Sauvegarde du graphique
ggsave("heatmap_contributions_SAB.png", width = 18, height = 10, dpi = 300)