04-bivarie.Rmd

# Analyse bivariée

Faire une analyse bivariée, c'est étudier la relation entre deux variables : sont-elles liées ? les valeurs de l'une influencent-elles les valeurs de l'autre ? ou sont-elles au contraire indépendantes ?

À noter qu'on va parler ici d'influence ou de lien, mais pas de relation de cause à effet : les outils présentés permettent de visualiser ou de déterminer une relation, mais des liens de causalité proprement dit sont plus difficiles à mettre en évidence. Il faut en effet vérifier que c'est bien telle variable qui influence telle autre et pas l'inverse, qu'il n'y a pas de "variable cachée", etc.

Là encore, le type d'analyse ou de visualisation est déterminé par la nature qualitative ou quantitative des deux variables.


## Croisement de deux variables qualitatives

### Tableaux croisés

On va continuer à travailler avec le jeu de données tiré de l'enquête *Histoire de vie* inclus dans l'extension `questionr`. On commence donc par charger l'extension, le jeu de données, et à le renommer en un nom plus court pour gagner un peu de temps de saisie au clavier :

```{r}
library(questionr)
data(hdv2003)
d <- hdv2003
```


Quand on veut croiser deux variables qualitatives, on fait un *tableau croisé*. Comme pour un tri à plat ceci s'obtient avec la fonction `table` de R, mais à laquelle on passe cette fois deux variables en argument. Par exemple, si on veut croiser la catégorie socio professionnelle et le sexe des enquêtés :

```{r}
table(d$qualif, d$sexe)
```

Pour pouvoir interpréter ce tableau on doit passer du tableau en effectifs au tableau en pourcentages ligne ou colonne. Pour cela, on peut utiliser les fonctions `lprop` et `cprop` de l'extension `questionr`, qu'on applique au tableau croisé précédent.

Pour calculer les pourcentages ligne :

```{r}
tab <- table(d$qualif, d$sexe)
lprop(tab)
```

Et pour les pourcentages colonne :

```{r}
cprop(tab)
```

```{block, type='rmdnote'}
Pour savoir si on doit faire des pourcentages ligne ou colonne, on pourra se référer à l'article suivant :

http://alain-leger.lescigales.org/textes/lignecolonne.pdf

En résumé, quand on fait un tableau croisé, celui-ci est parfaitement symétrique : on peut inverser les lignes et les colonnes, ça ne change pas son interprétation. Par contre, on a toujours en tête un "sens" de lecture dans le sens où on considère que l'une des variables *dépend* de l'autre. Par exemple, si on croise sexe et type de profession, on dira que le type de profession dépend du sexe, et non l'inverse : le type de profession est alors la variable *dépendante* (à expliquer), et le sexe la variable *indépendante* (explicative).

Pour faciliter la lecture d'un tableau croisé, il est recommandé de **faire les pourcentages sur la variable indépendante**. Dans notre exemple, la variable indépendante est le sexe, elle est en colonne, on calcule donc les pourcentages colonnes qui permettent de comparer directement, pour chaque sexe, la répartition des catégories socio-professionnelles.
```

### Test du χ²

Comme on travaille sur un échantillon et pas sur une population entière, on peut compléter ce tableau croisé par un test d'indépendance du χ². Celui-ci permet de rejeter l'hypothèse d'indépendance des lignes et des colonnes du tableau, c'est à dire de rejeter l'hypothèse que les écarts à l'indépendance observés seraient uniquement dus au biais d'échantillonnage (au fait qu'on n'a pas interrogé toute notre population).

Pour effectuer un test de ce type, on applique la fonction `chisq.test` au tableau croisé calculé précédemment :

```{r}
chisq.test(tab)
```

Le résultat nous indique trois valeurs : 

- `X-squared`, la valeur de la statistique du χ² pour notre tableau, c'est-à-dire une "distance" entre notre tableau observé et celui attendu si les deux variables étaient indépendantes.
- `df`, le nombre de degrés de libertés du test.
- `p-value`, le fameux *p*, qui indique la probabilité d'obtenir une valeur de la statistique du χ² au moins aussi extrême sous l'hypothèse d'indépendance.

Ici, le *p* est extrêmement petit (la notation `< 2.2e-16` indique qu'il est plus petit que la plus petite valeur calculable par R), donc certainement en-dessous du seuil de décision choisi préalablement au test (souvent 5%, soit 0.05). On peut donc rejeter l'hypothèse d'indépendance des lignes et des colonnes du tableau.

En complément du test du χ², on peut aussi regarder les *résidus* de ce test pour affiner la lecture du tableau. Ceux-ci s'obtiennent avec la fonction `chisq.residuals` de `questionr` :

```{r}
chisq.residuals(tab)
```

L'interprétation des résidus est la suivante :

- si la valeur du résidu pour une case est inférieure à -2, alors il y a une sous-représentation de cette case dans le tableau : les effectifs sont significativement plus faibles que ceux attendus sous l'hypothèse d'indépendance
- à l'inverse, si le résidu est supérieur à 2, il y a sur-représentatation de cette case
- si le résidu est compris entre -2 et 2, il n'y a pas d'écart à l'indépendance significatif

Les résidus peuvent être une aide utile à l'interprétation, notamment pour des tableaux de grande dimension.

### Représentation graphique

Il est possible de faire une représentation graphique d'un tableau croisé, par exemple avec la fonction `mosaicplot` :

```{r fig.height=6, fig.width=6}
mosaicplot(tab)
```



On peut améliorer ce graphique en colorant les cases selon les résidus du test du χ² (argument `shade = TRUE`) et en orientant verticalement les labels de colonnes (argument `las = 3`) :

```{r fig.height=6, fig.width=6}
mosaicplot(tab, las = 3, shade = TRUE)
```

Chaque rectangle de ce graphique
représente une case de tableau. Sa largeur correspond au pourcentage
des modalités en colonnes (il y'a beaucoup d'employés et d'ouvriers et
très peu d'"autres"). Sa hauteur correspond aux
pourcentages colonnes : la proportion d'hommes chez les cadres est
plus élevée que chez les employés. Enfin, la couleur de la case
correspond au résidu du test du χ² correspondant : les cases en
rouge sont sous-représentées, les cases en bleu sur-représentées, et
les cases blanches sont proches des effectifs attendus sous l'hypothèse
d'indépendance.




## Croisement d'une variable quantitative et d'une variable qualitative

### Représentation graphique

Croiser une variable quantitative et une variable qualitative, c'est essayer de voir si les valeurs de la variable quantitative se répartissent différemment selon la catégorie d'appartenance de la variable qualitative.

Pour cela, l'idéal est de commencer par une représentation graphique de type "boîte à moustache" à l'aide de la fonction `boxplot`. Par exemple, si on veut visualiser la répartition des âges selon la pratique ou non d'un sport, on va utiliser la syntaxe suivante :

```{r eval=FALSE}
boxplot(d$age ~ d$sport)
```

```{block type='rmdnote'}
Cette syntaxe de `boxplot` utilise une nouvelle notation de type "formule". Celle-ci est utilisée notamment pour la spécification des modèles de régression. Ici le `~` peut se lire comme "en fonction de" : on veut représenter le boxplot de l'âge en fonction du sport.
```

Ce qui va nous donner le résultat suivant :

```{r echo=FALSE}
boxplot(d$age ~ d$sport)
```





```{block type='rmdnote'}
L'interprétation d'un boxplot est la suivante : Les bords inférieurs et supérieurs du carré central représentent le premier et le troisième quartile de la variable représentée sur l'axe vertical. On a donc 50% de nos observations dans cet intervalle. Le trait horizontal dans le carré représente la médiane. Enfin, des "moustaches" s'étendent de chaque côté du carré, jusqu'aux valeurs minimales et maximales, avec une exception : si des valeurs sont éloignées du carré de plus de 1,5 fois l'écart interquartile (la hauteur du carré), alors on les représente sous forme de points (symbolisant des valeurs considérées comme "extrêmes").
```

Dans le graphique ci-dessus, on voit que ceux qui ont pratiqué un sport au cours des douze derniers mois ont l'air d'être sensiblement plus jeunes que les autres.


### Calculs d'indicateurs

On peut aussi vouloir comparer certains indicateurs (moyenne, médiane) d'une variable quantitative selon les modalités d'une variable qualitative. Si on reprend l'exemple précédent, on peut calculer la moyenne d'âge pour ceux qui pratiquent un sport et pour ceux qui n'en pratiquent pas.

Une première méthode pour cela est d'extraire de notre population autant de sous-populations qu'il y a de modalités dans la variable qualitative. On peut le faire notamment avec la fonction `filter` du package `dplyr` ^[Le package en question est présenté en détail dans la partie \@ref(dplyr).].

On commence par charger `dplyr` (en l'ayant préalablement installé) :

```{r, message=FALSE}
library(dplyr)
```

Puis on applique `filter` pour créer deux sous-populations, stockées dans deux nouveaux tableaux de données :

```{r}
d_sport <- filter(d, sport == "Oui")
d_nonsport <- filter(d, sport == "Non")
```

On peut ensuite utiliser ces deux nouveaux tableaux de données comme on en a l'habitude, et calculer les deux moyennes d'âge :

```{r}
mean(d_sport$age)
```
```{r}
mean(d_nonsport$age)
```

Une autre possibilité est d'utiliser la fonction `tapply`, qui prend en paramètre une variable quantitative, une variable qualitative et une fonction, puis applique automatiquement la fonction aux valeurs de la variables quantitative pour chaque niveau de la variable qualitative :

```{r}
tapply(d$age, d$sport, mean)
```

On verra dans la partie \@ref(dplyr) d'autres méthodes basées sur `dplyr` pour effectuer ce genre d'opérations.


### Tests statistiques

Un des tests les plus connus est le test du *t* de Student, qui permet de tester si les moyennes de deux sous-populations peuvent être considérées comme différentes (compte tenu des fluctuations aléatoires provenant du biais d'échantillonnage).

Un test *t* s'effectue à l'aide de la fonction `t.test`. Ainsi, on peut tester l'hypothèse d'égalité des âges moyens selon la pratique ou non d'un sport avec la commande suivante :

```{r}
t.test(d$age ~ d$sport)
```

Le résultat du test est significatif, avec un *p* extrêmement petit, et on peut rejeter l'hypothèse nulle d'égalité des moyennes des deux groupes. Le test nous donne même un intervalle de confiance à 95% pour la valeur de la différence entre les deux moyennes.

Nous sommes cependant allés un peu vite, et avons négligé le fait que le test *t* s'applique normalement à des distributions normales. On peut se faire un premier aperçu visuel de cette normalité en traçant les histogrammes des deux répartitions :

```{r}
hist(d_sport$age)
```

```{r}
hist(d_nonsport$age)
```

Si l'âge dans le groupe des non sportifs se rapproche d'une distribution normale, celui des sportifs en semble assez éloigné, notamment du fait de la limite d'âge à 18 ans imposée par construction de l'enquête.

On peut tester cette normalité à l'aide du test de Shapiro-Wilk et de la fonction `shapiro.test` : 


```{r}
shapiro.test(d_sport$age)
```

```{r}
shapiro.test(d_nonsport$age)
```

Le test est significatif dans les deux cas et rejette l'hypothèse d'une normalité des deux distributions.

Dans ce cas on peut faire appel à un test non-paramétrique, qui ne fait donc pas d'hypothèses sur les lois de distribution des variables testées, en l'occurrence le test des rangs de Wilcoxon, à l'aide de la fonction `wilcox.test` :

```{r}
wilcox.test(d$age ~ d$sport)
```

La valeur *p* étant à nouveau extrêmement petite, on peut rejeter l'hypothèse d'indépendance et considérer que les distributions des âges dans les deux sous-populations sont différentes.



## Croisement de deux variables quantitatives

Le jeu de données `hdv2003` comportant assez peu de variables quantitatives, on va s'intéresser maintenant à un autre jeu de données comportant des informations du recensement de la population de 2012. On le charge avec :

```{r}
data(rp2012)
```

Un nouveau tableau de données `rp2012` devrait apparaître dans votre environnement. Celui-ci comprend les 5170 communes de Métropole de plus de 2000 habitants, et une soixantaine de variables telles que le département, la population, le taux de chômage, etc. Pour une description plus complète et une liste des variables, voir section \@ref(rp2012).


### Représentation graphique

Quand on croise deux variables quantitatives, l'idéal est de faire une représentation graphique sous forme de nuage de points à l'aide de la fonction `plot`. On va représenter le croisement entre le pourcentage de cadres et le pourcentage de diplômés du supérieur dans la commune :

```{r}
plot(rp2012$cadres, rp2012$proprio)
```

Une représentation graphique est l'idéal pour visualiser l'existence d'un lien entre les deux variables. Voici quelques exemples d'interprétation :

```{r, echo = FALSE, fig.height=10, fig.width=6}
par(mfrow=c(3,2))

x <- rnorm(100)
y <- 2*x + 1 + rnorm(100, 0, 0.4)
plot(x, y, main = "Dépendance linéaire positive", 
     xlab = "", ylab = "", col="red")

x <- rnorm(100)
y <- -3*x + 15 + rnorm(100, 0, 0.4)
plot(x, y, main = "Dépendance linéaire négative", 
     xlab = "", ylab = "", col="red")

x <- rnorm(100)
y <- exp(x) + 1 + rnorm(100, 0, 0.4)
plot(x, y, main = "Dépendance non-linéaire monotone", 
     xlab = "", ylab = "", col="red")

x <- rnorm(100)
y <- 2*x^2 + 1 + rnorm(100, 0, 0.4)
plot(x, y, main = "Dépendance non-linéaire non monotone", 
     xlab = "", ylab = "", col="red")

x <- rnorm(100)
y <- rnorm(100)
plot(x, y, main = "Indépendance", 
     xlab = "", ylab = "", col="red")

x <- rnorm(100)
y <- rnorm(100, 0, 0.03)
plot(x, y, main = "Indépendance", 
     xlab = "", ylab = "", col="red", ylim = c(-2,2))

```


Dans ce premier graphique généré sur nos données, il semble difficile de mettre en évidence une relation de dépendance. Si par contre on croise le pourcentage de cadres et celui de diplômés du supérieur, on obtient une belle relation de dépendance linéaire.

```{r}
plot(rp2012$cadres, rp2012$dipl_sup)
```



### Calcul d'indicateurs 

En plus d'une représentation graphique, on peut calculer certains indicateurs permettant de mesurer le degré d'association de deux variables quantitatives.

#### Corrélation linéaire (Pearson)

La corrélation est une mesure du lien d'association *linéaire* entre deux variables quantitatives. Sa valeur varie entre -1 et 1. Si la corrélation vaut -1, il s'agit d'une association linéaire négative parfaite. Si elle vaut 1, il s'agit d'une association linéaire positive parfaite. Si elle vaut 0, il n'y a aucune association linéaire entre les variables.

On la calcule dans R à l'aide de la fonction `cor`.

Ainsi la corrélation entre le pourcentage de cadres et celui de diplômés du supérieur vaut :

```{r}
cor(rp2012$cadres, rp2012$dipl_sup)
```

Ce qui est extrêmement fort. Il y a donc un lien linéaire et positif entre les deux variables (quand la valeur de l'une augmente, la valeur de l'autre augmente également).

À l'inverse, la corrélation entre le pourcentage de cadres et le pourcentage de propriétaires vaut :

```{r}
cor(rp2012$cadres, rp2012$proprio)
```

Ce qui indique, pour nos données, une absence de liaison linéaire entre les deux variables.




#### Corrélation des rangs (Spearman)

Le coefficient de corrélation de Pearson ci-dessus fait une hypothèse forte sur les données : elles doivent être liées par une association linéaire. Quand ça n'est pas le cas mais qu'on est en présence d'une association monotone, on peut utiliser un autre coefficient, le coefficient de corrélation des rangs de Spearman.

Plutôt que de se baser sur les valeurs des variables, cette corrélation va se baser sur leurs rangs, c'est-à-dire sur leur position parmi les différentes valeurs prises par les variables.

Ainsi, si la valeur la plus basse de la première variable est associée à la valeur la plus basse de la deuxième, et ainsi de suite jusqu'à la valeur la plus haute, on obtiendra une corrélation de 1. Si la valeur la plus forte de la première variable est associée à la valeur la plus faible de la seconde, et ainsi de suite, et que la valeur la plus faible de la première est associée à la plus forte de la deuxième, on obtiendra une corrélation de -1. Si les rangs sont "mélangés", sans rapports entre eux, on obtiendra une corrélation autour de 0.


```{r, echo = FALSE, fig.height=10, fig.width=6}
par(mfrow=c(3,2))

x <- rnorm(100)
y <- 2*x + 1 + rnorm(100, 0, 0.4)
p <- round(cor(x, y), 2)
s <- round(cor(x, y, method = "spearman"), 2)
title <- paste0("Pearson : ", p, " - Spearman : ", s)  
plot(x, y, main = title, 
     xlab = "", ylab = "", col="red")

x <- rnorm(100)
y <- -3*x + 15 + rnorm(100, 0, 0.4)
p <- round(cor(x, y), 2)
s <- round(cor(x, y, method = "spearman"), 2)
title <- paste0("Pearson : ", p, " - Spearman : ", s)  
plot(x, y, main = title, 
     xlab = "", ylab = "", col="red")

x <- rnorm(100)
y <- exp(x) + 1 + rnorm(100, 0, 0.4)
p <- round(cor(x, y), 2)
s <- round(cor(x, y, method = "spearman"), 2)
title <- paste0("Pearson : ", p, " - Spearman : ", s)  
plot(x, y, main = title, 
     xlab = "", ylab = "", col="red")

x <- rnorm(100)
y <- 2*x^2 + 1 + rnorm(100, 0, 0.4)
p <- round(cor(x, y), 2)
s <- round(cor(x, y, method = "spearman"), 2)
title <- paste0("Pearson : ", p, " - Spearman : ", s)  
plot(x, y, main = title, 
     xlab = "", ylab = "", col="red")

x <- rnorm(100)
y <- rnorm(100)
p <- round(cor(x, y), 2)
s <- round(cor(x, y, method = "spearman"), 2)
title <- paste0("Pearson : ", p, " - Spearman : ", s)  
plot(x, y, main = title, 
     xlab = "", ylab = "", col="red")

x <- rnorm(20)
y <- x + 1 + rnorm(20, 0, 0.4)
x <- c(x, 2, 1.8)
y <- c(y, -2, -1.9)
p <- round(cor(x, y), 2)
s <- round(cor(x, y, method = "spearman"), 2)
title <- paste0("Pearson : ", p, " - Spearman : ", s)  
plot(x, y, main = title, 
     xlab = "", ylab = "", col="red", ylim = c(-2,2))



```

La corrélation des rangs a aussi pour avantage d'être moins sensibles aux valeurs extrêmes ou aux points isolés. On dit qu'elle est plus "robuste".

Pour calculer une corrélation de Spearman, on utilise la fonction `cor` mais avec l'argument `method = "spearman"` :

```{r}
cor(rp2012$cadres, rp2012$dipl_sup, method = "spearman")
```

### Régression linéaire

Quand on est en présence d'une association linéaire entre deux variables, on peut vouloir faire la régression linéaire d'une des variables sur l'autres. 

Une régression linéaire simple se fait à l'aide de la fonction `lm` :

```{r}
lm(rp2012$cadres ~ rp2012$dipl_sup)
```

```{block type='rmdnote'}
On retrouve avec `lm` la syntaxe "formule" déjà rencontrée avec `boxplot`. Elle permet ici de spécifier des modèles de régression : la variable dépendante se place à gauche du `~`, et la variable indépendante à droite. Si on souhaite faire une régression multiple avec plusieurs variables indépendantes, on aura une formule du type `dep ~ indep1 + indep2`. Il est également possible de spécifier des termes plus complexes, des interactions, etc.
```

`lm` nous renvoit par défaut les coefficients de la droite de régression :

- l'ordonnée à l'origine `(Intercept)` vaut 0.92
- le coefficient associé à `dipl_sup` vaut 1.08

Pour des résultats plus détaillés, on peut stocker le résultat de la régression dans un objet et utiliser la fonction `summary` :

```{r}
reg <- lm(rp2012$cadres ~ rp2012$dipl_sup)
summary(reg)
```

Ces résultats montrent notamment que les coefficients sont significativement
différents de 0. La part de cadres augmente donc bien avec celle de diplômés du supérieur. 

On peut enfin représenter la droite de régression sur notre nuage de points à l'aide de la fonction `abline` :

```{r}
plot(rp2012$dipl_sup, rp2012$cadres)
abline(reg, col="red")
```


## Exercices

**Exercice 1**

Dans le jeu de données `hdv2003`, faire le tableau croisé entre la catégorie socio-professionnelle (variable `qualif`) et le fait de croire ou non en l'existence des classes sociales (variable `clso`). Identifier la variable indépendante et la variable dépendante, et calculer les pourcentages ligne ou colonne. Interpréter le résultat.

```{r include=FALSE, eval=FALSE}
library(questionr)
data(hdv2003)
tab <- table(hdv2003$qualif, hdv2003$clso)
tab
lprop(tab)
chisq.test(tab)
mosaicplot(tab, shade=TRUE)
```


<div class="solution">
\iffalse
```{r eval=FALSE}
library(questionr)
data(hdv2003)
tab <- table(hdv2003$qualif, hdv2003$clso)

## Ici la variable indépendante est `qualif`, on calcule donc 
## les pourcentages lignes
lprop(tab)
```
\fi

</div>

Faire un test du χ². Peut-on rejeter l'hypothèse d'indépendance ?

<div class="solution">
\iffalse
```{r eval=FALSE}
chisq.test(tab)

## On peut rejeter l'hypothèse d'indépendance au seuil de 0.05, 
## et même au seuil de 0.001
```
\fi

</div>

Représenter ce tableau croisé sous la forme d'un `mosaicplot` en colorant les cases selon les résidus du test du χ².

<div class="solution">
\iffalse
```{r eval=FALSE}
mosaicplot(tab, shade = TRUE)
```
\fi

</div>


**Exercice 2**

Toujours sur le jeu de données `hdv2003`, faire le boxplot qui croise le nombre d'heures passées devant la télévision (variable `heures.tv`) avec le statut d'occupation (variable `occup`).

<div class="solution">
\iffalse
```{r eval=FALSE}
boxplot(hdv2003$heures.tv ~ hdv2003$occup)
```
\fi

</div>


Calculer la durée moyenne en fonction du statut d'occupation à l'aide de `tapply`.

<div class="solution">
\iffalse
```{r eval=FALSE}
tapply(hdv2003$heures.tv, hdv2003$occup, mean, na.rm = TRUE)

```
\fi

</div>



**Exercice 3**

Sur le jeu de données `rp2012`, représenter le nuage de points croisant le pourcentage de personnes sans diplôme (variable `dipl_aucun`) et le pourcentage de propriétaires (variable `proprio`).

<div class="solution">
\iffalse
```{r eval=FALSE}
library(questionr)
data(rp2012)
plot(rp2012$dipl_aucun, rp2012$proprio)
```
\fi

</div>

Calculer le coefficient de corrélation linéaire correspondant.

<div class="solution">
\iffalse
```{r eval=FALSE}
cor(rp2012$dipl_aucun, rp2012$proprio)
```
\fi

</div>