J’ai déjà publié, mais il y a longtemps, des articles sur la description des données. Ce que je vous propose dans ce nouvel article, ce sont les fonctions que j’utilise, actuellement, en routine pour :
Pour illustrer cet article, nous allons employer une version allégée du jeu de données heart_disease
qui est inclus dans le package funModeling
.
Pour cela, j’ai sélectionné uniquement quelques variables quantitatives et qualitatives. Et j’ai remplacé quelques valeurs par des données manquantes, afin de vérifier comment celles-ci sont prises en considération.
Vous pouvez créer ces données en employant le code ci-dessous :
#install.packages("funModeling")
library(funModeling)
library(tidyverse)
myHD <- heart_disease %>%
select(age:chest_pain,max_heart_rate , has_heart_disease)
str(myHD )
# remplacement de 30 valeurs d'âge par des NA
set.seed(1)
ind_age_na <- sample(303,30)
myHD$age[ind_age_na] <- NA
# remplacement de 30 valeurs de gender par des NA
set.seed(2)
ind_gender_na <- sample(303,10)
myHD$gender[ind_gender_na] <- NA
# remplacement de 30 valeurs de chest_pain par des NA
set.seed(3)
ind_chest_pain_na <- sample(303,50)
myHD$chest_pain[ind_chest_pain_na] <- NA
# remplacement de 30 valeurs de max_heart_rate par des NA
set.seed(4)
ind_max_heart_rate_na <- sample(303,30)
myHD$max_heart_rate[ind_max_heart_rate_na] <- NA
J’ai récemment découvert la fonction aggr()
du package VIM
. Cette fonction renvoie :
#install.packages("funModeling")
#install.packages("VIM")
library(funModeling)
library(VIM)
Etude_NA <- aggr(myHD,
col=c('navyblue','red'),
numbers=TRUE,
sortVars=TRUE,
labels=names(data),
cex.axis=.7, gap=3,
ylab=c("Histogram of missing data","Pattern"))
## Variables sorted by number of missings:
## Variable Count
## chest_pain 0.1650165
## age 0.0990099
## max_heart_rate 0.0990099
## gender 0.0330033
## has_heart_disease 0.0000000
Attention à laisser de la place à la fenêtre graphique pour obtenir les pourcentages sur la partie droite du graphique de droite.
Les colonnes représentent les variables, les lignes représentent les types de situations dans le data frame :
Si le sujet des données manquantes vous intéresse, vous pouvez aussi consulter mon article “Visualisation des données manquantes” :
La fonction qui a ma préférence, en ce moment, c’est la fonction descr()
du package summarytools
, parce qu’elle renvoie tous les paramètres réellement nécessaires (moyene, écart-type, min, quartile1, médiane, quartile3, max, MAD (Median Absolute Deviation), Skewness( pour évaluer si la distribution est symétrique ou pas), le nombre et le pourcentage de données non manquantes (N.Valid, et Pct.Valid)).
Il n’y a que le kurtosis (reflète l’aplatissement de la distribution) que je n’emploie pas. Et puis ce que je trouve très pratique, c’est qu’il n’est pas nécessaire de faire un subset contenant uniquement les variables quantitatives, la fonction fait le tri elle-même !
library(summarytools)
summarytools::descr(myHD)
## Descriptive Statistics
## myHD
## N: 303
##
## age max_heart_rate
## ----------------- -------- ----------------
## Mean 54.82 149.92
## Std.Dev 9.09 21.72
## Min 29.00 90.00
## Q1 49.00 136.00
## Median 56.00 153.00
## Q3 61.00 165.00
## Max 77.00 195.00
## MAD 8.90 22.24
## IQR 12.00 29.00
## CV 0.17 0.14
## Skewness -0.29 -0.46
## SE.Skewness 0.15 0.15
## Kurtosis -0.44 -0.36
## N.Valid 273.00 273.00
## Pct.Valid 90.10 90.10
On n’obtient pas directement le pourcentage de données manquantes, puisque c’est le pourcentage de données non manquantes qui est renvoyé. Mais dans tous les cas, celles-ci sont bien prises en considération.
Vous pouvez transposer la table de résultats pour avoir les variables en ligne et les paramètres statistiques en colonne, avec l’argument transpose=TRUE
, et supprimer les lignes d’entête (attention vous perdrez l’info sur le nombre total des données) avec l’argumentheadings=FALSE
, comme ceci :
summarytools::descr(myHD,
heading=FALSE,
transpose =TRUE)
##
## Mean Std.Dev Min Q1 Median Q3 Max MAD IQR
## -------------------- -------- --------- ------- -------- -------- -------- -------- ------- -------
## age 54.82 9.09 29.00 49.00 56.00 61.00 77.00 8.90 12.00
## max_heart_rate 149.92 21.72 90.00 136.00 153.00 165.00 195.00 22.24 29.00
##
## Table: Table continues below
##
##
##
## CV Skewness SE.Skewness Kurtosis N.Valid Pct.Valid
## -------------------- ------ ---------- ------------- ---------- --------- -----------
## age 0.17 -0.29 0.15 -0.44 273.00 90.10
## max_heart_rate 0.14 -0.46 0.15 -0.36 273.00 90.10
Si vous ne parvenez pas à utiliser, je vous recommande de mettre à jour le package summaytools et / ou R.
Et si vous voulez une description par groupe, vous pouvez le faire à l’aide d’une fonction by()
. Par exemple, pour faire une description par modalité de la variable gender, c’est-à-dire pour les hommes et les femmes séparément :
by(myHD,heart_disease$gender, summarytools::descr,heading=FALSE,transpose =TRUE)
## heart_disease$gender: female
## Group: heart_disease$gender = female
##
## Mean Std.Dev Min Q1 Median Q3 Max MAD IQR
## -------------------- -------- --------- ------- -------- -------- -------- -------- ------- -------
## age 55.69 9.49 34.00 50.00 57.00 63.00 76.00 8.90 13.00
## max_heart_rate 151.64 19.39 97.00 142.00 157.00 165.00 192.00 17.79 23.00
##
## Table: Table continues below
##
##
##
## CV Skewness SE.Skewness Kurtosis N.Valid Pct.Valid
## -------------------- ------ ---------- ------------- ---------- --------- -----------
## age 0.17 -0.28 0.25 -0.64 90.00 92.78
## max_heart_rate 0.13 -0.71 0.26 -0.12 85.00 87.63
## ------------------------------------------------------------
## heart_disease$gender: male
## Group: heart_disease$gender = male
##
## Mean Std.Dev Min Q1 Median Q3 Max MAD IQR
## -------------------- -------- --------- ------- -------- -------- -------- -------- ------- -------
## age 54.39 8.88 29.00 48.00 56.00 60.00 77.00 7.41 12.00
## max_heart_rate 149.15 22.71 90.00 132.00 151.00 165.50 195.00 25.20 33.25
##
## Table: Table continues below
##
##
##
## CV Skewness SE.Skewness Kurtosis N.Valid Pct.Valid
## -------------------- ------ ---------- ------------- ---------- --------- -----------
## age 0.16 -0.32 0.18 -0.37 183.00 88.83
## max_heart_rate 0.15 -0.36 0.18 -0.48 188.00 91.26
La fonction by()
, renvoie les paramètres descriptifs dans une structure de données de type list. Il est alors très facile d’exporter la description pour chaque groupe, dans une feuille séparée d’un même classeur Excel, en employant la fonction write_xlsx
du package writexl
, comme ceci :
library(writexl)
# création d'un dossier resultats à la racine du projet R
dir.create("resultats")
# écriture des description dans des feuilles séparées d'un même classeur
myHD_QT_by_gender_ls <- by(myHD,heart_disease$gender, summarytools::descr,heading=FALSE,transpose =TRUE)
myHD_QT_by_gender_ls %>%
writexl::write_xlsx(path = "resultats/myHD_descri_QT_gender .xlsx")
Pour les variables qualitatives, la fonction qui a ma préférence est freq()
du package summarytools,
également
. La fonction freq()
du package funModeling
est très similaire, et elle propose en plus la réalisation de graphiques, ce qui peut être très pratique.
Néanmoins, l’affichage de la fonction du package funModeling
est un peu étrange et rend difficile la distinction des groupes. Par contre, dans les deux cas, l’exportation ne peut, à ma connaissance, se faire uniquement dans un fichier word :
# summarytools
summarytools::freq(myHD)
## Frequencies
## myHD$gender
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ------------ ------ --------- -------------- --------- --------------
## female 93 31.74 31.74 30.69 30.69
## male 200 68.26 100.00 66.01 96.70
## <NA> 10 3.30 100.00
## Total 303 100.00 100.00 100.00 100.00
##
## myHD$chest_pain
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 1 18 7.11 7.11 5.94 5.94
## 2 43 17.00 24.11 14.19 20.13
## 3 70 27.67 51.78 23.10 43.23
## 4 122 48.22 100.00 40.26 83.50
## <NA> 50 16.50 100.00
## Total 303 100.00 100.00 100.00 100.00
##
## myHD$has_heart_disease
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## no 164 54.13 54.13 54.13 54.13
## yes 139 45.87 100.00 45.87 100.00
## <NA> 0 0.00 100.00
## Total 303 100.00 100.00 100.00 100.00
On peut voir que les données manquantes sont bien prises en considération, dans le calcul des pourcentages et pourcentages cumulés.
Remarque : Si vous utilisez la fonction freq()
du package funModeling
, vous pouvez faire défiler tous les plots en employant la flèche de gauche de la fenêtre de plot. Pour ne pas afficher les plots vous pouvez employer l’argument plot=FALSE
dans la fonction freq()
# avec funModeling
funModeling::freq(myHD, plot=TRUE)
gender frequency percentage cumulative_perc
1 male 200 66.01 66.01
2 female 93 30.69 96.70
3 <NA> 10 3.30 100.00
chest_pain frequency percentage cumulative_perc
1 4 122 40.26 40.26
2 3 70 23.10 63.36
3 <NA> 50 16.50 79.86
4 2 43 14.19 94.05
5 1 18 5.94 100.00
has_heart_disease frequency percentage
1 no 164 54.13
2 yes 139 45.87
cumulative_perc
1 54.13
2 100.00
Si vous souhaitez réaliser une description par groupe, vous pouvez, comme précédemment, employer une fonction by() :
by(myHD,myHD$gender, summarytools::freq)
## myHD$gender: female
## Frequencies
## myHD$gender
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ------------ ------ --------- -------------- --------- --------------
## female 93 100.00 100.00 100.00 100.00
## male 0 0.00 100.00 0.00 100.00
## <NA> 0 0.00 100.00
## Total 93 100.00 100.00 100.00 100.00
##
## myHD$chest_pain
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 1 3 3.90 3.90 3.23 3.23
## 2 12 15.58 19.48 12.90 16.13
## 3 29 37.66 57.14 31.18 47.31
## 4 33 42.86 100.00 35.48 82.80
## <NA> 16 17.20 100.00
## Total 93 100.00 100.00 100.00 100.00
##
## myHD$has_heart_disease
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## no 69 74.19 74.19 74.19 74.19
## yes 24 25.81 100.00 25.81 100.00
## <NA> 0 0.00 100.00
## Total 93 100.00 100.00 100.00 100.00
## ------------------------------------------------------------
## myHD$gender: male
## Frequencies
## myHD$gender
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ------------ ------ --------- -------------- --------- --------------
## female 0 0.00 0.00 0.00 0.00
## male 200 100.00 100.00 100.00 100.00
## <NA> 0 0.00 100.00
## Total 200 100.00 100.00 100.00 100.00
##
## myHD$chest_pain
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## 1 15 8.82 8.82 7.50 7.50
## 2 30 17.65 26.47 15.00 22.50
## 3 40 23.53 50.00 20.00 42.50
## 4 85 50.00 100.00 42.50 85.00
## <NA> 30 15.00 100.00
## Total 200 100.00 100.00 100.00 100.00
##
## myHD$has_heart_disease
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## no 91 45.50 45.50 45.50 45.50
## yes 109 54.50 100.00 54.50 100.00
## <NA> 0 0.00 100.00
## Total 200 100.00 100.00 100.00 100.00
Et pour l’exporter :
myHD_QL_by_gender <- by(myHD,myHD$gender, summarytools::freq)
capture.output(myHD_QL_by_gender, file="resultats/myHD_descri_QL_gender.doc")
J’ai déjà publié plusieurs articles concernant la réalisation de table 1 avec les packages tableone et Table1 :
Ici, je vais vous parler du package gtsummary()
qui est aussi très performant et facile à employer. Les sorties de ce package sont, par défaut, au format html. Pour les exporter au format word, il est nécessaire d’installer le package flextable
:
#install.packages("gtsummary")
library(gtsummary)
#install.packages("flextable")
library(flextable)
myHD %>%
tbl_summary()
Par défaut, les variables quantitatives, sont décrites par la médiane et l’intervalle interquartile (IQR).
Une règle parfois utilisée consiste à employer la médiane et l’ IQR lorsque la distribution n’est pas symétrique (skewness <-1 ou >1), et la moyenne et l’écart type si la distribution est symétrique.
Il est possible de spécifier d’autres paramètres (moyenne et écart type, par exemple), en utilisant l’argument statistic
, comme ceci :
myHD%>%
tbl_summary(
statistic = list(all_continuous() ~ "{mean} ({sd})"))
myHD%>%
tbl_summary(
statistic = list(c(age, max_heart_rate) ~ "{mean} ({sd})"))
Il est encore possible de choisir spécifiquement les variables décrites à l’aide de la moyenne et de l’écart type :
Le nombre de décimales peut être modifié à l’aide de l’argument digits :
myHD%>%
tbl_summary(
statistic = list(all_continuous() ~ "{mean} ({sd})"),
digits=all_continuous()~2)
Pour réaliser une table stratifiée, on emploie l’argument by()
:
myHD%>%
tbl_summary(by=gender,
statistic = list(all_continuous() ~ "{mean} ({sd})"),
digits=all_continuous()~2)
Pour ajouter un test statistique, il suffit d’ajouter la fonction add_p()
en employant le pipe (%>%) au préalable :
myHD%>%
tbl_summary(by=gender,
statistic = list(all_continuous() ~ "{mean} ({sd})"),
digits=all_continuous()~2) %>%
add_p()
Par défaut les tests relatifs aux variables quantitatives sont des tests non paramétriques. Il est possible d’utiliser des tests paramétriques. Comme ceci, par exemple :
myHD%>%
tbl_summary(by=gender,
statistic = list(c(age, max_heart_rate) ~ "{mean} ({sd})"),
digits=all_continuous()~2) %>%
add_p(
# réalisation de tests t pour toutes les variables
test = everything() ~ "t.test",
# hypothèse d'égalité des variances pour obtenir un test t (sinon test de Welch)
test.args = all_tests("t.test") ~ list(var.equal = TRUE)
)
Pour exporter, on utilise en combinaison les fonctions as_flex_table()
et save_as_doxc() :
myHD %>%
tbl_summary(by=gender) %>%
add_p() %>%
as_flex_table() %>%
flextable::save_as_docx(path="resultats/myHD_table_descr.docx")
gtsummary
et la création de tables descriptives synthétiques dans sa vignette.ggplot2
est indispensable. Vous trouverez une initiation dans mon article “Introduction à la visualisation sous R avec le package ggplot2“.ggplot2
, je vous recommande d’employer le package esquisse
qui vous permettra de réaliser vos premiers graphiques ultra facilement en drag and drop et de les exporter en jpeg ou en powerpoint. POur en savoir plus sur l’utilisation du package esquisse vous puvez consulter l’article suivant J’espère que cet article vous aidera à réaliser plus facilement vos analyses descriptives et à les exporter.
Je pense que ça vaut le coup d’ajouter le code pour exporter la description des variables quantitatives dans votre boite à outils.
Si vous utilisez d’autres fonctions en routine qui vous paraissent particulièrement efficaces et dont vous êtes très satisfait, s’il vous plait, partagez les en commentaire.
Et si cet article vous a plu, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee.Si cet article vous a plu, ou vous a été utile, et si vous le souhaitez, vous pouvez soutenir ce blog en faisant un don sur sa page Tipeee
Enregistrez vous pour recevoir gratuitement mes fiches “aide mémoire” (ou cheat sheets) qui vous permettront de réaliser facilement les principales analyses biostatistiques avec le logiciel R et pour être informés des mises à jour du site.
9 Responses
Merci Claire pour ce passionant article, je prends plaisir à te lire pour découvrir chaque fois de nouveaux packages et de nouvelles solutions.
La boite à outil est une très bonne idée, car je ne compte pas le nombre de fois ou j’ai répété les commandes de prise en main d’une base. Cependant, ne serait-ce pas plus intéressant de coder nos propres fonctions ?
Grâce aux fonctions : eval(parse(text = x)), il est possible d’évaluer du texte (en x) pour en faire du code, ce que permettrait malgré l’appel à une fonction de modifier le code dans la fonction (exemple : ajout d’une colonne Overall dans une table1, même si la fonction n’est initialement pas codée en ce sens).
Je serais inétressé d’avoir ton avis Claire, et pourquoi pas celui de la communauté.
Bonjour Pierre Etienne,
je programme très très peu, je n’utilise donc jamais ce genre de chose. Pouvez-vous nous faire un exemple reproductible ?
Merci
Merci encore une fois de plus Claire, je l’adore ton article, il n’est pas long mais très riche. Découverte de nouvelles fonctions encore plus intéressante
Super interessant . merci Claire
Merci Omar !
Merci pour Claire. c’est vraiment passionnant de lire votre article qui est synthétique mais très claire et riche. félicitation.
Bonjour Claire et merci pour tes aides.
Dites, vous n’avez pas de document ou lien de vos recherche me permettant de faire une verification de Données ?
en resumé verification de la qualité de donnée (Heigh Frequencies Checks).
merci encore.
Bonjour Alassane,
je ne connais pas du tout ce domaine, mais voici ce que j’ai trouvé :
“Supposons que vous ayez une série chronologique de données financières et que vous souhaitiez vérifier si elle contient des hautes fréquences.
Tout d’abord, vous pouvez utiliser la fonction spec.pgram() pour estimer la densité spectrale de puissance de votre série chronologique. Cela peut être fait en utilisant le code suivant :
library(tseries)
data <- read.csv("finance_data.csv") # importer vos données spec <- spec.pgram(data$returns) # calculer la densité spectrale de puissance des rendements plot(spec$freq, spec$spec, type="l", xlab="Frequency", ylab="Power")
Ce code importe d'abord vos données financières à partir d'un fichier CSV, puis utilise la fonction spec.pgram() pour estimer la densité spectrale de puissance des rendements de la série chronologique. La fonction plot() est ensuite utilisée pour tracer le résultat.
La figure résultante montrera une courbe de densité spectrale de puissance avec une fréquence sur l'axe x et une puissance sur l'axe y. Si la courbe montre une augmentation de la puissance à des fréquences élevées, cela peut indiquer la présence de hautes fréquences dans la série chronologique.
En outre, vous pouvez également utiliser la fonction acf() pour tracer la fonction d'autocorrélation de la série chronologique, comme suit :
acf(data$returns, lag.max=100, main="ACF of Returns")
La fonction acf() calcule la fonction d'autocorrélation de la série chronologique et la fonction lag.max spécifie le nombre maximal de décalages de temps à inclure dans le graphique. Si la fonction d'autocorrélation décroit rapidement, cela peut également indiquer la présence de hautes fréquences.
Ces deux approches sont couramment utilisées pour vérifier la présence de hautes fréquences dans les séries chronologiques financières."
J'espère que cela vous aide.
Bonne continuation