5.3 Mise en œuvre dans R
Pour calculer le khi-deux entre deux variables qualitatives, nous utilisons la fonction de base :
chisq.test(x = ..., y = ...) qui renvoie le nombre de degrés de liberté, les valeurs du khi-deux et de p.
# Importation du csv
dataHLM <- read.csv("data/bivariee/hlm.csv")
# Calcul du khi-deux avec la fonction de base chisq.test
chisq.test(x = dataHLM$Taille, y = dataHLM$Periode)##
## Pearson's Chi-squared test
##
## data: dataHLM$Taille and dataHLM$Periode
## X-squared = 63.543, df = 12, p-value = 5.063e-09Pour la construction du tableau de contingence, deux options sont possibles dépendamment de la structure de votre tableau de données. Premier cas de figure : votre tableau comprend une ligne par observation avec les différentes modalités dans deux colonnes (ici Periode et Taille). Dans la syntaxe ci-dessous, pour chacune des deux variables qualitatives, nous créons un facteur afin de spécifier un intitulé à chaque modalité (factor(levels =c(....), labels = c(..)). Puis, nous utilisons la fonction CrossTable du package gmodels. Pour obtenir les fréquences théoriques, les contributions locales au khi-deux et les déviations, nous spécifions les options suivantes : expected=TRUE, chisq=TRUE, resid=TRUE.
library("gmodels")
# Premiers enregistrements du tableau
head(dataHLM)## Periode Taille
## 1 5 1
## 2 5 1
## 3 5 2
## 4 5 1
## 5 5 1
## 6 5 2# La variable Periode comprend 5 modalités (de 1 à 5)
table(dataHLM$Periode)##
## 1 2 3 4 5
## 56 48 48 47 80# La variable Taille comprend 4 modalités (de 1 à 4)
table(dataHLM$Taille)##
## 1 2 3 4
## 80 53 99 47# Création d'un facteur pour les cinq modalités de la période de construction
dataHLM$Periode <- factor(dataHLM$Periode,
levels = c(1,2,3,4,5),
labels = c("<1975",
"1975-1979",
"1980-1984",
"1985-1989",
"1990-1994"))
# Création d'un facteur pour les quatre modalités de la taille des habitations
dataHLM$Taille <- factor(dataHLM$Taille,
levels = c(1,2,3,4),
labels = c("<25 log.",
"25-49",
"50-99",
"100 et +"))
# Pour construire un tableau de contingence, nous utilisons
# la fonction CrossTable du package gmodels.
# Les deux lignes ci-dessous sont mises en commentaire pour ne pas répéter le tableau.
# CrossTable(x=dataHLM$Taille, y=dataHLM$Periode, digits = 2,
# expected=TRUE, chisq = TRUE, resid = TRUE, format="SPSS")Deuxième cas de figure : vous disposez déjà d’un tableau de contingence, soit les fréquences observées (\(f_{ij}\)). Nous n’utilisons donc pas la fonction CrossTable, mais directement la fonction chisq.test.
# Importation des données
df1 <- read.csv("data/bivariee/data_transport.csv", stringsAsFactors = FALSE)
df1 # Visualisation du tableau## ModeTransport Homme Femme
## 1 Automobile, camion ou fourgonnette - conducteur 689400 561830
## 2 Automobile, camion ou fourgonnette - passager 21315 40010
## 3 Transport en commun 181435 238330
## 4 A pied 43715 54360
## 5 Bicyclette 24295 13765
## 6 Autre moyen 8395 6970Matrice <- as.matrix(df1[, c("Homme","Femme")])
dimnames(Matrice) <- list(unique(df1$ModeTransport), Sexe=c("Homme","Femme"))
# Notez que vous pouvez saisir vos données directement si vous avez peu d'observations
Femme <- c(689400, 21315, 181435, 43715, 24295, 8395) # Vecteur de valeurs pour les femmes
Homme <- c(561830, 40010, 238330, 54360, 13765, 6970) # Vecteur de valeurs pour les hommes
Matrice <- as.table(cbind(Femme, Homme)) # Création du tableau
# Nom des deux variables et de leurs modalités respectives
dimnames(Matrice) <- list(Transport=c("Automobile (conducteur)",
"Automobile (passager)",
"Transport en commun",
"À pied",
"Bicyclette",
"Autre moyen"),
Sexe=c("Homme","Femme"))
# Test du khi-deux
test <- chisq.test(Matrice)
print(test)##
## Pearson's Chi-squared test
##
## data: Matrice
## X-squared = 29134, df = 5, p-value < 2.2e-16# Fréquences observées (Fij)
test$observed## Sexe
## Transport Homme Femme
## Automobile (conducteur) 689400 561830
## Automobile (passager) 21315 40010
## Transport en commun 181435 238330
## À pied 43715 54360
## Bicyclette 24295 13765
## Autre moyen 8395 6970# Fréquences théoriques (FTij)
round(test$expected,0)## Sexe
## Transport Homme Femme
## Automobile (conducteur) 643313 607917
## Automobile (passager) 31530 29795
## Transport en commun 215820 203945
## À pied 50425 47650
## Bicyclette 19568 18492
## Autre moyen 7900 7465# Déviations (Fij - FTij)
round(test$observed-test$expected,0)## Sexe
## Transport Homme Femme
## Automobile (conducteur) 46087 -46087
## Automobile (passager) -10215 10215
## Transport en commun -34385 34385
## À pied -6710 6710
## Bicyclette 4727 -4727
## Autre moyen 495 -495# Contributions au khi-deux
round((test$observed-test$expected)^2/test$expected,2)## Sexe
## Transport Homme Femme
## Automobile (conducteur) 3301.74 3493.98
## Automobile (passager) 3309.37 3502.05
## Transport en commun 5478.22 5797.18
## À pied 892.81 944.80
## Bicyclette 1141.71 1208.19
## Autre moyen 31.04 32.85# Marges en ligne et en colonne
colSums(Matrice)## Homme Femme
## 968555 915265rowSums(Matrice)## Automobile (conducteur) Automobile (passager) Transport en commun À pied
## 1251230 61325 419765 98075
## Bicyclette Autre moyen
## 38060 15365# Grand total
sum(Matrice)## [1] 1883820# Pourcentages
round(Matrice/sum(Matrice)*100,2)## Sexe
## Transport Homme Femme
## Automobile (conducteur) 36.60 29.82
## Automobile (passager) 1.13 2.12
## Transport en commun 9.63 12.65
## À pied 2.32 2.89
## Bicyclette 1.29 0.73
## Autre moyen 0.45 0.37# Pourcentages en ligne
round(Matrice/rowSums(Matrice)*100,2)## Sexe
## Transport Homme Femme
## Automobile (conducteur) 55.10 44.90
## Automobile (passager) 34.76 65.24
## Transport en commun 43.22 56.78
## À pied 44.57 55.43
## Bicyclette 63.83 36.17
## Autre moyen 54.64 45.36# Pourcentages en colonne
round(Matrice/colSums(Matrice)*100,2)## Sexe
## Transport Homme Femme
## Automobile (conducteur) 71.18 58.01
## Automobile (passager) 2.33 4.37
## Transport en commun 18.73 24.61
## À pied 4.78 5.94
## Bicyclette 2.51 1.42
## Autre moyen 0.92 0.76Pour obtenir les autres mesures d’association (tableau 5.2), nous pourrons utiliser la syntaxe suivante :
df1 <- read.csv("data/bivariee/hlm.csv")
# Fonction pour calculer les autres mesures d'association
AutresMesuresKhi2 <- function(x, y){
testChi2 <- chisq.test(x, y) # Calcul du khi-deux
n <- sum(testChi2$observed) # Nombre d'observations
nc <- ncol(testChi2$observed) # Nombre de colonnes
l <- nrow(testChi2$observed) # Nombre de lignes
dl <- (nc-1)*(l-1) # Nombre de degrés de libertés
chi2 <- testChi2$statistic # Valeur du khi-deux
Pchi2 <- testChi2$p.value # P pour le khi-deux
#Ratio de vraisemblance du khi-deux
G <- 2*sum(testChi2$observed*log(testChi2$observed/testChi2$expected)) # G2
PG <- pchisq(G, df=dl, lower.tail = FALSE) # P pour le G22
# khi-deux de Mantel-Haenszel avec le package DescTools
MHTest <- DescTools::MHChisqTest(testChi2$observed)
MH <- MHTest$statistic
PMH <- MHTest$p.value
# Coefficient de correlation polychorique
df1 <- data.frame("x" = as.factor(x),
"y" = as.factor(y))
polychoricCorr <- correlation::cor_test(df1,"x","y",method = "polychoric")
polyR <- polychoricCorr$rho
polyP <- polychoricCorr$p
# Coefficient Phi et V de Cramer
phi <- sqrt(chi2/n)
vc <- sqrt(chi2/(n*min(nc-1,l-1)))
# Tableau pour les résultats
dfsortie <- data.frame(
Statistique = c("Khi-deux",
"Ratio de vraisemblance du khi-deux",
"Khi-deux de Mantel-Haenszel",
"Corrélation Polychorique",
"Coefficient de Phi",
"V de Cramer"),
Valeur = round(c(chi2, G, MH, polyR, phi, vc),3),
P = round(c(Pchi2, PG, PMH, polyP , NA, NA),10))
return(dfsortie)
}
dfkhi2 <- AutresMesuresKhi2(df1$Periode, df1$Taille)
# Impression du tableau avec le package stargazer
library(stargazer)
stargazer(dfkhi2, type="text", summary=FALSE, rownames=FALSE, align = FALSE, digits = 3,
title="Mesures d'association entre les deux variables qualitatives")| Statistique | Valeur | P |
|---|---|---|
| Khi-deux | 63,543 | 0 |
| Ratio de vraisemblance du khi-deux | 67,286 | 0 |
| Khi-deux de Mantel-Haenszel | 48,486 | 0 |
| Corrélation Polychorique | -0,479 | 0 |
| Coefficient de Phi | 0,477 | |
| V de Cramer | 0,276 |