diff --git a/M2/Data Visualisation/Exemple Projet/Application projet.Rmd b/M2/Data Visualisation/Exemple Projet/Application projet.Rmd new file mode 100644 index 0000000..6d67bf4 --- /dev/null +++ b/M2/Data Visualisation/Exemple Projet/Application projet.Rmd @@ -0,0 +1,1345 @@ +--- +title: "Projet de visualisation des données" +author: "Victoire de Salaberry" +date: "Années 2023-2024" +institute: "Université Paris-Dauphine" +departement: "Master 2 ISF" +output: + rmdformats::robobook: + highlight: kate + use_bookdown: true + css: style.css + lightbox : true + gallery: true + code_folding: hide +editor_options: + markdown: + wrap: 72 +--- + +```{r, include=FALSE, echo = FALSE} +# Installation des librairies nécesssaires + +# Pour les graphiques +library(ggplot2) # créer des graphiques +library(gridExtra) # organiser et agencer plusieurs graphiques +library(tidyr) # remodeler ou restructurer les données +library(plotly) # créer des graphiques interactifs +library(gmodels) # créer des graphiques spécifiques à la modélisation statistique + +# Pour les tableaux +library(knitr) +library(kableExtra) # créer des tableaux + +# Pour les corrélations +library(GGally) # extension à ggplot2 +library(corrplot) # visualiser la matrice de corrélations +library(magrittr) # pour les pipes (%>%) + +# Pour les modèles +library(glm2) # créer des modèles linéaires généralisés +library(caret) # ajuster des modèles et évaluer les performances +library(broom) # convertir les résultats de modèles statistiques en cadres de données "tidy" +library(randomForest) # créer des modèles Random Forest +library(rpart) # créer des arbres de décisions + +# pour les performances +library(pROC) # pour la courbe ROC +library(ROCR) # pour la courbe ROC +library(plotROC) # pour la courbe ROC +library(MLmetrics) # pour des mesures d'évaluation de la performance des modèles +library(Metrics) # pour des mesures d'évaluation de la performance des modèles +library(effects) # visualiser et interpréter les effets des variables dans des modèles statistiques +library(broom.helpers) # manipuler des résultats de modèles statistique +library(car) # effectuer des tests Anova +``` + +# Introduction + +Les maladies cardiovasculaires sont une cause majeure de mortalité dans le monde. La détection et l'intervention précoces sont essentielles pour prévenir celles-ci. Certains facteurs sont bien connus comme étant des symptômes de maladies cardiaques. + +Nous allons analyser un jeu de données pour confirmer ces symptômes, c'est-à-dire vérifier s'ils ont un réel impact sur la santé cardiaque d'un individu. Puis, à partir de ce même jeu de données, nous essaierons de prédire au mieux la variable cible grâce à des modèles statistiques dont nous étudierons les performances. + +La variable cible est la variable nommée *target*. Elle représente l'état du patient. Cette variable vaut 0 si l'individu est en bonne santé et vaut 1 si l'individu souffre d'une maladie cardiaque. + +# Observation et préparation des données + +## Observation des données + +Dans un premier temps, on importe les données et on affiche un aperçu de celles-ci pour les observer rapidement. +```{r} +# Téléchargement du jeu de données +df <- read.csv("heart.csv") + +# Affichage du jeu de données sous forme de tableau +rmarkdown::paged_table(df) +``` + + +On affiche aussi la dimension du jeu de données : +```{r} +# Création d'un dataframe contenant les dimensions +dim <- data.frame(dim(df)) +rownames(dim) <- c("Nombre de lignes","Nombre de colonnes") # Nom des lignes du tableau +colnames(dim) <- c("") # Nom de la colonne du tableau + +# Affichage des dimensions du jeu de données dans un tableau +dim %>% + kbl(caption = "Dimension du jeu de données") %>% + kable_styling() +``` + +Le jeu de données comporte 303 lignes, c'est-à-dire 303 observations (individus) et 14 colonnes, c'est-à-dire 14 variables. 13 d'entre elles sont les variables explicatives et la dernière est la variable cible (la variable *target* comme vu en introduction). + +On regarde ensuite le type des variables, pour savoir si elles sont quantitatives ou qualitatives ; numériques ou catégorielles, bien qu'on en ait déjà une idée grâce à l'affichage précédent. +Le tableau suivant présente une description récapitulative des variables. On pourra trouver une description plus détaillée de celles-ci dans la notice du projet. + +```{r, include = FALSE} +# Affichage des informations sur le type des données +str(df) +``` + +```{r} +# Affichage de la définition des variables dans un tableau +kable(caption = "Définition des variables", # Titre du tableau + # Création d'un dataframe avec les variables, leur description, leur type et leur support + data.frame( + Variable = c( + "age", "sex", "cp", "trestbps", "chol", "fbs", "restecg", "thalach", + "exang", "oldpeak", "slope", "ca", "thal"), + Description = c( + "Age de l'individu", + "Sexe de l'individu", + "Douleur thoracique", + "Pression artérielle au repos en mmHg", + "Niveau de cholestérol sérique du patient en mg/dL", + "Taux de glycémie à jeun", + "Résultats électrocardiographiques au repos", + "Fréquence cardiaque maximale atteinte", + "Angine de poitrine induite par l'exercice", + "Dépression du segment ST induite par l'exercice par rapport au repos", + "Pente du segment ST de pointe de l'exercice", + "Nombre de vaisseaux principaux colorés par la flourosopie", + "Flux sanguin vers le muscle cardiaque "), + Type = c( + rep("Entier", 9), + "Réel", + rep("Entier", 3)), + Support = c( + "⟦29;77⟧", "{0, 1}", "{0, 1, 2, 3}", "⟦94;200⟧", "⟦126; 564⟧", "{0, 1}", + "{0, 1, 2}", "⟦71;202⟧", "{0, 1}", "[0, 6.2]", "{0, 1, 2}", "⟦0;4⟧", "{0, 1, 2, 3}")), + booktabs=TRUE) # Mise en page et style du tableau +``` +Les variables sont toutes quantitatives et de type entier sauf *oldpeak* qui est un réel. +On a des variables + +* catégorielles : *sex, cp, fbs, restecg, exang, slope, ca* et *thal* ; +* numériques : *age, trestbps, chol, thalach* et *oldpeak*. + +**Remarque** : les variables *sex, cp, fbs, restecg, exang, slope, ca* et *thal* sont des variables numériques qu'on a classifié de catégorielles. Ces variables sont en fait des variables catégorielles, dont les catégories ont été présentées dans la notice, qui ont été "transformées" en variables numériques pour faciliter leur utilisation. + +## Valeurs manquantes + +On vérifie s'il y a des valeurs manquantes. Si c'est le cas, il faudra y remédier. +Le tableau ci-dessous affiche le nombre de valeur manquante par variable. +```{r} +# Vérification des valeurs manquantes dans le dataframe +val_manquantes <- sapply(df, function(x) sum(is.na(x))) + +# Affichage du nombre de valeurs manquantes par variable dans un tableau +val_manq <- data.frame(val_manquantes) # Création d'un dataframe avec le nombre de valeurs manquantes +colnames(val_manq) <- c("Nombre de valeur manquante") # Nom de la colonne du tableau +val_manq %>% + kbl(caption = "Nombre de valeur manquante par variable") %>% # Titre + kable_styling() # Mise en page et style du tableau +``` +Il n'y a donc pas de valeurs manquantes dans notre jeu de données. + +## Distribution des variables + +Après avoir bien observé le jeu de données, passons à l'étude graphique de la distribution des variables. + +### Distribution de la variable cible + +On étudie dans un premier temps la distribution de la variable cible puis celle des variables explicatives. + +```{r} +# Création d'un camembert avec les pourcentages +pie(table(df$target), # Table des effectifs pour la variable target + main = "Distribution de la variable cible (target)", # Titre du graphique + labels = c("Individus en bonne santé", + "Individus atteints \nd'une maladie cardiaque"), # Légende + col = c("lightblue", "lightcoral"), # Couleur des secteurs + cex = 0.8, # Taille de la légende + cex.main = 1) # Taille du titre + +# Création d'une fonction qui ajoute les valeurs en tant qu'étiquettes au centre (fonction créée par chatgpt) +text_pie = function(vector, labels = c(), cex = 0.8) { + vector = vector / sum(vector) * 2 * pi + temp = c() + j = 0 + l = 0 + for (i in 1:length(vector)) { + k = vector[i] / 2 + j = j + l + k + l = k + text(cos(j) / 2, sin(j) / 2, labels[i], cex = cex) + } + vector = temp +} + +# Ajout des étiquettes +percent_values <- round(100 * prop.table(table(df$target)), 1) # Calcul des pourcentages arrondis +percent <- paste(percent_values, "%") # Ajout du sigle % +text_pie(percent_values, + c(percent[1], percent[2]), + cex = 0.8) +``` +Le jeu de données est à peu près équilibré : on a environ 50% (54.5%) d'individus malades et 50% (45.5%) d'individus sains dans le jeu de données. + +### Distribution des variables explicatives + +On affiche maintenant la distribution des variables catégorielles. Cela va nous permettre d'observer les catégories de variable qui prédominent par rapport aux autres catégories de cette même variable. +```{r} +# Création d'un vecteur avec les variables catégorielles +categorical_var <- c("sex", "cp", "fbs", "restecg", "exang", "slope", "ca", "thal") + +# Organisation des graphiques en une grille 2x4 +par(mfrow = c(2, 4)) + +# Boucle qui crée un histogramme en baton par chaque variable catégorielle +for (var in categorical_var) { + barplot(table(df[[var]]), + col = "lightblue", # Couleur de l'histogramme + main = var, # Titre de l'histogramme + ylab = "Nombre d'individus") # Nom de l'axe des ordonnées +} +``` +Les variables catégorielles ne sont pas équilibrées dans le jeu de données : + +* *sex* : il y a deux fois plus d'hommes que de femmes. +* *cp* : beaucoup d'individus ne souffrent pas de douleurs thoraciques (0). Moitié moins souffrent de douleurs modérées (2). Seuls une vingtaine d'individus sur 300 ont de fortes douleurs (3). +* *fbs* : il y a 6 fois moins d'individus ayant un taux de glycémie anormal (> 120 mg/dL) que de personnes ayant un taux de glycémie normal (< 120 mg/dL). +* *restecg* : il y a presque autant d'individus dont les résultats électrocardiographiques sont normaux (0) que de personnes ayant des résultats avec des changements mineurs (1). Très peu (environ 5%) ont des résultats avec des anomalies importantes (2). +* *exang* : il y a 2 fois moins d'individus qui ont des angines de poitrine induites par l'exercice (1) que d'individus qui en n'ont pas (0). +* *slope* : la plupart des individus ont une dépression du segment ST (environ 85%). +* *ca* : plus de la moitié des individus n'ont pas de vaisseaux principaux colorés par la flourosopie. +* *thal* : la plupart des individus ont des défauts modérés (2) voir graves (3) de perfusion myocardique. + +On regarde maintenant la distribution des variables numériques. On trace un histogramme et le boxplot asssocié pour chaque variable. Le boxplot nous permet d'avoir la moyenne de chaque variable et d'observer les éventuelles valeurs aberrantes. +```{r} +# Création d'un vecteur avec les variables numériques +numerical_var <- c("age", "trestbps", "chol", "thalach", "oldpeak") + +# Création de listes pour les histogrammes et les boxplots +hist_plots <- list() +box_plots <- list() + +# Boucle pour créer les histogrammes et les boxplots pour chaque variable numérique +for (var in numerical_var) { + # Histogramme + histo <- ggplot(df, aes(x = .data[[var]])) + + geom_histogram(fill = "lightblue", # Couleur de remplissage + color = "black", # Couleur des contours + bins = 20, # Nombre de catégories + alpha = 0.7, # Transparence de la couleur de remplissage + linewidth = 0.2) + # Épaisseur des contours + labs(title = paste("Distribution de", var), # Titre des histogrammes + x = var, # Nom de l'axe des abscisses + y = "Nombre") + # Nom de l'axe des ordonnées + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 7, # Taille du titre + face = "bold"), # Titre en gras + axis.title.x = element_text(size = 7), # Taille du nom de l'axe des abscisses + axis.title.y = element_text(size = 7)) # Taille du nom de l'axe des ordonnée + + # Boxplot + boxplot <- ggplot(df, aes(x = 1, y = .data[[var]])) + + geom_boxplot(linewidth = 0.3, # Épaisseur des contours + outlier.size = 0.2) + # Taille des outliers + labs(title = paste("Boxplot de", var)) + # Titre des boxplots + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 7, # Taille du titre + face = "bold"), # Titre en gras + axis.title.x = element_text(size = 7), # Taille du nom de l'axe des abscisses + axis.title.y = element_text(size = 7)) # Taille du nom de l'axe des ordonnée + + # Stockage des histogrammes et des boxplots dans les listes + hist_plots[[var]] <- histo + box_plots[[var]] <- boxplot +} + +# Organisation des graphiques dans une grille 5x2 +grid.arrange(grobs = c(hist_plots, box_plots), + ncol = 5, + top = "Histogrammes et boxplots de chaque variable numérique") # Titre +``` +Grâce à ces graphiques ont peut observer que : + +* *age* : les individus ont entre 29 et 77 ans avec une moyenne d'âge de 54 ans. +* *trestbps* : la pression artérielle au repos de la plupart des individus est supérieure à 120 mmHg, ce qui est un niveau élevé voir à risque. +* *chol* : la plupart des individus ont un niveau de cholestérol sérique entre 200 et 300 mg/dL, ce qui est très élevé. +* *thalach* : la fréquence cardiaque maximale atteinte est entre 100 et 200 bpm, ce qui semble tout à fait normal car celle-ci diminue avec l'âge et dans notre jeu de données, il y a toutes les catégories d'âge. +* *oldpeak* : environ 1/3 des individus ont une dépression ST proche de 0. Les autres individus ont des valeurs plus élevées de dépression ST. + +On peut de plus observer certaines valeurs atypiques pour les variables *trestbps, chol, thalach* et *oldpeak*. Ce sont des valeurs extrêmes mais pas aberrantes (ce ne sont probablement pas des erreurs de saisie) donc on les garde. + +## Corrélations entre les variables + +### Corrélations entre les variables explicatives + +On étudie maintenant s'il y a des corrélations entre les variables explicatives. +On affiche pour cela la matrice des corrélations sous forme de graphique corrplot. + +```{r} +# Plot des corrélations entre les variables explicatives +corrplot(cor(df[, 1:13]), method = "color") +``` +D'après ce graphique, les variables ne semblent pas très corrélées. Nous pouvons cependant porter notre attention sur les 5 couples de variables les plus corrélés. + +Le premier couple de variables est celui composé des variables *oldpeak* et *slope*. Cette corrélation semble tout à fait intuitive car on rappelle que oldpeak représente la dépression ST induite par l'exercice par rapport au repos et que slope est la pente du segment ST de pointe de l'exercice. + +De plus, les variables *age* et *thalach* semblent avoir un coefficient de corrélation plus élevé que les autres variables ; ce qui est cohérent car, comme vu dans la description détaillée des variables (dans la notice), la fréquence cardiaque maximale diminue avec l'âge (on le vérifie plus bas). + +Les variables *exang* et *cp* ont aussi un coefficient de corrélation de couleur foncée. Là encore, cette corrélation n'est pas aberrante car ces variables concernent toutes les 2 des douleurs thoraciques. + +Le couple de variables *slope/thalach* est lui aussi un couple qui semble avoir un coefficient de corrélation plus élevé que la plupart des autres couples de variables. + +Enfin le couple *thalach* et *exang* est le dernier couple qui se détache des autres par la valeur, plutôt élevée, de son coefficient de corrélation. + +On affiche la matrice des corrélations pour voir plus précisément les coefficients de corrélations. +Les coefficients dont la valeur absolue est supérieure à 0.35 sont affichés en rouge. + +```{r} +# Calcul des corrélations entre les variables explicatives +correlation_matrix <- cor(df[, 1:13]) + +# Fonction pour formater les coefficients de corrélation (pour colorer les coefficients dont la valeur absolue est supérieure à 0.35 et arrondir les coefficients au centième près) +format_correlation <- function(x) { + if (abs(x) > 0.35) { + paste0("", round(x, 2), "") + } else { + round(x, 2) + } +} + +# Application de la fonction de formatage à la matrice de corrélation +formatted_correlation_matrix <- apply(correlation_matrix, 1, function(row) { + sapply(row, format_correlation) +}) + +# Affichage de la matrice de corrélation avec la mise en forme personnalisée +kable(formatted_correlation_matrix, + format = "html", + escape = FALSE) %>% + kable_styling() + +``` + +Naturellement, la matrice nous mène à la même conclusion que le graphique : il ne semble pas y avoir de corrélations significatives entre les variables. +Les variables les plus corrélées sont + +* *oldpeak* et *slope* avec un coefficient de corrélation de -0.58, +* *age* et *thalach* avec un coefficient de corrélation de -0.4, +* *exang* et *cp* avec un coefficient de corrélation de -0.39, +* *slope* et *thalach* avec un coefficient de corrélation de -0.39, +* *thalach* et *exang* avec un coefficient de corrélation de -0.38. + +On étudie graphiquement la relation entre les variables des couples les plus corrélés, pour avoir un représentation plus visuelle de celle-ci. + + +On trace ci-dessous la distribution de la la dépression ST induite par l'exercice (*oldpeak*) pour chaque niveau de la variable *slope*. + +```{r} +# Distribution de oldpeak en fonction de slope +ggplot(df, aes(x = oldpeak, fill = factor(slope))) + + geom_density(alpha = 0.5, # Transparence de la couleur de remplissage + linewidth = 0.3) + # Épaissseur des contours + + # Noms des axes et titre + labs(title = "Graphique de la distribution de la la dépression ST induite par l'exercice (oldpeak) \n en fonction de la pente du segment ST de pointe de l'exercice (slope)", # Titre + x = "Changements dans l'électrocardiogramme (oldpeak)", # Nom de l'axe des abscisses + y = "Nombre de patients", # Nom de l'axe des ordonnées + fill = "Pente du segment \nST de pointe (slope)") + # Titre de la légende + + # Couleur et texte de la légende + scale_fill_manual(values = c("0" = "#FFDB6D", + "1" = "#4E84C4", + "2" = "#C3D7A4")) + + + # Mise en forme du titre et de la légende + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 10, # Taille du titre + face = "bold"), # Titre en gras + legend.title = element_text(size = 9), # Taille du titre de la légende + legend.text = element_text(size = 8)) # Taille du texte de la légende +``` + +On remarque que lorsque la dépression ST induite par l'exercice (*oldpeak*) est proche de 0, la pente du segment ST de pointe de l'exercice (*slope*) a tendance à être de 2. Inversement, il n'y a pas de dépression du segment ST observée par rapport à la ligne de base lorsqu'il y a beaucoup de changements dans l'électrocardiogramme (ECG) qui se produisent en réponse à l'exercice physique. + +On confirme donc qu'à priori, il faut enlever une des 2 variables. Dans la suite, on va faire d'autres tests pour confirmer cette décision. + +On regarde ensuite plus précisement la relation entre les variables *age* et *thalach*. +```{r} +# Diagramme de dispersion entre age et thalach +ggplot(df, aes(x = age, y = thalach)) + + + # Ajout de la courbe de tendance + geom_smooth(method = "lm", + formula = y ~ x, # Lien linéaire (polynôme de degré 1) + se = FALSE, # Ne pas afficher l'intervalle de confiance + colour="black", # Couleur de la droite + linewidth = 0.3) + # Épaisseur du trait + + # Ajout des points + geom_point(size = 1, # Taille des points + color = "#4E84C4") + # Couleur des points + + # Nom des axes et titre + labs(title = "Graphique de la fréquence cardiaque maximale atteinte (thalach) \n en fonction de l'âge (age)", # Titre du graphique + x = "Âge (age)", # Nom de l'axe des abscisses + y = "Fréquence cardiaque maximale atteinte (thalach)") + # Nom de l'axe des ordonnées + + # Mise en forme du titre et de la légende + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 12, # Taille du titre + face = "bold")) # Titre en gras +``` +Cette figure fait clairement apparaître une relation décroissante et linéaire entre les deux variables. On en déduit que la fréquence cardiaque maximale diminue avec l'âge. + +On étudie maintenant graphiquement le lien entre la douleur thoracique (*cp*) et la présence / l'absence d'angine de poitrine induite par l'exercice (*exang*). +```{r} +# Histogramme entre cp et exang +ggplot(df, aes(x = cp, fill = factor(exang))) + + geom_bar(color = "black", # Couleur des contours + linewidth = 0.3, # Épaisseur des contours + alpha = 0.5) + # Transparence de la couleur de remplissage + + # Nom des axes et titre + labs(title = "Graphique de la distribution de la douleur thoracique (cp) en fonction de \n la présence / l'absence d'angine de poitrine induite par l'exercice (exang)", # Titre du graphique + x = "Douleur thoracique (cp)", # Nom de l'axe des abscisses + y = "Nombre d'individus", # Nom de l'axe des ordonnées + fill = "Angine thoracique \ninduite par l'exercice") + # Titre de la légende + + # Texte de la légende et couleur + scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4"), + labels = c("Non", "Oui")) + + + # Mise en forme du titre et de la légende + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 10, # Taille du titre + face = "bold"), # Titre en gras + legend.title = element_text(size = 9), # Taille du titre de la légende + legend.text = element_text(size = 8)) # Taille du texte de la légende +``` +On remarque que les individus n'ayant pas de douleur thoracique ont des angines thoraciques induites par l'exercice et inversement : les individus ayant des douleurs thoraciques n'ont pas d'angine thoracique induite par l'exercice. + +La relation entre la fréquence cardiaque maximale atteinte (*thalach*) et la pente du segment ST de pointe de l'exercice (*slope*) est représentée dans le graphique ci-dessous. +```{r} +# Distribution de thalach en fonction de slope +ggplot(df, aes(x = thalach, fill = factor(slope))) + + geom_histogram(color = "black", # Couleur des contours + linewidth = 0.3, # Épaisseur des contours + bins = 25, # Nombre de catégories + alpha = 0.5) + # Transparence de la couleur de remplsisage + + # Nom des axes et titre + labs(title = "Graphique de la distribution de la fréquence cardiaque maximale atteinte (thalach) \n en fonction de la pente du segment ST de pointe de l'exercice (slope)", # Titre du graphique + x ="Fréquence cardiaque maximale atteinte (thalach)", # Nom de l'axe des abscisses + y = "Nombre d'individus", # Nom de l'axe des ordonnées + fill = "Pente du segment \nST (slope)") + # Titre de la légende + + # Couleur de la légende + scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4", "2" = "#C3D7A4")) + + + # Mise en forme du titre et de la légende + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 10, # Taille du titre + face = "bold"), # Titre en gras + legend.title = element_text(size = 9), # Taille du titre de la légende + legend.text = element_text(size = 8)) # Taille du texte de la légende +``` +Ce graphique montre que plus la fréquence cardiaque maximale atteinte est élevée (entre 160 et 200), plus la pente du segment ST est importante (niveau 2). + +Enfin, on trace la distribution de la fréquence cardiaque maximale atteinte (*thalach*) en fonction de la présence / l'absence d'angine de poitrine induite par l'exercice (*exang*) . +```{r} +# Distribution de thalach en fonction de exang +ggplot(df, aes(x = thalach, fill = factor(exang))) + + geom_histogram(color = "black", # Couleur des contours + linewidth = 0.3, # Épaisseur des contours + bins = 30, # Nombre de catégories + alpha = 0.5) + # Transparence de la couleur de remplsisage + + # Noms des axes et titre + labs(title = "Graphique de la distribution de la fréquence cardiaque maximale atteinte (thalach) \n en fonction de la présence / l'absence d'angine de poitrine induite par l'exercice (exang) ", # Titre du graphique + x = "Fréquence cardiaque maximale atteinte (thalach)", # Nom de l'axe des abscisses + y = "Nombre d'individus", # Nom de l'axe des orodnnées + fill = "Angine thoracique induite \npar l'exercice (exang)") + # Titre de la légende + + # Couleur et texte de la légende + scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4"), + labels = c("Non", "Oui")) + + + # Mise en forme du titre et de la légende + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 10, # Taille du titre + face = "bold"), # Titre en gras + legend.title = element_text(size = 9), # Taille du titre de la légende + legend.text = element_text(size = 8)) # Taille du texte de la légende +``` +On observe que les individus ayant une fréquence cardiaque maximale atteinte plus faible (entre 50 et 150) ont, en général, des angines induites par l'exercice et inversement. + +On décide pour l'instant de garder les variables de chaque couple. On fera une étude plus approfondie par la suite. + +### Relation entre les variables explicatives et la variable cible + +#### Aperçu global + +Pour avoir un aperçu global et visuel des variables qui ont un impact sur la variable cible, on affiche un graphe des coefficients de corrélation entre la variable cible et chacune des autres variables. + +```{r} +# Matrice des corrélations de toutes les variables +cor_matrix <- cor(df[, 1:14]) + +# Récupérer les coefficients de corrélation entre la variable cible et les autres variables +target_correlations <- cor_matrix[, 14] # La 14ème colonne correspond à la variable cible + +# Création d'un dataframe avec les noms des variables et leurs corrélations avec la variable cible +correlation_data <- data.frame(variable = names(target_correlations)[-14], + correlation = unlist(target_correlations[-14])) + +# Création d'un graphique (barplot) des coefficients de corrélation +ggplot(correlation_data, aes(x = variable, y = correlation)) + + geom_bar(stat = "identity", # Hauteur des barres = coefficient de corrélation + fill = "#4E84C4", # Couleur de remplissage + alpha = 0.7, # Transparence de la couleur de remplissage + color = "black", # Couleur des contours + linewidth = 0.3) + # Épaisseur des contours + labs(title = 'Corrélations entre la variable cible et les autres variables', # Titre du graphique + x = 'Variables', # Nom de l'axe des abscisses + y = 'Coefficient de corrélation') + # Nom de l'axe des ordonnées + + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 12, # Taille du titre + face = "bold")) # Titre en gras +``` +Les variables qui semblent avoir un impact (un coefficient de corrélation proche de 1 en valeur absolue) sur la target sont : *ca, cp, exang, oldpeak, sex, slope, thal* et *thalach*. + +#### Relation entre les variables catégorielles et la variable cible + +On va maintenant regarder plus précisemment les relations entre les variables explicatives et la target. On commence par étudier les variables catégorielles puis on étudiera les variables numériques. + +Grâce à l'aperçu global, on a vu que seules les variables catégorielles *sex, cp, exang, slope, ca* et *thal* avaient potentiellement une relation avec la variable cible. On affiche alors un histogramme de chacune de ces variables en fonction de la target pour voir graphiquement cette relation. + +```{r} +# Création d'un vecteur avec les variables catégorielles qui semblent avoir un impact sur la variable cible +categorical_var_impact <- c("sex", "cp", "exang", "slope", "ca", "thal") + +# Création d'une liste pour stocker les graphiques +hist <- list() + +# Boucle qui crée un histogramme pour chaque variable du vecteur categorical_var_impactet et qui l'affiche +for (var in categorical_var_impact) { + h <- ggplot(df, aes(x = .data[[var]], fill = factor(target))) + + geom_bar(color = "black", # Couleur des contours + linewidth = 0.3, # Épaisseur des contours + position = "fill", # Représenter des proportions + alpha = 0.7) + # Transparence de la couleur de remplissage + + # Noms des axes et titre + labs(title = paste("Distribution de", var, "\nen fonction de la variable cible"), # Titre des graphiques + x = var, # Nom de l'axe des abscisses + y = "Nombre d'individus", # Nom de l'axe des ordonnées + fill = "Légende") + # Titre de la légende + + # Couleur et texte de la légende + scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4"), + labels = c("Sain", "Malade")) + + + # Mise en forme du titre et de la légende + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 8, # Taille du titre + face = "bold"), # Titre en gras + legend.title = element_text(size = 7), # Taille du titre de la légende + legend.text = element_text(size = 7), # Taille du texte de la légende + legend.key.size = unit(0.3, "cm"), # Taille des carrés de couleur de la légende + axis.title.x = element_text(size = 7), # Taille du nom de l'axe des abscisses + axis.title.y = element_text(size = 7)) # Taille du nom de l'axe des ordonnée + + hist[[var]] <- h +} + +# Organisation des graphiques dans une grille à 3 colonnes +grid.arrange(grobs = hist, ncol = 3) +``` + +On peut observer sur ces graphiques que les variables catégorielles selectionnées semblent effectivement avoir un impact sur la variable cible. En effet, les individus qui semblent plus enclin à avoir une maladie cardiaque sont : + +* les femmes, +* les individus ayant des douleurs thoraciques (de niveau 1, 2 ou 3), +* les individus ayant des anomalies mineures dans l’ECG (étonnant), +* les individus n'ayant pas d'angine de poitrine induite par l'exercice (là encore c'est étonnant), +* les individus ayant des changements dans le segment ST graves (de niveau 2), +* les individus n'ayant pas de vaisseaux principaux colorés par la flourosopie, +* les individus ayant des défauts modérés de perfusion myocardique (niveau 2). + +**Remarque** : des précisions sur les observations et conclusions étonnantes seront apportées dans la conclusion. + +#### Relation entre les variables numériques et la variable cible + +Enfin, grâce à l'aperçu global, on a vu que seules les variables numériques *thalach* et *oldpeak* avaient potentiellement une relation avec la variable cible. On affiche alors la densité de chacune de ces variables en fonction de la target pour voir graphiquement cette relation. +```{r} +# Graphique de target en fonction de thalach + +ggplot(df, aes(x = thalach, fill = factor(target))) + + geom_histogram(color = "black", # Couleur des contours + linewidth = 0.3, # Épaisseur des contours + bins = 20, # 20 classes pour plus de lisibilité + alpha = 0.5) + # Transparence de la couleur de remplissage + + # Noms des axes et titre + labs(title = "Graphique de la distribution de la fréquence cardiaque maximale atteinte (thalach) \n en fonction de la variable cible (target)", # Titre du grahique + x = "Fréquence cardiaque maximale atteinte (thalach)", # Nom de l'axe des abscisses + y = "Nombre d'individus", # Nom de l'axe des ordonnées + fill = "Légende") + # Nom du titre de la légende + + # Couleur et texte de la légende + scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4"), + labels = c("Sain", "Malade")) + + + # Mise en forme du titre et de la légende + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 12, # Taille du titre + face = "bold"), # Titre en gras + legend.title = element_text(size = 10), # Taille du titre de la légende + legend.text = element_text(size = 10), # Taille du texte de la légende + legend.key.size = unit(0.3, "cm"), # Taille des carrés de couleur de la légende + axis.title.x = element_text(size = 10), # Taille du nom de l'axe des abscisses + axis.title.y = element_text(size = 10)) # Taille du nom de l'axe des ordonnée +``` + +On peut observer sur ce graphique que la variable *thalach* semble effectivement avoir un impact sur la variable cible. En effet, les individus qui semblent plus enclin à avoir une maladie cardiaque sont les individus ayant une fréquence cardiaque maximale entre 150 et 200 bpm. + +Enfin, on trace ci-dessous la distribution de la la dépression ST induite par l'exercice (oldpeak) en fonction de la variable cible. + +```{r} +# Distribution de oldpeak en fonction de target +ggplot(df, aes(x = oldpeak, fill = factor(target))) + + geom_density(alpha = 0.5, # Transparence de la couleur de remplissage + linewidth = 0.3) + # Épaisseur des contours + + # Noms des axes et titre + labs(title = "Graphique de la distribution de la la dépression ST induite par l'exercice (oldpeak) \n en fonction de la variable cible (target)", # Titre du graphique + x = "Changements dans l'électrocardiogramme (oldpeak)", # Nom de l'axe des abscisses + y = "Nombre d'individus", # Nom de l'axe des ordonnées + fill = "Légende") + # Nom du titre de la légende + + # Couleur et texte de la légende + scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4"), + labels = c("Sain", "Malade")) + + + # Mise en forme du titre et de la légende + theme(plot.title = element_text(hjust = 0.5, # Position du titre + size = 12, # Taille du titre + face = "bold"), # Titre en gras + legend.title = element_text(size = 9), # Taille du titre de la légende + legend.text = element_text(size = 8)) # Taille du texte de la légende +``` +On conclut que les individus n'ayant pas de changement dans l'electrocardiogramme ont plus tendance à avoir une maladie cardiaque et inversement. + +## Sampling : division du jeu de données + +On divise le jeu de données en 2 ensembles : + +* un ensemble d'entrainement qui va permettre de construire le modèle, +* un ensemble de test pour évaluer les performances du modèle. + +On choisit une porportion souvent utilisée pour un jeu de données de petite taille : 80% des données sont utilisées pour former l'ensemble d'entrainement et 20% pour celui de test. +On affiche ci-dessous les dimensions des deux ensembles et celles du jeu de données de base. +```{r} +# Division des données en ensembles d'entraînement et de test +set.seed(123) # pour la reproductibilité +split <- rsample::initial_split(df, prop = 0.8) +train <- rsample::training(split) +test <- rsample::testing(split) + +# Création d'un dataframe avec les dimensions des différents ensembles +dimensions <- data.frame(dim(df), dim(train), dim(test)) +rownames(dimensions) <- row_names <- c("Nombre de lignes", + "Nombre de colonnes") # nom des lignes +colnames(dimensions) <- c("Dimension du jeu de données initial", + "Dimensions de l'ensemble d'entrainement", + "Dimensions de l'ensemble de test") # nom des colonnes + +# Affichage des dimensions des ensembles d'entraînement et de test +dimensions %>% + kbl() %>% + kable_styling() +``` +On regarde la distribution de la variable cible dans les 2 ensembles pour s'assurer que la répartition des classes dans l'ensemble de test et d'entrainement est représentative de l'ensemble complet. + +```{r} +# Organisation des graphiques en une grille 1x2 +par(mfrow = c(1, 2)) + +# Création du camembert pour l'ensemble train avec les pourcentages et légendes +pie(table(train$target), # Table des effectifs pour la variable target + main = "Distribution de la variable cible (target) \nsur l'ensemble d'entrainement", # Titre + labels = c("Individus en bonne santé", + "Individus atteints \nd'une maladie cardiaque"), # Légende + col = c("lightblue", "lightcoral"), # Définition de la couleur des secteurs + cex = 0.8, # Ajustement de la taille de la légende + cex.main = 1) # Ajustement de la taille du titre + +# Ajout des étiquettes +percent_values_train <- round(100 * prop.table(table(train$target)), 1) # Calcul des pourcentages arrondis +percent_train <- paste(percent_values_train, "%") # Ajout su sigle "%" +text_pie(percent_values_train, c(percent_train[1], percent_train[2]), cex=0.8) # Utilisation de la fonction text_pie définie plus haut + + +# Création du camembert pour l'ensemble test avec les pourcentages et légendes +pie(table(test$target), # Table des effectifs pour la variable target + main = "Distribution de la variable cible (target) \nsur l'ensemble de test", # Titre + labels = c("Individus en bonne santé", + "Individus atteints \nd'une maladie cardiaque"), # Légende + col = c("lightblue", "lightcoral"), # Définition de la couleur des secteurs + cex = 0.8, # Ajustement de la taille de la légende + cex.main = 1) # Ajustement de la taille du titre + +# Ajout des étiquettes +percent_values_test <- round(100 * prop.table(table(test$target)), 1) # Calcul des pourcentages arrondis +percent_test <- paste(percent_values_test, "%") # Ajout su sigle "%" +text_pie(percent_values_test, c(percent_test[1], percent_test[2]), cex=0.8) # Utilisation de la fonction text_pie définie plus haut + +``` + +# Fitting : création et entraînement du modèle + +## Déclaration du modèle + +On implémente maintenant une régression logistique à l'aide de la fonction "glm". On pourra trouver une explication du choix du modèle dans la notice. +On déclare aussi le modèle réduit à l'intercept car il sera utile dans la suite, pour des comparaisons de perfromance par exemple. + +```{r} +# Déclaration du modèle de regression logistique +mod <- glm(target ~ ., data = train, family = binomial) + +# Déclaration du modèle réduit à l'intercept (utile dans la suite) +mod_int <- glm(target~1, data = train, family = "binomial") +``` + +On affiche les statistqiues récapitulatives du modèle de regression logistique construit dans la sous-section section précédente. Les p-valeurs inférieures à 0.05 sont affichées en rouge car cela signifie que la variable associée à probablement de l'importance dans le modèle. + +```{r} +# Création d'un tableau qui contient les coefficients du modèle +tab_stat <- tidy(mod, conf.int = TRUE) + +# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05 +tab_stat <- tab_stat %>% + mutate(p.value = ifelse(p.value < 0.05, + cell_spec(round(p.value, 3), "html", color = "red"), + cell_spec(round(p.value, 3), "html"))) + +# Affichage du tableau avec mise en forme +tab_stat %>% + kbl(digits = c(0, 2, 2, 2, 5, 3, 3), # Arrondis + escape = FALSE, + caption = "Statistiques récapitulatives du modèle de regression logistique") %>% + kable_styling() +``` + +D'après cet affichage, les variables qui ont une importance sont : *sex, cp, exang, oldpeak, ca* et *thal*. + +On peut noter que la déviance du modèle est inférieure à celle du modèle réduit à l'intercept : + +```{r} +# Création d'un dataframe avec les déviances +deviance <- data.frame(mod$null.deviance, mod$deviance) +colnames(deviance) <- c("Déviance du modèle réduit à l'intercept", + "Déviance du modèle créé") # Nom des colonnes + +# Affichage des dimensions des ensembles d'entraînement et de test +deviance %>% + kbl(digits = c(2, 2), # Arrondis + caption = "Déviance des modèles") %>% + kable_styling() +``` +Le modèle est donc meilleur que le modèle réduit à l'inetrcept. + +## Sélection de modèles : tests et méthodes pas à pas + +On utilise le test Anova pour enlever les variables qui n'apportent pas d'informations complémentaires sur la variable à prédire. On affiche le résultat du test. Les p-valeurs inférieurs à 0.05 sont affichées en rouge. +```{r} +# Test Anova +Ano <- Anova(mod, + type = 3, # type III de l'Anova + test.statistic = "LR") # test de rapport de vraisemblance + +# Création d'un tableau qui contient les coefficients du modèle +tab_stat <- tidy(Ano, conf.int = TRUE) + +# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05 +tab_stat <- tab_stat %>% + mutate(p.value = ifelse(p.value < 0.05, + cell_spec(round(p.value, 3), "html", color = "red"), + cell_spec(round(p.value, 3), "html"))) + +# Affichage du tableau avec mise en forme +tab_stat %>% + kbl(digits = c(0, 3, 5), + escape = FALSE, + caption = "Résultats du test Anva") %>% + kable_styling() +``` +Seules les variables *sex, cp, exang, oldpeak, ca* et *thal* semblent avoir un impact sur la variable cible, ce qu'on avait déjà observé dans la partie "graphiques". On avait aussi remarqué que les variables *exang, slope* et *age* avaient une importance. On va continuer notre étude pour savoir si on les supprime ou non. + +On effectue maintenent des méthodes pas à pas pour sélectionner un modèle à l'aide d'une procédure basée sur la minimisation du critère AIC (Akaike Information Criterion). + +On réalise d'abord une méthode *backward* et on affiche les statistqiues descriptives du modèle obtenu. Les p_valeurs inférieures à 0.05 sont affichées en rouge. +```{r} +# Méthode backward +mod_back <- step(mod, direction = "backward", trace = FALSE) + +# Création d'un tableau qui contient les coefficients du modèle +tab_stat <- tidy(mod_back, conf.int = TRUE) + +# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05 +tab_stat <- tab_stat %>% + mutate(p.value = ifelse(p.value < 0.05, + cell_spec(round(p.value, 3), "html", color = "red"), + cell_spec(round(p.value, 3), "html"))) + +# Affichage du tableau avec mise en forme +tab_stat %>% + kbl(digits = c(0, 3, 3, 3, 4), # Arrondis + escape = FALSE, + caption = "Statistiques récapitulatives du modèle construit par la méthode backward") %>% + kable_styling() +``` +Les variables importantes sont celles du test Anova ainsi que *restecg* et *thalach*. + +On refait alors un test Anova sur le modèle obtenu par la méthode *backward* pour vérifier si les variables *restecg* et *thalach* sont réellement importantes. +```{r, include = FALSE} +# Test Anova +Ano_back <- Anova(mod_back, type = 3,test.statistic = "LR") + +# Création d'un tableau qui contient les coefficients du modèle +tab_stat <- tidy(Ano_back, conf.int = TRUE) + +# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05 +tab_stat <- tab_stat %>% + mutate(p.value = ifelse(p.value < 0.05, + cell_spec(round(p.value, 3), "html", color = "red"), + cell_spec(round(p.value, 3), "html"))) + +# Affichage du tableau avec mise en forme +tab_stat %>% + kbl(digits = c(0, 3, 3, 3, 4), # Arrondis + escape = FALSE, + caption = "Résultats du test Anova effectué sur le modèle construit à partir de la méthode backward") %>% + kable_styling() +``` +Le test donne le même résultat que la méthode *backward*. + +On réalise maintenant une méthode *forward* de sélection de modèle et on affiche là encore, les statistqiues descriptives du modèle obtenu. Les p_valeurs inférieures à 0.05 sont affichées en rouge. +```{r} +# Méthode forward +mod_for <- step(mod_int, + target ~ age + sex + trestbps + chol + fbs + restecg + thalach + + exang + oldpeak + slope + ca + thal, + data = train, + trace = FALSE, + direction = c("forward")) + +# Création d'un tableau qui contient les coefficients du modèle +tab_stat <- tidy(mod_for, conf.int = TRUE) + +# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05 +tab_stat <- tab_stat %>% + mutate(p.value = ifelse(p.value < 0.05, + cell_spec(round(p.value, 3), "html", color = "red"), + cell_spec(round(p.value, 3), "html"))) + +# Affichage du tableau avec mise en forme +tab_stat %>% + kbl(digits = c(0, 3, 3, 3, 4), # Arrondis + escape = FALSE, + caption = "Statistiques récapitulatives du modèle construit par la méthode forward") %>% + kable_styling() +``` + +Les variables importantes sont celles du test Anova sans la variable *cp*. + +On réalise une dernière méthode : la méthode *both* en affichant toujours les statistqiues descriptives du modèle obtenu. Les p_valeurs inférieures à 0.05 sont affichées en rouge. + +```{r} +# Méthode both +mod_both <- step(mod_int, + target ~ age + sex + trestbps + chol + fbs + restecg + thalach + + exang + oldpeak + slope + ca + thal, + data = train, + trace = F, + direction = c("both")) + +# Création d'un tableau qui contient les coefficients du modèle +tab_stat <- tidy(mod_both, conf.int = TRUE) + +# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05 +tab_stat <- tab_stat %>% + mutate(p.value = ifelse(p.value < 0.05, + cell_spec(round(p.value, 3), "html", color = "red"), + cell_spec(round(p.value, 3), "html"))) + +# Affichage du tableau avec mise en forme +tab_stat %>% + kbl(digits = c(0, 3, 3, 3, 4), # Arrondis + escape = FALSE, + caption = "Statistiques récapitulatives du modèle construit par la méthode both") %>% + kable_styling() + +``` +Le modèle garde les mêmes variables que celles obtenues avec la méthode *forward*. + +Enfin, on compare l'AIC des modèles obtenus par les 3 méthodes. + +```{r} +# Création d'un dataframe avec les AIC +deviance <- data.frame(AIC(mod_back), AIC(mod_for), AIC(mod_both)) +colnames(deviance) <- c("Méthode backward", + "Méthode forward", + "Méthode both") # Nom des colonnes + +# Affichage des AIC des modèles +deviance %>% + kbl(digits = c(1, 1, 1), # pour les arrondis + caption = "AIC des modèles") %>% # Titre du tableau + kable_styling() +``` + +## Construction du modèle final + +La méthode backward est la meilleure car elle a le plus petit AIC. On garde donc les variables du modèle construit à partir de cette méthode. +On définit alors notre modèle final avec les variables *restecg, cp, sex, exang, oldpeak, ca, thal* et *thalach*. On retire les variables *age, trestbps, chol, fbs* et *slope*. + +```{r} +# Construction du modèle final +mod_final <- glm(target ~ sex + exang + oldpeak + ca + thal + restecg + thalach + cp, + data = train, + family = "binomial") +``` + +# Validation du modèle + +Après avoir obtenu un modèle, il faut diagnostiquer la régression afin de valider ou non le modèle. +On passe alors à la validation de notre modèle ainsi construit. On regarde dans un premier temps la déviance du modèle puis on étudie les résidus et les *outliers*. + +## Déviance + +```{r} +# Création d'un dataframe avec les déviances +deviance <- data.frame(mod$null.deviance, mod_final$deviance) +colnames(deviance) <- c("Déviance du modèle réduit à l'intercept", + "Déviance du modèle final") # nom des colonnes + +# Affichage des dimensions des ensembles d'entraînement et de test +deviance %>% + kbl(digits = c(2, 2), caption = "Déviance des modèles") %>% # pour les arrondis + kable_styling() +``` + +La déviance de notre modèle étant plus petite que celle du modèle réduit à l'intercept, on le privilegie. + +## Résidus studentisés + +On affiche ci-dessous un nuage de points montrant les résidus studentisés du modèle final pour chaque observation. +On a aussi tracé sur ce même graphe des lignes de seuil pour repérer les observations ayant des résidus considérés comme atypiques en dehors de ces bornes. + +```{r} +par(mfrow = c(1, 1)) +plot(rstudent(mod_final), type = "p", cex = 0.5, ylab = "Résidus studentisés ", + col = "#0066CC", ylim = c(-3, 3)) +abline(h = c(-2, 2), col = "#CC3333") +``` +On remarque que les résidus suivent un schéma aléatoire, ce qui valide le modèle. +Les quelques points au-delà des lignes rouges sont les valeurs atypiques que nous avions déjà remarqué dans les boxplots. Nous avions décidé de garder ces valeurs. + +## Outliers + +Enfin, on affiche le graphe des *outliers*. +```{r} +# Plot des outliers +plot(mod_final,5) +``` +Il ne semble pas y avoir d'outliers car aucun point n'a une distance de Cook supérieure à 1. + +# Optimisation du seuil de décision + +On va maintenant optimsier le seuil à partir duquel on considère qu'un individu à une maladie cardiaque. Dans notre cas, on préfère naturellemnt avoir un petit nombre de faux négatifs (FN) c'est-à-dire un petit nombre de personnes qui n'ont pas été diagnostiquées comme positives alors qu'elles ont une maladie cardiaque. + +Pour optimiser le seuil, on va donc maximsier l'indice de Youden qui est un critère couramment utilisé pour minimiser le nombre de faux négatif. + +On affiche ici un graphe qui montre l'évolution de l'indice de Youden en fonction des seuils. On affiche sur ce même graphe le seuil optimal qui maximise cet indice. + +```{r, include = FALSE} +# Prédire les probabilités sur l'ensemble d'entraînement +predicted_probs_train <- predict(mod_final, newdata = train, type = "response") +``` + +```{r} +# Calcul des valeurs de seuil +val_seuil <- seq(0, 1, length.out = 100) + +# Calcul de l'indice de Youden pour chaque seuil +youden_index <- sapply(val_seuil, function(threshold) { + pred_pos <- predicted_probs_train > threshold + sens <- sum(pred_pos & (train$target == 1)) / sum(train$target == 1) + spec <- sum(!pred_pos & (train$target == 0)) / sum(train$target == 0) + return(sens + spec - 1) +}) + +# Trouver le seuil optimal +optimal <- val_seuil[which.max(youden_index)] + +# Créer le graphique de l'indice de Youden en fonction des seuils +plot(val_seuil, + youden_index, + type = "l", + col = "#0066CC", # Couleur de la courbe + xlab = "Seuils", # Nom de l'axe des abscisses + ylab = "Indice de Youden",# Nom de l'axe des ordonnées + main = "Graphique de l'évolution de l'ndice de Youden en fonction des seuils") # Titre du graphique + +# Afficher le seuil optimal +abline(v = optimal, + col = "#CC3333", # Couleur de la droite + lty = 2) # Épaisseur de la droite + +# Afficher la légende +text(optimal, + max(youden_index), + paste("Seuil optimal =", round(optimal, 3)), + pos = 4, + col = "#CC3333") +``` +Le seuil qui maximise l'indice de Youden vaut 0.556. On va donc garder ce seuil et prédire à patir de celui-ci. +On affiche ci-dessous l'évolution de la sensibilité et de la spécificité en fonction de la valeur du seuil. On affiche sur ce même-graphique la valeur du seuil optimal. + +```{r} +# Initialisation des vecteurs de sensibilité et spécificité +sensitivities <- c() +specificities <- c() + +for (threshold in val_seuil) { + pred_pos <- predicted_probs_train > threshold + sens <- sum(pred_pos & (train$target == 1)) / sum(train$target == 1) + spec <- sum(!pred_pos & (train$target == 0)) / sum(train$target == 0) + + sensitivities <- c(sensitivities, sens) + specificities <- c(specificities, spec) +} + +# Création du graphique avec la courbe de l'évolution de la sensibilité +plot(val_seuil, + sensitivities, + type = "l", + col = "#0066CC", # Couleur de la courbe + xlab = "Seuil", # Nom de l'axe des abscisses + ylab = "Valeur des paramètres", # Nom de l'axe des ordonnées + main = "Évolution de la sensibilité et de la spécificité en fonction des seuils") # Titre du graphique + +# Ajout de la courbe de l'évolution de la spécificité +lines(val_seuil, + specificities, + type = "l", + col = "#339966") # Couleur de la droite + +# Ajout de la légende +legend("bottomright", + legend = c("Sensibilité", "Spécificité"), + col = c("#0066CC", "#339966"), + lty = 1) + +# Ajout de la valeur du seuil optimal +abline(v = optimal, + col = "#CC3333", + lty = 2) + +# Ajout de texte +text(optimal, + max(youden_index), + paste("Seuil optimal =", round(optimal, 3)), + pos = 4, + col = "#CC3333") +``` +On remarque donc que la valeur optimale pour l'indice de Youden ne l'est pas pour la sensibilité ni la spécificité. Cela est du au fait qu'on a du faire un choix en selectionnant un indicateur à optimiser parmi plusieurs autres, ici, le nombre de faux négatifs. Cependant, au seuil optimal, la spécificité et la sensibilité sont satisfaisantes. On garde donc ce seuil. + +# Prédiction et performances + +## Prédiction + +On passe maintenant à la prédiction. On prédit de telle manière que si la probabilité prédite dépasse le seuil optimal, l'individu est considéré comme étant susceptible de développer une maladie cardiaque. +```{r, include = FALSE} +# Prédire les probabilités sur l'ensemble d'entraînement +predicted_probs_test <- predict(mod_final, newdata = test, type = "response") + +# Prédiction des valeurs sur l'ensemble de test +predicted <- ifelse(predicted_probs_test > optimal, 1, 0) +summary(predicted_probs_test) +``` + +On affiche ci-dessous les prédictions à côté des valeurs réelles dans un tableau. +```{r} +# Création d'un dataframe avec les observations et les prédictions +df_pred <- data.frame(Observed = test$target, Predicted = predicted) +rmarkdown::paged_table(df_pred) +``` + + +## Matrice de confusion + +Pour étudier les performances de prédiction, on affiche la matrice de confusion qui permet d'observer le nombre de vrais positifs, de faux positifs, de vrais négatifs et de faux négatifs. + +```{r} +# Création de la matrice de confusion +confusion_matrix <- table(factor(predicted), factor(test$target)) +# Affichage de la matrice de confusion +par(mfrow=c(1, 1)) +ctable <- as.table(confusion_matrix, + nrow = 2, + byrow = TRUE) +# Définir les étiquettes "malade" et "sain" pour les lignes et colonnes +rownames(ctable) <- colnames(ctable) <- c("malade (1)", "sain (0)") + +# Définir les étiquettes "Réel" et "Prédit" +dimnames(ctable) <- list(Réel = rownames(ctable), Prédit = colnames(ctable)) + +fourfoldplot(ctable, + color = c("#CC6666", "#CCCCCC"), + conf.level = 0, + margin = 1, + main = "Matrice de confusion \ndu modèle de regression logistique") + +``` +On obtient qu’on a prédit 3 faux négatifs, 7 faux positifs, 21 vrais positifs et 30 vrais négatifs. Ce résultat est tout à fait satisfaisait car le nombre de faux négatifs est assez faible, ce qui était notre objectif. + +## Indicateurs et métriques de performance + +La sensibilité, le F1 Score et la valeur prédictive négative sont des indicateurs importants pour évaluer la capacité du modèle à minimiser les faux négatifs, chacun apportant une compréhension particulière sur cette problématique. On s'intéresse donc ici plus particulièrement à ces indicateurs ainsi que l'accuracy (la précision globale). +```{r} +# Créer des facteurs avec des niveaux correspondants +predicted_factor <- factor(as.numeric(predicted), levels = c(0, 1)) +test_target_factor <- factor(as.numeric(test$target), levels = c(0, 1)) + +# Simulation de la sortie confusionMatrix +results <- confusionMatrix(predicted_factor, test_target_factor, mode = "everything", positive = "1") + +# Extraction de l'accuracy +acc <- data.frame(Valeur = results$overall["Accuracy"]) + +# Extraction des valeurs spécifiques +specific_metrics <- c("Sensitivity", "Neg Pred Value", "F1") +specific_values <- unlist(results$byClass[specific_metrics]) + +# Création du tableau avec les valeurs spécifiques +results_df_spe <- data.frame(Valeur = specific_values) + +# Affichage du tableau +rbind(acc, results_df_spe) %>% # Fusion des deux dataframe + kable(digits = 3, # Arrondi + caption = "Indicateurs et métriques de performance.") %>% # Titre + kable_styling() +``` + +Grâce au tableau, on observe que la valeur des indicateurs est plutôt élevée et que donc notre modèle n'est pas trop mauvais ! + +## Courbe ROC et AUC + +On trace enfin la courbe ROC de notre modèle afin de déterminer l'AUC de notre modèle. +```{r} +# Courbe ROC +pred <- prediction(predicted_probs_test, test$target) +roc <- performance(pred, measure = "tpr", x.measure = "fpr") + +# Calcul de l'AUC +auc <- performance(pred, measure = "auc") +auc_value <- round(unlist(slot(auc, "y.values")), 3) + +# Plot de la courbe ROC +plot(roc, xlab = "Taux de faux positifs", ylab = "Taux de vrais positifs", col = "#CC3333", main = "Courbe ROC de l'ensemble de test") +abline(a = 0, b = 1) + +# Ajout de l'AUC au graphique +text(0.5, 0.8, labels = paste("AUC =", auc_value), cex = 1.2, col = "black") +``` +Dans notre cas, l'AUC est proche de 1 (environ 0.9) donc la prédiction semble être bonne. + +# Comparaison avec d'autres modèles + +Pour terminer le projet, on compare le modèle de régression logistique à 2 autres modèles : un modèle de RadomForest et un modèle d’arbre de décision. + +## Déclaration des modèles + +On construit d'abord un modèle de forêt aléatoire (*Random Forest*) qui est populaire en raison de sa capacité à produire des prédictions précises, à gérer des ensembles de données complexes et bruités, et à éviter le surajustement. + +```{r, include = FALSE} + +# Créer un modèle de Random Forest +mod_rf <- randomForest(target ~., data = train, ntree = 500) + +# Prédire les probabilités sur l'ensemble d'entraînement +predicted_probs_train_rf <- predict(mod_rf, newdata = train, type = "response") + +# Calcul des valeurs de seuil +val_seuil <- seq(0, 1, length.out = 100) + +# Calcul de l'indice de Youden pour chaque seuil +youden_index <- sapply(val_seuil, function(threshold) { + pred_pos <- predicted_probs_train_rf > threshold + sens <- sum(pred_pos & (train$target == 1)) / sum(train$target == 1) + spec <- sum(!pred_pos & (train$target == 0)) / sum(train$target == 0) + return(sens + spec - 1) +}) + +# Trouver le seuil optimal +optimal_rf <- val_seuil[which.max(youden_index)] + +# Prédire les probabilités sur l'ensemble d'entraînement +predicted_probs_test_rf <- predict(mod_rf, newdata = test, type = "response") + +# Prédiction des valeurs sur l'ensemble de test +predicted_rf <- ifelse(predicted_probs_test_rf > optimal_rf, 1, 0) +summary(predicted_probs_test_rf) + +# Créer des facteurs avec des niveaux correspondants +predicted_factor_rf <- factor(as.numeric(predicted_rf), levels = c(0, 1)) +test_target_factor_rf <- factor(as.numeric(test$target), levels = c(0, 1)) + +# Simulation de la sortie confusionMatrix +results_rf <- confusionMatrix(predicted_factor_rf, test_target_factor_rf, mode = "everything", positive = "1") +``` + + +On construit ensuite un modèle d'arbre de décision, qui est plus approprié pour capturer des relations non linéaires entre les variables explicatives et la variable cible et des interactions complexes. +```{r, include = FALSE} +# Créer un modèle d'arbre de décision +mod_tree <- rpart(target ~ ., data = train, method = "class") + +# Prédire les probabilités sur l'ensemble d'entraînement +predicted_probs_train_tree <- predict(mod_tree, newdata = train)[, "1"] + +# Calcul des valeurs de seuil +val_seuil <- seq(0, 1, length.out = 100) + +# Calcul de l'indice de Youden pour chaque seuil +youden_index <- sapply(val_seuil, function(threshold) { + pred_pos <- predicted_probs_train_tree > threshold + sens <- sum(pred_pos & (train$target == 1)) / sum(train$target == 1) + spec <- sum(!pred_pos & (train$target == 0)) / sum(train$target == 0) + return(sens + spec - 1) +}) + +# Trouver le seuil optimal +optimal_tree <- val_seuil[which.max(youden_index)] + +# Prédire les probabilités sur l'ensemble d'entraînement +predicted_probs_test_tree <- predict(mod_tree, newdata = test)[, "1"] + +# Prédiction des valeurs sur l'ensemble de test +predicted_tree <- ifelse(predicted_probs_test_tree > optimal_tree, 1, 0) +summary(predicted_probs_test_tree) + +# Créer des facteurs avec des niveaux correspondants +predicted_factor_tree <- factor(as.numeric(predicted_tree), levels = c(0, 1)) +test_target_factor_tree <- factor(as.numeric(test$target), levels = c(0, 1)) + +# Simulation de la sortie confusionMatrix +results_tree <- confusionMatrix(predicted_factor_tree, test_target_factor_tree, mode = "everything", positive = "1") +``` + +On peut afficher les matrices de confusions des modèles construits. + +```{r} +par(mfrow=c(1, 2)) +# Affichage de la matrice de confusion de la forêt aléatoire +ctable <- as.table(results_rf, + nrow = 2, + byrow = TRUE) +# Définir les étiquettes "malade" et "sain" pour les lignes et colonnes +rownames(ctable) <- colnames(ctable) <- c("malade (1)", "sain (0)") + +# Définir les étiquettes "Réel" et "Prédit" +dimnames(ctable) <- list(Réel = rownames(ctable), Prédit = colnames(ctable)) + +fourfoldplot(ctable, + color = c("#CC6666", "#CCCCCC"), + conf.level = 0, + margin = 1, + main = "Matrice de confusion \ndu modèle de forêt aléatoire") + +# Affichage de la matrice de confusion de l'arbre de décision +ctable <- as.table(results_tree, + nrow = 2, + byrow = TRUE) +# Définir les étiquettes "malade" et "sain" pour les lignes et colonnes +rownames(ctable) <- colnames(ctable) <- c("malade (1)", "sain (0)") + +# Définir les étiquettes "Réel" et "Prédit" +dimnames(ctable) <- list(Réel = rownames(ctable), Prédit = colnames(ctable)) + +fourfoldplot(ctable, + color = c("#CC6666", "#CCCCCC"), + conf.level = 0, + margin = 1, + main = "Matrice de confusion \ndu modèle d'arbre de décision") + +``` + +## Comparaison des performances + +On va comparer la sensibilité, le F1 Score et la valeur prédictive négative car, comme vu précédemment, ce sont des indicateurs importants pour évaluer la capacité du modèle à minimiser les faux négatifs. On affiche aussi le nombre de faux négatifs. Enfin, on compare l'accuracy et l'AUC des différents modèles car c'est un indicateur souvent comparé. + +```{r} +# Calcul de l'AUC RF +pred_rf <- prediction(predicted_probs_test_rf, test$target) +auc_rf <- performance(pred_rf, measure = "auc") +auc_value_rf <- round(unlist(slot(auc_rf, "y.values")), 3) + +# Calcul de l'AUC Tree +pred_tree <- prediction(predicted_probs_test_tree, test$target) +# Calcul de l'AUC +auc_tree <- performance(pred_tree, measure = "auc") +auc_value_tree <- round(unlist(slot(auc_tree, "y.values")), 3) + +# Extraction des valeurs spécifiques +specific_metrics <- c("Sensitivity", "Neg Pred Value", "Precision") + +# Création d'un dataframe avec les valeurs spécifiques +metrics <- c("Sensibilité", "F1 Score", "Valeur prédictive négative", "Accuracy", "Nombre de faux négatifs", "AUC") +models <- c("Regression logistique", "Random Forest", "Arbre de décision") +data <- matrix(NA, nrow = length(models), ncol = length(metrics)) +colnames(data) <- metrics +rownames(data) <- models + +# Insertion des valeurs spécifiques dans le dataframe +data["Regression logistique", ] <- c(unlist(results$byClass[specific_metrics ]), results$overall['Accuracy'], results$table[1, 2], auc_value) +data["Random Forest",] <- c(unlist(results_rf$byClass[specific_metrics]), results_rf$overall['Accuracy'], results_rf$table[1, 2], auc_value_rf) +data["Arbre de décision",] <- c(unlist(results_tree$byClass[specific_metrics]), results_tree$overall['Accuracy'], results_tree$table[1, 2], auc_value_tree) + +# Affichage du tableau +data %>% + kbl(digits = c(3, 3, 3, 3, 1, 2), caption = "Comparaison des performances des différents modèles ") %>% # pour les arrondis et titres + kable_styling() +``` + +On obtient des modèles qui ont des performances assez similaires. + +# Conclusion + +Dans le cadre de ce projet, axé sur l’explication et la prédiction des maladies cardiaques, nous avons développé et évalué un modèle de régression logistique, qui est un modèle facilement interprétable. Les performances de ce modèle se sont avérées prometteuses, démontrant une bonne capacité à expliquer et prédire la présence de maladies cardiaques. Ce modèle a été comparé à des approches alternatives telles que la forêt aléatoire et l'arbre de décision, qui ont montré des performances similaires. + +D'autres modèles, comme le SVM, la classification naïve bayésienne, le KNN, ou le Gradient Boosting, pourraient être explorés pour améliorer la prédiction. + +Des observations paradoxales ont été constatées dans l'application. +Ces résultats inattendus pourraient s'expliquer par des interactions complexes entre les caractéristiques. Introduire des interactions entre les variables pourrait améliorer la capacité prédictive des modèles. +De plus, des biais statistiques dus à la taille limitée de l'échantillon pourraient impacter certaines catégories sous-représentées. +Enfin, certains facteurs externes importants, tels que les comportements, les facteurs environnementaux ou les conditions médicales spécifiques, n'ont pas été inclus dans l'étude. Prendre en compte ces éléments dans de futures recherches pourrait affiner les prédictions des maladies cardiaques. diff --git a/M2/Data Visualisation/Exemple Projet/Application-projet.html b/M2/Data Visualisation/Exemple Projet/Application-projet.html new file mode 100644 index 0000000..bccb66e --- /dev/null +++ b/M2/Data Visualisation/Exemple Projet/Application-projet.html @@ -0,0 +1,6980 @@ + + + + + + + + + + + Projet de visualisation des données + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + +
+
+
+ + + +
+
+ +
+ + +
+ + +

Projet de visualisation des données

+ + + + + + + +
+

1 Introduction

+

Les maladies cardiovasculaires sont une cause majeure de mortalité dans le monde. La détection et l’intervention précoces sont essentielles pour prévenir celles-ci. Certains facteurs sont bien connus comme étant des symptômes de maladies cardiaques.

+

Nous allons analyser un jeu de données pour confirmer ces symptômes, c’est-à-dire vérifier s’ils ont un réel impact sur la santé cardiaque d’un individu. Puis, à partir de ce même jeu de données, nous essaierons de prédire au mieux la variable cible grâce à des modèles statistiques dont nous étudierons les performances.

+

La variable cible est la variable nommée target. Elle représente l’état du patient. Cette variable vaut 0 si l’individu est en bonne santé et vaut 1 si l’individu souffre d’une maladie cardiaque.

+
+
+

2 Observation et préparation des données

+
+

2.1 Observation des données

+

Dans un premier temps, on importe les données et on affiche un aperçu de celles-ci pour les observer rapidement.

+
# Téléchargement du jeu de données 
+df <- read.csv("heart.csv")
+
+# Affichage du jeu de données sous forme de tableau
+rmarkdown::paged_table(df)
+
+ +
+

On affiche aussi la dimension du jeu de données :

+
# Création d'un dataframe contenant les dimensions 
+dim <- data.frame(dim(df)) 
+rownames(dim) <- c("Nombre de lignes","Nombre de colonnes") # Nom des lignes du tableau
+colnames(dim) <- c("") # Nom de la colonne du tableau
+
+# Affichage des dimensions du jeu de données dans un tableau
+dim %>% 
+  kbl(caption = "Dimension du jeu de données") %>%
+  kable_styling()
+ + + + + + + + + + + + + + + + + + +
+Table 2.1: Dimension du jeu de données +
+ +
+Nombre de lignes + +303 +
+Nombre de colonnes + +14 +
+

Le jeu de données comporte 303 lignes, c’est-à-dire 303 observations (individus) et 14 colonnes, c’est-à-dire 14 variables. 13 d’entre elles sont les variables explicatives et la dernière est la variable cible (la variable target comme vu en introduction).

+

On regarde ensuite le type des variables, pour savoir si elles sont quantitatives ou qualitatives ; numériques ou catégorielles, bien qu’on en ait déjà une idée grâce à l’affichage précédent. +Le tableau suivant présente une description récapitulative des variables. On pourra trouver une description plus détaillée de celles-ci dans la notice du projet.

+
# Affichage de la définition des variables dans un tableau 
+kable(caption = "Définition des variables", # Titre du tableau
+  # Création d'un dataframe avec les variables, leur description, leur type et leur support
+  data.frame( 
+        Variable = c(
+          "age", "sex", "cp", "trestbps", "chol", "fbs", "restecg", "thalach", 
+          "exang", "oldpeak", "slope", "ca", "thal"),
+        Description = c(
+          "Age de l'individu",
+          "Sexe de l'individu",
+          "Douleur thoracique",
+          "Pression artérielle au repos en mmHg",
+          "Niveau de cholestérol sérique du patient en mg/dL",
+          "Taux de glycémie à jeun",
+          "Résultats électrocardiographiques au repos",
+          "Fréquence cardiaque maximale atteinte",
+          "Angine de poitrine induite par l'exercice",
+          "Dépression du segment ST induite par l'exercice par rapport au repos",
+          "Pente du segment ST de pointe de l'exercice",
+          "Nombre de vaisseaux principaux colorés par la flourosopie",
+          "Flux sanguin vers le muscle cardiaque "), 
+        Type = c(
+          rep("Entier", 9),
+          "Réel",
+          rep("Entier", 3)),
+        Support = c(
+          "⟦29;77⟧", "{0, 1}", "{0, 1, 2, 3}", "⟦94;200⟧", "⟦126; 564⟧", "{0, 1}",
+          "{0, 1, 2}", "⟦71;202⟧", "{0, 1}", "[0, 6.2]", "{0, 1, 2}", "⟦0;4⟧", "{0, 1, 2, 3}")),
+    booktabs=TRUE) # Mise en page et style du tableau
+ + ++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Table 2.2: Définition des variables
VariableDescriptionTypeSupport
ageAge de l’individuEntier⟦29;77⟧
sexSexe de l’individuEntier{0, 1}
cpDouleur thoraciqueEntier{0, 1, 2, 3}
trestbpsPression artérielle au repos en mmHgEntier⟦94;200⟧
cholNiveau de cholestérol sérique du patient en mg/dLEntier⟦126; 564⟧
fbsTaux de glycémie à jeunEntier{0, 1}
restecgRésultats électrocardiographiques au reposEntier{0, 1, 2}
thalachFréquence cardiaque maximale atteinteEntier⟦71;202⟧
exangAngine de poitrine induite par l’exerciceEntier{0, 1}
oldpeakDépression du segment ST induite par l’exercice par rapport au reposRéel[0, 6.2]
slopePente du segment ST de pointe de l’exerciceEntier{0, 1, 2}
caNombre de vaisseaux principaux colorés par la flourosopieEntier⟦0;4⟧
thalFlux sanguin vers le muscle cardiaqueEntier{0, 1, 2, 3}
+

Les variables sont toutes quantitatives et de type entier sauf oldpeak qui est un réel.
+On a des variables

+
    +
  • catégorielles : sex, cp, fbs, restecg, exang, slope, ca et thal ;
    • +
  • numériques : age, trestbps, chol, thalach et oldpeak.
    • +
+

Remarque : les variables sex, cp, fbs, restecg, exang, slope, ca et thal sont des variables numériques qu’on a classifié de catégorielles. Ces variables sont en fait des variables catégorielles, dont les catégories ont été présentées dans la notice, qui ont été “transformées” en variables numériques pour faciliter leur utilisation.

+
+
+

2.2 Valeurs manquantes

+

On vérifie s’il y a des valeurs manquantes. Si c’est le cas, il faudra y remédier. +Le tableau ci-dessous affiche le nombre de valeur manquante par variable.

+
# Vérification des valeurs manquantes dans le dataframe
+val_manquantes <- sapply(df, function(x) sum(is.na(x)))
+
+# Affichage du nombre de valeurs manquantes par variable dans un tableau
+val_manq <- data.frame(val_manquantes) # Création d'un dataframe avec le nombre de valeurs manquantes
+colnames(val_manq) <- c("Nombre de valeur manquante") # Nom de la colonne du tableau
+val_manq %>% 
+  kbl(caption = "Nombre de valeur manquante par variable") %>% # Titre
+  kable_styling() # Mise en page et style du tableau
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+Table 2.3: Nombre de valeur manquante par variable +
+ +Nombre de valeur manquante +
+age + +0 +
+sex + +0 +
+cp + +0 +
+trestbps + +0 +
+chol + +0 +
+fbs + +0 +
+restecg + +0 +
+thalach + +0 +
+exang + +0 +
+oldpeak + +0 +
+slope + +0 +
+ca + +0 +
+thal + +0 +
+target + +0 +
+

Il n’y a donc pas de valeurs manquantes dans notre jeu de données.

+
+
+

2.3 Distribution des variables

+

Après avoir bien observé le jeu de données, passons à l’étude graphique de la distribution des variables.

+
+

2.3.1 Distribution de la variable cible

+

On étudie dans un premier temps la distribution de la variable cible puis celle des variables explicatives.

+
# Création d'un camembert avec les pourcentages 
+pie(table(df$target), #  Table des effectifs pour la variable target
+    main = "Distribution de la variable cible (target)", # Titre du graphique
+    labels = c("Individus en bonne santé", 
+               "Individus atteints \nd'une maladie cardiaque"), # Légende 
+    col = c("lightblue", "lightcoral"), # Couleur des secteurs
+    cex = 0.8, # Taille de la légende
+    cex.main = 1) # Taille du titre
+
+# Création d'une fonction qui ajoute les valeurs en tant qu'étiquettes au centre (fonction créée par chatgpt)
+text_pie = function(vector, labels = c(), cex = 0.8) {
+     vector = vector / sum(vector) * 2 * pi
+     temp = c()
+     j = 0
+     l = 0
+     for (i in 1:length(vector)) {
+          k = vector[i] / 2        
+          j =  j + l + k
+          l = k
+          text(cos(j) / 2, sin(j) / 2, labels[i], cex = cex)
+     }
+     vector = temp
+}
+
+# Ajout des étiquettes
+percent_values <- round(100 * prop.table(table(df$target)), 1) # Calcul des pourcentages arrondis
+percent <- paste(percent_values, "%") # Ajout du sigle %
+text_pie(percent_values,  
+         c(percent[1], percent[2]), 
+         cex = 0.8)
+

+Le jeu de données est à peu près équilibré : on a environ 50% (54.5%) d’individus malades et 50% (45.5%) d’individus sains dans le jeu de données.

+
+
+

2.3.2 Distribution des variables explicatives

+

On affiche maintenant la distribution des variables catégorielles. Cela va nous permettre d’observer les catégories de variable qui prédominent par rapport aux autres catégories de cette même variable.

+
# Création d'un vecteur avec les variables catégorielles 
+categorical_var <- c("sex", "cp", "fbs", "restecg", "exang", "slope", "ca", "thal")
+
+# Organisation des graphiques en une grille 2x4
+par(mfrow = c(2, 4)) 
+
+# Boucle qui crée un histogramme en baton par chaque variable catégorielle 
+for (var in categorical_var) {
+  barplot(table(df[[var]]), 
+          col = "lightblue", # Couleur de l'histogramme
+          main = var, # Titre de l'histogramme
+          ylab = "Nombre d'individus") # Nom de l'axe des ordonnées 
+}
+

+Les variables catégorielles ne sont pas équilibrées dans le jeu de données :

+
    +
  • sex : il y a deux fois plus d’hommes que de femmes.
    • +
  • cp : beaucoup d’individus ne souffrent pas de douleurs thoraciques (0). Moitié moins souffrent de douleurs modérées (2). Seuls une vingtaine d’individus sur 300 ont de fortes douleurs (3).
    • +
  • fbs : il y a 6 fois moins d’individus ayant un taux de glycémie anormal (> 120 mg/dL) que de personnes ayant un taux de glycémie normal (< 120 mg/dL).
    • +
  • restecg : il y a presque autant d’individus dont les résultats électrocardiographiques sont normaux (0) que de personnes ayant des résultats avec des changements mineurs (1). Très peu (environ 5%) ont des résultats avec des anomalies importantes (2).
    • +
  • exang : il y a 2 fois moins d’individus qui ont des angines de poitrine induites par l’exercice (1) que d’individus qui en n’ont pas (0).
    • +
  • slope : la plupart des individus ont une dépression du segment ST (environ 85%).
    • +
  • ca : plus de la moitié des individus n’ont pas de vaisseaux principaux colorés par la flourosopie.
    • +
  • thal : la plupart des individus ont des défauts modérés (2) voir graves (3) de perfusion myocardique.
    • +
+

On regarde maintenant la distribution des variables numériques. On trace un histogramme et le boxplot asssocié pour chaque variable. Le boxplot nous permet d’avoir la moyenne de chaque variable et d’observer les éventuelles valeurs aberrantes.

+
# Création d'un vecteur avec les variables numériques 
+numerical_var <- c("age", "trestbps", "chol", "thalach", "oldpeak")
+
+# Création de listes pour les histogrammes et les boxplots
+hist_plots <- list()
+box_plots <- list()
+
+# Boucle pour créer les histogrammes et les boxplots pour chaque variable numérique
+for (var in numerical_var) {
+  # Histogramme
+  histo <- ggplot(df, aes(x = .data[[var]])) + 
+    geom_histogram(fill = "lightblue", # Couleur de remplissage
+                   color = "black", # Couleur des contours 
+                   bins = 20, # Nombre de catégories
+                   alpha = 0.7, # Transparence de la couleur de remplissage
+                   linewidth = 0.2) + # Épaisseur des contours
+    labs(title = paste("Distribution de", var), # Titre des histogrammes 
+         x = var, # Nom de l'axe des abscisses 
+         y = "Nombre") + # Nom de l'axe des ordonnées 
+    theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 7, # Taille du titre
+                                  face = "bold"), # Titre en gras
+        axis.title.x = element_text(size = 7), # Taille du nom de l'axe des abscisses
+        axis.title.y = element_text(size = 7))  # Taille du nom de l'axe des ordonnée 
+
+  # Boxplot
+  boxplot <- ggplot(df, aes(x = 1, y = .data[[var]])) +
+    geom_boxplot(linewidth = 0.3, # Épaisseur des contours
+                 outlier.size = 0.2) + # Taille des outliers 
+    labs(title = paste("Boxplot de", var)) + # Titre des boxplots
+    theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 7, # Taille du titre
+                                  face = "bold"), # Titre en gras
+        axis.title.x = element_text(size = 7), # Taille du nom de l'axe des abscisses
+        axis.title.y = element_text(size = 7)) # Taille du nom de l'axe des ordonnée 
+
+  # Stockage des histogrammes et des boxplots dans les listes 
+  hist_plots[[var]] <- histo
+  box_plots[[var]] <- boxplot
+}
+
+# Organisation des graphiques dans une grille 5x2
+grid.arrange(grobs = c(hist_plots, box_plots), 
+             ncol = 5, 
+             top = "Histogrammes et boxplots de chaque variable numérique") # Titre
+

+Grâce à ces graphiques ont peut observer que :

+
    +
  • age : les individus ont entre 29 et 77 ans avec une moyenne d’âge de 54 ans.
    • +
  • trestbps : la pression artérielle au repos de la plupart des individus est supérieure à 120 mmHg, ce qui est un niveau élevé voir à risque.
    • +
  • chol : la plupart des individus ont un niveau de cholestérol sérique entre 200 et 300 mg/dL, ce qui est très élevé.
    • +
  • thalach : la fréquence cardiaque maximale atteinte est entre 100 et 200 bpm, ce qui semble tout à fait normal car celle-ci diminue avec l’âge et dans notre jeu de données, il y a toutes les catégories d’âge.
    • +
  • oldpeak : environ 1/3 des individus ont une dépression ST proche de 0. Les autres individus ont des valeurs plus élevées de dépression ST.
    • +
+

On peut de plus observer certaines valeurs atypiques pour les variables trestbps, chol, thalach et oldpeak. Ce sont des valeurs extrêmes mais pas aberrantes (ce ne sont probablement pas des erreurs de saisie) donc on les garde.

+
+
+
+

2.4 Corrélations entre les variables

+
+

2.4.1 Corrélations entre les variables explicatives

+

On étudie maintenant s’il y a des corrélations entre les variables explicatives. +On affiche pour cela la matrice des corrélations sous forme de graphique corrplot.

+
# Plot des corrélations entre les variables explicatives 
+corrplot(cor(df[, 1:13]), method = "color")
+

+D’après ce graphique, les variables ne semblent pas très corrélées. Nous pouvons cependant porter notre attention sur les 5 couples de variables les plus corrélés.

+

Le premier couple de variables est celui composé des variables oldpeak et slope. Cette corrélation semble tout à fait intuitive car on rappelle que oldpeak représente la dépression ST induite par l’exercice par rapport au repos et que slope est la pente du segment ST de pointe de l’exercice.

+

De plus, les variables age et thalach semblent avoir un coefficient de corrélation plus élevé que les autres variables ; ce qui est cohérent car, comme vu dans la description détaillée des variables (dans la notice), la fréquence cardiaque maximale diminue avec l’âge (on le vérifie plus bas).

+

Les variables exang et cp ont aussi un coefficient de corrélation de couleur foncée. Là encore, cette corrélation n’est pas aberrante car ces variables concernent toutes les 2 des douleurs thoraciques.

+

Le couple de variables slope/thalach est lui aussi un couple qui semble avoir un coefficient de corrélation plus élevé que la plupart des autres couples de variables.

+

Enfin le couple thalach et exang est le dernier couple qui se détache des autres par la valeur, plutôt élevée, de son coefficient de corrélation.

+

On affiche la matrice des corrélations pour voir plus précisément les coefficients de corrélations. +Les coefficients dont la valeur absolue est supérieure à 0.35 sont affichés en rouge.

+
# Calcul des corrélations entre les variables explicatives
+correlation_matrix <- cor(df[, 1:13])
+
+# Fonction pour formater les coefficients de corrélation (pour colorer les coefficients dont la valeur absolue est supérieure à 0.35 et arrondir les coefficients au centième près)
+format_correlation <- function(x) {
+  if (abs(x) > 0.35) {
+    paste0("<span style='color: red; font-weight: bold;'>", round(x, 2), "</span>")
+  } else {
+    round(x, 2)
+  }
+}
+
+# Application de la fonction de formatage à la matrice de corrélation
+formatted_correlation_matrix <- apply(correlation_matrix, 1, function(row) {
+  sapply(row, format_correlation)
+})
+
+# Affichage de la matrice de corrélation avec la mise en forme personnalisée
+kable(formatted_correlation_matrix, 
+      format = "html", 
+      escape = FALSE) %>%
+  kable_styling()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ +age + +sex + +cp + +trestbps + +chol + +fbs + +restecg + +thalach + +exang + +oldpeak + +slope + +ca + +thal +
+age + +1 + +-0.1 + +-0.07 + +0.28 + +0.21 + +0.12 + +-0.12 + +-0.4 + +0.1 + +0.21 + +-0.17 + +0.28 + +0.07 +
+sex + +-0.1 + +1 + +-0.05 + +-0.06 + +-0.2 + +0.05 + +-0.06 + +-0.04 + +0.14 + +0.1 + +-0.03 + +0.12 + +0.21 +
+cp + +-0.07 + +-0.05 + +1 + +0.05 + +-0.08 + +0.09 + +0.04 + +0.3 + +-0.39 + +-0.15 + +0.12 + +-0.18 + +-0.16 +
+trestbps + +0.28 + +-0.06 + +0.05 + +1 + +0.12 + +0.18 + +-0.11 + +-0.05 + +0.07 + +0.19 + +-0.12 + +0.1 + +0.06 +
+chol + +0.21 + +-0.2 + +-0.08 + +0.12 + +1 + +0.01 + +-0.15 + +-0.01 + +0.07 + +0.05 + +0 + +0.07 + +0.1 +
+fbs + +0.12 + +0.05 + +0.09 + +0.18 + +0.01 + +1 + +-0.08 + +-0.01 + +0.03 + +0.01 + +-0.06 + +0.14 + +-0.03 +
+restecg + +-0.12 + +-0.06 + +0.04 + +-0.11 + +-0.15 + +-0.08 + +1 + +0.04 + +-0.07 + +-0.06 + +0.09 + +-0.07 + +-0.01 +
+thalach + +-0.4 + +-0.04 + +0.3 + +-0.05 + +-0.01 + +-0.01 + +0.04 + +1 + +-0.38 + +-0.34 + +0.39 + +-0.21 + +-0.1 +
+exang + +0.1 + +0.14 + +-0.39 + +0.07 + +0.07 + +0.03 + +-0.07 + +-0.38 + +1 + +0.29 + +-0.26 + +0.12 + +0.21 +
+oldpeak + +0.21 + +0.1 + +-0.15 + +0.19 + +0.05 + +0.01 + +-0.06 + +-0.34 + +0.29 + +1 + +-0.58 + +0.22 + +0.21 +
+slope + +-0.17 + +-0.03 + +0.12 + +-0.12 + +0 + +-0.06 + +0.09 + +0.39 + +-0.26 + +-0.58 + +1 + +-0.08 + +-0.1 +
+ca + +0.28 + +0.12 + +-0.18 + +0.1 + +0.07 + +0.14 + +-0.07 + +-0.21 + +0.12 + +0.22 + +-0.08 + +1 + +0.15 +
+thal + +0.07 + +0.21 + +-0.16 + +0.06 + +0.1 + +-0.03 + +-0.01 + +-0.1 + +0.21 + +0.21 + +-0.1 + +0.15 + +1 +
+

Naturellement, la matrice nous mène à la même conclusion que le graphique : il ne semble pas y avoir de corrélations significatives entre les variables. +Les variables les plus corrélées sont

+
    +
  • oldpeak et slope avec un coefficient de corrélation de -0.58,
    +
    • +
  • age et thalach avec un coefficient de corrélation de -0.4,
    • +
  • exang et cp avec un coefficient de corrélation de -0.39,
    • +
  • slope et thalach avec un coefficient de corrélation de -0.39,
    • +
  • thalach et exang avec un coefficient de corrélation de -0.38.
    • +
+

On étudie graphiquement la relation entre les variables des couples les plus corrélés, pour avoir un représentation plus visuelle de celle-ci.

+

On trace ci-dessous la distribution de la la dépression ST induite par l’exercice (oldpeak) pour chaque niveau de la variable slope.

+
# Distribution de oldpeak en fonction de slope
+ggplot(df, aes(x = oldpeak, fill = factor(slope))) + 
+      geom_density(alpha = 0.5, # Transparence de la couleur de remplissage
+                   linewidth = 0.3) + # Épaissseur des contours 
+  
+      # Noms des axes et titre 
+      labs(title = "Graphique de la distribution de la la dépression ST induite par l'exercice (oldpeak) \n en fonction de la pente du segment ST de pointe de l'exercice (slope)", # Titre
+       x = "Changements dans l'électrocardiogramme (oldpeak)", # Nom de l'axe des abscisses 
+       y = "Nombre de patients", # Nom de l'axe des ordonnées 
+       fill = "Pente du segment \nST de pointe (slope)")  + # Titre de la légende
+  
+      # Couleur et texte de la légende
+      scale_fill_manual(values = c("0" = "#FFDB6D", 
+                                   "1" = "#4E84C4", 
+                                   "2" = "#C3D7A4")) +
+  
+      # Mise en forme du titre et de la légende 
+      theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 10, # Taille du titre
+                                  face = "bold"), # Titre en gras
+        legend.title = element_text(size = 9), # Taille du titre de la légende
+        legend.text = element_text(size = 8)) # Taille du texte de la légende
+

+

On remarque que lorsque la dépression ST induite par l’exercice (oldpeak) est proche de 0, la pente du segment ST de pointe de l’exercice (slope) a tendance à être de 2. Inversement, il n’y a pas de dépression du segment ST observée par rapport à la ligne de base lorsqu’il y a beaucoup de changements dans l’électrocardiogramme (ECG) qui se produisent en réponse à l’exercice physique.

+

On confirme donc qu’à priori, il faut enlever une des 2 variables. Dans la suite, on va faire d’autres tests pour confirmer cette décision.

+

On regarde ensuite plus précisement la relation entre les variables age et thalach.

+
# Diagramme de dispersion entre age et thalach
+ggplot(df, aes(x = age, y = thalach)) +
+  
+  # Ajout de la courbe de tendance
+  geom_smooth(method = "lm", 
+              formula = y ~ x, # Lien linéaire (polynôme de degré 1)
+              se = FALSE, # Ne pas afficher l'intervalle de confiance
+              colour="black", # Couleur de la droite
+              linewidth = 0.3) + # Épaisseur du trait
+  
+  # Ajout des points
+  geom_point(size = 1, # Taille des points
+             color = "#4E84C4") + # Couleur des points
+  
+  # Nom des axes et titre
+  labs(title = "Graphique de la fréquence cardiaque maximale atteinte (thalach) \n en fonction de l'âge (age)", # Titre du graphique
+       x = "Âge (age)", # Nom de l'axe des abscisses 
+       y = "Fréquence cardiaque maximale atteinte (thalach)") + # Nom de l'axe des ordonnées
+  
+  # Mise en forme du titre et de la légende 
+  theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 12, # Taille du titre
+                                  face = "bold")) # Titre en gras
+

+Cette figure fait clairement apparaître une relation décroissante et linéaire entre les deux variables. On en déduit que la fréquence cardiaque maximale diminue avec l’âge.

+

On étudie maintenant graphiquement le lien entre la douleur thoracique (cp) et la présence / l’absence d’angine de poitrine induite par l’exercice (exang).

+
# Histogramme entre cp et exang
+ggplot(df, aes(x = cp, fill = factor(exang))) + 
+      geom_bar(color = "black", # Couleur des contours
+               linewidth = 0.3, # Épaisseur des contours
+               alpha = 0.5) + # Transparence de la couleur de remplissage
+  
+      # Nom des axes et titre
+      labs(title = "Graphique de la distribution de la douleur thoracique (cp) en fonction de \n la présence / l'absence d'angine de poitrine induite par l'exercice (exang)", # Titre du graphique 
+           x = "Douleur thoracique (cp)", # Nom de l'axe des abscisses 
+           y = "Nombre d'individus", # Nom de l'axe des ordonnées
+           fill = "Angine thoracique \ninduite par l'exercice") + # Titre de la légende
+  
+      # Texte de la légende et couleur
+      scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4"), 
+                        labels = c("Non", "Oui")) +
+  
+      # Mise en forme du titre et de la légende 
+      theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 10, # Taille du titre
+                                  face = "bold"), # Titre en gras
+        legend.title = element_text(size = 9), # Taille du titre de la légende
+        legend.text = element_text(size = 8)) # Taille du texte de la légende
+

+On remarque que les individus n’ayant pas de douleur thoracique ont des angines thoraciques induites par l’exercice et inversement : les individus ayant des douleurs thoraciques n’ont pas d’angine thoracique induite par l’exercice.

+

La relation entre la fréquence cardiaque maximale atteinte (thalach) et la pente du segment ST de pointe de l’exercice (slope) est représentée dans le graphique ci-dessous.

+
# Distribution de thalach en fonction de slope
+ggplot(df, aes(x = thalach, fill = factor(slope))) + 
+      geom_histogram(color = "black", # Couleur des contours 
+                     linewidth = 0.3, # Épaisseur des contours
+                     bins = 25, # Nombre de catégories
+                     alpha = 0.5) + # Transparence de la couleur de remplsisage 
+  
+      # Nom des axes et titre
+      labs(title = "Graphique de la distribution de la fréquence cardiaque maximale atteinte (thalach) \n en fonction de la pente du segment ST de pointe de l'exercice (slope)", # Titre du graphique
+           x ="Fréquence cardiaque maximale atteinte (thalach)", # Nom de l'axe des abscisses
+           y = "Nombre d'individus", # Nom de l'axe des ordonnées
+           fill = "Pente du segment \nST (slope)") + # Titre de la légende
+  
+      # Couleur de la légende 
+      scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4", "2" = "#C3D7A4")) +
+  
+      # Mise en forme du titre et de la légende 
+      theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 10, # Taille du titre
+                                  face = "bold"), # Titre en gras
+        legend.title = element_text(size = 9), # Taille du titre de la légende
+        legend.text = element_text(size = 8)) # Taille du texte de la légende
+

+Ce graphique montre que plus la fréquence cardiaque maximale atteinte est élevée (entre 160 et 200), plus la pente du segment ST est importante (niveau 2).

+

Enfin, on trace la distribution de la fréquence cardiaque maximale atteinte (thalach) en fonction de la présence / l’absence d’angine de poitrine induite par l’exercice (exang) .

+
# Distribution de thalach en fonction de exang
+ggplot(df, aes(x = thalach, fill = factor(exang))) + 
+      geom_histogram(color = "black", # Couleur des contours 
+                     linewidth = 0.3, # Épaisseur des contours
+                     bins = 30, # Nombre de catégories
+                     alpha = 0.5) + # Transparence de la couleur de remplsisage 
+  
+      # Noms des axes et titre 
+      labs(title = "Graphique de la distribution de la fréquence cardiaque maximale atteinte (thalach) \n en fonction de la présence / l'absence d'angine de poitrine induite par l'exercice (exang) ", # Titre du graphique
+           x = "Fréquence cardiaque maximale atteinte (thalach)", # Nom de l'axe des abscisses
+           y = "Nombre d'individus", # Nom de l'axe des orodnnées
+           fill = "Angine thoracique induite \npar l'exercice (exang)") + # Titre de la légende
+  
+      # Couleur et texte de la légende
+      scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4"), 
+                      labels = c("Non", "Oui")) +
+  
+      # Mise en forme du titre et de la légende 
+      theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 10, # Taille du titre
+                                  face = "bold"), # Titre en gras
+        legend.title = element_text(size = 9), # Taille du titre de la légende
+        legend.text = element_text(size = 8)) # Taille du texte de la légende
+

+On observe que les individus ayant une fréquence cardiaque maximale atteinte plus faible (entre 50 et 150) ont, en général, des angines induites par l’exercice et inversement.

+

On décide pour l’instant de garder les variables de chaque couple. On fera une étude plus approfondie par la suite.

+
+
+

2.4.2 Relation entre les variables explicatives et la variable cible

+
+

2.4.2.1 Aperçu global

+

Pour avoir un aperçu global et visuel des variables qui ont un impact sur la variable cible, on affiche un graphe des coefficients de corrélation entre la variable cible et chacune des autres variables.

+
# Matrice des corrélations de toutes les variables 
+cor_matrix <- cor(df[, 1:14])
+
+# Récupérer les coefficients de corrélation entre la variable cible et les autres variables
+target_correlations <- cor_matrix[, 14] # La 14ème colonne correspond à la variable cible
+
+# Création d'un dataframe avec les noms des variables et leurs corrélations avec la variable cible
+correlation_data <- data.frame(variable = names(target_correlations)[-14], 
+                               correlation = unlist(target_correlations[-14]))
+
+# Création d'un graphique (barplot) des coefficients de corrélation
+ggplot(correlation_data, aes(x = variable, y = correlation)) +
+  geom_bar(stat = "identity", # Hauteur des barres = coefficient de corrélation
+           fill = "#4E84C4", # Couleur de remplissage 
+           alpha = 0.7, # Transparence de la couleur de remplissage
+           color = "black", # Couleur des contours 
+           linewidth = 0.3) + # Épaisseur des contours
+  labs(title = 'Corrélations entre la variable cible et les autres variables', # Titre du graphique
+       x = 'Variables', # Nom de l'axe des abscisses
+       y = 'Coefficient de corrélation') + # Nom de l'axe des ordonnées
+  
+  theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 12, # Taille du titre
+                                  face = "bold")) # Titre en gras
+

+Les variables qui semblent avoir un impact (un coefficient de corrélation proche de 1 en valeur absolue) sur la target sont : ca, cp, exang, oldpeak, sex, slope, thal et thalach.

+
+
+

2.4.2.2 Relation entre les variables catégorielles et la variable cible

+

On va maintenant regarder plus précisemment les relations entre les variables explicatives et la target. On commence par étudier les variables catégorielles puis on étudiera les variables numériques.

+

Grâce à l’aperçu global, on a vu que seules les variables catégorielles sex, cp, exang, slope, ca et thal avaient potentiellement une relation avec la variable cible. On affiche alors un histogramme de chacune de ces variables en fonction de la target pour voir graphiquement cette relation.

+
# Création d'un vecteur avec les variables catégorielles qui semblent avoir un impact sur la variable cible
+categorical_var_impact <- c("sex", "cp", "exang", "slope", "ca", "thal")
+
+# Création d'une liste pour stocker les graphiques
+hist <- list()
+
+# Boucle qui crée un histogramme pour chaque variable du vecteur categorical_var_impactet et qui l'affiche
+for (var in categorical_var_impact) {
+    h <- ggplot(df, aes(x = .data[[var]], fill = factor(target))) + 
+      geom_bar(color = "black", # Couleur des contours 
+               linewidth = 0.3, # Épaisseur des contours
+               position = "fill", # Représenter des proportions
+               alpha = 0.7) + # Transparence de la couleur de remplissage
+      
+      # Noms des axes et titre
+      labs(title = paste("Distribution de", var, "\nen fonction de la variable cible"), # Titre des graphiques 
+           x = var, # Nom de l'axe des abscisses 
+           y = "Nombre d'individus", # Nom de l'axe des ordonnées 
+           fill = "Légende") + # Titre de la légende
+      
+      # Couleur et texte de la légende 
+      scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4"), 
+                      labels = c("Sain", "Malade")) +
+      
+      # Mise en forme du titre et de la légende 
+      theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 8, # Taille du titre
+                                  face = "bold"), # Titre en gras
+        legend.title = element_text(size = 7), # Taille du titre de la légende
+        legend.text = element_text(size = 7), # Taille du texte de la légende
+        legend.key.size = unit(0.3, "cm"), # Taille des carrés de couleur de la légende
+        axis.title.x = element_text(size = 7), # Taille du nom de l'axe des abscisses
+        axis.title.y = element_text(size = 7))  # Taille du nom de l'axe des ordonnée 
+    
+    hist[[var]] <- h
+}
+
+# Organisation des graphiques dans une grille à 3 colonnes
+grid.arrange(grobs = hist, ncol = 3)
+

+

On peut observer sur ces graphiques que les variables catégorielles selectionnées semblent effectivement avoir un impact sur la variable cible. En effet, les individus qui semblent plus enclin à avoir une maladie cardiaque sont :

+
    +
  • les femmes,
    • +
  • les individus ayant des douleurs thoraciques (de niveau 1, 2 ou 3),
    • +
  • les individus ayant des anomalies mineures dans l’ECG (étonnant),
    • +
  • les individus n’ayant pas d’angine de poitrine induite par l’exercice (là encore c’est étonnant),
    • +
  • les individus ayant des changements dans le segment ST graves (de niveau 2),
    • +
  • les individus n’ayant pas de vaisseaux principaux colorés par la flourosopie,
    • +
  • les individus ayant des défauts modérés de perfusion myocardique (niveau 2).
    • +
+

Remarque : des précisions sur les observations et conclusions étonnantes seront apportées dans la conclusion.

+
+
+

2.4.2.3 Relation entre les variables numériques et la variable cible

+

Enfin, grâce à l’aperçu global, on a vu que seules les variables numériques thalach et oldpeak avaient potentiellement une relation avec la variable cible. On affiche alors la densité de chacune de ces variables en fonction de la target pour voir graphiquement cette relation.

+
# Graphique de target en fonction de thalach
+
+ggplot(df, aes(x = thalach, fill = factor(target))) + 
+      geom_histogram(color = "black", # Couleur des contours 
+                     linewidth = 0.3, # Épaisseur des contours
+                     bins = 20, # 20 classes pour plus de lisibilité
+                     alpha = 0.5) + # Transparence de la couleur de remplissage 
+  
+      # Noms des axes et titre
+      labs(title = "Graphique de la distribution de la fréquence cardiaque maximale atteinte (thalach) \n en fonction de la variable cible (target)", # Titre du grahique 
+           x = "Fréquence cardiaque maximale atteinte (thalach)", # Nom de l'axe des abscisses 
+           y = "Nombre d'individus", # Nom de l'axe des ordonnées 
+           fill = "Légende") + # Nom du titre de la légende 
+  
+      # Couleur et texte de la légende 
+      scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4"), 
+                        labels = c("Sain", "Malade")) +
+  
+      # Mise en forme du titre et de la légende 
+      theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 12, # Taille du titre
+                                  face = "bold"), # Titre en gras
+        legend.title = element_text(size = 10), # Taille du titre de la légende
+        legend.text = element_text(size = 10), # Taille du texte de la légende
+        legend.key.size = unit(0.3, "cm"), # Taille des carrés de couleur de la légende
+        axis.title.x = element_text(size = 10), # Taille du nom de l'axe des abscisses
+        axis.title.y = element_text(size = 10))  # Taille du nom de l'axe des ordonnée
+

+

On peut observer sur ce graphique que la variable thalach semble effectivement avoir un impact sur la variable cible. En effet, les individus qui semblent plus enclin à avoir une maladie cardiaque sont les individus ayant une fréquence cardiaque maximale entre 150 et 200 bpm.

+

Enfin, on trace ci-dessous la distribution de la la dépression ST induite par l’exercice (oldpeak) en fonction de la variable cible.

+
# Distribution de oldpeak en fonction de target
+ggplot(df, aes(x = oldpeak, fill = factor(target))) + 
+      geom_density(alpha = 0.5, # Transparence de la couleur de remplissage 
+                   linewidth = 0.3) + # Épaisseur des contours 
+  
+      # Noms des axes et titre 
+      labs(title = "Graphique de la distribution de la la dépression ST induite par l'exercice (oldpeak) \n en fonction de la variable cible (target)", # Titre du graphique
+       x = "Changements dans l'électrocardiogramme (oldpeak)", # Nom de l'axe des abscisses 
+       y = "Nombre d'individus", # Nom de l'axe des ordonnées 
+       fill = "Légende")  + # Nom du titre de la légende 
+  
+      # Couleur et texte de la légende
+      scale_fill_manual(values = c("0" = "#FFDB6D", "1" = "#4E84C4"), 
+                      labels = c("Sain", "Malade")) +
+  
+      # Mise en forme du titre et de la légende 
+      theme(plot.title = element_text(hjust = 0.5, # Position du titre 
+                                  size = 12, # Taille du titre
+                                  face = "bold"), # Titre en gras
+        legend.title = element_text(size = 9), # Taille du titre de la légende
+        legend.text = element_text(size = 8)) # Taille du texte de la légende
+

+On conclut que les individus n’ayant pas de changement dans l’electrocardiogramme ont plus tendance à avoir une maladie cardiaque et inversement.

+
+
+
+
+

2.5 Sampling : division du jeu de données

+

On divise le jeu de données en 2 ensembles :

+
    +
  • un ensemble d’entrainement qui va permettre de construire le modèle,
    • +
  • un ensemble de test pour évaluer les performances du modèle.
    • +
+

On choisit une porportion souvent utilisée pour un jeu de données de petite taille : 80% des données sont utilisées pour former l’ensemble d’entrainement et 20% pour celui de test. +On affiche ci-dessous les dimensions des deux ensembles et celles du jeu de données de base.

+
# Division des données en ensembles d'entraînement et de test
+set.seed(123) # pour la reproductibilité
+split <- rsample::initial_split(df, prop = 0.8)
+train <- rsample::training(split)
+test <- rsample::testing(split)
+
+# Création d'un dataframe avec les dimensions des différents ensembles
+dimensions <- data.frame(dim(df), dim(train), dim(test))
+rownames(dimensions) <- row_names <- c("Nombre de lignes",
+                                       "Nombre de colonnes") # nom des lignes 
+colnames(dimensions) <- c("Dimension du jeu de données initial", 
+                          "Dimensions de l'ensemble d'entrainement", 
+                          "Dimensions de l'ensemble de test") # nom des colonnes 
+
+# Affichage des dimensions des ensembles d'entraînement et de test
+dimensions %>% 
+  kbl() %>%
+  kable_styling()
+ + + + + + + + + + + + + + + + + + + + + + + +
+ +Dimension du jeu de données initial + +Dimensions de l’ensemble d’entrainement + +Dimensions de l’ensemble de test +
+Nombre de lignes + +303 + +242 + +61 +
+Nombre de colonnes + +14 + +14 + +14 +
+

On regarde la distribution de la variable cible dans les 2 ensembles pour s’assurer que la répartition des classes dans l’ensemble de test et d’entrainement est représentative de l’ensemble complet.

+
# Organisation des graphiques en une grille 1x2
+par(mfrow = c(1, 2)) 
+
+# Création du camembert pour l'ensemble train avec les pourcentages et légendes
+pie(table(train$target), #  Table des effectifs pour la variable target
+    main = "Distribution de la variable cible (target) \nsur l'ensemble d'entrainement", # Titre
+    labels = c("Individus en bonne santé", 
+               "Individus atteints \nd'une maladie cardiaque"), # Légende 
+    col = c("lightblue", "lightcoral"), # Définition de la couleur des secteurs
+    cex = 0.8, # Ajustement de la taille de la légende
+    cex.main = 1) # Ajustement de la taille du titre
+
+# Ajout des étiquettes
+percent_values_train <- round(100 * prop.table(table(train$target)), 1) # Calcul des pourcentages arrondis
+percent_train <- paste(percent_values_train, "%") # Ajout su sigle "%"
+text_pie(percent_values_train,  c(percent_train[1], percent_train[2]), cex=0.8) # Utilisation de la fonction text_pie définie plus haut 
+
+
+# Création du camembert pour l'ensemble test avec les pourcentages et légendes
+pie(table(test$target), #  Table des effectifs pour la variable target
+    main = "Distribution de la variable cible (target) \nsur l'ensemble de test", # Titre
+    labels = c("Individus en bonne santé", 
+               "Individus atteints \nd'une maladie cardiaque"), # Légende 
+    col = c("lightblue", "lightcoral"), # Définition de la couleur des secteurs
+    cex = 0.8, # Ajustement de la taille de la légende
+    cex.main = 1) # Ajustement de la taille du titre
+
+# Ajout des étiquettes
+percent_values_test <- round(100 * prop.table(table(test$target)), 1) # Calcul des pourcentages arrondis
+percent_test <- paste(percent_values_test, "%") # Ajout su sigle "%"
+text_pie(percent_values_test,  c(percent_test[1], percent_test[2]), cex=0.8) # Utilisation de la fonction text_pie définie plus haut
+

+
+
+
+

3 Fitting : création et entraînement du modèle

+
+

3.1 Déclaration du modèle

+

On implémente maintenant une régression logistique à l’aide de la fonction “glm”. On pourra trouver une explication du choix du modèle dans la notice. +On déclare aussi le modèle réduit à l’intercept car il sera utile dans la suite, pour des comparaisons de perfromance par exemple.

+
# Déclaration du modèle de regression logistique
+mod <- glm(target ~ ., data = train, family = binomial)
+
+# Déclaration du modèle réduit à l'intercept (utile dans la suite)
+mod_int <- glm(target~1, data = train, family = "binomial")
+

On affiche les statistqiues récapitulatives du modèle de regression logistique construit dans la sous-section section précédente. Les p-valeurs inférieures à 0.05 sont affichées en rouge car cela signifie que la variable associée à probablement de l’importance dans le modèle.

+
# Création d'un tableau qui contient les coefficients du modèle 
+tab_stat <- tidy(mod, conf.int = TRUE)
+
+# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05
+tab_stat <- tab_stat %>%
+  mutate(p.value = ifelse(p.value < 0.05, 
+                                    cell_spec(round(p.value, 3), "html", color = "red"), 
+                                    cell_spec(round(p.value, 3), "html")))
+
+# Affichage du tableau avec mise en forme
+tab_stat %>%
+  kbl(digits = c(0, 2, 2, 2, 5, 3, 3), # Arrondis
+      escape = FALSE, 
+      caption = "Statistiques récapitulatives du modèle de regression logistique") %>%
+  kable_styling()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+Table 3.1: Statistiques récapitulatives du modèle de regression logistique +
+term + +estimate + +std.error + +statistic + +p.value + +conf.low + +conf.high +
+(Intercept) + +4.56 + +2.79 + +1.63 + +0.103 + +-0.824 + +10.204 +
+age + +-0.01 + +0.03 + +-0.53 + +0.599 + +-0.063 + +0.036 +
+sex + +-1.54 + +0.51 + +-3.02 + +0.003 + +-2.590 + +-0.578 +
+cp + +0.77 + +0.20 + +3.94 + +0 + +0.398 + +1.169 +
+trestbps + +-0.02 + +0.01 + +-1.63 + +0.103 + +-0.040 + +0.003 +
+chol + +0.00 + +0.00 + +-0.76 + +0.446 + +-0.011 + +0.005 +
+fbs + +0.13 + +0.57 + +0.23 + +0.815 + +-0.962 + +1.270 +
+restecg + +0.66 + +0.39 + +1.69 + +0.09 + +-0.101 + +1.434 +
+thalach + +0.02 + +0.01 + +1.73 + +0.084 + +-0.002 + +0.043 +
+exang + +-1.08 + +0.45 + +-2.40 + +0.016 + +-1.969 + +-0.199 +
+oldpeak + +-0.68 + +0.24 + +-2.82 + +0.005 + +-1.176 + +-0.224 +
+slope + +0.26 + +0.41 + +0.64 + +0.522 + +-0.563 + +1.065 +
+ca + +-0.65 + +0.20 + +-3.21 + +0.001 + +-1.050 + +-0.254 +
+thal + +-1.02 + +0.33 + +-3.08 + +0.002 + +-1.684 + +-0.382 +
+

D’après cet affichage, les variables qui ont une importance sont : sex, cp, exang, oldpeak, ca et thal.

+

On peut noter que la déviance du modèle est inférieure à celle du modèle réduit à l’intercept :

+
# Création d'un dataframe avec les déviances
+deviance <- data.frame(mod$null.deviance, mod$deviance)
+colnames(deviance) <- c("Déviance du modèle réduit à l'intercept", 
+                          "Déviance du modèle créé") # Nom des colonnes 
+
+# Affichage des dimensions des ensembles d'entraînement et de test
+deviance %>% 
+  kbl(digits = c(2, 2), # Arrondis
+      caption = "Déviance des modèles") %>% 
+  kable_styling()
+ + + + + + + + + + + + + + +
+Table 3.2: Déviance des modèles +
+Déviance du modèle réduit à l’intercept + +Déviance du modèle créé +
+333.48 + +173.07 +
+

Le modèle est donc meilleur que le modèle réduit à l’inetrcept.

+
+
+

3.2 Sélection de modèles : tests et méthodes pas à pas

+

On utilise le test Anova pour enlever les variables qui n’apportent pas d’informations complémentaires sur la variable à prédire. On affiche le résultat du test. Les p-valeurs inférieurs à 0.05 sont affichées en rouge.

+
# Test Anova 
+Ano <- Anova(mod, 
+             type = 3, #  type III de l'Anova 
+             test.statistic = "LR") # test de rapport de vraisemblance
+
+# Création d'un tableau qui contient les coefficients du modèle 
+tab_stat <- tidy(Ano, conf.int = TRUE)
+
+# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05
+tab_stat <- tab_stat %>%
+  mutate(p.value = ifelse(p.value < 0.05, 
+                                    cell_spec(round(p.value, 3), "html", color = "red"), 
+                                    cell_spec(round(p.value, 3), "html")))
+
+# Affichage du tableau avec mise en forme
+tab_stat %>%
+  kbl(digits = c(0, 3, 5), 
+      escape = FALSE, 
+      caption = "Résultats du test Anva") %>%
+  kable_styling()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+Table 3.3: Résultats du test Anva +
+term + +statistic + +df + +p.value +
+age + +0.276 + +1 + +0.599 +
+sex + +10.168 + +1 + +0.001 +
+cp + +17.275 + +1 + +0 +
+trestbps + +2.732 + +1 + +0.098 +
+chol + +0.571 + +1 + +0.45 +
+fbs + +0.055 + +1 + +0.814 +
+restecg + +2.892 + +1 + +0.089 +
+thalach + +3.069 + +1 + +0.08 +
+exang + +5.758 + +1 + +0.016 +
+oldpeak + +8.806 + +1 + +0.003 +
+slope + +0.405 + +1 + +0.524 +
+ca + +10.405 + +1 + +0.001 +
+thal + +9.976 + +1 + +0.002 +
+

Seules les variables sex, cp, exang, oldpeak, ca et thal semblent avoir un impact sur la variable cible, ce qu’on avait déjà observé dans la partie “graphiques”. On avait aussi remarqué que les variables exang, slope et age avaient une importance. On va continuer notre étude pour savoir si on les supprime ou non.

+

On effectue maintenent des méthodes pas à pas pour sélectionner un modèle à l’aide d’une procédure basée sur la minimisation du critère AIC (Akaike Information Criterion).

+

On réalise d’abord une méthode backward et on affiche les statistqiues descriptives du modèle obtenu. Les p_valeurs inférieures à 0.05 sont affichées en rouge.

+
# Méthode backward
+mod_back <- step(mod, direction = "backward", trace = FALSE)
+
+# Création d'un tableau qui contient les coefficients du modèle 
+tab_stat <- tidy(mod_back, conf.int = TRUE)
+
+# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05
+tab_stat <- tab_stat %>%
+  mutate(p.value = ifelse(p.value < 0.05, 
+                                    cell_spec(round(p.value, 3), "html", color = "red"), 
+                                    cell_spec(round(p.value, 3), "html")))
+
+# Affichage du tableau avec mise en forme
+tab_stat %>%
+  kbl(digits = c(0, 3, 3, 3, 4), # Arrondis
+      escape = FALSE, 
+      caption = "Statistiques récapitulatives du modèle construit par la méthode backward") %>%
+  kable_styling()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+Table 3.4: Statistiques récapitulatives du modèle construit par la méthode backward +
+term + +estimate + +std.error + +statistic + +p.value + +conf.low + +conf.high +
+(Intercept) + +3.196 + +2.194 + +1.457 + +0.145 + +-1 + +7.644 +
+sex + +-1.326 + +0.461 + +-2.878 + +0.004 + +-2 + +-0.452 +
+cp + +0.771 + +0.194 + +3.967 + +0 + +0 + +1.168 +
+trestbps + +-0.019 + +0.011 + +-1.837 + +0.066 + +0 + +0.001 +
+restecg + +0.766 + +0.379 + +2.023 + +0.043 + +0 + +1.522 +
+thalach + +0.022 + +0.010 + +2.164 + +0.03 + +0 + +0.043 +
+exang + +-1.067 + +0.442 + +-2.415 + +0.016 + +-2 + +-0.201 +
+oldpeak + +-0.762 + +0.211 + +-3.613 + +0 + +-1 + +-0.368 +
+ca + +-0.629 + +0.192 + +-3.282 + +0.001 + +-1 + +-0.258 +
+thal + +-1.057 + +0.329 + +-3.211 + +0.001 + +-2 + +-0.427 +
+

Les variables importantes sont celles du test Anova ainsi que restecg et thalach.

+

On refait alors un test Anova sur le modèle obtenu par la méthode backward pour vérifier si les variables restecg et thalach sont réellement importantes.

+

Le test donne le même résultat que la méthode backward.

+

On réalise maintenant une méthode forward de sélection de modèle et on affiche là encore, les statistqiues descriptives du modèle obtenu. Les p_valeurs inférieures à 0.05 sont affichées en rouge.

+
# Méthode forward
+mod_for <- step(mod_int, 
+                target ~ age + sex + trestbps + chol + fbs + restecg + thalach
+                 + exang + oldpeak + slope + ca + thal, 
+                data = train,
+                trace = FALSE, 
+                direction = c("forward"))
+
+# Création d'un tableau qui contient les coefficients du modèle 
+tab_stat <- tidy(mod_for, conf.int = TRUE)
+
+# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05
+tab_stat <- tab_stat %>%
+  mutate(p.value = ifelse(p.value < 0.05, 
+                                    cell_spec(round(p.value, 3), "html", color = "red"), 
+                                    cell_spec(round(p.value, 3), "html")))
+
+# Affichage du tableau avec mise en forme
+tab_stat %>%
+  kbl(digits = c(0, 3, 3, 3, 4), # Arrondis
+      escape = FALSE,  
+      caption = "Statistiques récapitulatives du modèle construit par la méthode forward") %>%
+  kable_styling()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+Table 3.5: Statistiques récapitulatives du modèle construit par la méthode forward +
+term + +estimate + +std.error + +statistic + +p.value + +conf.low + +conf.high +
+(Intercept) + +0.706 + +1.607 + +0.439 + +0.661 + +-2 + +3.875 +
+oldpeak + +-0.634 + +0.187 + +-3.398 + +0.001 + +-1 + +-0.283 +
+exang + +-1.531 + +0.404 + +-3.789 + +0 + +-2 + +-0.753 +
+ca + +-0.668 + +0.179 + +-3.742 + +0 + +-1 + +-0.326 +
+thal + +-1.020 + +0.305 + +-3.345 + +0.001 + +-2 + +-0.434 +
+thalach + +0.026 + +0.009 + +2.740 + +0.006 + +0 + +0.044 +
+sex + +-1.029 + +0.411 + +-2.504 + +0.012 + +-2 + +-0.238 +
+restecg + +0.784 + +0.350 + +2.240 + +0.025 + +0 + +1.482 +
+

Les variables importantes sont celles du test Anova sans la variable cp.

+

On réalise une dernière méthode : la méthode both en affichant toujours les statistqiues descriptives du modèle obtenu. Les p_valeurs inférieures à 0.05 sont affichées en rouge.

+
# Méthode both
+mod_both <- step(mod_int, 
+                 target ~ age + sex + trestbps + chol + fbs + restecg + thalach
+                 + exang + oldpeak + slope + ca + thal, 
+                 data = train, 
+                 trace = F, 
+                 direction = c("both"))
+
+# Création d'un tableau qui contient les coefficients du modèle 
+tab_stat <- tidy(mod_both, conf.int = TRUE)
+
+# Mise en forme de la colonne p.value pour afficher en rouge les les p-valeurs inférieures à 0.05
+tab_stat <- tab_stat %>%
+  mutate(p.value = ifelse(p.value < 0.05, 
+                                    cell_spec(round(p.value, 3), "html", color = "red"), 
+                                    cell_spec(round(p.value, 3), "html")))
+
+# Affichage du tableau avec mise en forme
+tab_stat %>%
+  kbl(digits = c(0, 3, 3, 3, 4), # Arrondis
+      escape = FALSE,  
+      caption = "Statistiques récapitulatives du modèle construit par la méthode both") %>%
+  kable_styling()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+Table 3.6: Statistiques récapitulatives du modèle construit par la méthode both +
+term + +estimate + +std.error + +statistic + +p.value + +conf.low + +conf.high +
+(Intercept) + +0.706 + +1.607 + +0.439 + +0.661 + +-2 + +3.875 +
+oldpeak + +-0.634 + +0.187 + +-3.398 + +0.001 + +-1 + +-0.283 +
+exang + +-1.531 + +0.404 + +-3.789 + +0 + +-2 + +-0.753 +
+ca + +-0.668 + +0.179 + +-3.742 + +0 + +-1 + +-0.326 +
+thal + +-1.020 + +0.305 + +-3.345 + +0.001 + +-2 + +-0.434 +
+thalach + +0.026 + +0.009 + +2.740 + +0.006 + +0 + +0.044 +
+sex + +-1.029 + +0.411 + +-2.504 + +0.012 + +-2 + +-0.238 +
+restecg + +0.784 + +0.350 + +2.240 + +0.025 + +0 + +1.482 +
+

Le modèle garde les mêmes variables que celles obtenues avec la méthode forward.

+

Enfin, on compare l’AIC des modèles obtenus par les 3 méthodes.

+
# Création d'un dataframe avec les AIC
+deviance <- data.frame(AIC(mod_back), AIC(mod_for), AIC(mod_both))
+colnames(deviance) <- c("Méthode backward",
+                        "Méthode forward",
+                        "Méthode both") # Nom des colonnes 
+
+# Affichage des AIC des modèles
+deviance %>% 
+  kbl(digits = c(1, 1, 1), # pour les arrondis
+      caption = "AIC des modèles") %>% # Titre du tableau
+  kable_styling()
+ + + + + + + + + + + + + + + + +
+Table 3.7: AIC des modèles +
+Méthode backward + +Méthode forward + +Méthode both +
+194.5 + +209.6 + +209.6 +
+
+
+

3.3 Construction du modèle final

+

La méthode backward est la meilleure car elle a le plus petit AIC. On garde donc les variables du modèle construit à partir de cette méthode. +On définit alors notre modèle final avec les variables restecg, cp, sex, exang, oldpeak, ca, thal et thalach. On retire les variables age, trestbps, chol, fbs et slope.

+
# Construction du modèle final
+mod_final <- glm(target ~ sex + exang + oldpeak + ca + thal + restecg + thalach + cp,
+                 data = train, 
+                 family = "binomial")
+
+
+
+

4 Validation du modèle

+

Après avoir obtenu un modèle, il faut diagnostiquer la régression afin de valider ou non le modèle. +On passe alors à la validation de notre modèle ainsi construit. On regarde dans un premier temps la déviance du modèle puis on étudie les résidus et les outliers.

+
+

4.1 Déviance

+
# Création d'un dataframe avec les déviances
+deviance <- data.frame(mod$null.deviance, mod_final$deviance)
+colnames(deviance) <- c("Déviance du modèle réduit à l'intercept", 
+                        "Déviance du modèle final") # nom des colonnes 
+
+# Affichage des dimensions des ensembles d'entraînement et de test
+deviance %>% 
+  kbl(digits = c(2, 2), caption = "Déviance des modèles") %>% # pour les arrondis
+  kable_styling()
+ + + + + + + + + + + + + + +
+Table 4.1: Déviance des modèles +
+Déviance du modèle réduit à l’intercept + +Déviance du modèle final +
+333.48 + +178.01 +
+

La déviance de notre modèle étant plus petite que celle du modèle réduit à l’intercept, on le privilegie.

+
+
+

4.2 Résidus studentisés

+

On affiche ci-dessous un nuage de points montrant les résidus studentisés du modèle final pour chaque observation. +On a aussi tracé sur ce même graphe des lignes de seuil pour repérer les observations ayant des résidus considérés comme atypiques en dehors de ces bornes.

+
par(mfrow = c(1, 1))
+plot(rstudent(mod_final), type = "p", cex = 0.5, ylab = "Résidus studentisés ", 
+    col = "#0066CC", ylim = c(-3, 3))
+abline(h = c(-2, 2), col = "#CC3333")
+

+On remarque que les résidus suivent un schéma aléatoire, ce qui valide le modèle. +Les quelques points au-delà des lignes rouges sont les valeurs atypiques que nous avions déjà remarqué dans les boxplots. Nous avions décidé de garder ces valeurs.

+
+
+

4.3 Outliers

+

Enfin, on affiche le graphe des outliers.

+
# Plot des outliers
+plot(mod_final,5)
+

+Il ne semble pas y avoir d’outliers car aucun point n’a une distance de Cook supérieure à 1.

+
+
+
+

5 Optimisation du seuil de décision

+

On va maintenant optimsier le seuil à partir duquel on considère qu’un individu à une maladie cardiaque. Dans notre cas, on préfère naturellemnt avoir un petit nombre de faux négatifs (FN) c’est-à-dire un petit nombre de personnes qui n’ont pas été diagnostiquées comme positives alors qu’elles ont une maladie cardiaque.

+

Pour optimiser le seuil, on va donc maximsier l’indice de Youden qui est un critère couramment utilisé pour minimiser le nombre de faux négatif.

+

On affiche ici un graphe qui montre l’évolution de l’indice de Youden en fonction des seuils. On affiche sur ce même graphe le seuil optimal qui maximise cet indice.

+
# Calcul des valeurs de seuil
+val_seuil <- seq(0, 1, length.out = 100)
+
+# Calcul de l'indice de Youden pour chaque seuil
+youden_index <- sapply(val_seuil, function(threshold) {
+  pred_pos <- predicted_probs_train > threshold
+  sens <- sum(pred_pos & (train$target == 1)) / sum(train$target == 1)
+  spec <- sum(!pred_pos & (train$target == 0)) / sum(train$target == 0)
+  return(sens + spec - 1)
+})
+
+# Trouver le seuil optimal
+optimal <- val_seuil[which.max(youden_index)]
+
+# Créer le graphique de l'indice de Youden en fonction des seuils
+plot(val_seuil, 
+     youden_index, 
+     type = "l", 
+     col = "#0066CC", # Couleur de la courbe
+     xlab = "Seuils", # Nom de l'axe des abscisses 
+     ylab = "Indice de Youden",# Nom de l'axe des ordonnées 
+     main = "Graphique de l'évolution de l'ndice de Youden en fonction des seuils") # Titre du graphique 
+
+# Afficher le seuil optimal
+abline(v = optimal, 
+       col = "#CC3333", # Couleur de la droite
+       lty = 2) # Épaisseur de la droite 
+
+# Afficher la légende 
+text(optimal, 
+     max(youden_index), 
+     paste("Seuil optimal =", round(optimal, 3)), 
+     pos = 4, 
+     col = "#CC3333")
+

+Le seuil qui maximise l’indice de Youden vaut 0.556. On va donc garder ce seuil et prédire à patir de celui-ci. +On affiche ci-dessous l’évolution de la sensibilité et de la spécificité en fonction de la valeur du seuil. On affiche sur ce même-graphique la valeur du seuil optimal.

+
# Initialisation des vecteurs de sensibilité et spécificité
+sensitivities <- c()
+specificities <- c()
+
+for (threshold in val_seuil) {
+  pred_pos <- predicted_probs_train > threshold
+  sens <- sum(pred_pos & (train$target == 1)) / sum(train$target == 1)
+  spec <- sum(!pred_pos & (train$target == 0)) / sum(train$target == 0)
+  
+  sensitivities <- c(sensitivities, sens)
+  specificities <- c(specificities, spec)
+}
+
+# Création du graphique avec la courbe de l'évolution de la sensibilité
+plot(val_seuil, 
+     sensitivities, 
+     type = "l", 
+     col = "#0066CC", # Couleur de la courbe 
+     xlab = "Seuil", # Nom de l'axe des abscisses 
+     ylab = "Valeur des paramètres", # Nom de l'axe des ordonnées
+     main = "Évolution de la sensibilité et de la spécificité en fonction des seuils") # Titre du graphique
+
+# Ajout de la courbe de l'évolution de la spécificité
+lines(val_seuil, 
+      specificities, 
+      type = "l", 
+      col = "#339966") # Couleur de la droite 
+
+# Ajout de la légende 
+legend("bottomright", 
+       legend = c("Sensibilité", "Spécificité"), 
+       col = c("#0066CC", "#339966"), 
+       lty = 1)
+
+# Ajout de la valeur du seuil optimal
+abline(v = optimal, 
+       col = "#CC3333", 
+       lty = 2)
+
+# Ajout de texte
+text(optimal, 
+     max(youden_index),
+     paste("Seuil optimal =", round(optimal, 3)), 
+     pos = 4, 
+     col = "#CC3333")
+

+On remarque donc que la valeur optimale pour l’indice de Youden ne l’est pas pour la sensibilité ni la spécificité. Cela est du au fait qu’on a du faire un choix en selectionnant un indicateur à optimiser parmi plusieurs autres, ici, le nombre de faux négatifs. Cependant, au seuil optimal, la spécificité et la sensibilité sont satisfaisantes. On garde donc ce seuil.

+
+
+

6 Prédiction et performances

+
+

6.1 Prédiction

+

On passe maintenant à la prédiction. On prédit de telle manière que si la probabilité prédite dépasse le seuil optimal, l’individu est considéré comme étant susceptible de développer une maladie cardiaque.

+

On affiche ci-dessous les prédictions à côté des valeurs réelles dans un tableau.

+
# Création d'un dataframe avec les observations et les prédictions
+df_pred <- data.frame(Observed = test$target, Predicted = predicted)
+rmarkdown::paged_table(df_pred)
+
+ +
+
+
+

6.2 Matrice de confusion

+

Pour étudier les performances de prédiction, on affiche la matrice de confusion qui permet d’observer le nombre de vrais positifs, de faux positifs, de vrais négatifs et de faux négatifs.

+
# Création de la matrice de confusion
+confusion_matrix <- table(factor(predicted), factor(test$target))
+# Affichage de la matrice de confusion
+par(mfrow=c(1, 1))
+ctable <- as.table(confusion_matrix, 
+                   nrow = 2,
+                   byrow = TRUE)
+# Définir les étiquettes "malade" et "sain" pour les lignes et colonnes
+rownames(ctable) <- colnames(ctable) <- c("malade (1)", "sain (0)")
+
+# Définir les étiquettes "Réel" et "Prédit"
+dimnames(ctable) <- list(Réel = rownames(ctable), Prédit = colnames(ctable))
+
+fourfoldplot(ctable, 
+             color = c("#CC6666", "#CCCCCC"),
+             conf.level = 0, 
+             margin = 1, 
+             main = "Matrice de confusion \ndu modèle de regression logistique")
+

+On obtient qu’on a prédit 3 faux négatifs, 7 faux positifs, 21 vrais positifs et 30 vrais négatifs. Ce résultat est tout à fait satisfaisait car le nombre de faux négatifs est assez faible, ce qui était notre objectif.

+
+
+

6.3 Indicateurs et métriques de performance

+

La sensibilité, le F1 Score et la valeur prédictive négative sont des indicateurs importants pour évaluer la capacité du modèle à minimiser les faux négatifs, chacun apportant une compréhension particulière sur cette problématique. On s’intéresse donc ici plus particulièrement à ces indicateurs ainsi que l’accuracy (la précision globale).

+
# Créer des facteurs avec des niveaux correspondants
+predicted_factor <- factor(as.numeric(predicted), levels = c(0, 1))
+test_target_factor <- factor(as.numeric(test$target), levels = c(0, 1))
+
+# Simulation de la sortie confusionMatrix
+results <- confusionMatrix(predicted_factor, test_target_factor, mode = "everything", positive = "1")
+
+# Extraction de l'accuracy
+acc <- data.frame(Valeur = results$overall["Accuracy"])
+
+# Extraction des valeurs spécifiques
+specific_metrics <- c("Sensitivity", "Neg Pred Value", "F1")
+specific_values <- unlist(results$byClass[specific_metrics])
+
+# Création du tableau avec les valeurs spécifiques
+results_df_spe <- data.frame(Valeur = specific_values)
+
+# Affichage du tableau 
+rbind(acc, results_df_spe) %>% # Fusion des deux dataframe
+  kable(digits = 3, # Arrondi
+        caption = "Indicateurs et métriques de performance.") %>% # Titre
+  kable_styling()
+ + + + + + + + + + + + + + + + + + + + + + + + + + +
+Table 6.1: Table 6.2: Indicateurs et métriques de performance. +
+ +Valeur +
+Accuracy + +0.836 +
+Sensitivity + +0.909 +
+Neg Pred Value + +0.875 +
+F1 + +0.857 +
+

Grâce au tableau, on observe que la valeur des indicateurs est plutôt élevée et que donc notre modèle n’est pas trop mauvais !

+
+
+

6.4 Courbe ROC et AUC

+

On trace enfin la courbe ROC de notre modèle afin de déterminer l’AUC de notre modèle.

+
# Courbe ROC
+pred <- prediction(predicted_probs_test, test$target)
+roc <- performance(pred, measure = "tpr", x.measure = "fpr")
+
+# Calcul de l'AUC
+auc <- performance(pred, measure = "auc")
+auc_value <- round(unlist(slot(auc, "y.values")), 3)
+
+# Plot de la courbe ROC
+plot(roc, xlab = "Taux de faux positifs", ylab = "Taux de vrais positifs", col = "#CC3333", main = "Courbe ROC de l'ensemble de test")
+abline(a = 0, b = 1)
+
+# Ajout de l'AUC au graphique
+text(0.5, 0.8, labels = paste("AUC =", auc_value), cex = 1.2, col = "black")
+

+Dans notre cas, l’AUC est proche de 1 (environ 0.9) donc la prédiction semble être bonne.

+
+
+
+

7 Comparaison avec d’autres modèles

+

Pour terminer le projet, on compare le modèle de régression logistique à 2 autres modèles : un modèle de RadomForest et un modèle d’arbre de décision.

+
+

7.1 Déclaration des modèles

+

On construit d’abord un modèle de forêt aléatoire (Random Forest) qui est populaire en raison de sa capacité à produire des prédictions précises, à gérer des ensembles de données complexes et bruités, et à éviter le surajustement.

+

On construit ensuite un modèle d’arbre de décision, qui est plus approprié pour capturer des relations non linéaires entre les variables explicatives et la variable cible et des interactions complexes.

+

On peut afficher les matrices de confusions des modèles construits.

+
par(mfrow=c(1, 2))
+# Affichage de la matrice de confusion de la forêt aléatoire
+ctable <- as.table(results_rf, 
+                   nrow = 2,
+                   byrow = TRUE)
+# Définir les étiquettes "malade" et "sain" pour les lignes et colonnes
+rownames(ctable) <- colnames(ctable) <- c("malade (1)", "sain (0)")
+
+# Définir les étiquettes "Réel" et "Prédit"
+dimnames(ctable) <- list(Réel = rownames(ctable), Prédit = colnames(ctable))
+
+fourfoldplot(ctable, 
+             color = c("#CC6666", "#CCCCCC"),
+             conf.level = 0, 
+             margin = 1, 
+             main = "Matrice de confusion \ndu modèle de forêt aléatoire")
+
+# Affichage de la matrice de confusion de l'arbre de décision
+ctable <- as.table(results_tree, 
+                   nrow = 2,
+                   byrow = TRUE)
+# Définir les étiquettes "malade" et "sain" pour les lignes et colonnes
+rownames(ctable) <- colnames(ctable) <- c("malade (1)", "sain (0)")
+
+# Définir les étiquettes "Réel" et "Prédit"
+dimnames(ctable) <- list(Réel = rownames(ctable), Prédit = colnames(ctable))
+
+fourfoldplot(ctable, 
+             color = c("#CC6666", "#CCCCCC"),
+             conf.level = 0, 
+             margin = 1, 
+             main = "Matrice de confusion \ndu modèle d'arbre de décision")
+

+
+
+

7.2 Comparaison des performances

+

On va comparer la sensibilité, le F1 Score et la valeur prédictive négative car, comme vu précédemment, ce sont des indicateurs importants pour évaluer la capacité du modèle à minimiser les faux négatifs. On affiche aussi le nombre de faux négatifs. Enfin, on compare l’accuracy et l’AUC des différents modèles car c’est un indicateur souvent comparé.

+
# Calcul de l'AUC RF
+pred_rf <- prediction(predicted_probs_test_rf, test$target)
+auc_rf <- performance(pred_rf, measure = "auc")
+auc_value_rf <- round(unlist(slot(auc_rf, "y.values")), 3)
+
+# Calcul de l'AUC Tree
+pred_tree <- prediction(predicted_probs_test_tree, test$target)
+# Calcul de l'AUC
+auc_tree <- performance(pred_tree, measure = "auc")
+auc_value_tree <- round(unlist(slot(auc_tree, "y.values")), 3)
+
+# Extraction des valeurs spécifiques
+specific_metrics <- c("Sensitivity", "Neg Pred Value", "Precision")
+
+# Création d'un dataframe avec les valeurs spécifiques
+metrics <- c("Sensibilité", "F1 Score", "Valeur prédictive négative", "Accuracy", "Nombre de faux négatifs", "AUC")
+models <- c("Regression logistique", "Random Forest", "Arbre de décision")
+data <- matrix(NA, nrow = length(models), ncol = length(metrics))
+colnames(data) <- metrics 
+rownames(data) <- models
+
+# Insertion des valeurs spécifiques dans le dataframe
+data["Regression logistique", ] <- c(unlist(results$byClass[specific_metrics ]), results$overall['Accuracy'], results$table[1, 2], auc_value)
+data["Random Forest",] <- c(unlist(results_rf$byClass[specific_metrics]), results_rf$overall['Accuracy'], results_rf$table[1, 2], auc_value_rf)
+data["Arbre de décision",] <- c(unlist(results_tree$byClass[specific_metrics]), results_tree$overall['Accuracy'], results_tree$table[1, 2], auc_value_tree)
+
+# Affichage du tableau
+data %>%
+  kbl(digits = c(3, 3, 3, 3, 1, 2), caption = "Comparaison des performances des différents modèles ") %>% # pour les arrondis et titres  
+  kable_styling()
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+Table 7.1: Comparaison des performances des différents modèles +
+ +Sensibilité + +F1 Score + +Valeur prédictive négative + +Accuracy + +Nombre de faux négatifs + +AUC +
+Regression logistique + +0.909 + +0.875 + +0.811 + +0.836 + +3 + +0.92 +
+Random Forest + +0.879 + +0.833 + +0.784 + +0.803 + +4 + +0.94 +
+Arbre de décision + +0.939 + +0.909 + +0.795 + +0.836 + +2 + +0.89 +
+

On obtient des modèles qui ont des performances assez similaires.

+
+
+
+

8 Conclusion

+

Dans le cadre de ce projet, axé sur l’explication et la prédiction des maladies cardiaques, nous avons développé et évalué un modèle de régression logistique, qui est un modèle facilement interprétable. Les performances de ce modèle se sont avérées prometteuses, démontrant une bonne capacité à expliquer et prédire la présence de maladies cardiaques. Ce modèle a été comparé à des approches alternatives telles que la forêt aléatoire et l’arbre de décision, qui ont montré des performances similaires.

+

D’autres modèles, comme le SVM, la classification naïve bayésienne, le KNN, ou le Gradient Boosting, pourraient être explorés pour améliorer la prédiction.

+

Des observations paradoxales ont été constatées dans l’application. +Ces résultats inattendus pourraient s’expliquer par des interactions complexes entre les caractéristiques. Introduire des interactions entre les variables pourrait améliorer la capacité prédictive des modèles. +De plus, des biais statistiques dus à la taille limitée de l’échantillon pourraient impacter certaines catégories sous-représentées. +Enfin, certains facteurs externes importants, tels que les comportements, les facteurs environnementaux ou les conditions médicales spécifiques, n’ont pas été inclus dans l’étude. Prendre en compte ces éléments dans de futures recherches pourrait affiner les prédictions des maladies cardiaques.

+
+ + + +
+
+
+
+
+ + + + + + + + + + + + diff --git a/M2/Data Visualisation/Exemple Projet/heart.csv b/M2/Data Visualisation/Exemple Projet/heart.csv new file mode 100644 index 0000000..1939fc6 --- /dev/null +++ b/M2/Data Visualisation/Exemple Projet/heart.csv @@ -0,0 +1,304 @@ +age,sex,cp,trestbps,chol,fbs,restecg,thalach,exang,oldpeak,slope,ca,thal,target +63,1,3,145,233,1,0,150,0,2.3,0,0,1,1 +37,1,2,130,250,0,1,187,0,3.5,0,0,2,1 +41,0,1,130,204,0,0,172,0,1.4,2,0,2,1 +56,1,1,120,236,0,1,178,0,0.8,2,0,2,1 +57,0,0,120,354,0,1,163,1,0.6,2,0,2,1 +57,1,0,140,192,0,1,148,0,0.4,1,0,1,1 +56,0,1,140,294,0,0,153,0,1.3,1,0,2,1 +44,1,1,120,263,0,1,173,0,0,2,0,3,1 +52,1,2,172,199,1,1,162,0,0.5,2,0,3,1 +57,1,2,150,168,0,1,174,0,1.6,2,0,2,1 +54,1,0,140,239,0,1,160,0,1.2,2,0,2,1 +48,0,2,130,275,0,1,139,0,0.2,2,0,2,1 +49,1,1,130,266,0,1,171,0,0.6,2,0,2,1 +64,1,3,110,211,0,0,144,1,1.8,1,0,2,1 +58,0,3,150,283,1,0,162,0,1,2,0,2,1 +50,0,2,120,219,0,1,158,0,1.6,1,0,2,1 +58,0,2,120,340,0,1,172,0,0,2,0,2,1 +66,0,3,150,226,0,1,114,0,2.6,0,0,2,1 +43,1,0,150,247,0,1,171,0,1.5,2,0,2,1 +69,0,3,140,239,0,1,151,0,1.8,2,2,2,1 +59,1,0,135,234,0,1,161,0,0.5,1,0,3,1 +44,1,2,130,233,0,1,179,1,0.4,2,0,2,1 +42,1,0,140,226,0,1,178,0,0,2,0,2,1 +61,1,2,150,243,1,1,137,1,1,1,0,2,1 +40,1,3,140,199,0,1,178,1,1.4,2,0,3,1 +71,0,1,160,302,0,1,162,0,0.4,2,2,2,1 +59,1,2,150,212,1,1,157,0,1.6,2,0,2,1 +51,1,2,110,175,0,1,123,0,0.6,2,0,2,1 +65,0,2,140,417,1,0,157,0,0.8,2,1,2,1 +53,1,2,130,197,1,0,152,0,1.2,0,0,2,1 +41,0,1,105,198,0,1,168,0,0,2,1,2,1 +65,1,0,120,177,0,1,140,0,0.4,2,0,3,1 +44,1,1,130,219,0,0,188,0,0,2,0,2,1 +54,1,2,125,273,0,0,152,0,0.5,0,1,2,1 +51,1,3,125,213,0,0,125,1,1.4,2,1,2,1 +46,0,2,142,177,0,0,160,1,1.4,0,0,2,1 +54,0,2,135,304,1,1,170,0,0,2,0,2,1 +54,1,2,150,232,0,0,165,0,1.6,2,0,3,1 +65,0,2,155,269,0,1,148,0,0.8,2,0,2,1 +65,0,2,160,360,0,0,151,0,0.8,2,0,2,1 +51,0,2,140,308,0,0,142,0,1.5,2,1,2,1 +48,1,1,130,245,0,0,180,0,0.2,1,0,2,1 +45,1,0,104,208,0,0,148,1,3,1,0,2,1 +53,0,0,130,264,0,0,143,0,0.4,1,0,2,1 +39,1,2,140,321,0,0,182,0,0,2,0,2,1 +52,1,1,120,325,0,1,172,0,0.2,2,0,2,1 +44,1,2,140,235,0,0,180,0,0,2,0,2,1 +47,1,2,138,257,0,0,156,0,0,2,0,2,1 +53,0,2,128,216,0,0,115,0,0,2,0,0,1 +53,0,0,138,234,0,0,160,0,0,2,0,2,1 +51,0,2,130,256,0,0,149,0,0.5,2,0,2,1 +66,1,0,120,302,0,0,151,0,0.4,1,0,2,1 +62,1,2,130,231,0,1,146,0,1.8,1,3,3,1 +44,0,2,108,141,0,1,175,0,0.6,1,0,2,1 +63,0,2,135,252,0,0,172,0,0,2,0,2,1 +52,1,1,134,201,0,1,158,0,0.8,2,1,2,1 +48,1,0,122,222,0,0,186,0,0,2,0,2,1 +45,1,0,115,260,0,0,185,0,0,2,0,2,1 +34,1,3,118,182,0,0,174,0,0,2,0,2,1 +57,0,0,128,303,0,0,159,0,0,2,1,2,1 +71,0,2,110,265,1,0,130,0,0,2,1,2,1 +54,1,1,108,309,0,1,156,0,0,2,0,3,1 +52,1,3,118,186,0,0,190,0,0,1,0,1,1 +41,1,1,135,203,0,1,132,0,0,1,0,1,1 +58,1,2,140,211,1,0,165,0,0,2,0,2,1 +35,0,0,138,183,0,1,182,0,1.4,2,0,2,1 +51,1,2,100,222,0,1,143,1,1.2,1,0,2,1 +45,0,1,130,234,0,0,175,0,0.6,1,0,2,1 +44,1,1,120,220,0,1,170,0,0,2,0,2,1 +62,0,0,124,209,0,1,163,0,0,2,0,2,1 +54,1,2,120,258,0,0,147,0,0.4,1,0,3,1 +51,1,2,94,227,0,1,154,1,0,2,1,3,1 +29,1,1,130,204,0,0,202,0,0,2,0,2,1 +51,1,0,140,261,0,0,186,1,0,2,0,2,1 +43,0,2,122,213,0,1,165,0,0.2,1,0,2,1 +55,0,1,135,250,0,0,161,0,1.4,1,0,2,1 +51,1,2,125,245,1,0,166,0,2.4,1,0,2,1 +59,1,1,140,221,0,1,164,1,0,2,0,2,1 +52,1,1,128,205,1,1,184,0,0,2,0,2,1 +58,1,2,105,240,0,0,154,1,0.6,1,0,3,1 +41,1,2,112,250,0,1,179,0,0,2,0,2,1 +45,1,1,128,308,0,0,170,0,0,2,0,2,1 +60,0,2,102,318,0,1,160,0,0,2,1,2,1 +52,1,3,152,298,1,1,178,0,1.2,1,0,3,1 +42,0,0,102,265,0,0,122,0,0.6,1,0,2,1 +67,0,2,115,564,0,0,160,0,1.6,1,0,3,1 +68,1,2,118,277,0,1,151,0,1,2,1,3,1 +46,1,1,101,197,1,1,156,0,0,2,0,3,1 +54,0,2,110,214,0,1,158,0,1.6,1,0,2,1 +58,0,0,100,248,0,0,122,0,1,1,0,2,1 +48,1,2,124,255,1,1,175,0,0,2,2,2,1 +57,1,0,132,207,0,1,168,1,0,2,0,3,1 +52,1,2,138,223,0,1,169,0,0,2,4,2,1 +54,0,1,132,288,1,0,159,1,0,2,1,2,1 +45,0,1,112,160,0,1,138,0,0,1,0,2,1 +53,1,0,142,226,0,0,111,1,0,2,0,3,1 +62,0,0,140,394,0,0,157,0,1.2,1,0,2,1 +52,1,0,108,233,1,1,147,0,0.1,2,3,3,1 +43,1,2,130,315,0,1,162,0,1.9,2,1,2,1 +53,1,2,130,246,1,0,173,0,0,2,3,2,1 +42,1,3,148,244,0,0,178,0,0.8,2,2,2,1 +59,1,3,178,270,0,0,145,0,4.2,0,0,3,1 +63,0,1,140,195,0,1,179,0,0,2,2,2,1 +42,1,2,120,240,1,1,194,0,0.8,0,0,3,1 +50,1,2,129,196,0,1,163,0,0,2,0,2,1 +68,0,2,120,211,0,0,115,0,1.5,1,0,2,1 +69,1,3,160,234,1,0,131,0,0.1,1,1,2,1 +45,0,0,138,236,0,0,152,1,0.2,1,0,2,1 +50,0,1,120,244,0,1,162,0,1.1,2,0,2,1 +50,0,0,110,254,0,0,159,0,0,2,0,2,1 +64,0,0,180,325,0,1,154,1,0,2,0,2,1 +57,1,2,150,126,1,1,173,0,0.2,2,1,3,1 +64,0,2,140,313,0,1,133,0,0.2,2,0,3,1 +43,1,0,110,211,0,1,161,0,0,2,0,3,1 +55,1,1,130,262,0,1,155,0,0,2,0,2,1 +37,0,2,120,215,0,1,170,0,0,2,0,2,1 +41,1,2,130,214,0,0,168,0,2,1,0,2,1 +56,1,3,120,193,0,0,162,0,1.9,1,0,3,1 +46,0,1,105,204,0,1,172,0,0,2,0,2,1 +46,0,0,138,243,0,0,152,1,0,1,0,2,1 +64,0,0,130,303,0,1,122,0,2,1,2,2,1 +59,1,0,138,271,0,0,182,0,0,2,0,2,1 +41,0,2,112,268,0,0,172,1,0,2,0,2,1 +54,0,2,108,267,0,0,167,0,0,2,0,2,1 +39,0,2,94,199,0,1,179,0,0,2,0,2,1 +34,0,1,118,210,0,1,192,0,0.7,2,0,2,1 +47,1,0,112,204,0,1,143,0,0.1,2,0,2,1 +67,0,2,152,277,0,1,172,0,0,2,1,2,1 +52,0,2,136,196,0,0,169,0,0.1,1,0,2,1 +74,0,1,120,269,0,0,121,1,0.2,2,1,2,1 +54,0,2,160,201,0,1,163,0,0,2,1,2,1 +49,0,1,134,271,0,1,162,0,0,1,0,2,1 +42,1,1,120,295,0,1,162,0,0,2,0,2,1 +41,1,1,110,235,0,1,153,0,0,2,0,2,1 +41,0,1,126,306,0,1,163,0,0,2,0,2,1 +49,0,0,130,269,0,1,163,0,0,2,0,2,1 +60,0,2,120,178,1,1,96,0,0,2,0,2,1 +62,1,1,128,208,1,0,140,0,0,2,0,2,1 +57,1,0,110,201,0,1,126,1,1.5,1,0,1,1 +64,1,0,128,263,0,1,105,1,0.2,1,1,3,1 +51,0,2,120,295,0,0,157,0,0.6,2,0,2,1 +43,1,0,115,303,0,1,181,0,1.2,1,0,2,1 +42,0,2,120,209,0,1,173,0,0,1,0,2,1 +67,0,0,106,223,0,1,142,0,0.3,2,2,2,1 +76,0,2,140,197,0,2,116,0,1.1,1,0,2,1 +70,1,1,156,245,0,0,143,0,0,2,0,2,1 +44,0,2,118,242,0,1,149,0,0.3,1,1,2,1 +60,0,3,150,240,0,1,171,0,0.9,2,0,2,1 +44,1,2,120,226,0,1,169,0,0,2,0,2,1 +42,1,2,130,180,0,1,150,0,0,2,0,2,1 +66,1,0,160,228,0,0,138,0,2.3,2,0,1,1 +71,0,0,112,149,0,1,125,0,1.6,1,0,2,1 +64,1,3,170,227,0,0,155,0,0.6,1,0,3,1 +66,0,2,146,278,0,0,152,0,0,1,1,2,1 +39,0,2,138,220,0,1,152,0,0,1,0,2,1 +58,0,0,130,197,0,1,131,0,0.6,1,0,2,1 +47,1,2,130,253,0,1,179,0,0,2,0,2,1 +35,1,1,122,192,0,1,174,0,0,2,0,2,1 +58,1,1,125,220,0,1,144,0,0.4,1,4,3,1 +56,1,1,130,221,0,0,163,0,0,2,0,3,1 +56,1,1,120,240,0,1,169,0,0,0,0,2,1 +55,0,1,132,342,0,1,166,0,1.2,2,0,2,1 +41,1,1,120,157,0,1,182,0,0,2,0,2,1 +38,1,2,138,175,0,1,173,0,0,2,4,2,1 +38,1,2,138,175,0,1,173,0,0,2,4,2,1 +67,1,0,160,286,0,0,108,1,1.5,1,3,2,0 +67,1,0,120,229,0,0,129,1,2.6,1,2,3,0 +62,0,0,140,268,0,0,160,0,3.6,0,2,2,0 +63,1,0,130,254,0,0,147,0,1.4,1,1,3,0 +53,1,0,140,203,1,0,155,1,3.1,0,0,3,0 +56,1,2,130,256,1,0,142,1,0.6,1,1,1,0 +48,1,1,110,229,0,1,168,0,1,0,0,3,0 +58,1,1,120,284,0,0,160,0,1.8,1,0,2,0 +58,1,2,132,224,0,0,173,0,3.2,2,2,3,0 +60,1,0,130,206,0,0,132,1,2.4,1,2,3,0 +40,1,0,110,167,0,0,114,1,2,1,0,3,0 +60,1,0,117,230,1,1,160,1,1.4,2,2,3,0 +64,1,2,140,335,0,1,158,0,0,2,0,2,0 +43,1,0,120,177,0,0,120,1,2.5,1,0,3,0 +57,1,0,150,276,0,0,112,1,0.6,1,1,1,0 +55,1,0,132,353,0,1,132,1,1.2,1,1,3,0 +65,0,0,150,225,0,0,114,0,1,1,3,3,0 +61,0,0,130,330,0,0,169,0,0,2,0,2,0 +58,1,2,112,230,0,0,165,0,2.5,1,1,3,0 +50,1,0,150,243,0,0,128,0,2.6,1,0,3,0 +44,1,0,112,290,0,0,153,0,0,2,1,2,0 +60,1,0,130,253,0,1,144,1,1.4,2,1,3,0 +54,1,0,124,266,0,0,109,1,2.2,1,1,3,0 +50,1,2,140,233,0,1,163,0,0.6,1,1,3,0 +41,1,0,110,172,0,0,158,0,0,2,0,3,0 +51,0,0,130,305,0,1,142,1,1.2,1,0,3,0 +58,1,0,128,216,0,0,131,1,2.2,1,3,3,0 +54,1,0,120,188,0,1,113,0,1.4,1,1,3,0 +60,1,0,145,282,0,0,142,1,2.8,1,2,3,0 +60,1,2,140,185,0,0,155,0,3,1,0,2,0 +59,1,0,170,326,0,0,140,1,3.4,0,0,3,0 +46,1,2,150,231,0,1,147,0,3.6,1,0,2,0 +67,1,0,125,254,1,1,163,0,0.2,1,2,3,0 +62,1,0,120,267,0,1,99,1,1.8,1,2,3,0 +65,1,0,110,248,0,0,158,0,0.6,2,2,1,0 +44,1,0,110,197,0,0,177,0,0,2,1,2,0 +60,1,0,125,258,0,0,141,1,2.8,1,1,3,0 +58,1,0,150,270,0,0,111,1,0.8,2,0,3,0 +68,1,2,180,274,1,0,150,1,1.6,1,0,3,0 +62,0,0,160,164,0,0,145,0,6.2,0,3,3,0 +52,1,0,128,255,0,1,161,1,0,2,1,3,0 +59,1,0,110,239,0,0,142,1,1.2,1,1,3,0 +60,0,0,150,258,0,0,157,0,2.6,1,2,3,0 +49,1,2,120,188,0,1,139,0,2,1,3,3,0 +59,1,0,140,177,0,1,162,1,0,2,1,3,0 +57,1,2,128,229,0,0,150,0,0.4,1,1,3,0 +61,1,0,120,260,0,1,140,1,3.6,1,1,3,0 +39,1,0,118,219,0,1,140,0,1.2,1,0,3,0 +61,0,0,145,307,0,0,146,1,1,1,0,3,0 +56,1,0,125,249,1,0,144,1,1.2,1,1,2,0 +43,0,0,132,341,1,0,136,1,3,1,0,3,0 +62,0,2,130,263,0,1,97,0,1.2,1,1,3,0 +63,1,0,130,330,1,0,132,1,1.8,2,3,3,0 +65,1,0,135,254,0,0,127,0,2.8,1,1,3,0 +48,1,0,130,256,1,0,150,1,0,2,2,3,0 +63,0,0,150,407,0,0,154,0,4,1,3,3,0 +55,1,0,140,217,0,1,111,1,5.6,0,0,3,0 +65,1,3,138,282,1,0,174,0,1.4,1,1,2,0 +56,0,0,200,288,1,0,133,1,4,0,2,3,0 +54,1,0,110,239,0,1,126,1,2.8,1,1,3,0 +70,1,0,145,174,0,1,125,1,2.6,0,0,3,0 +62,1,1,120,281,0,0,103,0,1.4,1,1,3,0 +35,1,0,120,198,0,1,130,1,1.6,1,0,3,0 +59,1,3,170,288,0,0,159,0,0.2,1,0,3,0 +64,1,2,125,309,0,1,131,1,1.8,1,0,3,0 +47,1,2,108,243,0,1,152,0,0,2,0,2,0 +57,1,0,165,289,1,0,124,0,1,1,3,3,0 +55,1,0,160,289,0,0,145,1,0.8,1,1,3,0 +64,1,0,120,246,0,0,96,1,2.2,0,1,2,0 +70,1,0,130,322,0,0,109,0,2.4,1,3,2,0 +51,1,0,140,299,0,1,173,1,1.6,2,0,3,0 +58,1,0,125,300,0,0,171,0,0,2,2,3,0 +60,1,0,140,293,0,0,170,0,1.2,1,2,3,0 +77,1,0,125,304,0,0,162,1,0,2,3,2,0 +35,1,0,126,282,0,0,156,1,0,2,0,3,0 +70,1,2,160,269,0,1,112,1,2.9,1,1,3,0 +59,0,0,174,249,0,1,143,1,0,1,0,2,0 +64,1,0,145,212,0,0,132,0,2,1,2,1,0 +57,1,0,152,274,0,1,88,1,1.2,1,1,3,0 +56,1,0,132,184,0,0,105,1,2.1,1,1,1,0 +48,1,0,124,274,0,0,166,0,0.5,1,0,3,0 +56,0,0,134,409,0,0,150,1,1.9,1,2,3,0 +66,1,1,160,246,0,1,120,1,0,1,3,1,0 +54,1,1,192,283,0,0,195,0,0,2,1,3,0 +69,1,2,140,254,0,0,146,0,2,1,3,3,0 +51,1,0,140,298,0,1,122,1,4.2,1,3,3,0 +43,1,0,132,247,1,0,143,1,0.1,1,4,3,0 +62,0,0,138,294,1,1,106,0,1.9,1,3,2,0 +67,1,0,100,299,0,0,125,1,0.9,1,2,2,0 +59,1,3,160,273,0,0,125,0,0,2,0,2,0 +45,1,0,142,309,0,0,147,1,0,1,3,3,0 +58,1,0,128,259,0,0,130,1,3,1,2,3,0 +50,1,0,144,200,0,0,126,1,0.9,1,0,3,0 +62,0,0,150,244,0,1,154,1,1.4,1,0,2,0 +38,1,3,120,231,0,1,182,1,3.8,1,0,3,0 +66,0,0,178,228,1,1,165,1,1,1,2,3,0 +52,1,0,112,230,0,1,160,0,0,2,1,2,0 +53,1,0,123,282,0,1,95,1,2,1,2,3,0 +63,0,0,108,269,0,1,169,1,1.8,1,2,2,0 +54,1,0,110,206,0,0,108,1,0,1,1,2,0 +66,1,0,112,212,0,0,132,1,0.1,2,1,2,0 +55,0,0,180,327,0,2,117,1,3.4,1,0,2,0 +49,1,2,118,149,0,0,126,0,0.8,2,3,2,0 +54,1,0,122,286,0,0,116,1,3.2,1,2,2,0 +56,1,0,130,283,1,0,103,1,1.6,0,0,3,0 +46,1,0,120,249,0,0,144,0,0.8,2,0,3,0 +61,1,3,134,234,0,1,145,0,2.6,1,2,2,0 +67,1,0,120,237,0,1,71,0,1,1,0,2,0 +58,1,0,100,234,0,1,156,0,0.1,2,1,3,0 +47,1,0,110,275,0,0,118,1,1,1,1,2,0 +52,1,0,125,212,0,1,168,0,1,2,2,3,0 +58,1,0,146,218,0,1,105,0,2,1,1,3,0 +57,1,1,124,261,0,1,141,0,0.3,2,0,3,0 +58,0,1,136,319,1,0,152,0,0,2,2,2,0 +61,1,0,138,166,0,0,125,1,3.6,1,1,2,0 +42,1,0,136,315,0,1,125,1,1.8,1,0,1,0 +52,1,0,128,204,1,1,156,1,1,1,0,0,0 +59,1,2,126,218,1,1,134,0,2.2,1,1,1,0 +40,1,0,152,223,0,1,181,0,0,2,0,3,0 +61,1,0,140,207,0,0,138,1,1.9,2,1,3,0 +46,1,0,140,311,0,1,120,1,1.8,1,2,3,0 +59,1,3,134,204,0,1,162,0,0.8,2,2,2,0 +57,1,1,154,232,0,0,164,0,0,2,1,2,0 +57,1,0,110,335,0,1,143,1,3,1,1,3,0 +55,0,0,128,205,0,2,130,1,2,1,1,3,0 +61,1,0,148,203,0,1,161,0,0,2,1,3,0 +58,1,0,114,318,0,2,140,0,4.4,0,3,1,0 +58,0,0,170,225,1,0,146,1,2.8,1,2,1,0 +67,1,2,152,212,0,0,150,0,0.8,1,0,3,0 +44,1,0,120,169,0,1,144,1,2.8,0,0,1,0 +63,1,0,140,187,0,0,144,1,4,2,2,3,0 +63,0,0,124,197,0,1,136,1,0,1,0,2,0 +59,1,0,164,176,1,0,90,0,1,1,2,1,0 +57,0,0,140,241,0,1,123,1,0.2,1,0,3,0 +45,1,3,110,264,0,1,132,0,1.2,1,0,3,0 +68,1,0,144,193,1,1,141,0,3.4,1,2,3,0 +57,1,0,130,131,0,1,115,1,1.2,1,1,3,0 +57,0,1,130,236,0,0,174,0,0,1,1,2,0 diff --git a/M2/Data Visualisation/tp1/3-td_ggplot2 - enonce.Rmd b/M2/Data Visualisation/tp1/3-td_ggplot2 - enonce.Rmd new file mode 100644 index 0000000..23ff695 --- /dev/null +++ b/M2/Data Visualisation/tp1/3-td_ggplot2 - enonce.Rmd @@ -0,0 +1,629 @@ +--- +title: "Prise en main de ggplot2" +author: "Quentin Guibert" +date: "Année 2025-2026" +institute: "Université Paris-Dauphine | Master ISF" +lang: fr +link-citations: true +output: + rmdformats::robobook: + highlight: kate + use_bookdown: true + css: style.css + lightbox : true + gallery: true + code_folding: show + theme: flatly + toc_float: + collapsed: no +editor_options: + markdown: + wrap: 72 +# bibliography: references.bib +--- + +```{r setup, include=FALSE} +## Global options +knitr::opts_chunk$set( + cache = FALSE, + warning = FALSE, + message = FALSE, + fig.retina = 2, + fig.height = 6, + fig.width = 12 +) +options(encoding = 'UTF-8') +``` + +```{r, echo = FALSE, fig.keep= 'none'} +# Chargement des librairies graphiques +library(lattice) +library(grid) +library(ggplot2) +require(gridExtra) +library(locfit) +library(scales) +library(formattable) +library(RColorBrewer) +library(plotly) +library(dplyr) +library(tidyr) +library(rmarkdown) +library(ggthemes) +library(cowplot) +library(kableExtra) +``` + +# Objectifs du TP + +L'objectif de ce TP vise à se familiariser avec le package **ggplot2** +de `R`. Il s'agit de faire des manipulations graphiques élémentaires et +d'interpréter les résultats de ces visualisations. + +Dans un premier temps, vous pouvez suivre l'exemple introductif en +répliquant le code fourni. Dans un deuxième temps, il convient de +réaliser l'exercice point par point. + +# Prérequis + +- Avoir installer `R` [ici](https://www.r-project.org/). +- Avoir installer un IDE, par exemple `RStudio` + [ici](https://posit.co/download/rstudio-desktop/). +- Créer un nouveau projet (`File`, puis `New Projet`) dans un dossier + sur votre ordinateur. +- Télécharger [ici](https://moodle.psl.eu/course/view.php?id=33799) + les fichiers nécessaires au TD. + +Vous pouvez ensuite écrire vos codes soit : + +- En ouvrant un nouveau script `.R` ; +- En ouvrant le ouvrant le rapport Rmarkdown `3-td_ggplot2 - enonce`. + Certains codes sont partiels et sont à compléter (indication `???`). + N'oubliez pas de modifier l'option `eval = TRUE` pour que les + calculs puissent être réalisés. + +# Exemple introductif + +Pour illustrer cette première partie, nous reprenons l'exemple +introductif fourni par @wickham2023 sur le jeu de données `penguins` du +package **palmerpenguins**. Ce jeu de données s'intèresse des mesures +réalisées sur des manchots sur 3 îles de l'archipelle Palmer. + +## Données + +Dans un premier temps, il faut installer le package et le charger. + +```{r} +# install.packages("palmerpenguins") +library(palmerpenguins) +``` + +Ce jeu de données contient 344 observations où chaque ligne correspond à +un individu. + +```{r} +paged_table(penguins, options = list(rows.print = 15)) +``` + +On se concentre plus particulièrement sur les variables suivantes : + +- `species` : l'espèce de manchot ; +- `flipper_length_mm` : la longueur de la nageoire en mm ; +- `body_mass_g` : la masse en gramme. + +Pour plus détails, voir l'aide `?penguins`. + +## But de la visualisation + +On s'intéresse au lien entre le masse et la taille des nageoires des +manchots : + +- ceux dont les nageoires sont les plus longues sont-ils plus lourds + que les manchots aux nageoires courtes ? +- si oui quelle est le type de relation (linéaire, croissante, + décroissante, ...) ? +- quels facteurs influencent également cette relation (lieu, l'espèce, + ... ) ? + +On cherche à recréer la figure suivante. + +![](images\manchot.png) + +## Création de la figure étape par étape + +### Etape 1 : Scatterplot {.unnumbered} + +On commence par créer un scatterplot pour examiner la relation entre la +masse et la taille de la nageoire. + +```{r} +ggplot( + data = penguins, + mapping = aes(x = flipper_length_mm, y = body_mass_g) +) + + geom_point() +``` + +Cette figure fait clairement apparaître une relation croissante et a +priori linéaire entre les deux variables. + +::: remark-box +Un message d'erreur apparaît pour deux individus avec des données +manquantes. Ils sont automatiquement exclus. +::: + +### Etape 2 : Ajout d'élements esthétiques {.unnumbered} + +On cherche à présent exhiber le rôle de l’espèce à partir d'une couleur. +Trois espèces sont présents, ainsi l'ajout de 3 couleurs à la figure ne +devrait pas surcharger le graphique. + +```{r} +ggplot( + data = penguins, + mapping = aes(x = flipper_length_mm, y = body_mass_g, color = species) +) + + geom_point() +``` + +Compte tenu du nombre important de points, nous pouvons renforcer les +différences par espèce en ajoutant une variation de forme aux points. + +```{r} +ggplot( + data = penguins, + mapping = aes(x = flipper_length_mm, y = body_mass_g) +) + + geom_point( + mapping = aes( + color = species, + shape = species + ) + ) +``` + +### Etape 3 : Ajout d'une géométrie {.unnumbered} + +Voyons à présent comment interpréter la nature de la relation entre +masse et longueur de la nageoire. Pour ce faire, nous essayons d'ajout +des courbes de tendance. Nous commençons par une tendance linéaire. + +```{r} +ggplot( + data = penguins, + mapping = aes(x = flipper_length_mm, y = body_mass_g) +) + + geom_point( + mapping = aes( + color = species, + shape = species + ) + ) + + geom_smooth(method = "lm") +``` + +La même figure peut être générée par espèce en déplaçant l'argument +`color = species`. + +```{r} +ggplot( + data = penguins, + mapping = aes(x = flipper_length_mm, y = body_mass_g, color = species) +) + + geom_point( + mapping = aes( + shape = species + ) + ) + + geom_smooth(method = "lm") +``` + +Les pentes entre les espèces ne sont pas si éloignées. Nous décidons que +conserver une relation commune pour toutes espèces. Pour tester si la +nature linéaire de la relation est a priori une bonne hypothèse, nous +considérons un lissage non-paramétrique. + +```{r} +ggplot( + data = penguins, + mapping = aes(x = flipper_length_mm, y = body_mass_g) +) + + geom_point( + mapping = aes( + color = species, + shape = species + ) + ) + + geom_smooth(method = "loess") +``` + +L'ajout d'un lissage non-paramétrique permet d'affiner l’adéquation aux +données, mais sans pour autant clairement remettre en cause la tendance +linéaire qui sera donc conservée. + +### Etape 4 : Ajout des titres et changement de thème {.unnumbered} + +Afin de finaliser la figure, nous ajouter : + +- un titre ; +- un sous-titre ; +- des titres aux axes ; +- un titre à la légende. + +Ces informations sont ajoutées avec `labs()`. + +De plus, nous modifions le thème avec la commande `theme_bw()`. + +```{r} +ggplot( + data = penguins, + mapping = aes(x = flipper_length_mm, y = body_mass_g) +) + + geom_point(aes(color = species, shape = species)) + + geom_smooth(method = "lm") + + labs( + title = "Masse et taille de la nageoire", + subtitle = "Manchots d'Adelie, a + jugulaire et de Gentoo", + x = "Longueur de la nageoire (mm)", + y = "Masse (g)", + color = "Espece", + shape = "Espece" + ) + + scale_color_colorblind() + + theme_bw() +``` + +------------------------------------------------------------------------ + +# Exercice + +## Données + +Nous travaillons avec les jeux de données `FreMTPL2freq` et +`FreMTPL2sev` du package **Casdatasets**. Ces données ont été +préalablement pré-formatées et regroupées. + +Ce jeux de données regroupent les caractéristiques de 677 991 polices de +responsabilité civile automobile, observées principalement sur une +année. Dans les données regroupées, on dispose des numéros de sinistre +par police, des montants de sinistre correspondants, des +caractéristiques du risque et du nombre de sinistres. + +On présente ci-dessous un aperçu des données. + +```{r begin} +# Folds +fold <- getwd() + +# Load data +# load(paste0(fold, "/data/datafreMPTL.RData")) +load(paste0(fold, "/M2/Data Visualisation/tp1", "/data/datafreMPTL.RData")) +paged_table(dat, options = list(rows.print = 15)) +``` + +Le tableau suivant présente une définition des variables. + +```{r} +kableExtra::kable( + data.frame( + Variable = c( + "IDpol", + "Exposure", + "VehPower", + "VehAge", + "DrivAge", + "BonusMalus", + "VehBrand", + "VehGas", + "Area", + "Density", + "Region", + "ClaimTotal", + "ClaimNb" + ), + Description = c( + "Identifiant de la police", + "Exposition au risque", + "Puissance du véhicule", + "Age du véhicule en année", + "Age du conducteur en année", + "Coefficient de bonus-malus", + "Marque du véhicule", + "Carburant du véhicule", + "Catégorie correspondant à la densité de la zone assurée", + "Densité de population", + "Region (selon la classication 1970-2015)", + "Montant total des sinistres", + "Nombre de sinistres sur la période" + ), + Type = c( + rep("Reel", 2), + rep("Entier", 4), + rep("Cat", 3), + "Entier", + "Cat", + rep("Reel", 2) + ) + ), + booktabs = TRUE +) +# Short summary +str(dat) +``` + +Pour plus de détails, consulter l'aide `?CASdatasets::freMTPL2freq`. + +------------------------------------------------------------------------ + +## But de la visualisation + +Nous effectuons une première analyse descriptive de données et cherchons +à étudier la relation entre : + +- la fréquence, calculée avec les variables `ClaimNb` et `Exposure` + (période d'exposition en année). +- les variables `Area` et `DrivAge`. + +Le but de la visualisation est de fait ressortir les liens entre la +fréquence et ces deux variables. + +### Etape 1 : Visualisation de la fréquence et de l'exposition {.unnumbered} + +::: exercise-box +A partir des données `dat` : + +- afficher les statistiques descriptives du nombre de sinistres + `ClaimNb` et de la variable `Exposure` ; +- afficher des histogrammes pour visualiser leur distribution ; +- afficher les figures côte a côte avec la fonction `plot_grid()`. + +Essayer de choisir un thème de couleur et un écartement des barres de +l'histogramme facilitant sa lisibilité. +::: + +::: indice-box +On pourra développer une fonction qui utilise `geom_histogram()` sous la +package **ggplot2**. +::: + +```{r, fig.height = 6, fig.width = 12} +# Descriptive statistics +summary(dat$ClaimNb) +summary(dat$Exposure) + +p1 <- ggplot(dat) + + geom_histogram( + aes(x = ClaimNb), + binwidth = 0.25, + fill = "lightblue", + color = "black" + ) + + labs( + title = "Distribution du nombre de sinistres", + x = "Nombre de sinistres", + y = "Effectif" + ) + + theme_bw() + +p2 <- ggplot(dat) + + geom_histogram( + aes(x = Exposure), + binwidth = 0.05, + fill = "lightblue", + color = "black" + ) + + labs(title = "Exposition", x = "Nombre de sinistres", y = "Effectif") + + theme_bw() + +plot_grid(p1, p2, ncol = 2) +``` + +### Etape 2 : Calculer la fréquence {.unnumbered} + +::: exercise-box +Construire un tableau présentant l’exposition cumulée et le nombre +d’observations avec 0 sinistre, 1 sinistre, … +::: + +```{r} +dat %>% + group_by(ClaimNb) %>% + summarise(n = n(), Exposure = round(sum(Exposure), 0)) %>% + kable( + col_names = c( + "Nombre de sinistres", + "Nombres d'observations", + "Exposition totale" + ) + ) %>% + kable_styling(full_width = F) +``` + +```{r} +pf_freq <- round(sum(dat$ClaimNb) / sum(dat$Exposure), 4) +pf_freq +`` +` +Ce calcul de fréquence sera ensuite utile pour l'affichage des +résultats. + +### Etape 3 : Calculer l'exposition et la fréquence par variable {.unnumbered} + +::: exercise-box +Pour la variable `DrivAge`, présenter : + +1. un histogramme de l'exposition en fonction de cette variable. +2. un histogramme de la fréquence moyenne de sinistres en fonction de + cette variable. + +Remplacer ensuite le second histogramme par un scatter plot avec une +courbe de tendance. Est-ce plus clair ? + +**Indice** + +On pourra développer une fonction qui utilise `geom_bar()` sous la +package **ggplot2**. +::: + +```{r, eval = FALSE} +# On regroupe selon les modalites de la DrivAge +# l'exposition, le nombre de sinistres et la frequence +df_plot <- dat %>% + group_by(DrivAge) %>% + summarize(exp = Exposure, nb_claims = ClaimNb, freq = nb_claims / exp) + +# Histogramme exposition +ggplot(df_plot) + + geom_bar( + aes(x = DrivAge, y = exp), + stat = "identity", + fill = "lightblue", + color = "blue" + ) + + labs( + title = "Exposition par âge du conducteur", + x = "Âge du conducteur", + y = "Exposition" + ) + + theme_minimal() + +# Histogramme frequence +ggplot(df_plot) + + geom_bar( + aes(x = DrivAge, y = freq), + stat = "identity", + fill = "lightblue", + color = "blue" + ) + + labs( + title = "Fréquence par âge du conducteur", + x = "Âge du conducteur", + y = "Fréquence" + ) + + theme_minimal() +``` + +```{r} + +# Scatter plot frequence + +# A compléter + +``` + +### Etape 4 : Examiner l'intéraction avec une autre variable {.unnumbered} + +::: exercise-box +A partir du scatter plot réalisé à l'étape précédente, distinguer les +évolutions de fréquence en fonction de `DrivAge` et de `BonusMalus`. + +Ce graphique vous paraît-il transmettre un message clair ? Proposez des +améliorations en modifiant les variables `DrivAge` et `BonusMalus`. +::: + +```{r} +# On regroupe selon les modalites de la DrivAge et de Area +# l'exposition, le nombre de sinistres et la frequence + +# A compléter + +``` + +On propose 4 ajustements : + +- Exclure les âges extrêmes au-delà de 85 ans pour lesquels + l'exposition est très faible. +- Faire des classes d'âges. +- Limiter le Bonus-Malus à 125. +- Faire des classes de Bonus-Malus. + +```{r} +# Classes d'âges pour Bonus-Malus +lim_classes <- c(50, 75, 100, 125, Inf) + +# Exclusion des donnees "extremes" et faire les regroupement +df_plot <- dat %>% + filter(DrivAge <= 85, BonusMalus <= 125) %>% + # regroupement en classes d'ages de 5 ans + mutate(DrivAge = ceiling(pmin(DrivAge, 85) / 5) * 5) %>% + mutate(BonusMalus = cut(BonusMalus, + breaks = lim_classes, include.lowest = TRUE)) + +# On regroupe selon les modalites de la DrivAge et de Area +# l'exposition, le nombre de sinistres et la frequence + +# A compléter + +``` + +### Conclure {.unnumbered} + +::: exercise-box +Comparer à présenter comment l'exposition se répartie entre âge et +bonus-malus. +::: + +```{r, fig.height = 6, fig.width = 12} + +# A compléter +``` + +### Bonus - Analyse des couples {.unnumbered} + +::: exercise-box +En traitant toutes les variables comme des variables catégorielles, +analyser graphiquement comment évolue la fréquence de sinistres selon +les couples de variables. + +Compléter pour cela la fonction suivante et appliquer la à différents +couples. + +```{r, eval = F} +# Fonction d'analyse bivariée +# df : nom du data.frame +# var1 : nom de la variable explicative 1 +# var2 : nom de la variable explicative 2 +plot_pairwise_disc <- function(df, var1, var2) +{ + df <- rename(df, "varx" = all_of(var1), "vary" = all_of(var2)) + +# replace variable vname by the binning variable + if(is.numeric(df$varx)) + { + df <- df %>% + mutate(varx = ntile(varx, 5)) + } + + if(is.numeric(df$vary)) + { + df <- df %>% + mutate(vary = ntile(vary, 5), + vary = factor(vary)) + } + + df %>% + group_by(??) %>% + summarize(exp = ??, + nb_claims = ??, + freq = ??, + .groups = "drop") %>% + ggplot(aes(x = ??, + y = ??, + colour = ??, + group = vary), + alpha = 0.3) + + geom_point() + geom_line() + theme_bw() + + labs(x = var1, y = "Frequence", colour = var2) +} +` +`` +::: + +# Informations de session {.unnumbered} + +```{r} +sessionInfo() +``` + +# Références \ No newline at end of file diff --git a/M2/Data Visualisation/tp1/3-td_ggplot2---enonce.html b/M2/Data Visualisation/tp1/3-td_ggplot2---enonce.html new file mode 100644 index 0000000..b614d23 --- /dev/null +++ b/M2/Data Visualisation/tp1/3-td_ggplot2---enonce.html @@ -0,0 +1,4910 @@ + + + + + + + + + + + Prise en main de ggplot2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+ + + +
+
+
+ + + +
+
+ +
+ + +
+ + +

Prise en main de ggplot2

+ + + + + + + +
+

1 Objectifs du TP

+

L’objectif de ce TP vise à se familiariser avec le package ggplot2 +de R. Il s’agit de faire des manipulations graphiques élémentaires et +d’interpréter les résultats de ces visualisations.

+

Dans un premier temps, vous pouvez suivre l’exemple introductif en +répliquant le code fourni. Dans un deuxième temps, il convient de +réaliser l’exercice point par point.

+
+
+

2 Prérequis

+
    +
  • Avoir installer R ici.
    • +
  • Avoir installer un IDE, par exemple RStudio +ici.
    • +
  • Créer un nouveau projet (File, puis New Projet) dans un dossier +sur votre ordinateur.
    • +
  • Télécharger ici +les fichiers nécessaires au TD.
    • +
+

Vous pouvez ensuite écrire vos codes soit :

+
    +
  • En ouvrant un nouveau script .R ;
    • +
  • En ouvrant le ouvrant le rapport Rmarkdown 3-td_ggplot2 - enonce. +Certains codes sont partiels et sont à compléter (indication ???). +N’oubliez pas de modifier l’option eval = TRUE pour que les +calculs puissent être réalisés.
    • +
+
+
+

3 Exemple introductif

+

Pour illustrer cette première partie, nous reprenons l’exemple +introductif fourni par @wickham2023 sur le jeu de données penguins du +package palmerpenguins. Ce jeu de données s’intèresse des mesures +réalisées sur des manchots sur 3 îles de l’archipelle Palmer.

+
+

3.1 Données

+

Dans un premier temps, il faut installer le package et le charger.

+
# install.packages("palmerpenguins")
+library(palmerpenguins)
+

Ce jeu de données contient 344 observations où chaque ligne correspond à +un individu.

+
paged_table(penguins, options = list(rows.print = 15))
+
+ +
+

On se concentre plus particulièrement sur les variables suivantes :

+
    +
  • species : l’espèce de manchot ;
    • +
  • flipper_length_mm : la longueur de la nageoire en mm ;
    • +
  • body_mass_g : la masse en gramme.
    • +
+

Pour plus détails, voir l’aide ?penguins.

+
+
+

3.2 But de la visualisation

+

On s’intéresse au lien entre le masse et la taille des nageoires des +manchots :

+
    +
  • ceux dont les nageoires sont les plus longues sont-ils plus lourds +que les manchots aux nageoires courtes ?
    • +
  • si oui quelle est le type de relation (linéaire, croissante, +décroissante, …) ?
    • +
  • quels facteurs influencent également cette relation (lieu, l’espèce, +… ) ?
    • +
+

On cherche à recréer la figure suivante.

+

+
+
+

3.3 Création de la figure étape par étape

+
+

Etape 1 : Scatterplot

+

On commence par créer un scatterplot pour examiner la relation entre la +masse et la taille de la nageoire.

+
ggplot(
+  data = penguins,
+  mapping = aes(x = flipper_length_mm, y = body_mass_g)) +
+  geom_point()
+

+

Cette figure fait clairement apparaître une relation croissante et a +priori linéaire entre les deux variables.

+
+

Un message d’erreur apparaît pour deux individus avec des données +manquantes. Ils sont automatiquement exclus.

+
+
+
+

Etape 2 : Ajout d’élements esthétiques

+

On cherche à présent exhiber le rôle de l’espèce à partir d’une couleur. +Trois espèces sont présents, ainsi l’ajout de 3 couleurs à la figure ne +devrait pas surcharger le graphique.

+
ggplot(
+  data = penguins,
+  mapping = aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
+  geom_point()
+

+

Compte tenu du nombre important de points, nous pouvons renforcer les +différences par espèce en ajoutant une variation de forme aux points.

+
ggplot(
+  data = penguins,
+  mapping = aes(x = flipper_length_mm, y = body_mass_g)) +
+  geom_point(mapping = aes(
+    color = species,
+    shape = species))
+

+
+
+

Etape 3 : Ajout d’une géométrie

+

Voyons à présent comment interpréter la nature de la relation entre +masse et longueur de la nageoire. Pour ce faire, nous essayons d’ajout +des courbes de tendance. Nous commençons par une tendance linéaire.

+
ggplot(
+  data = penguins,
+  mapping = aes(x = flipper_length_mm, y = body_mass_g)) +
+  geom_point(mapping = aes(
+    color = species,
+    shape = species)) +
+    geom_smooth(method = "lm")
+

+

La même figure peut être générée par espèce en déplaçant l’argument +color = species.

+
ggplot(
+  data = penguins,
+  mapping = aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
+  geom_point(mapping = aes(
+    shape = species)) +
+    geom_smooth(method = "lm")
+

+

Les pentes entre les espèces ne sont pas si éloignées. Nous décidons que +conserver une relation commune pour toutes espèces. Pour tester si la +nature linéaire de la relation est a priori une bonne hypothèse, nous +considérons un lissage non-paramétrique.

+
ggplot(
+  data = penguins,
+  mapping = aes(x = flipper_length_mm, y = body_mass_g)) +
+  geom_point(mapping = aes(
+    color = species,
+    shape = species)) +
+    geom_smooth(method = "loess")
+

+

L’ajout d’un lissage non-paramétrique permet d’affiner l’adéquation aux +données, mais sans pour autant clairement remettre en cause la tendance +linéaire qui sera donc conservée.

+
+
+

Etape 4 : Ajout des titres et changement de thème

+

Afin de finaliser la figure, nous ajouter :

+
    +
  • un titre ;
    • +
  • un sous-titre ;
    • +
  • des titres aux axes ;
    • +
  • un titre à la légende.
    • +
+

Ces informations sont ajoutées avec labs().

+

De plus, nous modifions le thème avec la commande theme_bw().

+
ggplot(
+  data = penguins,
+  mapping = aes(x = flipper_length_mm, y = body_mass_g)
+) +
+  geom_point(aes(color = species, shape = species)) +
+  geom_smooth(method = "lm") +
+  labs(
+    title = "Masse et taille de la nageoire",
+    subtitle = "Manchots d'Adelie, a
+    jugulaire et de Gentoo",
+    x = "Longueur de la nageoire (mm)", y = "Masse (g)",
+    color = "Espece", shape = "Espece"
+  ) +
+  scale_color_colorblind() +
+  theme_bw()
+

+
+
+
+
+
+

4 Exercice

+
+

4.1 Données

+

Nous travaillons avec les jeux de données FreMTPL2freq et +FreMTPL2sev du package Casdatasets. Ces données ont été +préalablement pré-formatées et regroupées.

+

Ce jeux de données regroupent les caractéristiques de 677 991 polices de +responsabilité civile automobile, observées principalement sur une +année. Dans les données regroupées, on dispose des numéros de sinistre +par police, des montants de sinistre correspondants, des +caractéristiques du risque et du nombre de sinistres.

+

On présente ci-dessous un aperçu des données.

+
# Folds
+fold <- getwd()
+
+# Load data
+load(paste0(fold, "/data/datafreMPTL.RData"))
+paged_table(dat, options = list(rows.print = 15))
+
+ +
+

Le tableau suivant présente une définition des variables.

+
kableExtra::kable(
+  data.frame(
+    Variable = c("IDpol", "Exposure", "VehPower", "VehAge",
+                 "DrivAge", "BonusMalus", "VehBrand", "VehGas",
+                 "Area", "Density", "Region", 
+                 "ClaimTotal", "ClaimNb"),
+    Description = c(
+      "Identifiant de la police",
+      "Exposition au risque",
+      "Puissance du véhicule",
+      "Age du véhicule en année",
+      "Age du conducteur en année",
+      "Coefficient de bonus-malus",
+      "Marque du véhicule",
+      "Carburant du véhicule",
+      "Catégorie correspondant à la densité de la zone assurée",
+      "Densité de population",
+      "Region (selon la classication 1970-2015)",
+      "Montant total des sinistres",
+      "Nombre de sinistres sur la période"
+    ), 
+    Type = c(
+      rep("Reel", 2),
+      rep("Entier", 4),
+      rep("Cat", 3),
+      "Entier",
+      "Cat",
+      rep("Reel", 2)
+    )
+  ),
+  booktabs=TRUE)
+ +++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
VariableDescriptionType
IDpolIdentifiant de la policeReel
ExposureExposition au risqueReel
VehPowerPuissance du véhiculeEntier
VehAgeAge du véhicule en annéeEntier
DrivAgeAge du conducteur en annéeEntier
BonusMalusCoefficient de bonus-malusEntier
VehBrandMarque du véhiculeCat
VehGasCarburant du véhiculeCat
AreaCatégorie correspondant à la densité de la zone assuréeCat
DensityDensité de populationEntier
RegionRegion (selon la classication 1970-2015)Cat
ClaimTotalMontant total des sinistresReel
ClaimNbNombre de sinistres sur la périodeReel
+
# Short summary
+str(dat)
+
## 'data.frame':    135601 obs. of  13 variables:
+##  $ IDpol     : num  5089457 2215917 2202141 1077071 3044896 ...
+##  $ Exposure  : num  0.21 0.25 1 0.07 0.09 0.08 1 1 0.49 0.8 ...
+##  $ VehPower  : int  7 5 7 7 7 7 13 4 7 4 ...
+##  $ VehAge    : int  10 2 6 4 7 4 3 14 8 2 ...
+##  $ DrivAge   : int  36 40 79 42 45 53 53 32 31 42 ...
+##  $ BonusMalus: int  50 85 54 57 100 53 50 60 71 50 ...
+##  $ VehBrand  : Factor w/ 11 levels "B1","B2","B3",..: 1 3 10 2 7 2 1 2 3 9 ...
+##  $ VehGas    : Factor w/ 2 levels "Diesel","Regular": 2 2 1 1 1 1 2 2 2 2 ...
+##  $ Area      : Factor w/ 6 levels "A","B","C","D",..: 5 5 3 5 6 3 5 4 5 2 ...
+##  $ Density   : int  8 8 6 8 10 5 8 7 9 4 ...
+##  $ Region    : Factor w/ 21 levels "Alsace","Aquitaine",..: 21 7 7 21 12 20 7 16 21 7 ...
+##  $ ClaimTotal: num  0 0 0 0 0 0 0 0 0 0 ...
+##  $ ClaimNb   : num  0 0 0 0 0 0 0 0 0 0 ...
+

Pour plus de détails, consulter l’aide ?CASdatasets::freMTPL2freq.

+
+
+
+

4.2 But de la visualisation

+

Nous effectuons une première analyse descriptive de données et cherchons +à étudier la relation entre :

+
    +
  • la fréquence, calculée avec les variables ClaimNb et Exposure +(période d’exposition en année).
    • +
  • les variables Area et DrivAge.
    • +
+

Le but de la visualisation est de fait ressortir les liens entre la +fréquence et ces deux variables.

+
+

Etape 1 : Visualisation de la fréquence et de l’exposition

+
+

A partir des données dat :

+
    +
  • afficher les statistiques descriptives du nombre de sinistres +ClaimNb et de la variable Exposure ;
    • +
  • afficher des histogrammes pour visualiser leur distribution ;
    • +
  • afficher les figures côte a côte avec la fonction plot_grid().
    • +
+

Essayer de choisir un thème de couleur et un écartement des barres de +l’histogramme facilitant sa lisibilité.

+
+
+

On pourra développer une fonction qui utilise geom_histogram() sous la +package ggplot2.

+
+
# Descriptive statistics
+summary(dat$ClaimNb)
+
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
+## 0.00000 0.00000 0.00000 0.03919 0.00000 3.00000
+
summary(dat$Exposure)
+
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
+## 0.002732 0.170000 0.490000 0.527782 0.990000 1.000000
+
p1 <- ggplot(dat) +
+  geom_histogram(aes(x = ClaimNb), binwidth = 1, fill = "lightblue", color = "black") +
+  labs(title = "Distribution du nombre de sinistres", x = "Nombre de sinistres", y = "Effectif") +
+  theme_minimal()
+
+p2 <- ggplot(dat) +
+  geom_histogram(aes(x = Exposure), binwidth = 0.05, fill = "lightblue", color = "black") +
+  labs(title = "Distribution du nombre de sinistres", x = "Nombre de sinistres", y = "Effectif") +
+  theme_minimal()
+
+plot_grid(p1, p2, ncol = 2)
+

+
+
+

Etape 2 : Calculer la fréquence

+
+

Construire un tableau présentant l’exposition cumulée et le nombre +d’observations avec 0 sinistre, 1 sinistre, …

+
+
freq <- dat %>% 
+  summarize(total_exposure = sum(Exposure),
+            total_claims = sum(ClaimNb),
+            frequency = total_claims / total_exposure)
+
+freq
+
##   total_exposure total_claims  frequency
+## 1       71567.74         5314 0.07425133
+

Ce calcul de fréquence sera ensuite utile pour l’affichage des +résultats.

+
+
+

Etape 3 : Calculer l’exposition et la fréquence par variable

+
+

Pour la variable DrivAge, présenter :

+
    +
  1. un histogramme de l’exposition en fonction de cette variable.
    2. +
  2. un histogramme de la fréquence moyenne de sinistres en fonction de +cette variable.
    3. +
+

Remplacer ensuite le second histogramme par un scatter plot avec une +courbe de tendance. Est-ce plus clair ?

+

Indice

+

On pourra développer une fonction qui utilise geom_bar() sous la +package ggplot2.

+
+
# On regroupe selon les modalites de la DrivAge 
+# l'exposition, le nombre de sinistres et la frequence 
+df_plot <- dat %>%
+    group_by(DrivAge) %>% 
+    summarize(exp = Exposure,
+              nb_claims = ClaimNb,
+              freq = nb_claims / exp)
+
+# Histogramme exposition
+ggplot(df_plot) +
+  geom_bar(aes(x = DrivAge, y = exp), stat = "identity", fill = "lightblue", color = "blue") +
+  labs(title = "Exposition par âge du conducteur", x = "Âge du conducteur", y = "Exposition") +
+  theme_minimal()
+
+# Histogramme frequence
+ggplot(df_plot) +
+  geom_bar(aes(x = DrivAge, y = freq), stat = "identity", fill = "lightblue", color = "blue") +
+  labs(title = "Fréquence par âge du conducteur", x = "Âge du conducteur", y = "Fréquence") +
+  theme_minimal()
+
# Scatter plot frequence
+
+# A compléter
+
+
+

Etape 4 : Examiner l’intéraction avec une autre variable

+
+

A partir du scatter plot réalisé à l’étape précédente, distinguer les +évolutions de fréquence en fonction de DrivAge et de BonusMalus.

+

Ce graphique vous paraît-il transmettre un message clair ? Proposez des +améliorations en modifiant les variables DrivAge et BonusMalus.

+
+
# On regroupe selon les modalites de la DrivAge et de Area
+# l'exposition, le nombre de sinistres et la frequence 
+
+# A compléter
+

On propose 4 ajustements :

+
    +
  • Exclure les âges extrêmes au-delà de 85 ans pour lesquels +l’exposition est très faible.
    • +
  • Faire des classes d’âges.
    • +
  • Limiter le Bonus-Malus à 125.
    • +
  • Faire des classes de Bonus-Malus.
    • +
+
# Classes d'âges pour Bonus-Malus
+lim_classes <- c(50, 75, 100, 125, Inf)
+
+# Exclusion des donnees "extremes" et faire les regroupement
+df_plot <- dat %>%
+  filter(DrivAge <= 85, BonusMalus <= 125) %>%
+  # regroupement en classes d'ages de 5 ans
+  mutate(DrivAge = ceiling(pmin(DrivAge, 85) / 5) * 5) %>%
+  mutate(BonusMalus = cut(BonusMalus, 
+                          breaks = lim_classes, include.lowest = TRUE))
+  
+# On regroupe selon les modalites de la DrivAge et de Area
+# l'exposition, le nombre de sinistres et la frequence 
+
+# A compléter
+
+
+

Conclure

+
+

Comparer à présenter comment l’exposition se répartie entre âge et +bonus-malus.

+
+
# A compléter
+
+
+

Bonus - Analyse des couples

+
+

En traitant toutes les variables comme des variables catégorielles, +analyser graphiquement comment évolue la fréquence de sinistres selon +les couples de variables.

+

Compléter pour cela la fonction suivante et appliquer la à différents +couples.

+
# Fonction d'analyse bivariée
+# df : nom du data.frame
+# var1 : nom de la variable explicative 1
+# var2 : nom de la variable explicative 2
+plot_pairwise_disc <- function(df, var1, var2)
+{
+ df <- rename(df, "varx" = all_of(var1), "vary" = all_of(var2))
+ 
+#  replace variable vname by the binning variable
+  if(is.numeric(df$varx))
+  {
+    df <- df %>%
+      mutate(varx = ntile(varx, 5))
+  }
+ 
+ if(is.numeric(df$vary))
+  {
+    df <- df %>%
+      mutate(vary = ntile(vary, 5),
+             vary = factor(vary))
+  }
+ 
+ df %>% 
+   group_by(??) %>%
+   summarize(exp = ??,
+             nb_claims = ??,
+             freq = ??,
+             .groups = "drop") %>%
+   ggplot(aes(x = ??, 
+              y = ??,
+              colour = ??,
+              group = vary),
+          alpha = 0.3) + 
+   geom_point() + geom_line() + theme_bw() +
+   labs(x = var1,  y = "Frequence", colour = var2)
+}
+
+
+
+
+
+

Informations de session

+
sessionInfo()
+
## R version 4.5.1 (2025-06-13)
+## Platform: aarch64-apple-darwin20
+## Running under: macOS Tahoe 26.0.1
+## 
+## Matrix products: default
+## BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
+## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
+## 
+## locale:
+## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
+## 
+## time zone: Europe/Paris
+## tzcode source: internal
+## 
+## attached base packages:
+## [1] grid      stats     graphics  grDevices utils     datasets  methods  
+## [8] base     
+## 
+## other attached packages:
+##  [1] palmerpenguins_0.1.1 kableExtra_1.4.0     cowplot_1.2.0       
+##  [4] ggthemes_5.1.0       rmarkdown_2.30       tidyr_1.3.1         
+##  [7] dplyr_1.1.4          plotly_4.11.0        RColorBrewer_1.1-3  
+## [10] formattable_0.2.1    scales_1.4.0         locfit_1.5-9.12     
+## [13] gridExtra_2.3        ggplot2_4.0.0        lattice_0.22-7      
+## 
+## loaded via a namespace (and not attached):
+##  [1] sass_0.4.10       generics_0.1.4    xml2_1.4.0        stringi_1.8.7    
+##  [5] digest_0.6.37     magrittr_2.0.4    evaluate_1.0.5    bookdown_0.45    
+##  [9] fastmap_1.2.0     Matrix_1.7-4      jsonlite_2.0.0    mgcv_1.9-3       
+## [13] httr_1.4.7        purrr_1.1.0       viridisLite_0.4.2 lazyeval_0.2.2   
+## [17] textshaping_1.0.3 jquerylib_0.1.4   cli_3.6.5         rlang_1.1.6      
+## [21] splines_4.5.1     withr_3.0.2       cachem_1.1.0      yaml_2.3.10      
+## [25] tools_4.5.1       vctrs_0.6.5       R6_2.6.1          lifecycle_1.0.4  
+## [29] stringr_1.5.2     htmlwidgets_1.6.4 pkgconfig_2.0.3   pillar_1.11.1    
+## [33] bslib_0.9.0       gtable_0.3.6      data.table_1.17.8 glue_1.8.0       
+## [37] rmdformats_1.0.4  systemfonts_1.3.1 xfun_0.53         tibble_3.3.0     
+## [41] tidyselect_1.2.1  rstudioapi_0.17.1 knitr_1.50        farver_2.1.2     
+## [45] nlme_3.1-168      htmltools_0.5.8.1 labeling_0.4.3    svglite_2.2.1    
+## [49] compiler_4.5.1    S7_0.2.0
+
+
+

5 Références

+
+ + + +
+
+
+
+
+ + + + + + + + + + + + diff --git a/M2/Data Visualisation/tp1/data/datafreMPTL.RData b/M2/Data Visualisation/tp1/data/datafreMPTL.RData new file mode 100644 index 0000000..02d0d1f Binary files /dev/null and b/M2/Data Visualisation/tp1/data/datafreMPTL.RData differ