L'idée était de vous faire mettre en forme un jeu de données vous-même, avant d'effectuer une ACP. En effet l'ACP directe n'est pas possible car il y a trop de lignes incomplètes, des colonnes a priori peu pertinentes, et des données de type séries temporelles (que l'on ramènera à une seule valeur).
# Import du jeu de données brut
data <- read.csv("https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/owid-covid-data.csv")
# Affichage d'un aperçu du jeu de données ; fonctions head() et summary() = bons réflexes :)
dim(data)
head(data)
| iso_code | continent | location | date | total_cases | new_cases | new_cases_smoothed | total_deaths | new_deaths | new_deaths_smoothed | ⋯ | female_smokers | male_smokers | handwashing_facilities | hospital_beds_per_thousand | life_expectancy | human_development_index | excess_mortality_cumulative_absolute | excess_mortality_cumulative | excess_mortality | excess_mortality_cumulative_per_million | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | ⋯ | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1 | AFG | Asia | Afghanistan | 2020-02-24 | 5 | 5 | NA | NA | NA | NA | ⋯ | NA | NA | 37.746 | 0.5 | 64.83 | 0.511 | NA | NA | NA | NA |
| 2 | AFG | Asia | Afghanistan | 2020-02-25 | 5 | 0 | NA | NA | NA | NA | ⋯ | NA | NA | 37.746 | 0.5 | 64.83 | 0.511 | NA | NA | NA | NA |
| 3 | AFG | Asia | Afghanistan | 2020-02-26 | 5 | 0 | NA | NA | NA | NA | ⋯ | NA | NA | 37.746 | 0.5 | 64.83 | 0.511 | NA | NA | NA | NA |
| 4 | AFG | Asia | Afghanistan | 2020-02-27 | 5 | 0 | NA | NA | NA | NA | ⋯ | NA | NA | 37.746 | 0.5 | 64.83 | 0.511 | NA | NA | NA | NA |
| 5 | AFG | Asia | Afghanistan | 2020-02-28 | 5 | 0 | NA | NA | NA | NA | ⋯ | NA | NA | 37.746 | 0.5 | 64.83 | 0.511 | NA | NA | NA | NA |
| 6 | AFG | Asia | Afghanistan | 2020-02-29 | 5 | 0 | 0.714 | NA | NA | 0 | ⋯ | NA | NA | 37.746 | 0.5 | 64.83 | 0.511 | NA | NA | NA | NA |
# On ne verra pas de méthodes s'appliquant sur des séries temporelles dans ce cours,
# donc on se limite ici à des statistiques globales prises au dernier jour :
dmax <- max(data$date) #28/11/2021 aujourd'hui, 29/11 demain etc...
data <- data[data$date == dmax,] #filtre sur les lignes, donc *avant* la virgule
dim(data)
summary(data)
iso_code continent location date
Length:209 Length:209 Length:209 Length:209
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
total_cases new_cases new_cases_smoothed total_deaths
Min. : 1 Min. : 0 Min. : 0.0 Min. : 1
1st Qu.: 24024 1st Qu.: 0 1st Qu.: 16.4 1st Qu.: 598
Median : 213982 Median : 58 Median : 149.3 Median : 3573
Mean : 5276354 Mean : 9348 Mean : 12071.7 Mean : 108749
3rd Qu.: 1020130 3rd Qu.: 1306 3rd Qu.: 1725.1 3rd Qu.: 18752
Max. :261504022 Max. :430360 Max. :560242.9 Max. :5199821
NA's :2 NA's :2 NA's :2 NA's :10
new_deaths new_deaths_smoothed total_cases_per_million
Min. : 0.00 Min. : -8.143 Min. : 8.6
1st Qu.: 0.00 1st Qu.: 0.143 1st Qu.: 4385.9
Median : 0.00 Median : 2.286 Median : 38622.5
Mean : 98.91 Mean : 142.700 Mean : 55584.1
3rd Qu.: 20.00 3rd Qu.: 24.929 3rd Qu.: 91788.6
Max. :4683.00 Max. :6934.429 Max. :249777.5
NA's :10 NA's :2 NA's :3
new_cases_per_million new_cases_smoothed_per_million total_deaths_per_million
Min. : 0.000 Min. : 0.000 Min. : 3.101
1st Qu.: 0.000 1st Qu.: 1.397 1st Qu.: 111.542
Median : 5.818 Median : 27.471 Median : 588.116
Mean : 122.415 Mean : 182.958 Mean : 946.956
3rd Qu.: 84.296 3rd Qu.: 185.441 3rd Qu.:1577.086
Max. :2061.251 Max. :2036.909 Max. :6028.523
NA's :3 NA's :3 NA's :11
new_deaths_per_million new_deaths_smoothed_per_million reproduction_rate
Min. : 0.0000 Min. :-0.2540 Min. : NA
1st Qu.: 0.0000 1st Qu.: 0.0095 1st Qu.: NA
Median : 0.0020 Median : 0.2840 Median : NA
Mean : 1.3115 Mean : 1.8690 Mean :NaN
3rd Qu.: 0.8818 3rd Qu.: 1.4585 3rd Qu.: NA
Max. :17.8140 Max. :17.5650 Max. : NA
NA's :11 NA's :3 NA's :209
icu_patients icu_patients_per_million hosp_patients
Min. : 20.0 Min. : 0.448 Min. : 172
1st Qu.: 145.5 1st Qu.:10.521 1st Qu.: 1290
Median : 217.0 Median :15.591 Median : 2606
Mean : 1836.3 Mean :18.046 Mean : 9121
3rd Qu.: 488.0 3rd Qu.:27.576 3rd Qu.: 4157
Max. :11350.0 Max. :34.093 Max. :43820
NA's :202 NA's :202 NA's :203
hosp_patients_per_million weekly_icu_admissions
Min. : 18.51 Min. :49
1st Qu.: 56.82 1st Qu.:49
Median :123.74 Median :49
Mean :177.57 Mean :49
3rd Qu.:140.30 3rd Qu.:49
Max. :619.12 Max. :49
NA's :203 NA's :208
weekly_icu_admissions_per_million weekly_hosp_admissions
Min. :5.274 Min. : 84
1st Qu.:5.274 1st Qu.: 207
Median :5.274 Median : 330
Mean :5.274 Mean :15227
3rd Qu.:5.274 3rd Qu.:22798
Max. :5.274 Max. :45266
NA's :208 NA's :206
weekly_hosp_admissions_per_million new_tests total_tests
Min. : 9.041 Min. : NA Min. : NA
1st Qu.: 23.453 1st Qu.: NA 1st Qu.: NA
Median : 37.864 Median : NA Median : NA
Mean : 60.958 Mean :NaN Mean :NaN
3rd Qu.: 86.916 3rd Qu.: NA 3rd Qu.: NA
Max. :135.969 Max. : NA Max. : NA
NA's :206 NA's :209 NA's :209
total_tests_per_thousand new_tests_per_thousand new_tests_smoothed
Min. : NA Min. : NA Min. : NA
1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
Median : NA Median : NA Median : NA
Mean :NaN Mean :NaN Mean :NaN
3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
Max. : NA Max. : NA Max. : NA
NA's :209 NA's :209 NA's :209
new_tests_smoothed_per_thousand positive_rate tests_per_case
Min. : NA Min. : NA Min. : NA
1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
Median : NA Median : NA Median : NA
Mean :NaN Mean :NaN Mean :NaN
3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
Max. : NA Max. : NA Max. : NA
NA's :209 NA's :209 NA's :209
tests_units total_vaccinations people_vaccinated
Length:209 Min. :1.006e+05 Min. :5.430e+04
Class :character 1st Qu.:5.125e+06 1st Qu.:3.389e+06
Mode :character Median :2.698e+07 Median :1.593e+07
Mean :4.863e+08 Mean :2.532e+08
3rd Qu.:1.248e+08 3rd Qu.:7.141e+07
Max. :7.943e+09 Max. :4.269e+09
NA's :147 NA's :150
people_fully_vaccinated total_boosters new_vaccinations
Min. :4.632e+04 Min. : 180081 Min. : 5
1st Qu.:2.358e+06 1st Qu.: 1117167 1st Qu.: 4953
Median :9.838e+06 Median : 2845397 Median : 31632
Mean :1.961e+08 Mean : 24568327 Mean : 2301120
3rd Qu.:5.631e+07 3rd Qu.: 14101808 3rd Qu.: 297380
Max. :3.382e+09 Max. :232082543 Max. :31551065
NA's :149 NA's :176 NA's :158
new_vaccinations_smoothed total_vaccinations_per_hundred
Min. : 62 Min. : 12.66
1st Qu.: 8121 1st Qu.: 86.50
Median : 72898 Median :125.30
Mean : 1818702 Mean :114.53
3rd Qu.: 557528 3rd Qu.:146.55
Max. :29928829 Max. :194.59
NA's :147 NA's :147
people_vaccinated_per_hundred people_fully_vaccinated_per_hundred
Min. : 8.21 Min. : 4.85
1st Qu.:45.85 1st Qu.:38.06
Median :63.53 Median :58.21
Mean :58.80 Mean :51.50
3rd Qu.:71.56 3rd Qu.:64.66
Max. :82.86 Max. :79.72
NA's :150 NA's :149
total_boosters_per_hundred new_vaccinations_smoothed_per_million
Min. : 0.030 Min. : 157
1st Qu.: 3.940 1st Qu.: 1281
Median : 7.280 Median : 2654
Mean : 9.736 Mean : 3123
3rd Qu.:10.930 3rd Qu.: 4613
Max. :43.850 Max. :12001
NA's :176 NA's :147
new_people_vaccinated_smoothed new_people_vaccinated_smoothed_per_hundred
Min. : 26 Min. :0.0050
1st Qu.: 3169 1st Qu.:0.0445
Median : 23272 Median :0.0690
Mean : 456827 Mean :0.0898
3rd Qu.: 138212 3rd Qu.:0.1225
Max. :7565461 Max. :0.3900
NA's :150 NA's :150
stringency_index population population_density median_age
Min. :16.67 Min. :8.120e+02 Min. : 1.98 Min. :15.10
1st Qu.:35.19 1st Qu.:2.140e+06 1st Qu.: 37.73 1st Qu.:22.00
Median :43.28 Median :9.870e+06 Median : 83.48 Median :29.40
Mean :45.80 Mean :1.527e+08 Mean : 441.97 Mean :30.23
3rd Qu.:63.43 3rd Qu.:3.584e+07 3rd Qu.: 208.35 3rd Qu.:38.70
Max. :72.22 Max. :7.875e+09 Max. :20546.77 Max. :48.20
NA's :195 NA's :1 NA's :16 NA's :22
aged_65_older aged_70_older gdp_per_capita extreme_poverty
Min. : 1.144 Min. : 0.526 Min. : 661.2 Min. : 0.10
1st Qu.: 3.466 1st Qu.: 2.048 1st Qu.: 3823.2 1st Qu.: 0.60
Median : 6.211 Median : 3.699 Median : 11840.9 Median : 2.50
Mean : 8.586 Mean : 5.420 Mean : 18800.1 Mean :13.93
3rd Qu.:13.914 3rd Qu.: 8.607 3rd Qu.: 26777.6 3rd Qu.:21.40
Max. :27.049 Max. :18.493 Max. :116935.6 Max. :77.60
NA's :24 NA's :23 NA's :20 NA's :84
cardiovasc_death_rate diabetes_prevalence female_smokers male_smokers
Min. : 79.37 Min. : 0.990 Min. : 0.10 Min. : 7.70
1st Qu.:173.49 1st Qu.: 5.310 1st Qu.: 1.90 1st Qu.:22.35
Median :244.66 Median : 7.155 Median : 6.30 Median :32.25
Mean :264.32 Mean : 8.098 Mean :10.57 Mean :32.88
3rd Qu.:334.39 3rd Qu.:10.232 3rd Qu.:19.02 3rd Qu.:41.35
Max. :724.42 Max. :30.530 Max. :44.00 Max. :78.10
NA's :22 NA's :17 NA's :63 NA's :65
handwashing_facilities hospital_beds_per_thousand life_expectancy
Min. : 1.188 Min. : 0.100 Min. :53.28
1st Qu.:20.105 1st Qu.: 1.300 1st Qu.:67.77
Median :49.542 Median : 2.398 Median :74.17
Mean :50.271 Mean : 3.001 Mean :72.92
3rd Qu.:81.569 3rd Qu.: 3.856 3rd Qu.:78.10
Max. :98.999 Max. :13.800 Max. :86.75
NA's :114 NA's :39 NA's :13
human_development_index excess_mortality_cumulative_absolute
Min. :0.3940 Min. : NA
1st Qu.:0.6020 1st Qu.: NA
Median :0.7400 Median : NA
Mean :0.7225 Mean :NaN
3rd Qu.:0.8290 3rd Qu.: NA
Max. :0.9570 Max. : NA
NA's :20 NA's :209
excess_mortality_cumulative excess_mortality
Min. : NA Min. : NA
1st Qu.: NA 1st Qu.: NA
Median : NA Median : NA
Mean :NaN Mean :NaN
3rd Qu.: NA 3rd Qu.: NA
Max. : NA Max. : NA
NA's :209 NA's :209
excess_mortality_cumulative_per_million
Min. : NA
1st Qu.: NA
Median : NA
Mean :NaN
3rd Qu.: NA
Max. : NA
NA's :209
# Le jeu de données comprend quelques "lignes résumé" (continents, catégories de revenus...).
# En observant un peu, on remarque que leur code ISO démarre par "OWID_" :
owid_lines <- startsWith(data$iso_code, "OWID_")
data[owid_lines,]
| iso_code | continent | location | date | total_cases | new_cases | new_cases_smoothed | total_deaths | new_deaths | new_deaths_smoothed | ⋯ | female_smokers | male_smokers | handwashing_facilities | hospital_beds_per_thousand | life_expectancy | human_development_index | excess_mortality_cumulative_absolute | excess_mortality_cumulative | excess_mortality | excess_mortality_cumulative_per_million | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | ⋯ | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 1299 | OWID_AFR | Africa | 2021-11-28 | 8637754 | 5439 | 5127.143 | 222573 | 72 | 139.000 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 7014 | OWID_ASI | Asia | 2021-11-28 | 81941193 | 77981 | 90983.143 | 1216288 | 1226 | 1468.714 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 40764 | OWID_EUR | Europe | 2021-11-28 | 73954046 | 297363 | 364354.571 | 1411846 | 2731 | 3893.000 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 41440 | OWID_EUN | European Union | 2021-11-28 | 46193217 | 213516 | 257871.714 | 841808 | 951 | 1801.143 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 53290 | OWID_HIC | High income | 2021-11-28 | 114125485 | 288779 | 390332.000 | 1823157 | 1060 | 2495.714 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 57839 | OWID_INT | International | 2021-11-28 | 721 | 0 | 0.000 | 15 | 0 | 0.000 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 65670 | OWID_KOS | Europe | Kosovo | 2021-11-28 | 161081 | 6 | 10.714 | 2983 | 0 | 0.143 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| 72584 | OWID_LIC | Low income | 2021-11-28 | 1377792 | 422 | 704.143 | 37259 | 20 | 16.429 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 73261 | OWID_LMC | Lower middle income | 2021-11-28 | 63654375 | 41241 | 52539.286 | 1152335 | 1237 | 1626.714 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 90945 | OWID_NAM | North America | 2021-11-28 | 57752998 | 35075 | 78981.429 | 1165076 | 395 | 1021.571 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 93220 | OWID_OCE | Oceania | 2021-11-28 | 308152 | 1462 | 1626.571 | 3282 | 6 | 16.000 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 113670 | OWID_SAM | South America | 2021-11-28 | 38908293 | 13040 | 19170.000 | 1180741 | 253 | 396.143 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 130135 | OWID_UMC | Upper middle income | 2021-11-28 | 82344757 | 99918 | 116667.429 | 2187055 | 2366 | 2795.571 | ⋯ | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | |
| 134625 | OWID_WRL | World | 2021-11-28 | 261504022 | 430360 | 560242.857 | 5199821 | 4683 | 6934.429 | ⋯ | 6.434 | 34.635 | 60.13 | 2.705 | 72.58 | 0.737 | NA | NA | NA | NA |
# Hormis le Kosovo, les autres lignes sont à supprimer car on retrouvera ces indicateurs (et plus encore)
# via l'analyse ACP. Renommons aussi OWID_KOS après avoir vérifier que ce code n'est pas utilisé :
any(data$iso_code == "KOS")
kosovo_index = which(data$iso_code == "OWID_KOS")
data[kosovo_index, "iso_code"] <- "KOS"
# Supprimons alors les lignes "aggrégées" :
owid_lines <- startsWith(data$iso_code, "OWID_")
data <- data[!owid_lines,]
dim(data)
# Pour supprimer les lignes avec valeurs manquantes, il y a plusieurs options.
# J'en indique ici 3, de la plus complexe à la plus simple :
data1 <- data[apply(data, 1, function(row) all(!is.na(row))),] #méthode 1
data2 <- data[complete.cases(data),] #méthode 2
data3 <- na.omit(data) #méthode 3
# On remarque alors qu'il n'y a plus aucune ligne ==> il faut restreindre les colonnes
# (ou "deviner" les valeurs manquantes d'une manière ou d'une autre : voir le package missMDA)
# (ici on se contente de la version simple : pas de données manquantes en entrée).
nrow(data3)
# Que représentent les variables ?
colnames(data)
# Variables "new_*" : instantané journalier d'un certain indicateur. On ne s'y intéressera pas ici (cf. plus bas).
# De même pour les variables "weekly_*" (indicateurs hebdomadaires, j'imagine). Reste :
selection <- colnames(data)[!startsWith(colnames(data), "new_") & !startsWith(colnames(data), "weekly_")]
# Colonnes avec +50% de valeurs renseignées
selection <- selection[ apply(data[,selection], 2, function(col) sum(!is.na(col)) > nrow(data)/2) ]
selection
On y voit enfin plus clair :
data$tests_units #???
selection <- selection[selection != "tests_units"]
newData <- na.omit(data[,selection])
nrow(newData)
92 lignes est raisonnable (proche de 50% de la taille du jeu de données initial). Cependant, pour être cohérent il faut en plus choisir un type de variable : absolu, ou relatif ? Je préfère les indicateurs relatifs (*_per_million, *_density) :
selection <- selection[!selection %in% c("total_cases", "total_deaths", "population")]
newData <- na.omit(data[,selection])
nrow(newData) #92 encore
rownames(newData) <- newData$iso_code #pour l'affichage des individus
newData <- newData[,-c(1,4)] #suppression des colonnes "code ISO" et "date", désormais inutiles
Note : toute l'analyse jusqu'ici aurait pu se faire aussi facilement avec un autre langage, Python par exemple.
À partir d'ici cependant, le package R FactoMineR est très pratique (pas d'équivalent Python (?!))
# ...On est enfin prêt pour l'ACP !
library(FactoMineR)
res.pca <- PCA(newData, quali.sup=1:2, ncp=6, graph=FALSE)
options(repr.plot.width=15, repr.plot.height=10)
plotInd <- plot(res.pca, choix="ind", invisible="quali")
plotVar <- plot(res.pca, choix="var")
library(gridExtra)
grid.arrange(plotInd, plotVar, ncol=2)
Warning message: “ggrepel: 2 unlabeled data points (too many overlaps). Consider increasing max.overlaps”
Bon, il semble qu'aucun individu ne se détache. En fait certaines lignes ont presque toutes leurs valeurs renseignées, et une fois complétées à l'aide de Google on trouve des individus extrêmes (Monaco, Singapour).
which(data$location %in% c("Monaco", "Singapore"))
data[158,selection]
| iso_code | continent | location | date | total_cases_per_million | total_deaths_per_million | population_density | median_age | aged_65_older | aged_70_older | gdp_per_capita | extreme_poverty | cardiovasc_death_rate | diabetes_prevalence | female_smokers | male_smokers | hospital_beds_per_thousand | life_expectancy | human_development_index | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| 109831 | SGP | Asia | Singapore | 2021-11-28 | 48111.89 | 128.539 | 7915.731 | 42.4 | 12.922 | 7.049 | 85535.38 | NA | 92.243 | 10.99 | 5.2 | 28.3 | 2.4 | 83.62 | 0.938 |
Côté variables, aged_65_older et aged_70_older apparaissent très corrélées (en fait même confondues). C'est logique, on gardera donc seulement aged_70_older après vérification numérique :
cor(newData[,c("aged_65_older", "aged_70_older")])
| aged_65_older | aged_70_older | |
|---|---|---|
| aged_65_older | 1.000000 | 0.994226 |
| aged_70_older | 0.994226 | 1.000000 |
newData <- subset(newData, select=-aged_65_older)
res.pca <- PCA(newData, quali.sup=1:2, ncp=6, graph=FALSE)
plotInd <- plot(res.pca, choix="ind", invisible="quali", habillage=1)
plotVar <- plot(res.pca, choix="var")
library(gridExtra)
grid.arrange(plotInd, plotVar, ncol=2)
Environ 65% de l'inertie expliquée dans ce premier plan ACP (49.3 + 14.2). Le cercle des corrélations oppose logiquement "richesse" à droite (HDI, PIB/hab), avec "pauvreté" à gauche. Sur le nuage des individus cela correspond grossièrement à l'opposition Europe occidentale / Afrique (quelques exceptions : Tunisie, Seychelles, ...).
Il est intéressant de constater que les indicateurs de richesse sont très corrélés positivement au nombre de morts par millions, lui-même anti-corrélé avec extreme_poverty : le COVID frapperait plutôt les riches ? Mais pourquoi donc, puisque le virus est partout ? Et bien un élément de réponse se trouve dans ce même plan ACP : aged_70_older => on y vit plus vieux, et dans une moindre mesure diabetes_prevalence => plus de cas de diabète (à vérifier numériquement).
cor(newData[,c("aged_70_older", "total_deaths_per_million", "human_development_index",
"diabetes_prevalence", "extreme_poverty")])
| aged_70_older | total_deaths_per_million | human_development_index | diabetes_prevalence | extreme_poverty | |
|---|---|---|---|---|---|
| aged_70_older | 1.000000000 | 0.6264610 | 0.8196467 | -0.008746564 | -0.5752816 |
| total_deaths_per_million | 0.626461040 | 1.0000000 | 0.5441791 | 0.117930733 | -0.4798006 |
| human_development_index | 0.819646681 | 0.5441791 | 1.0000000 | 0.222428371 | -0.7786659 |
| diabetes_prevalence | -0.008746564 | 0.1179307 | 0.2224284 | 1.000000000 | -0.4010230 |
| extreme_poverty | -0.575281608 | -0.4798006 | -0.7786659 | -0.401022956 | 1.0000000 |
La corrélation (resp. anti-corrélation) aged_70_older avec HDI (resp. extreme_poverty) et total_deaths est vérifiée. De même, on observe une légère corrélation positive (resp. négative) entre diabetes_prevalence et HDI (resp. extreme_poverty).
Ensuite, le taux de fumeuses semble très corrélé à l'âge médian : les femmes auraient plus tendance à fumer dans les pays où l'on vit plus vieux (donc en général plus riches). Ce n'est pas le cas de male_smokers : le taux de fumeurs n'indique quant à lui pas grand chose. De même, et plus étonnament, le taux de mortalité par maladies cardiovasculaires (infarctus j'imagine) ne paraît corrélé à rien - si ce n'est justement et assez logiquement, la proportion de fumeurs : "According to the American Heart Association, cardiovascular disease accounts for about 800,000 U.S. deaths every year,5 making it the leading cause of all deaths in the United States. Of those, nearly 20 percent are due to cigarette smoking." [https://www.fda.gov/tobacco-products/health-effects-tobacco-use/how-smoking-affects-heart-health#]
La coloration par continents montre une opposition haut/bas entre Europe de l'est et Europe de l'ouest + USA/Canada/Israel/Corée/Australie. Il semble y avoir relativement plus de fumeurs en Géorgie/Ukraine/Russie. les pays d'Amérique centrale et du sud sont plus bas, donc a priori moins touchés par les décès par infarctus et comportant moins de fumeurs. Il n'y a pas assez de pays d'Océanie pour en dire grand chose, et l'Asie est répartie un peu partout, montrant une grande inhomogénéité en comparaison aux autres continents.
Vérifions notre analyse en regardant de plus près quelques individus :
indivs_indices <- rownames(newData) %in% c("LUX", "UKR", "NER", "ECU")
newData[indivs_indices,c("location", "total_deaths_per_million", "aged_70_older", "male_smokers",
"cardiovasc_death_rate", "human_development_index")]
| location | total_deaths_per_million | aged_70_older | male_smokers | cardiovasc_death_rate | human_development_index | |
|---|---|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| ECU | Ecuador | 1851.919 | 4.458 | 12.3 | 140.448 | 0.759 |
| LUX | Luxembourg | 1364.179 | 9.842 | 26.0 | 128.275 | 0.916 |
| NER | Niger | 10.107 | 1.378 | 15.4 | 238.339 | 0.394 |
| UKR | Ukraine | 2078.482 | 11.133 | 47.4 | 539.849 | 0.779 |
Luxembourg : population âgée, HDI élevé, 2x moins de fumeurs qu'en Ukraine mais 2x + qu'en Equateur.
Niger : population jeune, HDI bas, peu de fumeurs, très peu de morts du COVID.
Bref, passons au second plan ACP :
plotInd <- plot(res.pca, choix="ind", invisible="quali", habillage=1, axes=3:4)
plotVar <- plot(res.pca, choix="var", axes=3:4)
library(gridExtra)
grid.arrange(plotInd, plotVar, ncol=2)
Warning message: “ggrepel: 5 unlabeled data points (too many overlaps). Consider increasing max.overlaps”
Peu d'inertie expliquée dans ce plan (à peine 18%), mais une observation intéressante : anti-corrélation population_density et total_deaths_per_million ? À vérifier numériquement bien sûr car cette dernière flèche est loin du bord. Ce serait cependant cohérent : densément peuplé => contaminations plus faciles => plus de cas => plus de personnes très fragiles touchées => plus de morts.
On note aussi l'anti-corrélation entre diabetes_prevalence et extreme_poverty, déjà un peu observée dans le premier plan. Vérification numérique :
indivs_indices <- rownames(newData) %in% c("MLT", "EGY", "MNE", "MWI")
newData[indivs_indices,c("location", "total_deaths_per_million", "population_density",
"diabetes_prevalence", "extreme_poverty")]
| location | total_deaths_per_million | population_density | diabetes_prevalence | extreme_poverty | |
|---|---|---|---|---|---|
| <chr> | <dbl> | <dbl> | <dbl> | <dbl> | |
| EGY | Egypt | 195.159 | 97.999 | 17.31 | 1.3 |
| MWI | Malawi | 117.317 | 197.519 | 3.94 | 71.4 |
| MLT | Malta | 906.801 | 1454.037 | 8.83 | 0.2 |
| MNE | Montenegro | 3638.240 | 46.280 | 10.08 | 1.0 |
Opposition Égypte / malawi vérifiée sur l'axe diabète/pauvreté, ainsi que malte/Montenegro sur l'axe morts_par_million/densité.
À l'URL indiquée, on lit clairement ce que représente chaque variable :
Note : gill = branchie, stalk = tige, veil = voile (merci Google Translate).
# Récupération du jeu de données
data <- read.csv("https://archive.ics.uci.edu/ml/machine-learning-databases/mushroom/agaricus-lepiota.data",
header=FALSE) #attention: pas de header (on s'en rend vite compte)
dim(data)
head(data)
| V1 | V2 | V3 | V4 | V5 | V6 | V7 | V8 | V9 | V10 | ⋯ | V14 | V15 | V16 | V17 | V18 | V19 | V20 | V21 | V22 | V23 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | ⋯ | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | |
| 1 | p | x | s | n | t | p | f | c | n | k | ⋯ | s | w | w | p | w | o | p | k | s | u |
| 2 | e | x | s | y | t | a | f | c | b | k | ⋯ | s | w | w | p | w | o | p | n | n | g |
| 3 | e | b | s | w | t | l | f | c | b | n | ⋯ | s | w | w | p | w | o | p | n | n | m |
| 4 | p | x | y | w | t | p | f | c | n | n | ⋯ | s | w | w | p | w | o | p | k | s | u |
| 5 | e | x | s | g | f | n | f | w | b | k | ⋯ | s | w | w | p | w | o | e | n | a | g |
| 6 | e | x | y | y | t | a | f | c | b | n | ⋯ | s | w | w | p | w | o | p | k | n | g |
23 variables : la première indique si le champignon est comestible (edible) ou non (poisonous).
On peut se poser diverses questions passionnantes sur les correspondances entre modalités, par exemple :
# Commençons par donner des noms aux colonnes :
colnames(data) <- c("edible", "cap-shape", "cap-surface", "cap-color", "bruises", "odor", "gill-attachment",
"gill-spacing", "gill-size", "gill-color", "stalk-shape", "stalk-root",
"stalk-surface-above-ring", "stalk-surface-below-ring", "stalk-color-above-ring",
"stalk-color-below-ring", "veil-type", "veil-color", "ring-number", "ring-type",
"spore-print-color", "population", "habitat")
# Essayons odeur / habitat :
eff_obs <- table(subset(data, select=c("odor", "habitat")))
eff_obs
habitat odor d g l m p u w a 48 176 0 128 48 0 0 c 192 0 0 0 0 0 0 f 624 576 192 0 624 144 0 l 48 176 0 128 48 0 0 m 36 0 0 0 0 0 0 n 1816 1092 256 36 40 96 192 p 0 128 0 0 0 128 0 s 192 0 192 0 192 0 0 y 192 0 192 0 192 0 0
Certaines lignes sont identiques : odeurs a/l et s/y + m/c (proportionnelles : c'est pareil). Elles peuvent d'ailleurs être regroupées sans changer le résultat de l'analyse : faisons ça.
eff_obs[1,] <- 2 * eff_obs[1,]
eff_obs[2,] <- eff_obs[2,] + eff_obs[5,]
eff_obs[8,] <- 2 * eff_obs[8,]
rownames(eff_obs)[c(1,2,8)] <- c("a-l", "c-m", "s-y")
eff_obs <- eff_obs[c(1:3, 6:8),]
eff_obs
habitat odor d g l m p u w a-l 96 352 0 256 96 0 0 c-m 228 0 0 0 0 0 0 f 624 576 192 0 624 144 0 n 1816 1092 256 36 40 96 192 p 0 128 0 0 0 128 0 s-y 384 0 384 0 384 0 0
# Effectifs théoriques, pour comparaison :
eff_theo <- as.matrix(rowSums(eff_obs)) %*% colSums(eff_obs) / sum(eff_obs)
colnames(eff_theo) <- colnames(eff_obs)
round(eff_theo)
| d | g | l | m | p | u | w | |
|---|---|---|---|---|---|---|---|
| a-l | 310 | 212 | 82 | 29 | 113 | 36 | 19 |
| c-m | 88 | 60 | 23 | 8 | 32 | 10 | 5 |
| f | 837 | 571 | 221 | 78 | 304 | 98 | 51 |
| n | 1367 | 933 | 361 | 127 | 497 | 160 | 83 |
| p | 99 | 68 | 26 | 9 | 36 | 12 | 6 |
| s-y | 446 | 305 | 118 | 41 | 162 | 52 | 27 |
res.ca <- CA(eff_obs, graph=FALSE)
plot(res.ca)
a-l (odeur) et m (habitat) sont du même côté : les prés contiennent surtout des champignons qui sentent l'amande ou l'anis (le tableau numérique - effectifs observés - indique qu'il faut lire l'association ainsi ; on peut aussi s'aider des effectifs théoriques : "m" est significativement plus associé avec a-l qu'en cas d'indépendance).
f, s-y sont du même côté que l et p : les champignons d'odeur "foul", "spicy" et "fishy" (épicé, nauséabond, poisson : bref des champis qu'on n'a pas très envie de manger a priori) se trouvent significativement plus dans les feuilles et sur les chemins, en comparaison aux effectifs théoriques. Pourquoi, ne me demandez pas.
p (odeur) / u (habitat) semble ressortir : odeur âcre (merci G...) en milieu urbain. On n'aurait pas très envie de manger des champignons poussant au milieu du béton, donc disons que c'est conforme à l'intuition.
Enfin, reste n/w (odeur/habitat) et c-m/d : champignons inodore dans les déchets (heu... ? admettons ^^), et champis sentant la "créosote" (charbon ? ...apparemment c'est un mot français) et le moisi dans les bois : assez cohérent disons.
J'ai un peu la flemme d'écrire autant de détails pour les deux autres cas, alors j'irai à l'essentiel :
eff_obs <- table(subset(data, select=c("cap-color", "edible")))
eff_obs
eff_theo <- as.matrix(rowSums(eff_obs)) %*% colSums(eff_obs) / sum(eff_obs)
colnames(eff_theo) <- colnames(eff_obs)
round(eff_theo)
res.ca <- CA(eff_obs, graph=FALSE)
#plot(res.ca) #erreur: seulement deux modalités => unidimensionnel
# => "plot" manuel.. (?!)
xrange <- range(res.ca$row$coord, res.ca$col$coord)
plot(NULL, xlim=xrange, ylim=c(-1,1))
points(res.ca$row$coord, rep(0, length(res.ca$row$coord)), col="blue", pch=19, cex=2)
points(res.ca$col$coord, rep(0, length(res.ca$col$coord)), col="red", pch=19, cex=2)
text(res.ca$row$coord, rep(-0.2, length(res.ca$row$coord)), rownames(eff_obs), col="blue", cex=2)
text(res.ca$col$coord, rep(0.2, length(res.ca$col$coord)), colnames(eff_obs), col="red", cex=2)
edible
cap-color e p
b 48 120
c 32 12
e 624 876
g 1032 808
n 1264 1020
p 56 88
r 16 0
u 16 0
w 720 320
y 400 672
| e | p | |
|---|---|---|
| b | 87 | 81 |
| c | 23 | 21 |
| e | 777 | 723 |
| g | 953 | 887 |
| n | 1183 | 1101 |
| p | 75 | 69 |
| r | 8 | 8 |
| u | 8 | 8 |
| w | 539 | 501 |
| y | 555 | 517 |
b, y, p, e : plutôt poison - jaune/rouge/rose/"chamois"
n, g : on sait pas trop, faire attention - marron/gris
w, c : plutôt comestibles (3 chances sur 4 disons) - blanc/"cannelle" (heu... cannelle != marron ?! c'est subtil !)
r, u : tout est bon ! - vert/pourpre
Ok pour jaune/rouge/rose = méfiance, mais à part ça, je n'avais aucune intuition sur le sujet et cette étude ne m'inspirera pas spécialement confiance en présence d'un champignon violet =)
Note : deux lignes identiques, r = u.
Finalement, couleur des branchies vs. couleur de la tige "sous anneau" (c'est précis !) ?
eff_obs <- table(subset(data, select=c("gill-color", "stalk-color-below-ring")))
eff_obs
eff_theo <- as.matrix(rowSums(eff_obs)) %*% colSums(eff_obs) / sum(eff_obs)
colnames(eff_theo) <- colnames(eff_obs)
round(eff_theo)
res.ca <- CA(eff_obs, graph=FALSE)
plot(res.ca)
stalk-color-below-ring
gill-color b c e g n o p w y
b 0 0 0 0 0 0 864 864 0
e 0 0 48 0 0 0 0 48 0
g 144 0 0 0 144 0 144 320 0
h 144 0 0 0 144 0 144 300 0
k 0 0 0 0 0 0 0 408 0
n 0 0 0 144 0 64 144 696 0
o 0 0 0 0 0 64 0 0 0
p 144 0 0 144 144 0 288 772 0
r 0 0 0 0 0 0 0 24 0
u 0 0 0 144 0 0 144 204 0
w 0 18 48 144 80 0 144 748 20
y 0 18 0 0 0 64 0 0 4
| b | c | e | g | n | o | p | w | y | |
|---|---|---|---|---|---|---|---|---|---|
| b | 92 | 8 | 20 | 123 | 109 | 41 | 398 | 932 | 5 |
| e | 5 | 0 | 1 | 7 | 6 | 2 | 22 | 52 | 0 |
| g | 40 | 3 | 9 | 53 | 47 | 18 | 173 | 406 | 2 |
| h | 39 | 3 | 9 | 52 | 46 | 17 | 169 | 395 | 2 |
| k | 22 | 2 | 5 | 29 | 26 | 10 | 94 | 220 | 1 |
| n | 56 | 5 | 12 | 74 | 66 | 25 | 241 | 566 | 3 |
| o | 3 | 0 | 1 | 5 | 4 | 2 | 15 | 35 | 0 |
| p | 79 | 7 | 18 | 106 | 94 | 35 | 344 | 805 | 4 |
| r | 1 | 0 | 0 | 2 | 2 | 1 | 6 | 13 | 0 |
| u | 26 | 2 | 6 | 35 | 31 | 12 | 113 | 266 | 1 |
| w | 64 | 5 | 14 | 85 | 76 | 28 | 277 | 649 | 4 |
| y | 5 | 0 | 1 | 6 | 5 | 2 | 20 | 46 | 0 |
Warning message: “ggrepel: 8 unlabeled data points (too many overlaps). Consider increasing max.overlaps”
e/e : rouge/rouge. Ok.
y/o : jaune (branchies) / orange (tige). Ok, c'est assorti.
o/c : orange (branchies) / cannelle (tige). De même, ce n'est pas choquant.
Le reste est proche de l'indépendance d'après le graphe.