#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)
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)
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)
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()
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)
# 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>
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")
# 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
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.
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.
#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
# 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")
# 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
# 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)
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")
# 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
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.
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.
#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
# 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")
# 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
# 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)
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")
# 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
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.
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.
#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
# 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")
# 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
# 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)