INTRODUCTION

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 ?


IMPORTATION DES DONNÉES

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

DESCRIPTION STATISTIQUE DES DONNEES

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.

NETTOYAGE DES DONNÉES

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

CREATION DES NOUVELLES VARIABLES (feature engineering)

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)
  )

EXTRACTION DU JEU DE DONNÉES SPÉCIFIQUE POUR LA CLASSIFICATION

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)")

Choix du nombre de clusters (k)

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.

Choix : k = 4

Une coupe du dendrogramme autour d’une hauteur d’environ 55 conduit à la formation de quatre groupes distincts :

  • La branche principale de gauche se subdivise en deux sous-groupes homogènes.
  • La branche principale de droite se subdivise également en deux sous-groupes.
  • On obtient ainsi 4 clusters bien séparés.

Ce choix représente un bon compromis entre : - homogénéité intra-cluster, - séparation inter-clusters, - interprétabilité des profils environnementaux.


Visualisation des coupes

Classification par k-means

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.

Application de k-means (k = 4)

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>

Cluster 1 – Versants actifs (contamination modérée)

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.

Cluster 2 – Relief exposé (contamination intermédiaire)

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.

Cluster 3 – Zones hydrologiquement actives (faible contamination)

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.

Cluster 4 – Zones d’accumulation (contamination élevée)

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"
    )
  )

Visualisation : boxplot du taux de chlordécone par cluster

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")