Ce document a été généré directement depuis RStudio en utilisant l’outil Markdown.

Un certain nombre de packages ont été utilisés pendant cette session : il faudra les installer soit depuis le site du CRAN quand le package est disponible:

install.packages(c(
  "devtools", # permet d'utiliser la fonction install_github()
  "essurvey", # permet d'importer directement les données ESS
  "foreign",   # permet d'importer des formats de données Stata, SPSS
  "ggridges", # permet de représenter les densités non paramétriques 
  "gridExtra", # permet d'afficher plusieurs graphiques ggplot2 ensemble
  "labelled", # recherche de mots-clés
  "tidyverse"  # ensemble de packages pour le data scientist
  )
)

Soit on les télécharge depuis Github en utilisant la fonction install_github() inclus dans le package devtools lorsque le package n’est pas encore publié :

devtools::install_github('rensa/ggflags')

Voici quelques références :

1 Importation de données

En général, il existe autant de fonctions d’importation qu’il existe de formats de fichiers. Par exemple :

Selon le type de formats, la fonction n’est pas nécessairement disponible par défaut, il faudra charger un package

1.1 Solution 1

  1. Télécharger la base depuis : https://www.europeansocialsurvey.org/download.html?file=ESS8e02_1&y=2016

  2. Enregistrer la base dans un répertoire de travail nommé progedo

  3. De-zipper le fichier : cela a pour effet de créer un sous-repertoire ESS8e02_1.stata

  4. Copier le fichier session_R.R dans le répertoire progedo

  5. Dans le répertoire progedo, créer un sous-répertoire figures

  6. Ouvrir le fichier session_R.R depuis RStudio.

La commande getwd() permet de vérifier que le répertoire courant est progedo

getwd()
## [1] "/media/thibault/My Passport/confinement/progedo"

Pour importer un jeu de données qui a été obtenu depuis Stata, il faut utiliser la fonction read.dta() du package foreign. Si le répertoire courant est progedo, on fait :

my_data_8 <- foreign::read.dta("ESS8e02_1.stata/ESS8e02_1.dta")

Sinon, on peut écrire le chemin d’accès en entier :

my_data_8 <- foreign::read.dta("/media/thibault/My Passport/confinement/progedo/ESS8e02_1.stata/ESS8e02_1.dta")

1.2 Solution 2

La base de données peut être téléchargée directement avec la fonction import_rounds() du package essurvey

my_data_8_tb <- essurvey::import_rounds(8, 
                     ess_email = "thibault.laurent@tse-fr.eu",
                     format = 'stata')

2 Apprentissage de la base

Dans R, un jeu de données est stocké sous forme d’un data.frame ou tibble:

class(my_data_8)
## [1] "data.frame"
class(my_data_8_tb)
## [1] "tbl_df"     "tbl"        "data.frame"

Pour connaître un certain nombre d’informations sur le jeu de données, on peut utiliser la commande str() :

str(my_data_8)

Cette commande retourne:

Le jeu de données est volumineux et on peut vouloir interroger la base sur le nom des variables. Par exemple, pour trouver toutes les variables qui font appel à “poor” dans la table ou dans les meta-données, on utilise la commande look_for() du package labelled

library(labelled) # install from CRAN first
labelled::look_for(my_data_8, 'poor')
## pos   variable label                                col_type values             
## <chr> <chr>    <chr>                                <chr>    <chr>              
## 100   impcntr  Allow many/few immigrants from poor… fct      Allow many to come…
## ​      ​         ​                                     ​         Allow some         
## ​      ​         ​                                     ​         Allow a few        
## ​      ​         ​                                     ​         Allow none         
## ​      ​         ​                                     ​         Refusal            
## ​      ​         ​                                     ​         Don't know         
## ​      ​         ​                                     ​         No answer

Un data.frame peut être vu comme une liste contenant \(p\) vecteurs dont la taille est la même.

Pour visualiser la table en entier, on peut faire :

View(my_data_8)

Pour afficher les premières lignes, on utilise la fonction head():

head(my_data_8)

Pour extraire des éléments d’un data.frame, il existe plusieurs syntaxes possibles. Ici, les commandes suivantes permettent d’afficher les 5 premières observations de la variable cntry :

my_data_8[1:5, "cntry"]
## [1] "AT" "AT" "AT" "AT" "AT"
my_data_8[1:5, 6]
## [1] "AT" "AT" "AT" "AT" "AT"
my_data_8$cntry[1:5]
## [1] "AT" "AT" "AT" "AT" "AT"

Les types principaux dans R sont :

my_data_8[1:5, "idno"]
## [1]  1  2  4  6 10
my_data_8[1:5, "pweight"]
## [1] 0.3703929 0.3703929 0.3703929 0.3703929 0.3703929
my_data_8[1:5, "cntry"]
## [1] "AT" "AT" "AT" "AT" "AT"
my_data_8[1:5, "trstplc"]
## [1] No trust at all 3               9               2              
## [5] Complete trust 
## 14 Levels: No trust at all 1 2 3 4 5 6 7 8 9 Complete trust ... No answer

Pour appliquer une fonction (ici class()) sur chaque élément (ou variable) du data.frame, on utilise la fonction sapply() de la manière suivante :

class_var <- sapply(my_data_8, class)
table(class_var)
## class_var
## character    factor   integer   numeric 
##        11       479        41         3

On constate que la plupart des variables sont de type character ou factor

3 Statistique univariée

3.1 Variable qualitative (ou catégorielle)

Dans ce jeu de données, la majorité des variables sont catégorielles. Sur ce type de variable, on s’intéresse à la table des effecifs (ou table de contingence).

Pour calculer cette table, on a deux solutions :

  • on utilise la fonction table()
cntry_tab <- table(my_data_8$cntry)

Sur une table de contingence, on peut utiliser la fonction prop.table() pour obtenir les pourcentages

prop.table(cntry_tab)
## 
##         AT         BE         CH         CZ         DE         EE         ES 
## 0.04528353 0.03978642 0.03435691 0.05111857 0.06425305 0.04548629 0.04411201 
##         FI         FR         GB         HU         IE         IL         IS 
## 0.04336855 0.04663528 0.04413454 0.03636200 0.06211278 0.05760696 0.01982562 
##         IT         LT         NL         NO         PL         PT         RU 
## 0.05916147 0.04780679 0.03787145 0.03480749 0.03816433 0.02861198 0.05474576 
##         SE         SI 
## 0.03494266 0.02944556

Pour afficher dans un document la table, on utilise la fonction kable() du package knitr :

# d'abord, on crée un objet qui contient les effectifs et les proportions
my_tab <- rbind(round(cntry_tab, 0),
        round(prop.table(cntry_tab) * 100, 2))
row.names(my_tab) <- c("effectifs", "proportions")
knitr::kable(t(my_tab))
effectifs proportions
AT 2010 4.53
BE 1766 3.98
CH 1525 3.44
CZ 2269 5.11
DE 2852 6.43
EE 2019 4.55
ES 1958 4.41
FI 1925 4.34
FR 2070 4.66
GB 1959 4.41
HU 1614 3.64
IE 2757 6.21
IL 2557 5.76
IS 880 1.98
IT 2626 5.92
LT 2122 4.78
NL 1681 3.79
NO 1545 3.48
PL 1694 3.82
PT 1270 2.86
RU 2430 5.47
SE 1551 3.49
SI 1307 2.94 
library(tidyverse)
my_data_8 %>%
  count(cntry) %>%
  mutate(pct = round(100 * n / sum(n), 2)) %>%
  arrange(-n)
##    cntry    n  pct
## 1     DE 2852 6.43
## 2     IE 2757 6.21
## 3     IT 2626 5.92
## 4     IL 2557 5.76
## 5     RU 2430 5.47
## 6     CZ 2269 5.11
## 7     LT 2122 4.78
## 8     FR 2070 4.66
## 9     EE 2019 4.55
## 10    AT 2010 4.53
## 11    GB 1959 4.41
## 12    ES 1958 4.41
## 13    FI 1925 4.34
## 14    BE 1766 3.98
## 15    PL 1694 3.82
## 16    NL 1681 3.79
## 17    HU 1614 3.64
## 18    SE 1551 3.49
## 19    NO 1545 3.48
## 20    CH 1525 3.44
## 21    SI 1307 2.94
## 22    PT 1270 2.86
## 23    IS  880 1.98

Pour représenter un variable catégorielle, on utilise le diagramme en barres. On utilise ici la syntaxe à la ggplot qui permet d’obtenir plus simplement des graphiques prêts à être diffusés. La syntaxe est la suivante :

p <- ggplot(my_data_8, aes(x = cntry)) +
  geom_bar() 
p

A partir des commandes précédentes, on peut customiser le graphique en modifiant par exemple les légendes :

p + labs(title = "Répartition des individus par pays",
         caption = "Source: ESS 2016",
         x = "Pays", y = "Effectifs")

On peut re-ordonner les barres du plus grand effectif au plus petit :

my_data_8 %>%
  count(cntry) %>%
  mutate(cntry = fct_reorder(cntry, n, .desc = TRUE)) %>%
  ggplot(aes(x = cntry, y = n)) +
  geom_bar(stat = 'identity') + 
  labs(title = "Répartition des individus par pays",
                    caption = "Source: ESS 2016",
                    x = "Pays", y = "Effectifs")

Enfin, on peut même ajouter des drapeaux pour faciliter la lecture du graphique :

# on utilise un autre codage pour les pays
my_data_8$country <- tolower(my_data_8$cntry) 
# On crée une base de données qui va nous permettre d'afficher aux points (x, y) le drapeau : 
sum_country <- my_data_8 %>%
  group_by(country) %>%
  summarise(y = n()) %>%
  mutate(x = 1:length(unique(my_data_8$cntry)))
# on affiche les drapeaux
library(ggflags)
my_data_8 %>%
  count(country) %>%
  mutate(country = fct_reorder(country, n, .desc = TRUE)) %>%
  ggplot(aes(x = country, y = n, country = as.character(country))) +
  geom_bar(stat = 'identity') + 
  labs(title = "Répartition des individus par pays",
       caption = "Source: ESS 2016",
       x = "Pays", y = "Effectifs") + 
  geom_flag(aes(x = country, y = n))

Pour sauvegarder un graphique ggplot, il suffit d’utiliser la fonction ggsave()

ggsave("Figures/my_plot.pdf")
## Saving 7 x 5 in image

3.2 Variable quantitative

La fonction summary() renvoit des statistiques de base : minimum, maximum, 1er, 2ème, 3ème quartile, moyenne :

summary(my_data_8$agea)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   15.00   34.00   49.00   49.14   64.00  100.00     155

Pour représenter une variable quantitative :

  • on utilise le diagramme en bâton (si la variable est discrète, i.e. le nombre de valeurs possibles de la variable est fini)
plot(table(my_data_8$agea), xlab = "Age", ylab = "Effectifs")

  • l’histogramme (si la variable est continue, i.e. les valeurs possibles de la variable sont dans l’espace des réels)
hist(my_data_8$agea, xlab = "Age", ylab = "Effectifs", 
     main = "Histogramme avec classe d amplitudes égales")

3.3 Transformer une variable quantitative en variable qualitative et inversement

On va transformer la variable âge (quantitative) en 4 classes d’âge. Pour cela, on procède de la façon suivante :

my_data_8 <- my_data_8 %>%
  mutate(agea_quali = ifelse(agea <= 20, "<=20", 
                        ifelse(agea > 20 & agea <= 40, "]20-40]", 
                          ifelse(agea > 40 & agea <= 60, "]40-60]", "> 60"))))

Beacoup de variables qualitatives sont en réalité des notes comprises entre 0 et 10. On peut donc les convertir en variable quantitative discrète. Ici, on transforme la variable “confiance en la police” en integer (on fait la même chose pour les variables “confiance à la justice”, “confiance aux partis politiques” et “confiance aux hommes/femmes politiques”) :

my_data_8 <- my_data_8 %>%
  mutate(trstplc_quanti = as.integer(as.character(ifelse(
           trstplc == "No trust at all", 0, 
             ifelse(trstplc == "Complete trust", 10, trstplc)))),
         trstprt_quanti = as.integer(as.character(ifelse(
           trstprt == "No trust at all", 0, 
             ifelse(trstprt == "Complete trust", 10, trstprt)))),
         trstplt_quanti = as.integer(as.character(ifelse(
           trstplt == "No trust at all", 0, 
             ifelse(trstplt == "Complete trust", 10, trstplt)))),
         trstlgl_quanti = as.integer(as.character(ifelse(
           trstlgl == "No trust at all", 0, 
             ifelse(trstlgl == "Complete trust", 10, trstlgl)))) 
)

Dans la suite de l’étude, on va utiliser les variables suivantes :

  • avis sur l’immigration : impcntr
  • confiance aux partis politiques : trstprt et trstplc_quanti
  • confiance aux hommes femmes politiques : trstplt et trstplt_quanti
  • confiance dans la police : trstplc et trstplc_quanti
  • confiance dans la justice : trstlgl et trstlgl_quanti
  • pays
  • age

On va créer un sous jeu de données qui ne contient que ces variables :

extract <- my_data_8 %>%
  select(impcntr, trstplc, trstprt, trstplt, trstlgl,
         trstplc_quanti, trstprt_quanti, trstplt_quanti, trstlgl_quanti,
         country, agea, agea_quali, happy)

4 Gestion des valeurs manquantes

Pour connaître le nombre de valeurs manquantes, on applique la fonction sapply() vue précédemment à la fonction is.na() qui indique quelles sont les observations d’un vecteur qui ont des valeurs manquantes :

sapply(extract, function(x) length(which(is.na(x))))
##        impcntr        trstplc        trstprt        trstplt        trstlgl 
##           1450            320            855            646            849 
## trstplc_quanti trstprt_quanti trstplt_quanti trstlgl_quanti        country 
##            320            855            646            849              0 
##           agea     agea_quali          happy 
##            155            155            215

Pour traiter ces valeurs manquantes, il existe plusieurs alternatives :

  1. Supprimer ces valeurs en utilisant la commande na.omit()
extract_1 <- na.omit(extract)
  1. Remplacer par la moyenne pour les variables quantitatives et le mode pour les variables qualitatives
replace_na <- function(x) {
  if(class(x) %in% c("numeric", "integer")) {
    x[is.na(x)] <- mean(x, na.rm = T)
  } else {
    x[is.na(x)] <- names(which.max(table(x)))
  }
  x
}
extract_2 <- as.data.frame(lapply(extract, replace_na))
  1. Utiliser des méthodes plus sophistiquées : imputation par régression linéaire, méthode du plus proche voisin, forêts aléatoires (voir https://cran.r-project.org/web/views/MissingData.html)

5 Etude bivariée ou croisement de variables

5.1 Croisement de deux variables qualitatives

Quand on croise deux variables qualitatives, on peut représenter un tableau croisé, toujours avec la fonction table() :

tab_2 <- table(extract_2$country, extract_2$impcntr)

On peut ensuite calculer les proportions par ligne:

knitr::kable(round(prop.table(tab_2, margin = 1) * 100, 2))
Allow many to come and live here Allow some Allow a few Allow none Refusal Don’t know No answer
at 7.81 29.00 40.80 22.39 0 0 0
be 17.44 48.24 25.14 9.17 0 0 0
ch 14.56 49.31 29.38 6.75 0 0 0
cz 1.90 22.21 44.51 31.38 0 0 0
de 20.51 46.63 25.91 6.94 0 0 0
ee 4.26 31.10 37.94 26.70 0 0 0
es 27.53 39.58 26.30 6.59 0 0 0
fi 10.75 35.84 46.03 7.38 0 0 0
fr 14.11 48.12 24.35 13.43 0 0 0
gb 12.91 47.93 28.89 10.26 0 0 0
hu 0.37 9.54 31.10 58.98 0 0 0
ie 14.00 44.32 29.63 12.04 0 0 0
il 6.77 30.86 26.40 35.98 0 0 0
is 40.45 41.25 16.70 1.59 0 0 0
it 10.89 36.63 34.96 17.52 0 0 0
lt 6.22 37.46 32.19 24.13 0 0 0
nl 13.62 49.49 28.38 8.51 0 0 0
no 23.17 52.49 22.07 2.27 0 0 0
pl 7.20 41.44 37.78 13.58 0 0 0
pt 14.57 57.09 20.24 8.11 0 0 0
ru 5.02 28.81 32.10 34.07 0 0 0
se 32.75 54.48 10.64 2.13 0 0 0
si 8.95 45.22 32.59 13.24 0 0 0

Le graphique de base est le diagramme en barre avec les effectifs. La version la plus simple du diagramme en barre bi-dimensionnel est la suivante :

p <- ggplot(extract_2, aes(x = country, fill = impcntr)) +
  geom_bar() 
p

Pour mieux voir le lien entre les deux variables, il est plus intéressant d’afficher les proportions (ou profils-lignes). Pour cela, il faut créer une table qui contient une colonne :

  • avec le nom du pays
  • avec les modalités de impcntr
  • les proportions correspondant au croisement pays x impcntr

Pour faire cela :

graph_1 <- extract_2 %>%
  select(country, impcntr) %>%
  group_by(country, impcntr) %>%
  summarise(n = n()) %>%
  group_by(country) %>%
  mutate(prop = n / sum(n))

On représente ensuite le diagramme :

p <- ggplot(graph_1, aes(x = country, y = prop, fill = impcntr)) +
  geom_bar(stat="identity") 
p

Pour ordonner les pays en fonction de ceux qui sont le plus contre : pour cela, il faut changer l’ordre des niveaux dans le factor :

graph_1$country <- factor(graph_1$ country, levels = graph_1 %>%
  filter(impcntr == "Allow none") %>%
  arrange(-prop) %>%
  select(country) %>%
  pull())

Enfin, pour avoir un graphique prêt à la publication, on ajoute les légendes et les drapeaux :

p <- ggplot(graph_1, aes(x = country, y = prop, fill = impcntr, 
                         country = as.character(country))) +
  geom_bar(stat="identity") 
p + labs(title = "Opinion sur l'accueil des immigrants en Europe",
         caption = "Source: ESS 2016",
         x = "Pays", y = "Proportions") + 
  scale_fill_discrete(name = "Opinion", 
        labels = c("Favorable", "Moyennement", "Peu", "Pas du tout")) + 
  geom_flag(aes(x = country, y = 0))

5.2 Croisement d’une variable quantitative et une variable qualitative

On peut utiliser la boîte à moustache parallèle :

p <- ggplot(data = extract_2, aes(x = happy, y = agea)) +
  geom_boxplot()
p

Une alternative est de représenter les distributions non paramétriques de la variable quantitative pour chaque modalité de la variable qualitative :

library(ggridges)
ggplot(extract_2, 
       aes(x = agea, 
           y = happy, 
           fill = happy)) +
  geom_density_ridges() + 
  theme_ridges() +
  labs("Age by levels happy") +
  theme(legend.position = "none")
## Picking joint bandwidth of 3.75

Souvent, on cherche à comparer les moyennes dans chaque groupe. Dans un premier temps, on peut représenter les moyennes sous forme d’un dot chart ou graphique de Cleveland. Pour cela, on calcule la moyenne de la confiance en la police par pays (on fait la même chose pour les 3 autres variables de confiances pour la suite de la session) et on trie les observations:

plotdata <- extract_2 %>%
  group_by(country) %>%
  summarize(trstplc = mean(trstplc_quanti),
            trstlgl = mean(trstlgl_quanti),
            trstprt = mean(trstprt_quanti),
            trstplt = mean(trstplt_quanti))

On représente ensuite le dotchart de la façon suivante :

ggplot(plotdata, 
       aes(x = trstplc, 
           y = reorder(country, 
                       trstplc),
           country = country)) +
  geom_point() + 
  geom_flag(aes(x = 0, y = country)) + 
  labs(title = "Confiance en la police par pays",
       caption = "Source: ESS 2016",
       x = "Note moyenne", y = "Pays")

Si on veut également les moyennes ainsi que les barres d’erreurs (qui représentent ici un intervalle de confiance en supposant une distribution gaussienne des données dans chaque groupe avec des variances égales), on calcule d’abord les moyennes et écarts-types dans chaque groupe groupe. Ici, on s’intéresse au niveau du bohneur en fonction de l’âge :

plot_mean <- extract_2 %>%
  group_by(happy) %>%
  summarize(n = n(),
            mean = mean(agea),
            sd = sd(agea),
            se = sd / sqrt(n),
            ci = qt(0.975, df = n - 1) * sd / sqrt(n))

On représente ensuite les moyennes avec leurs intervalles de confiance

ggplot(plot_mean, aes(x = happy,
                     y = mean)) +
  geom_point(size = 3, col = "red") +
  geom_line(size = 1) +
  geom_errorbar(aes(ymin = mean - ci, 
                    ymax = mean + ci), 
                width = .1, col = "blue")

On peut filter uniquement les habitants d’Allemagne et de France et regarder si l’avis est le même selon le pays :

plot_mean <- extract_2 %>%
  filter(country %in% c("fr", "de")) %>%
  group_by(happy, country) %>%
  summarize(n = n(),
            mean = mean(agea),
            sd = sd(agea),
            se = sd / sqrt(n),
            ci = qt(0.975, df = n - 1) * sd / sqrt(n))
ggplot(plot_mean, aes(x = happy,
                     y = mean,
                     color = country,
                     group = country)) +
  geom_point(size = 3) +
  geom_line(size = 1) +
  geom_errorbar(aes(ymin = mean - ci, 
                    ymax = mean + ci), 
                width = .1)

5.3 Lien entre deux variables quantitatives

Pour représenter le lien entre deux variables quantitatives, on utilise le nuage de points. Ici, on s’intéresse aux données de confiance aggrégées par pays.

p1 <- ggplot(plotdata, 
        aes(x = trstprt, y = trstplt)) + 
  # geom_flag() +
  geom_point() + 
  geom_smooth(method = "lm", col = "blue") + 
  xlab("Confiance aux partis politiques") +                 
  ylab("Confiance aux hommes politiques") +         
  ggtitle("Confiance aux institutions") +
  geom_flag(data = plotdata, aes(x = trstprt, 
                                 y = trstplt, country = country))

p2 <- ggplot(plotdata, 
             aes(x = trstlgl, y = trstplc)) + 
  # geom_flag() +
  geom_point() + 
  geom_smooth(method = "lm", col = "blue") + 
  xlab("Confiance à la justice") +                  
  ylab("Confiance à la police") +         
  ggtitle("Confiance aux institutions") +
  geom_flag(data = plotdata, aes(x = trstlgl, y = trstplc, country = country))
gridExtra::grid.arrange(p1, p2, ncol = 2)