La chlordécone est un pesticide organochloré utilisé aux Antilles françaises jusqu’au début des années 1990. En raison de sa forte persistance dans les sols, de sa mobilité dans les eaux et de sa bioaccumulation dans les chaînes alimentaires, elle constitue aujourd’hui un enjeu sanitaire, environnemental et socio-économique majeur.
Problématique :
Quels facteurs environnementaux et agricoles expliquent les niveaux
de chlordécone observés, et quelles zones/parcelles doivent être
priorisées pour la surveillance ?
#install.packages("VIM") # si pas encore installé
library(VIM)
## Warning: le package 'VIM' a été compilé avec la version R 4.3.3
## Le chargement a nécessité le package : colorspace
## Warning: le package 'colorspace' a été compilé avec la version R 4.3.3
## Le chargement a nécessité le package : grid
## VIM is ready to use.
## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues
##
## Attachement du package : 'VIM'
## L'objet suivant est masqué depuis 'package:datasets':
##
## sleep
library(readxl)
## Warning: le package 'readxl' a été compilé avec la version R 4.3.3
library(dplyr)
## Warning: le package 'dplyr' a été compilé avec la version R 4.3.3
##
## Attachement du package : 'dplyr'
## Les objets suivants sont masqués depuis 'package:stats':
##
## filter, lag
## Les objets suivants sont masqués depuis 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(stringr)
library(lubridate)
## Warning: le package 'lubridate' a été compilé avec la version R 4.3.3
##
## Attachement du package : 'lubridate'
## Les objets suivants sont masqués depuis 'package:base':
##
## date, intersect, setdiff, union
BaseCLD2026 <- read_excel("D:/Downloads/BaseCLD2026.xlsx")
summary(BaseCLD2026)
## ID ANNEE COMMU_LAB RAIN
## Min. :20003 Min. :2010 Length:31126 Length:31126
## 1st Qu.:21868 1st Qu.:2013 Class :character Class :character
## Median :23412 Median :2016 Mode :character Mode :character
## Mean :38451 Mean :2015
## 3rd Qu.:62784 3rd Qu.:2018
## Max. :63714 Max. :2019
##
## Sol_simple type_sol Date_prelevement Date_enregistrement
## Length:31126 Length:31126 Length:31126 Length:31126
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Date_analyse Operateur_chld Taux_Chlordecone Operateur_5b
## Length:31126 Length:31126 Length:31126 Length:31126
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Taux_5b_hydro histoBanane_Histo_ban mnt_tpi_mean mnt_tri_mean
## Length:31126 Min. :1.000 Length:31126 Length:31126
## Class :character 1st Qu.:1.000 Class :character Class :character
## Mode :character Median :2.000 Mode :character Mode :character
## Mean :1.673
## 3rd Qu.:2.000
## Max. :3.000
## NA's :17983
## mnt_rugosite_mean mnt_ombrage_mean mnt_exposition_mean mnt_pente_mean
## Length:31126 Length:31126 Length:31126 Length:31126
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## X Y
## Length:31126 Length:31126
## Class :character Class :character
## Mode :character Mode :character
##
##
##
##
Observation : Le résumé statistique montre que plusieurs variables supposées quantitatives sont importées au format character. On observe également un volume important de valeurs manquantes et la présence potentielle de valeurs aberrantes. Un nettoyage des données est donc nécessaire avant toute analyse.
1) Conversion des variables numériques Certaines variables (Taux_Chlordecone,Taux_5b_hydro, pluviométrie, indices topographiques, pente, coordonnées) doivent être converties en valeurs numériques afin de permettre les analyses statistiques et spatiales.
vars_num <- c("Taux_Chlordecone","Taux_5b_hydro","mnt_tpi_mean", "mnt_tri_mean",
"mnt_rugosite_mean", "mnt_ombrage_mean",
"mnt_exposition_mean", "mnt_pente_mean","X", "Y")
BaseCLD2026 <- BaseCLD2026 %>%
mutate(across(all_of(vars_num), ~ as.numeric(as.character(.))))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(all_of(vars_num), ~as.numeric(as.character(.)))`.
## Caused by warning:
## ! NAs introduits lors de la conversion automatique
2) Quantification des valeurs manquantes
na_count <-colSums(is.na(BaseCLD2026))
na_count[na_count>0]
## COMMU_LAB Sol_simple type_sol
## 298 74 2609
## Taux_5b_hydro histoBanane_Histo_ban mnt_tpi_mean
## 23 17983 28
## mnt_tri_mean mnt_rugosite_mean mnt_ombrage_mean
## 28 28 28
## mnt_exposition_mean mnt_pente_mean
## 28 28
summary(BaseCLD2026)
## ID ANNEE COMMU_LAB RAIN
## Min. :20003 Min. :2010 Length:31126 Length:31126
## 1st Qu.:21868 1st Qu.:2013 Class :character Class :character
## Median :23412 Median :2016 Mode :character Mode :character
## Mean :38451 Mean :2015
## 3rd Qu.:62784 3rd Qu.:2018
## Max. :63714 Max. :2019
##
## Sol_simple type_sol Date_prelevement Date_enregistrement
## Length:31126 Length:31126 Length:31126 Length:31126
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Date_analyse Operateur_chld Taux_Chlordecone Operateur_5b
## Length:31126 Length:31126 Min. : 0.0010 Length:31126
## Class :character Class :character 1st Qu.: 0.0024 Class :character
## Mode :character Mode :character Median : 0.0033 Mode :character
## Mean : 0.6677
## 3rd Qu.: 0.4100
## Max. :17.3500
##
## Taux_5b_hydro histoBanane_Histo_ban mnt_tpi_mean mnt_tri_mean
## Min. :0.0000 Min. :1.000 Min. :-29.0658 Min. : 0.000
## 1st Qu.:0.0010 1st Qu.:1.000 1st Qu.: -0.8750 1st Qu.: 2.375
## Median :0.0033 Median :2.000 Median : 0.1200 Median : 3.822
## Mean : Inf Mean :1.673 Mean : 0.2118 Mean : 4.313
## 3rd Qu.:0.0150 3rd Qu.:2.000 3rd Qu.: 1.2669 3rd Qu.: 5.625
## Max. : Inf Max. :3.000 Max. : 26.0119 Max. :29.358
## NA's :23 NA's :17983 NA's :28 NA's :28
## mnt_rugosite_mean mnt_ombrage_mean mnt_exposition_mean mnt_pente_mean
## Min. : 0.00 Min. : 24.79 Min. : 0.00 Min. : 0.00
## 1st Qu.: 7.50 1st Qu.:159.39 1st Qu.: 89.76 1st Qu.: 10.80
## Median :12.00 Median :177.10 Median :169.31 Median : 17.90
## Mean :13.44 Mean :175.03 Mean :170.17 Mean : 20.27
## 3rd Qu.:17.76 3rd Qu.:193.29 3rd Qu.:249.27 3rd Qu.: 27.10
## Max. :85.40 Max. :255.00 Max. :359.14 Max. :124.95
## NA's :28 NA's :28 NA's :28 NA's :28
## X Y
## Min. :690703 Min. :1594393
## 1st Qu.:703247 1st Qu.:1614066
## Median :713492 Median :1624496
## Mean :713040 Mean :1624037
## 3rd Qu.:721006 3rd Qu.:1633394
## Max. :734339 Max. :1645435
##
3) Détection et traitement des valeurs aberrantes Les valeurs Inf ne doivent pas être remplacées par la moyenne car elles resultent d’une dividion par 0. Elles doivent être traitées comme des valeurs manquantes (NA).
sum(is.infinite(BaseCLD2026$Taux_5b_hydro))
## [1] 586
BaseCLD2026 <- BaseCLD2026 %>%
mutate(
Taux_5b_hydro = ifelse(is.infinite(Taux_5b_hydro),
NA,
Taux_5b_hydro)
)
Gestion des formats de dates
BaseCLD2026 <- BaseCLD2026 %>%
mutate(Date_prelevement = dmy(Date_prelevement),Date_enregistrement = dmy(Date_enregistrement), Date_analyse=dmy(Date_analyse))
Resultat des données nettoyées
str(BaseCLD2026)
## tibble [31,126 × 22] (S3: tbl_df/tbl/data.frame)
## $ ID : num [1:31126] 20143 20143 20143 20143 20143 ...
## $ ANNEE : num [1:31126] 2010 2010 2010 2010 2010 ...
## $ COMMU_LAB : chr [1:31126] "GROS-MORNE" "GROS-MORNE" "GROS-MORNE" "GROS-MORNE" ...
## $ RAIN : chr [1:31126] "2000-3000" "2000-3000" "2000-3000" "2000-3000" ...
## $ Sol_simple : chr [1:31126] "Andosol" "Andosol" "Andosol" "Andosol" ...
## $ type_sol : chr [1:31126] "Intergrades Sols … allophane relativement ‚volu‚s/Ferrisols compacts" "Intergrades Sols … allophane relativement ‚volu‚s/Ferrisols compacts" "Intergrades Sols … allophane relativement ‚volu‚s/Ferrisols compacts" "Intergrades Sols … allophane relativement ‚volu‚s/Ferrisols compacts" ...
## $ Date_prelevement : Date[1:31126], format: "2007-05-24" "2007-05-24" ...
## $ Date_enregistrement : Date[1:31126], format: "2007-05-24" "2007-05-24" ...
## $ Date_analyse : Date[1:31126], format: "2007-05-24" "2007-05-24" ...
## $ Operateur_chld : chr [1:31126] "=" "=" "=" "=" ...
## $ Taux_Chlordecone : num [1:31126] 4.6 4.6 4.6 4.6 4.6 4.6 4.6 0.0033 0.0033 0.0033 ...
## $ Operateur_5b : chr [1:31126] "=" "=" "=" "=" ...
## $ Taux_5b_hydro : num [1:31126] 0.07 0.07 0.07 0.07 0.07 0.07 0.07 NA 0.015 0.015 ...
## $ histoBanane_Histo_ban: num [1:31126] 2 2 3 1 2 3 2 1 1 1 ...
## $ mnt_tpi_mean : num [1:31126] 5.806 5.684 2.239 4.038 0.597 ...
## $ mnt_tri_mean : num [1:31126] 8.03 7.92 7.11 7.53 6.64 ...
## $ mnt_rugosite_mean : num [1:31126] 21.6 20.9 20.1 23.4 20.2 ...
## $ mnt_ombrage_mean : num [1:31126] 131 135 139 122 134 ...
## $ mnt_exposition_mean : num [1:31126] 79.4 77 76.1 92.4 83.9 ...
## $ mnt_pente_mean : num [1:31126] 39 38.1 35 38.3 33.9 ...
## $ X : num [1:31126] 714301 714304 714309 714294 714304 ...
## $ Y : num [1:31126] 1626344 1626354 1626360 1626321 1626341 ...
summarise(BaseCLD2026)
## # A tibble: 1 × 0
Afin enrichir l’information avant analyse, on va créer des nouvelles variables À partir de la date :ANNEE,MOIS,TRIMESTRE,SAISON (sèche / humide),et recoder les variables categorielles et traitement des NA avant la classification
BaseCLD2026 <- BaseCLD2026 %>%
mutate(
ANNEE = year(Date_prelevement),
MOIS = month(Date_prelevement),
SAISON = case_when(
MOIS %in% c(12,1,2,3,4) ~ "Carême (saison sèche)",
MOIS %in% c(5,6,7,8,9,10,11) ~ "Hivernage (saison humide)"
),
TRIMESTRE = quarter(Date_prelevement)
)
Les méthodes de classification non supervisée (CAH et k-means) nécessitant l’absence de valeurs manquantes et des variables numériques comparables, une table dédiée va etre construite. Les variables environnementales continues ont été sélectionnées, les valeurs manquantes vont etre imputées par la médiane et les données standardisées. Selection des variables numériques
vars_clust <- c(
"Taux_Chlordecone",
"Taux_5b_hydro",
"mnt_tpi_mean",
"mnt_tri_mean",
"mnt_rugosite_mean",
"mnt_ombrage_mean",
"mnt_exposition_mean",
"mnt_pente_mean",
"X",
"Y"
)
Base_clust <- BaseCLD2026 %>%
select(all_of(vars_clust)) %>%
mutate(across(everything(),
~ ifelse(is.na(.), median(., na.rm = TRUE), .)))
Imputation KNN Les valeurs manquantes des variables environnementales ont été imputées par la méthode des K plus proches voisins (KNN, k = 5). Cette méthode non paramétrique estime les valeurs manquantes à partir des observations les plus similaires selon la distance euclidienne dans l’espace des variables standardisées. Elle permet de préserver la structure multivariée des données avant application des méthodes de classification non supervisée (CAH et k-means).
Base_clust <- kNN(Base_clust, k = 5, imp_var = FALSE)
## Warning in kNN(Base_clust, k = 5, imp_var = FALSE): Nothing to impute, because
## no NA are present (also after using makeNA)
# 1) Standardiser (les distances KNN sont sensibles aux échelles)
Base_clust <- as.data.frame(scale(Base_clust))
#CLASSIFICATION ASCENDANTE HIÉRARCHIQUE (CAH - Ward.D2)
set.seed(123)
# Taille d’échantillon (3000 )
n_sample <- min(3000, nrow(Base_clust))
idx <- sample(seq_len(nrow(Base_clust)), n_sample)
Base_sample <- Base_clust[idx, ]
dist_mat <- dist(Base_sample, method = "euclidean")
cah <- hclust(dist_mat, method = "ward.D2")
plot(cah, labels = FALSE, hang = -1,
main = paste("Dendrogramme - CAH (Ward.D2) sur échantillon n =", n_sample),
xlab = "", ylab = "Hauteur (inertie)")
L’analyse du dendrogramme issu de la Classification Ascendante Hiérarchique (méthode de Ward.D2) permet d’identifier les ruptures majeures dans la structure des données.
Sur le dendrogramme, on observe des sauts importants de hauteur autour de :
Ces sauts correspondent à des fusions entre groupes relativement éloignés. Plus la hauteur est élevée, plus les groupes fusionnés sont dissemblables.
Une coupe du dendrogramme autour d’une hauteur d’environ 55 conduit à la formation de quatre groupes distincts :
Ce choix représente un bon compromis entre : - homogénéité intra-cluster, - séparation inter-clusters, - interprétabilité des profils environnementaux.
Après avoir exploré la structure des données avec la CAH et retenu un
nombre de classes k = 4, on applique l’algorithme
k-means à l’ensemble des observations.
Le k-means partitionne les données en k groupes en minimisant
l’inertie intra-classe (distance aux centroïdes).
Comme la méthode dépend d’une initialisation aléatoire, on fixe une
graine (set.seed) et on utilise plusieurs initialisations
(nstart) afin d’obtenir une solution plus stable.
k <- 4
set.seed(123)
kmeans_res <- kmeans(Base_clust, centers = k, nstart = 25)
# Effectifs par cluster
table(kmeans_res$cluster)
##
## 1 2 3 4
## 5290 4085 12484 9267
#Ajout de l’étiquette de cluster à la base
BaseCLD2026 <- BaseCLD2026 %>%
mutate(cluster_kmeans = factor(kmeans_res$cluster))
##Interprétation : profils moyens par cluster
profil_kmeans <- BaseCLD2026 %>%
group_by(cluster_kmeans) %>%
summarise(across(all_of(vars_clust), ~ mean(.x, na.rm = TRUE))) %>%
arrange(cluster_kmeans)
profil_kmeans
## # A tibble: 4 × 11
## cluster_kmeans Taux_Chlordecone Taux_5b_hydro mnt_tpi_mean mnt_tri_mean
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 1 0.650 5.21 0.641 7.70
## 2 2 0.600 3.49 0.293 7.00
## 3 3 0.196 7.21 -0.00604 2.91
## 4 4 1.34 1.02 0.224 3.08
## # ℹ 6 more variables: mnt_rugosite_mean <dbl>, mnt_ombrage_mean <dbl>,
## # mnt_exposition_mean <dbl>, mnt_pente_mean <dbl>, X <dbl>, Y <dbl>
Ce cluster présente une forte pente, une rugosité élevée et une
complexité topographique marquée.
Les concentrations en chlordécone sont intermédiaires.
Ces parcelles correspondent à des zones de versant où le ruissellement
peut favoriser le transport des polluants.
Les parcelles de ce groupe présentent une pente élevée et une forte
exposition du terrain.
Les concentrations restent modérées.
L’orientation et les conditions microclimatiques peuvent influencer la
dynamique de dégradation et de redistribution.
Ce cluster se caractérise par une faible pente mais un indice
hydrologique élevé.
Les niveaux de contamination sont les plus faibles.
Ces parcelles pourraient être associées à une dynamique de drainage ou
de dilution.
Ce groupe présente les concentrations les plus élevées en
chlordécone.
Il correspond à des parcelles relativement planes et peu
accidentées.
Ces conditions topographiques favorisent l’accumulation et la
persistance du contaminant.
Ce cluster constitue une priorité pour la surveillance environnementale.
Categorisation des clusters
BaseCLD2026 <- BaseCLD2026 %>%
mutate(
cluster_label = case_when(
cluster_kmeans == 1 ~ "Versants actifs - contamination modérée",
cluster_kmeans == 2 ~ "Relief exposé - contamination intermédiaire",
cluster_kmeans == 3 ~ "Zones hydrologiques - faible contamination",
cluster_kmeans == 4 ~ "Zones d'accumulation - contamination élevée"
)
)
ggplot(BaseCLD2026, aes(x = cluster_kmeans, y = Taux_Chlordecone)) +
geom_boxplot() +
labs(title = "Taux de chlordécone par cluster (k-means)",
x = "Cluster", y = "Taux_Chlordecone")
Visualisation : projection PCA des clusters
pca <- prcomp(Base_clust, center = TRUE, scale. = FALSE)
pca_df <- data.frame(
PC1 = pca$x[, 1],
PC2 = pca$x[, 2],
cluster = BaseCLD2026$cluster_kmeans
)
ggplot(pca_df, aes(PC1, PC2, color = cluster)) +
geom_point(alpha = 0.5) +
labs(title = "Projection PCA - clusters k-means")
Analyse croisée : saison, commune, sol
# Saison vs cluster
if ("SAISON" %in% names(BaseCLD2026)) {
table(BaseCLD2026$cluster_kmeans, BaseCLD2026$SAISON, useNA = "ifany")
}
##
## Carême (saison sèche) Hivernage (saison humide)
## 1 2032 3258
## 2 1507 2578
## 3 4650 7834
## 4 3787 5480
# Commune vs cluster (si COMMU_LAB existe)
if ("COMMU_LAB" %in% names(BaseCLD2026)) {
BaseCLD2026 %>%
filter(!is.na(COMMU_LAB)) %>%
count(cluster_kmeans, COMMU_LAB, sort = TRUE) %>%
group_by(cluster_kmeans) %>%
slice_head(n = 5)
}
## # A tibble: 20 × 3
## # Groups: cluster_kmeans [4]
## cluster_kmeans COMMU_LAB n
## <fct> <chr> <int>
## 1 1 MORNE-ROUGE(LE) 607
## 2 1 SAINT-JOSEPH 516
## 3 1 LORRAIN(LE) 424
## 4 1 GROS-MORNE 414
## 5 1 ROBERT(LE) 310
## 6 2 MORNE-ROUGE(LE) 466
## 7 2 LORRAIN(LE) 352
## 8 2 ROBERT(LE) 279
## 9 2 SAINT-ESPRIT 246
## 10 2 GROS-MORNE 212
## 11 3 LAMENTIN(LE) 1672
## 12 3 DUCOS 1423
## 13 3 FRANCOIS(LE) 1177
## 14 3 ROBERT(LE) 1136
## 15 3 SAINT-JOSEPH 1133
## 16 4 MORNE-ROUGE(LE) 2652
## 17 4 LORRAIN(LE) 987
## 18 4 SAINTE-MARIE 888
## 19 4 BASSE-POINTE 836
## 20 4 GROS-MORNE 793
ggplot(BaseCLD2026, aes(x = X, y = Y, color = cluster_kmeans)) +
geom_point(alpha = 0.5, size = 1) +
labs(title = "Répartition spatiale des clusters (k-means)",
x = "X", y = "Y")