super projet

---

La majorité des graphiques présentés sont interactifs, n’hésitez pas à passer votre souris dessus pour avoir les valeurs exactes !🙂

Introduction

Notre recherche se concentre sur l’identification des facteurs prédictifs de la démence. Notre variable dépendante, le “CDR” (Clinical Dementia Rate), est un indicateur neurologique établi en 1993 par Morris. Cette variable est ici binaire, avec 0 pour les individus non atteints de démence et 1 pour ceux qui le sont. Nous examinons également plusieurs autres variables potentiellement pertinentes, notamment le sexe, l’âge, le niveau d’éducation, le niveau socio-culturel, le volume cérébral total, le volume intercrânien, le niveau d’atrophie, ainsi que le Mini Mental State Examination (MMSE). Bien que le MMSE soit un test neurologique largement utilisé en raison de sa rapidité, il est toutefois souvent critiqué dans la littérature pour sa moindre précision et sa fiabilité variable. Nous avons l’opportunité de vérifier ces points en l’incluant dans notre étude parallèlement au CDR.

I. Analyse descriptive

Justification

L’analyse descriptive permet de comprendre la distribution et les caractéristiques générales des données, offrant ainsi un aperçu préliminaire des tendances et des associations potentielles entre les variables étudiées.

Il est important de s’assurer que les niveaux des variables sont assez équilibré, car sinon ils pourraient être sous-représentés dans les analyses à venir, si c’est le cas il sera ainsi important de le prendre en compte lors de l’interprétation.

Graphique des variables qualitatives

Pour simplifier la création des graphiques, nous pouvons écrire une fonction qui inclut tous les paramètres nécessaires. Cela nous permettra de générer des ggplot plus rapidement.

fastggplot <- function(data, x_var, title, x_label, y_label) {
  ggplot(data, aes(x = factor(data[[x_var]]), fill = factor(data[[x_var]]))) +
    geom_bar(color = "black", alpha = 0.8) +
    
    labs(title = title, x = x_label, y = y_label, fill = x_var) +
    theme_minimal() +
    theme(axis.text = element_text(size = 12),
          axis.title = element_text(size = 14, face = "bold"),
          plot.title = element_text(size = 16, face = "bold"))
}

CDR <- fastggplot(demence2, "CDR", "Répartition par score CDR", "Score CDR", "Nombre de cas")
CDR %>% plotly::ggplotly()

MMSE <- fastggplot(demence2, "MMSE", "Répartition par score MMSE", "Score MMSE", "Nombre de cas")

MMSE <- MMSE + theme(axis.text.x = element_text(angle = 45, hjust = 1))

MMSE %>% plotly::ggplotly()

EDUC <- fastggplot(demence2, "EDUC", "Répartition par niveau d'éducation", "Niveau d'éducation", "Nombre de cas")

EDUC <- EDUC + theme(axis.text.x = element_text(angle = 45, hjust = 1))

EDUC %>% plotly::ggplotly()

genre = fastggplot(demence2, "sexe", "Répartition par sexe", "Genre", "Nombre de cas")
genre %>% plotly::ggplotly()

Graphique des variables quantitatives

hist_age <- ggplot(demence2, aes(x = Age)) +
  geom_histogram(binwidth = 0.01, fill = "green", color = "black", alpha = 0.8) +
  labs(title = "Distribution de l'âge", x = "", y = "Fréquence") +
  theme_minimal()
hist_age %>% plotly::ggplotly()

hist_nwbv <- ggplot(demence2, aes(x = Volume_cérébral_total)) +
  geom_histogram(binwidth = 0.01, fill = "green", color = "black", alpha = 0.8) +
  labs(title = "Distribution du volume cérébral total", x = "Volume cérébral total", y = "Fréquence") +
  theme_minimal()
hist_nwbv %>% plotly::ggplotly()

hist_etiv <- ggplot(demence2, aes(x = Volume_intercranien)) +
  geom_histogram(binwidth = 100, fill = "orange", color = "black", alpha = 0.8) +
  labs(title = "Distribution du volume intercranien", x = "eTIV", y = "Fréquence") +
  theme_minimal()
hist_etiv %>% plotly::ggplotly()

hist_asf <- ggplot(demence2, aes(x = Indice_atrophie)) +
  geom_histogram(binwidth = 0.1, fill = "blue", color = "black", alpha = 0.8) +
  labs(title = "Distribution de l'indice d'atrophie", x = "Indice_atrophie", y = "Fréquence") +
  theme_minimal()
hist_asf %>% plotly::ggplotly()

Pour aller plus loin dans l’analyse descriptive, on peut également faire des graphiques qui associent notre variable dépendante, le CDR, aux autres.

Graphique associant la variable dépendante (CDR) et les variables indépendantes

CDRSEXE = ggplot(data = demence2, aes(x = CDR, fill = sexe)) +
  geom_bar(position = "dodge") +
  labs(title = "Comparaison de CDR par sexe",
       x = "CDR",
       y = "Nombre de participants",
       fill = "Sexe")
CDRSEXE %>% plotly::ggplotly()

CDREDUC = ggplot(data = demence2, aes(x = CDR, fill = EDUC)) +
  geom_bar(position = "dodge") +
  labs(title = "Comparaison de CDR par niveau d'éducation",
       x = "CDR",
       y = "Nombre de participants",
       fill = "EDUC")

CDREDUC %>% plotly::ggplotly()

CDRSES = ggplot(data = demence2, aes(x = CDR, fill = SES)) +
  geom_bar(position = "dodge") +
  labs(title = "Comparaison de CDR par niveau socio-économique",
       x = "CDR",
       y = "Nombre de participants",
       fill = "SES")

CDRSES %>% plotly::ggplotly()

CDRSES = ggplot(data = demence2, aes(x = CDR, fill = MMSE)) +
  geom_bar(position = "dodge") +
  labs(title = "Comparaison de CDR par MMSE",
       x = "CDR",
       y = "Nombre de participants",
       fill = "MMSE")

CDRSES %>% plotly::ggplotly()

CDReTIV = ggplot(demence2, aes(x = CDR, y = Volume_cérébral_total)) +
  geom_boxplot(fill = "lightgreen", color = "black") +
  labs(x = "CDR", y = "Volume cérébral total", title = "Distribution du volume cérébral total par CDR")

CDReTIV %>% plotly::ggplotly()

CDRASF = ggplot(demence2, aes(x = CDR, y = Indice_atrophie)) +
  geom_boxplot(fill = "darkred", color = "black") +
  labs(x = "CDR", y = "ASF", title = "Distribution de l'indice d'atrophie par CDR")
CDRASF %>% plotly::ggplotly()

Corrélation de point-bisériale, CDR + X

La corrélation point-bisériale mesure la relation entre une variable continue et une variable binaire. Un coefficient proche de +1 ou -1 indique une corrélation forte, tandis qu’une valeur proche de zéro suggère une faible corrélation. L’utilisation de cette méthode peut nous permettre d’avoir un premier aperçu des relations possibles entre notre variable dépendante et nos variables indépendantes.

    FORcorb <- demence2 
    # créer une variable binaire à partir de CDR
    FORcorb$CDR_binaire <- ifelse(FORcorb$CDR == "Dément", 1, 0)
    
    FORcorb$CDR_numeric <- FORcorb$CDR
    
    # sélection des variables pertinentes
    variables_quantitatives <- FORcorb[, c("Age", "Volume_intercranien", "Volume_cérébral_total", "Indice_atrophie")]
    
    variable_binaire <- as.numeric(FORcorb$CDR_numeric)
    
    # créer une liste pour stocker les résultats
    results <- list()
    
    # on calcule les coefficient
    for (variable in names(variables_quantitatives)) {
      correlation <- cor.test(variable_binaire, variables_quantitatives[[variable]], method = "pearson")
      results[[variable]] <- correlation$estimate
    }
    
# transformer en data frame
results_df <- data.frame(
  Variable = c("Age", "Volume_intercranien", "Volume_cérébral_total", "Indice_atrophie"),
  Correlation = unlist(results)
)

    corbis <- ggplot(results_df, aes(x = Variable, y = Correlation, fill = Correlation)) +
      geom_bar(stat = "identity", color = "black") +
      labs(title = "Coefficient de corrélation point-bisériale avec CDR",
           x = "Variable",
           y = "Coefficient de corrélation") +
      ylim(-1, 1) +  # Définir la limite de l'axe y de -1 à 1
      theme_minimal() +  # Utiliser un thème minimal
      theme(axis.text.x = element_text(angle = 45, hjust = 1)) +  # rotation des étiquettes de l'axe x
      theme(panel.grid.major.y = element_line(color = "gray", linetype = "NULL")) +  # Ajouter une grille en pointillés
      scale_fill_gradient2(low = "blue", mid = "white", high = "red", 
                           midpoint = 0, limits = c(-1, 1), name = "Corrélation")  # Échelle de couleur personnalisée
    
   suppressMessages({
  corbis_plotly <- ggplotly(corbis)
})
    # affichage du graphique Plotly
     suppressMessages(corbis_plotly)

Coefficient de Cramér, CDR + X

Le coefficient de Cramér mesure l’association entre deux variables catégorielles. Calculé à partir des tables de contingence, il varie de 0 (aucune association) à 1 (association parfaite). Une valeur proche de 0 indique une faible relation, tandis que près de 1 indique une forte association. Il sert à évaluer la dépendance entre les variables.

# on crée d'abord des tables de contingence pour chaque paire de variables

t1 = table(demence2$sexe, demence2$CDR)
t2 = table(demence2$MMSE,demence2$CDR)
t3 = table(demence2$EDUC,demence2$CDR)
t4 = table(demence2$SES,demence2$CDR)

# calcul des les statistiques d'association pour chaque table
ma1 = assocstats(t1)
ma2 = assocstats(t2)
ma3 = assocstats(t3)
ma4 = assocstats(t4)

data <- data.frame(
  Variable = c("Sexe", "MMSE", "EDUC", "SES"),
  CoefficientdeCramer = c(ma1$cramer, ma2$cramer, ma3$cramer, ma4$cramer)
)
Cramer_V =  c(ma1$cramer, ma2$cramer, ma3$cramer, ma4$cramer)

# créer le graphique à barres
GGCram = ggplot(data, aes(x = Variable, y = CoefficientdeCramer, fill = CoefficientdeCramer)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.7, color = "black") +
  scale_fill_gradient2(low = "darkblue", mid = "white", high = "red", 
                       midpoint = 0.5, limits = c(0, 1), 
                       name = "Indice de Cramér") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Cramér en fonction de CDR",
       x = "Variables",
       y = "Cramér's V") +
  theme_bw()

GGCram = GGCram + theme(axis.text.x = element_text(angle = 45, hjust = 1))

GGCram %>% plotly::ggplotly()

II. CART

Justification

L’arbre de classification et de régression (CART) est un modèle d’apprentissage automatique robuste qui segmente de manière récursive les données en sous-groupes, facilitant ainsi la prise de décisions. Son avantage réside dans sa facilité de lecture et d’interprétation.

Pour chacun des modèles prédictifs créés dans cette étude, une validation croisée a été effectuée en divisant l’échantillon en un ensemble de test (regroupant 20% des données) et un ensemble d’apprentissage (regroupant 80% des données).

En plus des modèles traditionnels, des modèles de type dummy ont été mis en place pour servir de référence lors de l’évaluation des performances des modèles prédictifs, permettant ainsi de comparer leur efficacité.

Il a également été pris en compte des métriques telles que la précision, le rappel, la spécificité et l’AUC (Area Under the Curve) pour évaluer la performance des modèles prédictifs, fournissant ainsi une analyse complète de leur capacité à classifier correctement les données

Finalement, une boucle a été créée pour automatiser le processus d’évaluation des modèles sur différents ensembles de données, garantissant ainsi une analyse exhaustive et reproductible de leur performance. Cette boucle, répétée 100 fois, vise à évaluer les métriques d’évaluation telles que la précision, le rappel, la spécificité et l’AUC, permettant ainsi d’éliminer les variations aléatoires et d’obtenir des résultats robustes.

    nrun <- 100
    res <- c()
    performance_metrics <- list()
    
    # Stocker les noms des métriques
    metric_names <- c("Précision", "Rappel", "Spécificité", "AUC")
    
    for (irun in 1:nrun){
      df <- demence2[sample(nrow(demence2)), ]
      train <- df[1:280, ]
      test <- df[281:354, ]
      
      tree <- rpart(CDR ~ ., train)
      
      preds_dummy <- rep(names(table(train$CDR))[which.max(table(train$CDR))], nrow(test))
      predictions <- predict(tree, test, type="class")
      
      res <- rbind(res, c(
        mean(preds_dummy == test$CDR),
        mean(predictions == test$CDR)
      ))
      
      # calcul des métriques d'évaluation
      conf_matrix <- table(predictions = predictions, actual = test$CDR)
      
      precision <- confusionMatrix(predictions, test$CDR)$byClass["Pos Pred Value"]
      rappel <- confusionMatrix(predictions, test$CDR)$byClass["Sensitivity"]
      specificite <- confusionMatrix(predictions, test$CDR)$byClass["Specificity"]
      roc_curve <- roc(test$CDR, as.numeric(predictions))
      auc_value <- auc(roc_curve)

      performance_metrics[[irun]] <- c(precision, rappel, specificite, auc_value)
    }
    
    # on calcule la moyenne des métriques pour toutes les itérations
    average_metrics <- colMeans(do.call(rbind, performance_metrics), na.rm = TRUE)
    
    # créer un data frame avec les noms des métriques
    metric_summary <- data.frame(
      Métrique = metric_names,
      Moyenne = average_metrics
    )
    
    colnames(res) <- c("dummy", "tree")
    
    means = colMeans(res)
    
    df2 <- data.frame(Model = names(means), Accuracy = means)
    
    
    # créer le ggplot
    CompaCART =   ggplot(df2, aes(x = Model, y = Accuracy, fill = Model)) +
      geom_col() +
      labs(title = "Comparaison entre Dummy et Tree",
           x = "Modèle", y = "Précision") +
      theme_minimal()
    
    CompaCART %>% plotly::ggplotly()

Le dummy a obtenu 0.58, contre 0.73 pour notre modèle, nous pouvons donc continuer notre CART !

Robustesse du modèle

 # créer un ggplot
    recapCART = ggplot(metric_summary, aes(x = Métrique, y = Moyenne, fill = Métrique)) +
      geom_col() +
      labs(title = "CART (100 itérations)",
           x = "", y = "Moyenne") +
        theme_minimal()
    
recapCART <- recapCART + theme(axis.text.x = element_text(angle = 45, hjust = 1))

    recapCART %>% plotly::ggplotly()

Notre modèle a obtenu de très bonnes métriques d’évaluations, avec 0.80 pour l’AUC, 0.83 pour la précision, 0.71 pour le rappel et 0.90 pour la spécificité, cela montre quele modèle a une capacité élevée à discriminer entre les classes positives et négatives, tout en maintenant un bon équilibre entre la précision et le rappel, et en minimisant à la fois les faux positifs et les faux négatifs.

Importance des variables

On peut regarder l’importance des variables dans notre modèles

    tree_global <- rpart(CDR ~ ., df)
# Les importances des variables
importances <- c(Age = 14.44873,
                 EDUC = 33.85687,
                 Indice_atrophie = 16.66886,
                 MMSE = 37.60806,
                 SES = 34.07162,
                 sexe = 22.16193,
                 Volume_intercranien = 22.69176,
                 Volume_cérébral_total = 52.96523)

# Créer un data frame pour les importances
importances_df <- data.frame(Variables = names(importances), Importance = importances)

# Créer le ggplot
ggplot(importances_df, aes(x = Variables, y = Importance)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Importance des variables")

Visualisation de l’arbre de décision

rpart.plot(tree_global)

Interprétation de l’arble :

1. MMSE: Le premier nœud divise les individus en fonction de leur score au Mini-Mental State Examination (MMSE). Un score inférieur à la limite indiquée (non spécifiée ici) mène directement à un diagnostic de démence (probabilité de 1,00 ou 100%).

2. Si le score MMSE est supérieur à la limite, l’arbre se divise en deux branches, une pour les hommes et une pour les femmes (dénoté par “sexe = M”).

3. Pour les Femmes:

   • Le nœud suivant regarde le ‘Volume_intracrânien’. Si c’est moins de 1609, il se divise encore selon l’âge.
     • Si l’âge est >= 79, la probabilité de démence augmente (0,37 ou 37%).
     • Pour ceux avec un ‘Volume_intracrânien’ >= 1609, la division se fait selon le ‘Socioeconomic Status’ (SES).
       • Un SES <= 1,4 conduit à une division suivante basée sur le niveau d’éducation (EDUC).
         • Pour un EDUC entre 14 et 16, la probabilité de démence est plus basse (0,10 ou 10%).

4. Pour les Hommes:

   • Si le ‘Volume_intracrânien’ est moins de 1486, cela amène à une autre division en fonction du ‘Volume_cérébral_total’.
     • Si le ‘Volume_cérébral_total’ est < 0,7, alors la probabilité de démence est de 0,46 (46%).
     • Si le ‘Volume_cérébral_total’ est >= 0,7, et l’âge < 77, la probabilité de démence est plus faible (0,19 ou 19%).
     • Pour un ‘Volume_cérébral_total’ >= 0,77 et un SES entre 3 et 4,5, les probabilités de démence varient de 0,16 (16%) à 0,39 (39%) selon l’âge et le SES.

III. Random Forest

Justification

La forêt aléatoire (Random Forest) est une technique puissante en apprentissage automatique. Elle construit plusieurs arbres de décision lors de l’entraînement et les combine pour obtenir des prédictions plus robustes et précises. Cette méthode segmente les données de manière récursive comme l’arbre de classification et de régression (CART), mais au lieu de s’appuyer sur un seul arbre, elle assemble les résultats de multiples arbres pour produire une prédiction finale.

Nous avons suivi une démarche similaire au CART. Validation croisée avec division en ensembles de test (20%) et d’apprentissage (80%). Des modèles de type dummy ont été utilisés comme référence. Les performances ont été mesurées par précision, rappel, spécificité et AUC. Une boucle a été mise en place, répétée 100 fois pour des résultats robustes.

nrun <- 100
res <- c()
performance_metrics <- list()

# Stocker les noms des métriques
metric_names <- c("Précision", "Rappel", "Spécificité", "AUC")

for (irun in 1:nrun){
  dff <- demence2[sample(nrow(demence2)), ]
  train <- dff[1:280, ]
  test <- dff[281:354, ]
  
  # Utiliser un modèle Random Forest au lieu de CART
  forest <- randomForest(CDR ~ ., data = train)
  
  preds_dummy <- rep(names(table(train$CDR))[which.max(table(train$CDR))], nrow(test))
  predictions <- predict(forest, test)
  
  res <- rbind(res, c(
    mean(preds_dummy == test$CDR),
    mean(predictions == test$CDR)
  ))
  
  # Calcul des métriques d'évaluation
  conf_matrix <- table(predictions = predictions, actual = test$CDR)
  
  precision <- confusionMatrix(predictions, test$CDR)$byClass["Pos Pred Value"]
  rappel <- confusionMatrix(predictions, test$CDR)$byClass["Sensitivity"]
  specificite <- confusionMatrix(predictions, test$CDR)$byClass["Specificity"]
  roc_curve <- roc(test$CDR, as.numeric(predictions))
  auc_value <- auc(roc_curve)
  
  # Stocker les métriques pour chaque itération
  performance_metrics[[irun]] <- c(precision, rappel, specificite, auc_value)}

Robustesse du modèle

colnames(res) <- c("dummy", "forest")


means = colMeans(res)

df2 <- data.frame(Model = names(means), Accuracy = means)

# Créer le ggplot
CompaRF <- ggplot(df2, aes(x = Model, y = Accuracy, fill = Model)) +
  geom_col() +
  labs(title = "Comparaison entre Dummy et Random Forest",
       x = "Modèle", y = "Précision") +
  theme_minimal()

CompaRF %>% plotly::ggplotly()

Le dummy a obtenu 0.58, contre 0.82 pour notre modèle, nous pouvons donc continuer notre analyse !

average_metrics <- colMeans(do.call(rbind, performance_metrics), na.rm = TRUE)

# Créer un data frame avec les noms des métriques
metric_summary <- data.frame(
  Métrique = metric_names,
  Moyenne = average_metrics
)

# Créer un ggplot
recapRF <- ggplot(metric_summary, aes(x = Métrique, y = Moyenne, fill = Métrique)) +
  geom_col() +
  labs(title = "Random Forest (100 itérations)",
       x = "", y = "Moyenne") +
     theme_minimal()

recapRF <- recapRF + theme(axis.text.x = element_text(angle = 45, hjust = 1))

recapRF %>% plotly::ggplotly()

Notre modèle a obtenu de très bonnes métriques d’évaluations, avec 0.80 pour l’AUC, 0.83 pour la précision, 0.71 pour le rappel et 0.90 pour la spécificité, cela montre quele modèle a une capacité élevée à discriminer entre les classes positives et négatives, tout en maintenant un bon équilibre entre la précision et le rappel, et en minimisant à la fois les faux positifs et les faux négatifs.

Optimisation du modèle

On peut essayer d’optimiser notre modèle pour qu’il ait d’encore meilleur résultat :

forest_global <- randomForest(CDR ~ ., data = dff, importance = TRUE, ntree = 5000)


plot(forest_global$err.rate[, 1], type = "l", xlab = "nombre d'arbres", ylab = "erreur OOB")

Ici, on voit le niveau d’erreur en ordonnée, pour le nombre d’arbres donnés en abscisse. On peut remarquer que le niveau d’erreur se stabilise autour de 1000 arbres. Nous allons donc choisir ce chiffre pour notre modèle.

ntree <- 1000

# selection du mtry optimal

tuned_rf <- tuneRF(dff[, -which(names(dff) == "CDR")], dff$CDR, ntree = ntree)

mtry = 2  OOB error = 16.1% 
Searching left ...
mtry = 1    OOB error = 22.6% 
-0.4035088 0.05 
Searching right ...
mtry = 4    OOB error = 14.97% 
0.07017544 0.05 
mtry = 8    OOB error = 17.8% 
-0.1886792 0.05 

plot(tuned_rf)

Le paramètre “mtry” contrôle le nombre de variables à considérer à chaque division d’un nœud de l’arbre de décision. Sur ce graphique, on a en ordonnée le niveau d’erreur et en abscisse le mtry. On voit que le mtry optimal est de 4, car c’est celui qui a le niveau d’erreur le plus bas.

On peut donc refaire notre modèle en prenant en compte ces deux paramètres.

Importance des variables

# on peut maintenant spécifier nos paramétres
ntree <- 1000
mtry <- 4

# on crée le modèle avec les paramètres spécifiés
forest_global <- randomForest(CDR ~ ., data = dff, importance = TRUE, ntree = ntree, mtry = mtry)

data <- data.frame(
  Variable = c("Sexe", "Âge", "Éducation", "SES", "MMSE", "Volume_intercranien", "Volume_cérébral_total", "Indice atrophie"),
  MeanDecreaseAccuracy = c(0.04451271, 0.03256549, 0.03458807, 0.03844448, 0.08953097, 0.04785418, 0.07876710, 0.05066611),
  MeanDecreaseGini = c(8.341651, 22.764785, 12.559660, 12.870613, 30.298598, 24.410822, 35.275160, 25.007553)
)

# Graphique pour MeanDecreaseAccuracy
plot_accuracy <- ggplot(data, aes(x = reorder(Variable, MeanDecreaseAccuracy), y = MeanDecreaseAccuracy)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  labs(x = "Variable", y = "Mean Decrease Accuracy", title = "Mean Decrease Accuracy") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Graphique pour MeanDecreaseGini
plot_gini <- ggplot(data, aes(x = reorder(Variable, MeanDecreaseGini), y = MeanDecreaseGini)) +
  geom_bar(stat = "identity", fill = "salmon") +
  labs(x = "Variable", y = "Gini", title = "Gini") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

grid.arrange(plot_accuracy, plot_gini, ncol = 2)

En Random Forest, “Mean Decrease Accuracy” (MDA) évalue comment l’exactitude du modèle diminue lorsque les valeurs d’une variable sont permutées, mesurant ainsi son importance pour la prédiction.

En revanche, “Gini” mesure comment une variable réduit l’impureté des nœuds lorsqu’elle est utilisée pour diviser les données, évaluant ainsi son efficacité à séparer les classes. Bien que les deux mesures évaluent l’importance des variables, elles le font à travers des perspectives différentes, elles sont complémentaires et permettent de solidifier les résultats.

Ici, on peut voir que, pour la MDA, l’âge est la variable la moins importante de notre modèle, on pourrait faire deux groupes, avec :

• Age, Education, SES, Sexe, Volume intercrânien, Indice d’atrophie, qui ont chacun eu des résultats assez homogènes.
• Volume cérébral total et MMSE, qui survole plus notre modèle.

Pour Gini, on a plutôt trois groupes, avec :

• Sexe, Education, SES
• Age, Volume intercrânien et indice d’atrophie
• MMSE et Volume cérébral total, qui sont une nouvelle fois au dessus en terme d’importance

Maintenant que nous avons une idée de l’importance de nos variables, nous pouvons faire une régression logistique, qui pourrait permettre une analyse plus détaillée.

IV. Régression logistique binaire

Justification

La régression logistique binaire est une extension du modèle linéaire adaptée aux cas où la variable dépendante est binaire. Elle modélise la relation entre les variables indépendantes et la probabilité d’occurrence de l’événement binaire, permettant d’estimer les probabilités de succès pour chaque niveau des variables indépendantes. En utilisant cette méthode, on peut calculer les p-values pour évaluer la significativité des effets des variables indépendantes sur la variable dépendante.

Nous avons adopté une approche semblable aux méthodes précédentes, avec une validation croisée effectuée en divisant les ensembles en test (20%) et en apprentissage (80%). Des modèles de type dummy ont servi de référence. Les performances ont été évaluées en termes de précision, de rappel, de spécificité et d’AUC. Une boucle répétitive, réitérée 100 fois, a été instaurée pour garantir des résultats robustes.

demence2$CDR_binaire <- ifelse(demence2$CDR == "Dément", 1, 0)
  
  demence2 <- subset(demence2, select = -CDR)
  demence2$CDR_binaire <- as.factor(demence2$CDR_binaire)

  nrun <- 1000
  res_lr <- c()
  res_precision <- c()
  res_recall <- c()
  res_specificity <- c()
  res_auc <- c()
  
  for (irun in 1:nrun) {
    dfflr <- demence2[sample(nrow(demence2)), ]
    trainlog <- dfflr[1:280, ]
    testlog <- dfflr[281:354, ]
    colnames(demence2)
    
      suppressWarnings({
    lr <- glm(CDR_binaire ~ sexe + Age + MMSE + Indice_atrophie + Volume_cérébral_total + 
                Indice_atrophie + EDUC + SES, data = trainlog, family = binomial())
  })
      
    
    dummylr <- rep(names(table(trainlog$CDR_binaire))[which.max(table(trainlog$CDR_binaire))], nrow(testlog))
    predictionslr <- predict(lr, testlog, type = "response")
    predictionslr <- ifelse(predictionslr > 0.5, 1, 0)
    
    res_lr <- rbind(res_lr, c(
      mean(dummylr == testlog$CDR_binaire),
      mean(predictionslr == testlog$CDR_binaire)
    ))
    
    # calcul des métriques
    predictions_lr <- predict(lr, testlog, type = "response")
    predictions_lr <- ifelse(predictions_lr > 0.5, 1, 0)
    predictions_lr <- factor(predictions_lr)
    
    confusion_matrix <- confusionMatrix(predictions_lr, testlog$CDR_binaire)
    
    precision <- confusion_matrix$byClass["Pos Pred Value"]
    recall <- confusion_matrix$byClass["Sensitivity"]
    specificity <- confusion_matrix$byClass["Specificity"]
    
    # calcul du ROC
    roc_curve_lr <- roc(testlog$CDR_binaire, as.numeric(predictions_lr))
    auc_value_lr <- auc(roc_curve_lr)
    
    # on store les résultat
    res_precision <- c(res_precision, precision)
    res_recall <- c(res_recall, recall)
    res_specificity <- c(res_specificity, specificity)
    res_auc <- c(res_auc, auc_value_lr)
  }
  
  
  colnames(res_lr) <- c("dummy", "lr")

  
  means_lr <- colMeans(res_lr)
  
  df_lr <- data.frame(Model = names(means_lr), Accuracy = means_lr)

Robustesse du modèle

   # on crée le ggplot pour la comparaison entre Dummy et Logistic Regression
 ggLR = ggplot(df_lr, aes(x = Model, y = Accuracy, fill = Model)) +
    geom_col() +
    labs(title = "Comparaison entre Dummy et LR",
         x = "Modèle", y = "Précision") +
    theme_minimal()
  ggLR

Le dummy a de moins bons résultat que le modèle, on peut donc continuer notre analyse !

 # on prend les moyenne des métriques
  mean_precision <- mean(res_precision)
  mean_recall <- mean(res_recall)
  mean_specificity <- mean(res_specificity)
  mean_auc <- mean(res_auc)
  
  # data frame des métrique
  df_metrics <- data.frame(Métrique = c("Precision", "Rappel", "Specificité", "AUC"),
                           Mean_Value = c(mean_precision, mean_recall, mean_specificity, mean_auc))
  
  # graph des métriques
  RecapRLB = ggplot(df_metrics, aes(x = Métrique, y = Mean_Value, fill = Métrique)) +
    geom_bar(stat = "identity") +
    labs(title = "Régression Logistique (100 itérations)",
         x = "", y = "Moyenne") +
        theme_minimal()
  
RecapRLB <- RecapRLB + theme(axis.text.x = element_text(angle = 45, hjust = 1))
RecapRLB%>% plotly::ggplotly()

Le modèle a de très bons résultats pour l’AUC, la précision, et le rappel, avec respectivement 0.72, 0.74 et 0.87. Toutefois, la sensbilité est plus basse, avec 0.58, cela indique que le modèle pourrait manquer certaines observations positives.

{
  lr_global <- glm(CDR_binaire ~ sexe + Age + MMSE + EDUC + 
                   Volume_intercranien + Indice_atrophie + SES + Volume_cérébral_total, 
                   data = dfflr, maxit = 1000, family = binomial())
}

Avis : glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

Ce message d’avis apparait car R nous prévient que des variables sont certainement multicolinéaires (la multicolinéarité des variables rend ininterprétable les résultats en régression logistique), nous prendrons cela en compte après avoir affiner notre modèle.

suppressWarnings({
  Anova(lr_global)
})

Analysis of Deviance Table (Type II tests)

Response: CDR_binaire
                      LR Chisq Df Pr(>Chisq)    
sexe                    21.662  1  3.252e-06 ***
Age                     10.340  1   0.001301 ** 
MMSE                    54.956  2  1.166e-12 ***
EDUC                    11.614  4   0.020466 *  
Volume_intercranien      0.018  1   0.892697    
Indice_atrophie          0.272  1   0.602225    
SES                     12.709  4   0.012787 *  
Volume_cérébral_total   21.446  1  3.639e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sélection de variables par étapes

Pour assurer l’interprétabilité des résultats de la régression logistique, il est impératif que toutes les variables de notre modèle soient significatives (ici, nous prendrons un seuil alpha égale à 5%). Ainsi, pour atteindre cet objectif, on retire une à une les variables non significatives (car le retrait d’une variable peut modifier le niveau de significativité des autres, impactant ainsi l’analyse globale du modèle), en commençant par éliminer celle qui présente la p-value la moins significative. Ce processus garantit que seules les variables les plus pertinentes demeurent dans le modèle final.

Pour commencer, on retire donc le volume intercranien :

lr_global <- glm(CDR_binaire ~ sexe + Age + MMSE + EDUC + 
                    Indice_atrophie + SES + Volume_cérébral_total, 
                   data = dfflr, maxit = 1000, family = binomial())

Avis : glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

{
  lr_global <- glm(CDR_binaire ~ sexe + Age + MMSE + EDUC + 
                    Indice_atrophie + SES + Volume_cérébral_total, 
                   data = dfflr, maxit = 1000, family = binomial())
}

Avis : glm.fit: des probabilités ont été ajustées numériquement à 0 ou 1

suppressWarnings({
  Anova(lr_global)
})

Analysis of Deviance Table (Type II tests)

Response: CDR_binaire
                      LR Chisq Df Pr(>Chisq)    
sexe                    21.794  1  3.036e-06 ***
Age                     10.570  1   0.001149 ** 
MMSE                    55.413  2  9.273e-13 ***
EDUC                    11.769  4   0.019153 *  
Indice_atrophie         10.519  1   0.001182 ** 
SES                     12.838  4   0.012097 *  
Volume_cérébral_total   21.791  1  3.040e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

On peut voir que les p-values de nos variables ont évolué, toutes sont devenues significatives, ce qui est une bonne nouvelle. Toutefois, on ne peut toujours pas interpréter notre modèle de régression logistique. Avant cela, il est nécessaire de prendre en compte la possible multicollinéarité des variables.

Prise en compte du VIF

Le VIF (Variance Inflation Factor) en régression est une mesure de l’ampleur à laquelle la variance d’un coefficient de régression est augmentée en raison de la corrélation entre les prédicteurs dans le modèle. Un VIF élevé suggère une forte multicollinéarité entre les variables explicatives, ce qui peut rendre les estimations de coefficients peu fiables. On peut ainsi regarder les VIFs de nos variables :

vif(lr_global)

                              GVIF Df GVIF^(1/(2*Df))
sexe                  2.019870e+00  1        1.421221
Age                   1.812819e+00  1        1.346410
MMSE                  1.000000e+00  2        1.000000
EDUC                  1.067767e+07  4        7.560657
Indice_atrophie       2.164453e+00  1        1.471208
SES                   1.181484e+07  4        7.656909
Volume_cérébral_total 2.015482e+00  1        1.419677

On peut créer un GGplot pour que ce soit plus parlant :

# créer un ggplot pour mieux visualiser
data <- data.frame(
  Variable = c("sexe", "Age", "MMSE", "EDUC", "Volume_cérébral_t", "Indice_atrophie", "SES"),
  GVIF = c(1.421221, 1.346410, 1.000000, 7.560657, 1.419677, 1.471208, 7.656909)
)

# Graphique
ggvif = ggplot(data, aes(x = Variable, y = GVIF, fill = Variable)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_hline(yintercept = 5, color = "red", linetype = "solid") + # Ligne rouge à y = 5
  labs(title = "GVIF pour chaque variable",
       x = "Variable",
       y = "GVIF Raised") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggvif %>% plotly::ggplotly()

Dans la littérature, il est couramment admis que le VIF peut poser problème lorsque sa valeur dépasse 5. Cette limite est représentée par la ligne rouge sur le graphique. Ici, on peut voir que le niveau d’éducation et le niveau socio-économique posent problème, nous allons retirer la variable qui a le VIF le plus élevé, à savoir le niveau socio-économique, qui est de 7.65.

lr_global2 <- glm(CDR_binaire ~ sexe + Age + MMSE + EDUC + 
                   Volume_cérébral_total + Indice_atrophie, 
                   data = dfflr, maxit = 1000, family = binomial())

vif(lr_global2)

                          GVIF Df GVIF^(1/(2*Df))
sexe                  1.754377  1        1.324529
Age                   1.759665  1        1.326524
MMSE                  1.000000  2        1.000000
EDUC                  1.201903  4        1.023255
Volume_cérébral_total 1.902207  1        1.379205
Indice_atrophie       1.967118  1        1.402540

# on va créer un ggplot pour mieux visualiser
data <- data.frame(
  Variable = c("sexe", "Indice_atrophie", "MMSE", "Volume_substance_blanche", "Indice_atrophie", "EDUC"),
  GVIF = c(1.324529, 1.402540, 1.000000, 1.379205, 1.419677, 1.023255)
)

# Graphique
ggvif2 = ggplot(data, aes(x = Variable, y = GVIF, fill = Variable)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_hline(yintercept = 5, color = "red", linetype = "solid") + # Ligne rouge à y = 5
  labs(title = "GVIF pour chaque variable",
       x = "Variable",
       y = "Niveau") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggvif2 %>% plotly::ggplotly()

Sans cette variable, tous les VIFS sont inférieurs à 1,5 et ne posent donc plus de problème. On peut donc continuer l’analyse et recalculer nos p-values.

Modèle final

Anova(lr_global2)

Analysis of Deviance Table (Type II tests)

Response: CDR_binaire
                      LR Chisq Df Pr(>Chisq)    
sexe                    22.386  1  2.230e-06 ***
Age                      7.946  1   0.004820 ** 
MMSE                    50.322  2  1.182e-11 ***
EDUC                    12.490  4   0.014058 *  
Volume_cérébral_total   20.794  1  5.115e-06 ***
Indice_atrophie         10.406  1   0.001256 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Maintenant que nous avons vérifié toutes nos conditions, nous pouvons interpréter notre modèle final. Nous avons donc retiré le volume intercranien ainsi que le niveau socio-économique, nous laissant :

“sexe”, “MMSE”, “Volume_cérébral_total”, qui ont des p-values très faibles (inférieures à 0,001), ce qui est indiqué par trois étoiles (***), suggérant une forte association avec la variable dépendante.

Les variables “Age” et “Indice_atrophie” ont des p-valeurs de 0.004820 et 0.001256 respectivement, ce qui est indiqué par deux étoiles, suggérant une association significative, mais avec un niveau de confiance légèrement moins élevé que celles avec trois étoiles.

Ce qui nous laisse la variable “EDUC”, qui, bien que présentant une p-valeur légèrement plus élevée que les autres variables dans notre analyse (0,014058), reste statistiquement significative au niveau de confiance de 95%.

Vérification des conditions d’utilisation

Dans la régression logistique, un modèle linéaire généralisé, l’adéquation est évaluée en examinant les résidus de Pearson et de déviance. Les résidus de Pearson mesurent la différence entre les valeurs observées et attendues, tandis que les résidus de déviance évaluent l’ajustement du modèle aux données. Ces résidus sont souvent testés à l’aide de distributions du chi-deux pour détecter les points atypiques ou mal ajustés dans les données.

binnedplot(fitted(lr_global2), residuals(lr_global2, type = "response"))

Le graphique produit par binnedplot indique un intervalle de prédiction à 95% (lignes grises) pour les résidus moyens. Ainsi, si le modèle binomial est bon, environ 95% des résidus devraient se trouver à l’intérieur de cet intervalle. Ici, c’est le cas de 13 résidus sur 16, soit 81.25%, ce qui est moins performant de ce qu’on attendrait dans l’idéal !

L’estimation des paramètres des modèles linéaires généralisés n’utilise pas la méthode des moindres carrés, ce qui fait que le R², couramment utilisé en modèle linaire générale pour juger du degré d’efficacité du modèle, n’a pas de sens ici.

Pour pallier cela, on peut utiliser un pseudo R², plusieurs existe, on a par exemple le R² de McFadden, qui se calcule tel que la fraction de la déviance du modèle nul expliquée par le modèle incluant les prédicteurs.

pseudo_R2 <- 1 - lr_global2$deviance/lr_global2$null.deviance
pseudo_R2

[1] 0.3249663

On obtient un R² ajusté de 0.32. Le R² ajusté va de 0 à 1, 0 veut dire que le modèle ne parvient pas à expliquer quoi que ce soit en fonction de la variable dépendante, et à 1 qu’il y arrive parfaitement, ce qui indique un ajustement modéré du modèle aux données.

Pour être sûr, on peut utiliser une autre méthode qui calcule le R² ajusté, quantifie la différence entre la valeur moyenne prédite pour les observations où la variable dépendante est égale à 1 et celle pour les observations où la variable dépendante est égale à 0. Cela permet d’évaluer la capacité du modèle à différencier les deux catégories de réponse. ainsi, il mesure la force de prédiction du modèle en termes de différenciation entre les deux niveaux de la variable dépendante.

r2_tjur <- mean(lr_global2$fitted.values[lr_global2$y == 1]) - mean(lr_global2$fitted.values[lr_global2$y == 0])
r2_tjur

[1] 0.36894

Avec cette méthode, on obtient 0.36, ce qui est assez semblable au R² de McFadden calculé précedemment et consolide ce résultat.

En conclusion, notre analyse révèle des associations significatives entre le sexe, l’âge, le MMSE, le crébral total, l’indice d’atrophie et le niveau d’éducation avec la variable dépendante, suggérant leur importance dans notre modèle.

La vérification des conditions d’applications montre un résultat en dessous de ce qu’on aurait pu attendre pour Les résidus de pearson, et un R² ajusté montrant une efficacité modéré, les résultats de notre régression logistique binomiale, bien que très intéressants, sont donc à prendre avec des pincettes.

V. Analyse discriminante linéaire (LDA)

Justification

L’analyse discriminante linéaire offre une méthode efficace pour simplifier les modèles, améliorer la précision de la classification et fournir des résultats interprétables. Pour effectuer une LDA, il est nécessaire que toutes les variables soient quantitatives, et il ne doit pas y avoir de colinéarité entre les variables.

Comme nous avons des variables qualitatives, nous allons d’abord faire une Analyse Factorielle des Correspondances (AFC).

AFC

Nous avons vu lors de la régression logistique qu’il y avait une multicolinéarité entre la variable représentant le niveau d’éducation (EDUC) et la variable représentant le niveau socio-culturel (SES). Comme précedemment, nous allons donc retirer la variable qui avait le VIF le plus haut, c’est à dire le niveau socio-culturel.

Cela dit, on peut faire notre AFC :

 data_afc <- demence2[, c("sexe", "MMSE", "EDUC")]
  
  afc_result <- MCA(data_afc, graph = FALSE)
 
  head(afc_result$ind$coord)

       Dim 1      Dim 2       Dim 3      Dim 4       Dim 5
1  0.2145599  0.5027782 -0.07376176 -0.6477649 -0.19594701
2  0.2145599  0.5027782 -0.07376176 -0.6477649 -0.19594701
6 -0.8155078 -0.1621775  0.50210218 -0.1128201  0.62699299
7 -0.8155078 -0.1621775  0.50210218 -0.1128201  0.62699299
8  0.1340800  0.2572015 -0.16953925  0.6477649  0.08310553
9  0.1340800  0.2572015 -0.16953925  0.6477649  0.08310553

On peut voir qu’on a cinq dimensions dans notre AFC, ceux-ci représentent nos variables qualitatives, on va donc les inclure dans notre jeu de données de base, puis retirer les variables qualitatives afin de pouvoir faire notre LDA.

 lda_input <- as.data.frame(afc_result$ind$coord[, paste0("Dim ", 1:5)]) 
  demence2 <- cbind(demence2, lda_input)
  demence2 <- demence2[, -c(1, 3, 4, 5)]

Préparation des données pour la LDA

Cela fait, on peut préparer nos données pour faire la LDA :

  gp = as.factor(demence2$CDR)
  matrice = as.matrix(demence2)
  data = data.frame(group=demence2$CDR, matrice)
  colnames(data)

 [1] "group"                 "Age"                   "Volume_intercranien"  
 [4] "Volume_cérébral_total" "Indice_atrophie"       "CDR_binaire"          
 [7] "Dim.1"                 "Dim.2"                 "Dim.3"                
[10] "Dim.4"                 "Dim.5"                

  data$Age <- as.numeric(data$Age)
  
  data$Volume_intercranien<- as.numeric(data$Volume_intercranien)
  data$Volume_cérébral_total <- as.numeric(data$Volume_cérébral_total)
  data$Indice_atrophie <- as.numeric(data$Indice_atrophie)
  data$Dim.1 <- as.numeric(data$Dim.1)
  data$Dim.2 <- as.numeric(data$Dim.2)
  data$Dim.3 <- as.numeric(data$Dim.3)
  data$Dim.4 <- as.numeric(data$Dim.4)
  data$Dim.5 <- as.numeric(data$Dim.5)
  
  colonnes_numeriques <- sapply(data, is.numeric)
  DataCR<- scale(data[,colonnes_numeriques])

  data2 <- as.data.frame(DataCR)

  data2$CDR <- demence2$CDR_binaire

  data2$CDR <- as.factor(data2$CDR)


  gp = as.factor(data2$CDR)
  Datalda = data.frame(group=gp,data2)

Mise en place de la LDA

Comme pour les méthodes précédentes, on fait de la validation croisée, avec 80 % pour l’entraînement et 20 % pour le test, puis on compare notre modèle à un modèle dummy et on calcule les métriques d’évaluation sur l’échantillon test, avec une boucle de 100 répétitions.

  # créer une liste vide pour stocker les performances
  performances_metrics <- list()
  
    # boucle sur les itérations
  set.seed(100)
  nrun <- 100
  for (irun in 1:nrun){
    df<- data2[sample(nrow(data2)), ]
    
    
    train <- df[1:280, ]
    test <- df[281:354, ]
    reslda <- lda(CDR ~ . , data=train)
    predlda <- predict(reslda,test)
    # matrice de confusion
    confusion_matrix <- table(test$CDR, predlda$class)
    # métriques pour chaque itération de validation croisée
    confusion_matrices <- confusionMatrix(test$CDR, predlda$class)
    
    # métriques de performance
    
    Sensibilite <- confusion_matrices$byClass['Sensitivity']
    Specificite <- confusion_matrices$byClass['Specificity']
    Precision <- confusion_matrices$byClass['Pos Pred Value']
    # calcul l'AUC
    AUC <- roc(test$CDR, as.numeric(predlda$x))$auc
    #stocker les metriques pour chaque iteration
    performances_metrics[[irun]] = c(Sensibilite, Specificite, Precision, AUC)
  }
  metric_noms=c("Sensibilite", "Specificite", "Precision", "AUC")
  #calculer la moyenne des metriques pour toutes les iterations
  avr_metrics=colMeans(do.call(rbind,performances_metrics))
  metric_summary <- data.frame(Metrique = metric_noms,Moyenne = avr_metrics)
 
  ##créer un ggplot
  recapLDA = ggplot(metric_summary, aes(x = Metrique, y = Moyenne, fill = Metrique)) +
    geom_col() +
    labs(title = "LDA (100 itérations)",
         x = "", y = "Moyenne") +
    theme_minimal()
  recapLDA %>% plotly::ggplotly()

  recapLDA <- recapLDA + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
  
  
    recapLDA %>% plotly::ggplotly()

Notre LDA démontre des performances solides, avec une AUC de 0,80 indiquant une bonne capacité de discrimination entre les classes positives et négatives. De plus, une précision élevée de 0,90 suggère que la majorité des prédictions positives sont correctes. Bien que la sensibilité soit légèrement plus basse à 0,70, elle indique tout de même une capacité à capturer une proportion significative des observations positives. Enfin, une spécificité de 0,78 souligne la capacité du modèle à identifier correctement les observations négatives.

reslda2 <- lda(CDR ~ . , data=df)
predlda2 <- predict(reslda,df)

boxplot(predlda2$x~gp)

 # extraire les coefficients du modèle LDA
  coefficients <- coef(reslda2)
  
  # Créer un dataframe à partir des coefficients
  coefficients_df <- data.frame(
    Variable = rownames(coefficients),
    Coefficient = coefficients[, 1]  # Extraire uniquement la première colonne de coefficients
  )
  
  # Créer le graphique avec les coefficients extraits
  ggplot(coefficients_df, aes(x = Variable, y = Coefficient)) +
    geom_bar(stat = "identity", fill = "skyblue", width = 0.5) +
    labs(title = "Coefficients de la dimension discriminante (LD1)",
         x = "", y = "Coefficient") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

V. Conclusion

grob_recapCART <- ggplotGrob(recapCART)
grob_recapRF <- ggplotGrob(recapRF)
grob_RecapRLB <- ggplotGrob(RecapRLB)
grob_recapLDA <- ggplotGrob(recapLDA)

# VI. Bibliographie