La version originale de ce chapitre a été écrite par Julien Barnier dans le cadre du support de cours Introduction à R.

On entend par statistique bivariée l’étude des relations entre deux variables, celles-ci pouvant être quantitatives ou qualitatives. La statistique bivariée fait partie de la statistique descriptive.

La statistique univariée a quant à elle déjà été abordée dans un chapitre dédié.

Comme dans la partie précédente, on travaillera sur les jeux de données fournis avec l’extension questionr et tiré de l’enquête Histoire de vie et du recensement 1999 :

library(questionr)
data(hdv2003)
d <- hdv2003
data(rp99)

Deux variables quantitatives

La comparaison de deux variables quantitatives se fait en premier lieu graphiquement, en représentant l’ensemble des couples de valeurs. On peut ainsi représenter les valeurs du nombre d’heures passées devant la télévision selon l’âge.

plot(d$age, d$heures.tv)
Nombre d’heures de télévision selon l’âge

Le fait que des points sont superposés ne facilite pas la lecture du graphique. On peut utiliser une représentation avec des points semi-transparents.

plot(d$age, d$heures.tv, pch = 19, col = rgb(1, 0, 0, 0.1))
Nombre d’heures de télévision selon l’âge avec semi-transparence

Plus sophistiqué, on peut faire une estimation locale de densité et représenter le résultat sous forme de « carte ». Pour cela on commence par isoler les deux variables, supprimer les observations ayant au moins une valeur manquante à l’aide de la fonction complete.cases, estimer la densité locale à l’aide de la fonction kde2d de l’extension MASS1 et représenter le tout à l’aide d’une des fonctions image, contour ou filled.contour

library(MASS)
tmp <- d[, c("age", "heures.tv")]
tmp <- tmp[complete.cases(tmp), ]
filled.contour(kde2d(tmp$age, tmp$heures.tv), color = terrain.colors)
Représentation de l’estimation de densité locale

Une représentation alternative de la densité locale peut être obtenue avec la fonction smoothScatter.

smoothScatter(d[, c("age", "heures.tv")])
Représentation alternative de l’estimation de densité locale

Dans tous les cas, il n’y a pas de structure très nette qui semble se dégager. On peut tester ceci mathématiquement en calculant le coefficient de corrélation entre les deux variables à l’aide de la fonction cor :

cor(d$age, d$heures.tv, use = "complete.obs")
[1] 0.1776249

L’option use permet d’éliminer les observations pour lesquelles l’une des deux valeurs est manquante. Le coefficient de corrélation est très faible.

On va donc s’intéresser plutôt à deux variables présentes dans le jeu de données rp99, la part de diplômés du supérieur et la proportion de cadres dans les communes du Rhône en 1999.

À nouveau, commençons par représenter les deux variables.

plot(rp99$dipl.sup, rp99$cadres, ylab = "Part des cadres", xlab = "Part des diplomês du supérieur")
Proportion de cadres et proportion de diplômés du supérieur

Ça ressemble déjà beaucoup plus à une relation de type linéaire.

Calculons le coefficient de corrélation :

cor(rp99$dipl.sup, rp99$cadres)
[1] 0.8975282

C’est beaucoup plus proche de 1. On peut alors effectuer une régression linéaire complète en utilisant la fonction lm :

reg <- lm(cadres ~ dipl.sup, data = rp99)
summary(reg)

Call:
lm(formula = cadres ~ dipl.sup, data = rp99)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.6905 -1.9010 -0.1823  1.4913 17.0866 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.24088    0.32988   3.762 0.000203 ***
dipl.sup     1.38352    0.03931  35.196  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.281 on 299 degrees of freedom
Multiple R-squared:  0.8056,    Adjusted R-squared:  0.8049 
F-statistic:  1239 on 1 and 299 DF,  p-value: < 2.2e-16

Le résultat montre que les coefficients sont significativement différents de 0. La part de cadres augmente donc avec celle de diplômés du supérieur (ô surprise). On peut très facilement représenter la droite de régression à l’aide de la fonction abline.

plot(rp99$dipl.sup, rp99$cadres, ylab = "Part des cadres", xlab = "Part des diplômés du supérieur")
abline(reg, col = "red")
Régression de la proportion de cadres par celle de diplômés du supérieur

On remarquera que le premier argument passé à la fonction lm a une syntaxe un peu particulière. Il s’agit d’une formule, utilisée de manière générale dans les modèles statistiques. On indique la variable d’intérêt à gauche et la variable explicative à droite, les deux étant séparées par un tilde (obtenu sous Windows en appuyant simultanément sur les touches Alt Gr et 2). On remarquera que les noms des colonnes de notre tableau de données ont été écrites sans guillemets.

Dans le cas présent, nous avons calculé une régression linéaire simple entre deux variables, d’où l’écriture cadres ∼ dipl.sup. Si nous avions voulu expliquer une variable z par deux variables x et y, nous aurions écrit z ∼ x + y. Il est possible de spécifier des modèles encore plus complexes.

Pour un aperçu de la syntaxe des formules sous R, voir http://ww2.coastal.edu/kingw/statistics/R-tutorials/formulae.html.

Trois variables ou plus

Lorsque l’on souhaite représenter trois variables quantitatives simultanément, il est possible de réaliser un nuage de points représentant les deux premières variables sur l’axe horizontal et l’axe vertical et en faisant varier la taille des points selon la troisième variable, en utilisant l’argument cex de la fonction plot.

plot(rp99$dipl.aucun, rp99$tx.chom, cex = rp99$pop.tot/10^4)
Nuage de points avec taille des points proportionnels à une troisième variable

Lorsque l’on étudie un plus grand nombres de variables quantitatives, il est peut être utile de réaliser une matrice de nuages de points, qui compare chaque variable deux à deux et qui s’obtient facilement avec la fonction pairs.

pairs(rp99[, c("proprio", "hlm", "locataire", "maison")])
Matrice de nuages de points

Une variable quantitative et une variable qualitative

Représentations graphiques

Quand on parle de comparaison entre une variable quantitative et une variable qualitative, on veut en général savoir si la distribution des valeurs de la variable quantitative est la même selon les modalités de la variable qualitative. En clair : est ce que l’âge de ceux qui écoutent du hard rock est différent de l’âge de ceux qui n’en écoutent pas ?

Là encore, l’idéal est de commencer par une représentation graphique. Les boîtes à moustaches (boxplot en anglais) sont parfaitement adaptées pour cela.

Si on a construit des sous-populations d’individus écoutant ou non du hard rock, on peut utiliser la fonction boxplot.

d.hard <- subset(d, hard.rock == "Oui")
d.non.hard <- subset(d, hard.rock == "Non")
boxplot(d.hard$age, d.non.hard$age)
Boxplot de la répartition des âges (sous-populations)

Mais construire les sous-populations n’est pas nécessaire. On peut utiliser directement la version de boxplot prenant une formule en argument.

boxplot(age ~ hard.rock, data = d)
Boxplot de la répartition des âges (formule)

À première vue, ô surprise, la population écoutant du hard rock a l’air sensiblement plus jeune. Peut-on le tester mathématiquement ?

Tests statistiques

On peut calculer la moyenne d’âge des deux groupes en utilisant la fonction tapply2 :

tapply(d$age, d$hard.rock, mean)
     Non      Oui 
48.30211 27.57143 

L’écart est important. Est-il statistiquement significatif ? Pour cela on peut faire un test t de Student comparaison de moyennes à l’aide de la fonction t.test :

t.test(d$age ~ d$hard.rock)

    Welch Two Sample t-test

data:  d$age by d$hard.rock
t = 9.6404, df = 13.848, p-value = 1.611e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 16.11379 25.34758
sample estimates:
mean in group Non mean in group Oui 
         48.30211          27.57143 

Le test est extrêmement significatif. L’intervalle de confiance à 95 % de la différence entre les deux moyennes va de 14,5 ans à 21,8 ans.

La valeur affichée pour p est de 1.611e-07. Cette valeur peut paraître étrange pour les non avertis. Cela signifie tout simplement 1,611 multiplié par 10 à la puissance -7, autrement dit 0,0000001611. Cette manière de représenter un nombre est couramment appelée notation scientifique.

Pour plus de détails, voir http://fr.wikipedia.org/wiki/Notation_scientifique.

Nous sommes cependant allés un peu vite en besogne, car nous avons négligé une hypothèse fondamentale du test t : les ensembles de valeur comparés doivent suivre approximativement une loi normale et être de même variance3. Comment le vérifier ?

D’abord avec un petit graphique composés de deux histogrammes :

par(mfrow = c(1, 2))
hist(d$age[d$hard.rock == "Oui"], main = "Hard rock", col = "red")
hist(d$age[d$hard.rock == "Non"], main = "Sans hard rock", col = "red")
Distribution des âges pour appréciation de la normalité

La fonction par permet de modifier de nombreux paramètres graphiques. par(mfrow = c(1, 2)) sert à indiquer que l’on souhaite afficher deux graphiques sur une même fenêtre, plus précisément que la fenêtre doit comporter une ligne et deux colonnes.

Ça a l’air à peu près bon pour les « Sans hard rock », mais un peu plus limite pour les fans de Metallica, dont les effectifs sont d’ailleurs assez faibles. Si on veut en avoir le coeur net on peut utiliser le test de normalité de Shapiro-Wilk avec la fonction shapiro.test :

shapiro.test(d$age[d$hard.rock == "Oui"])

    Shapiro-Wilk normality test

data:  d$age[d$hard.rock == "Oui"]
W = 0.86931, p-value = 0.04104
shapiro.test(d$age[d$hard.rock == "Non"])

    Shapiro-Wilk normality test

data:  d$age[d$hard.rock == "Non"]
W = 0.98141, p-value = 2.079e-15

Visiblement, le test estime que les distributions ne sont pas suffisamment proches de la normalité dans les deux cas.

Et concernant l’égalité des variances ?

tapply(d$age, d$hard.rock, var)
      Non       Oui 
285.62858  62.72527 

L’écart n’a pas l’air négligeable. On peut le vérifier avec le test d’égalité des variances fourni par la fonction var.test :

var.test(d$age ~ d$hard.rock)

    F test to compare two variances

data:  d$age by d$hard.rock
F = 4.5536, num df = 1985, denom df = 13, p-value = 0.003217
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 1.751826 8.694405
sample estimates:
ratio of variances 
          4.553644 

La différence est très significative. En toute rigueur le test t n’aurait donc pas pu être utilisé.

Damned ! Ces maudits tests statistiques vont-ils nous empêcher de faire connaître au monde entier notre fabuleuse découverte sur l’âge des fans de Sepultura ? Non ! Car voici qu’approche à l’horizon un nouveau test, connu sous le nom de Wilcoxon/Mann-Whitney. Celui-ci a l’avantage d’être non-paramétrique, c’est à dire de ne faire aucune hypothèse sur la distribution des échantillons comparés. Par contre il ne compare pas des différences de moyennes mais des différences de médianes :

wilcox.test(d$age ~ d$hard.rock)

    Wilcoxon rank sum test with continuity correction

data:  d$age by d$hard.rock
W = 23980, p-value = 2.856e-06
alternative hypothesis: true location shift is not equal to 0

Ouf ! La différence est hautement significative4. Nous allons donc pouvoir entamer la rédaction de notre article pour la Revue française de sociologie.

Deux variables qualitatives

La comparaison de deux variables qualitatives s’appelle en général un tableau croisé. C’est sans doute l’une des analyses les plus fréquentes lors du traitement d’enquêtes en sciences sociales.

Tableau croisé

La manière la plus simple d’obtenir un tableau croisé est d’utiliser la fonction table en lui donnant en paramètres les deux variables à croiser. En l’occurrence nous allons croiser un recodage du niveau de qualification regroupé avec le fait de pratiquer un sport.

On commence par calculer la variable recodée et par afficher le tri à plat des deux variables :

d$qualif2 <- as.character(d$qualif)
d$qualif2[d$qualif %in% c("Ouvrier specialise", "Ouvrier qualifie")] <- "Ouvrier"
d$qualif2[d$qualif %in% c("Profession intermediaire", "Technicien")] <- "Intermediaire"
table(d$qualif2)

        Autre         Cadre       Employe Intermediaire       Ouvrier 
           58           260           594           246           495 

Le tableau croisé des deux variables s’obtient de la manière suivante :

table(d$sport, d$qualif2)
     
      Autre Cadre Employe Intermediaire Ouvrier
  Non    38   117     401           127     381
  Oui    20   143     193           119     114

Il est tout à fait possible de croiser trois variables ou plus. Par exemple :

table(d$sport, d$cuisine, d$sexe)
, ,  = Homme

     
      Non Oui
  Non 401 129
  Oui 228 141

, ,  = Femme

     
      Non Oui
  Non 358 389
  Oui 132 222

Une alternative à la fonction table est la fonction xtabs. On indiquera à cette dernière le croisement à effectuer à l’aide d’une formule puis l’objet contenant nos données. Comme il ne s’agit pas d’un modèle avec une variable à expliquer, toutes les variables seront indiquées à la droite du symbole et séparées par +.

xtabs(~sport, d)
sport
 Non  Oui 
1277  723 
xtabs(~sport + cuisine, d)
     cuisine
sport Non Oui
  Non 759 518
  Oui 360 363
xtabs(~sport + cuisine + sexe, d)
, , sexe = Homme

     cuisine
sport Non Oui
  Non 401 129
  Oui 228 141

, , sexe = Femme

     cuisine
sport Non Oui
  Non 358 389
  Oui 132 222

On remarquera que le rendu par défaut est en général plus lisible car le nom des variables est indiqué, permettant de savoir quelle variable est affichée en colonnes et laquelle en lignes.

Si l’on utilise des données labellisées, la fonction xtabs ne prendra pas en compte les étiquettes de valeur.

data(fecondite)
xtabs(~educ + region, femmes)
    region
educ   1   2   3   4
   0 387 213 282 256
   1 179  53  86 142
   2 123  57  37 131
   3  18   1   2  33

On pourra alors utiliser la fonction ltabs de l’extension question, qui fonctionne exactement comme xtabs, à ceci près qu’elle prendra en compte les étiquettes de variable et de valeur quand elles existent.

ltabs(~educ + region, femmes)
                        region: Région de résidence
educ: Niveau d'éducation [1] Nord [2] Est [3] Sud [4] Ouest
          [0] aucun           387     213     282       256
          [1] primaire        179      53      86       142
          [2] secondaire      123      57      37       131
          [3] supérieur        18       1       2        33

Pourcentages en ligne et en colonne

On n’a cependant que les effectifs, ce qui rend difficile les comparaisons. L’extension questionr fournit des fonctions permettant de calculer facilement les pourcentages lignes, colonnes et totaux d’un tableau croisé.

Les pourcentages lignes s’obtiennent avec la fonction lprop5. Celle-ci s’applique au tableau croisé généré par table ou xtabs  :

tab <- table(d$sport, d$qualif2)
lprop(tab)
          
           Autre Cadre Employe Intermediaire Ouvrier Total
  Non        3.6  11.0  37.7    11.9          35.8   100.0
  Oui        3.4  24.3  32.8    20.2          19.4   100.0
  Ensemble   3.5  15.7  35.9    14.9          29.9   100.0
tab <- xtabs(~sport + qualif2, d)
lprop(tab)
          qualif2
sport      Autre Cadre Employe Intermediaire Ouvrier Total
  Non        3.6  11.0  37.7    11.9          35.8   100.0
  Oui        3.4  24.3  32.8    20.2          19.4   100.0
  Ensemble   3.5  15.7  35.9    14.9          29.9   100.0

Les pourcentages ligne ne nous intéressent guère ici. On ne cherche pas à voir quelle est la proportion de cadres parmi ceux qui pratiquent un sport, mais plutôt quelle est la proportion de sportifs chez les cadres. Il nous faut donc des pourcentages colonnes, que l’on obtient avec la fonction cprop :

cprop(tab)
       qualif2
sport   Autre Cadre Employe Intermediaire Ouvrier Ensemble
  Non    65.5  45.0  67.5    51.6          77.0    64.4   
  Oui    34.5  55.0  32.5    48.4          23.0    35.6   
  Total 100.0 100.0 100.0   100.0         100.0   100.0   

Dans l’ensemble, le pourcentage de personnes ayant pratiqué un sport est de 35,6 %. Mais cette proportion varie fortement d’une catégorie professionnelle à l’autre : 55,0 % chez les cadres contre 23,0 % chez les ouvriers.

Enfin, les pourcentage totaux s’obtiennent avec la fonction prop :

prop(tab)
       qualif2
sport   Autre Cadre Employe Intermediaire Ouvrier Total
  Non     2.3   7.1  24.3     7.7          23.0    64.4
  Oui     1.2   8.7  11.7     7.2           6.9    35.6
  Total   3.5  15.7  35.9    14.9          29.9   100.0

À noter qu’on peut personnaliser l’affichage de ces tableaux de pourcentages à l’aide de différentes options, dont digits qui règle le nombre de décimales à afficher et percent qui indique si on souhaite ou non rajouter un symbole % dans chaque case du tableau. Cette personnalisation peut se faire directement au moment de la génération du tableau et dans ce cas elle sera utilisée par défaut :

ctab <- cprop(tab, digits = 2, percent = TRUE)
ctab
       qualif2
sport   Autre   Cadre   Employe Intermediaire Ouvrier Ensemble
  Non    65.52%  45.00%  67.51%  51.63%        76.97%  64.37% 
  Oui    34.48%  55.00%  32.49%  48.37%        23.03%  35.63% 
  Total 100.00% 100.00% 100.00% 100.00%       100.00% 100.00% 

ou bien ponctuellement en passant les mêmes arguments à la fonction print :

ctab <- cprop(tab)
print(ctab, percent = TRUE)
       qualif2
sport   Autre  Cadre  Employe Intermediaire Ouvrier Ensemble
  Non    65.5%  45.0%  67.5%   51.6%         77.0%   64.4%  
  Oui    34.5%  55.0%  32.5%   48.4%         23.0%   35.6%  
  Total 100.0% 100.0% 100.0%  100.0%        100.0%  100.0%  

Représentation graphique

On peut obtenir une représentation graphique synthétisant l’ensemble des résultats obtenus sous la forme d’un graphique en mosaïque grâce à la fonction mosaicplot.

mosaicplot(qualif2 ~ sport, data = d, shade = TRUE, main = "Graphe en mosaïque")
Exemple de graphe en mosaïque

Comment interpréter ce graphique haut en couleurs6 ? Chaque rectangle représente une case de tableau. Sa largeur correspond aux pourcentages en colonnes (il y a beaucoup d’employés et d’ouvriers et très peu d’« Autre »). Sa hauteur correspond aux pourcentages en lignes : la proportion de sportifs 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 statistiquement proches de l’hypothèse d’indépendance.

Les graphiques en mosaïque permettent notamment de représenter des tableaux croisés à 3 ou 4 dimensions, voire plus.

L’extension vcd fournie une fonction mosaic fournissant plus d’options pour la création d’un graphique en mosaïque, permettant par exemple d’indiquer quelles variables doivent être affichées horizontalement ou verticalement, ou encore de colorier le contenu des rectangles en fonction d’une variable donnée, …

library(vcd)
mosaic(~sport + cuisine + sexe, d, highlighting = "sexe", main = "Exemple de graphique en mosaïque à 3 dimensions")

Lorsque l’on s’intéresse principalement aux variations d’une variable selon une autre, par exemple ici à la pratique du sport selon le niveau de qualification, il peut être intéressant de présenter les pourcentages en colonne sous la forme de barres cumulées.

barplot(cprop(tab, total = FALSE), main = "Pratique du sport selon le niveau de qualification")
Exemple de barres cumulées

Tests statistiques

Pour un test de comparaison de proportions, un test du Chi² ou encore un test exact de Fisher, on pourra se référer au chapitre dédié aux test statistiques de comparaison.


  1. MASS est installée par défaut avec la version de base de R.

  2. La fonction tapply est présentée plus en détails dans le chapitre Manipulation de données.

  3. Concernant cette seconde condition, t.test propose une option nommée var.equal qui permet d’utiliser une approximation dans le cas où les variances ne sont pas égales.

  4. Ce test peut également fournir un intervalle de confiance avec l’option conf.int=TRUE.

  5. Il s’agit en fait d’un alias pour les francophones de la fonction rprop.

  6. Sauf s’il est imprimé en noir et blanc…