Le test du chi2 : un couteau suisse pour statisticien

Commençons par charger quelques modules utiles.

# comme toujours, nous utiliserons dans ce tutoriel la syntaxe du "tidyverse"
library(tidyverse)
# nous utiliserons "cowplot" pour faire des figures multiples facilement
library(cowplot)

Le test du \(\chi^2\) (lire “ki deux”) est central dans la culture statistique du fait de ses très nombreuses utilisations. Il peut en effet servir :

1. à tester l’ajustement d’une distribution discrète à une loi donnée
2. à tester l’ajustement d’une distribution discrète à une famille de lois donnée (paramétrée)
3. à tester l’indépendance de deux distributions appariées
4. à tester l’homogénéité de plusieurs échantillons de lois discrètes.

Par ailleurs, bien que ces 4 points d’utilisation concernent normalement des variables prenant des valeurs dans un ensemble discret et fini, des tests du \(\chi^2\) sont aussi régulièrement effectués pour répondre aux même questions sur des variables prenant des valeurs dans un ensemble continu, auquel cas une première étape de discrétisation des données sera effectuée.

L’objectif de ce cours est de comprendre les grands principes théoriques sous-jacents au test du \(\chi^2\), et d’être capable de réaliser des tests avec R dans les 4 cas qui nous intéressent. Commençons par le plus simple : le test d’ajustement à une distribution discrète.

Ajustement à une distribution discrète

Contexte formalisé

On dispose d’un échantillon de \(n\) valeurs observées \((X_1, X_2, ..., X_n)\) qui prennent leurs valeurs dans un ensemble discret de \(d\) valeurs \(\mathcal{X} = \{ x_1, x_2, ..., x_d \}\).

On considère que nos \(n\) observations proviennent de variables aléatoires indépendantes et identiquement distribuées dans une même loi qu’on note \(p = (p_1, p_2, ..., p_d)\), où \(p_i\) est la probabilité que la variable aléatoire prenne la valeur \(x_i\).

Par ailleurs, on dispose d’une loi de référence sur \(\mathcal{X}\) qu’on appelle \(p^\text{ref}\), à laquelle on souhaiterait comparer notre échantillon.
Plus précisément, on souhaite tester :

• H0 : \(p = p^\text{ref}\)
• H1 : \(p \neq p^\text{ref}\).

Contexte pratique

On simule un “toy dataset” très simple inspiré d’essais cliniques, contenant deux variables discrètes uniquement. On maîtrise ainsi la loi génératrice de nos variables :

• le sexe est tiré dans une loi uniforme sur l’ensemble \(\{ F, M \}\),
• le groupe est tiré dans une loi à valeurs dans \(\{ \text{Control, Treatment1, Treatment2} \}\), qui place la même probabilité de 0.3 sur les deux groupes traitement et une masse de probabilité de 0.4 sur le groupe contrôle.

n <- 100
d <- data.frame(sex = sample(size = n, 
                             x = c("F", "M"), 
                             prob = c(0.5, 0.5), 
                             replace = T),
                group = sample(size = n, 
                               x = c("Control", "Treatment1", "Treatment2"), 
                               prob = c(0.4, 0.3, 0.3), 
                               replace = T))

Plusieurs questions peuvent se poser autour de ce dataset :

• est-ce que les sexes sont équilibrés ? autrement dit, est-ce que la loi génératrice des sexes est \(p^\text{ref} = (0.5, 0.5)\) ?
• est-ce que les groupes sont équilibrés ? i.e. est-ce que la loi génératrice est uniforme sur nos trois groupes ?

Ce que dit la théorie

On peut estimer chacun des \(p_i\) par la fréquence empirique \(\hat{p}_i\) de l’observation \(x_i\), c’est à dire : \[ \hat{p}_i = \frac{1}{n} \sum_{k=1}^n \delta_{X_k = x_i} \] où on note \(\delta_{X_k = x_i}\) la variable indicatrice de l’événement \(X_k = x_i\), c’est à dire la variable valant 1 lorsque \(X_k = x_i\) et 0 sinon.

On souhaite ensuite mesurer une certaine idée de “distance” entre la distribution empiriquement observée \(\hat{p}\) et celle de référence \(p^\text{ref}\). Cette distance est appelée “pseudo-distance du \(\chi^2\)”, et elle s’exprime comme : \[ D^2_n( \hat{p}, p^\text{ref} ) = n \sum_{i=1}^d \frac{\left(\hat{p}_i - p^\text{ref}_i \right)^2}{p^\text{ref}_i} \]

Un théorème nous assure alors que

• sous H0, \(D^2_n ( \hat{p}, p^\text{ref} )\) tend en loi vers une distribution du \(\chi^2\) à \(d-1\) degrés de liberté.
• sous H1, \(D^2_n ( \hat{p}, p^\text{ref} )\) tend vers l’infini presque sûrement.

Construction du test

Ce théorème est un théorème asymptotique : on ne peut donc l’appliquer pour effectuer notre test sur de vraies données que pour de grands effectifs \(n\) ! On demande aussi généralement à ce que \(n p^\text{ref}_i\) soit au moins égal à 5 pour tout \(i\).

Si ces conditions d’effectifs suffisants sont respectées, le test fonctionne ainsi :

• la statistique de test \(D^2_n( \hat{p}, p^\text{ref} )\) est calculée,
• on fixe l’erreur de premier ordre \(\alpha\) acceptable,
• on calcule \(q_{1-\alpha}\), le quantile de la loi du \(\chi^2\) à \(d-1\) degrés de liberté de niveau \(1-\alpha\),
• on rejette H0 si \(D^2_n( \hat{p}, p^\text{ref} ) > q_{1-\alpha}\).

Intuitivement, on mesure une pseudo-distance entre notre distribution empirique et notre distribution de référence, et on rejette l’hypothèse nulle d’égalité des deux distributions lorsque cette pseudo-distance devient trop grande.

Et avec R ? version courte

En une ligne très simple, est-ce que le sexe est distribué de façon uniforme ?

# nous affichons d'abord les effectifs des différentes catégories
table(d$sex)

## 
##  F  M 
## 58 42

# puis nous affichons le résultat du test du chi2
chisq.test(table(d$sex), p=c(0.5,0.5))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(d$sex)
## X-squared = 2.56, df = 1, p-value = 0.1096

La lecture du compte-rendu de R présente :

• la valeur de la statistique de test \(D^2_n\) qui est appelée X-squared ici,
• le nombre de degrés de liberté (df pour degrees of freedom) de la loi du \(\chi^2\) que suit notre statistique de test sous H0,
• et, très important, la p-value, qui est la probabilité que la statistique de test soit au moins aussi extrême qu’observée, sous H0.

Si la p-value est inférieure au risque de premier espèce \(\alpha\) qu’on s’est fixé, on s’autorise à rejeter H0.

Ici, la p-value est grande, on conserve donc H0 par défaut.

Et nos groupes alors, sont-ils distribués selon une loi uniforme, associant une probabilité \((\frac{1}{3}, \frac{1}{3}, \frac{1}{3})\) à chacun des trois groupes ? Pour gagner encore du temps, lorsqu’on compare notre distribution à une loi de référence uniforme, il n’est pas nécessaire de spécifier à R la loi de référence : implicitement, elle est uniforme.

table(d$group)

## 
##    Control Treatment1 Treatment2 
##         30         29         41

chisq.test(table(d$group))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(d$group)
## X-squared = 2.66, df = 2, p-value = 0.2645

La p-value est ici faible. Selon les domaines scientifiques, l’erreur de première espèce qu’on s’autorise à commettre peut varier grandement. Dans le domaine des essais cliniques, les valeurs classiquement utilisées peuvent être 0.05, 0.025 ou encore 0.001, selon les conséquences qui sont associées au fait de rejeter H0 alors qu’on était sous H0.

Pour finir, peut-être était-il écrit au protocole que 40% des inclusions concerneraient le groupe contrôle, et que le reste serait réparti de façon uniforme entre les deux groupes traitement. Si on souhaite vérifier l’adéquation avec cette loi de référence (non-uniforme cette fois-ci), vous l’aurez deviné, la syntaxe est tout simplement :

chisq.test(table(d$group), p=c(0.4, 0.3, 0.3))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(d$group)
## X-squared = 6.5667, df = 2, p-value = 0.0375

Cette fois-ci, la p-value est forte, et on conserve H0. Ces trois tests étaient donc bien cohérents avec la façon dont les données ont été simulées.

Et avec R ? version longue

Une bonne façon de se convaincre que tout ça fonctionne bien consiste à calculer manuellement les différentes quantités d’intérêt, et surtout à représenter graphiquement la statistique de test et la loi du \(\chi^2\) que notre statistique de test est censée suivre sous H0.

# on peut calculer et afficher les fréquences empiriques
t <- table(d$group)
hatp <- t/sum(t)
hatp

## 
##    Control Treatment1 Treatment2 
##       0.30       0.29       0.41

# qu'on souhaite comparer aux probabilité de référence
pref <- rep(1/3, 3)

# et les utiliser pour calculer la statistique du Chi2 : 
D <- n * sum((hatp - pref) ^ 2 / pref)
print(paste("La statistique de test du chi2 est", D))

## [1] "La statistique de test du chi2 est 2.66"

# La p-value est la probabilité d'avoir une valeur de statistique de test supérieure à la valeur observée, sous H0
pval <- 1 - pchisq(D, df=2)

# On représente la densité de la loi du chi2 qui nous intéresse
x <- seq(0, 10, 0.01)
data_chisq <- data.frame(x = x,
                         dchisq = dchisq(x, df=2),
                         pchisq = pchisq(x, df=2))

plot_densite <- data_chisq %>% ggplot(aes(x = x, y = dchisq)) +
    geom_line() +
    geom_vline(xintercept = D, lty="dashed") +
    geom_area(data=data_chisq[x>D,], fill="red") +
    ylab("Densité de la statistique de test D sous H0") +
    xlab("d") +
    theme_bw()

plot_proba <- data_chisq %>% ggplot() +
    geom_line(aes(x = x, y = 1-pchisq)) +
    geom_hline(yintercept = pval, lty="dashed", color="red") +
    geom_vline(xintercept = D, lty="dashed") +
    ylab("Probabilité de voir une statistique de test D supérieure à la valeur en abscisse") +
    xlab("d") +
    theme_bw()

plot_grid(plot_densite, plot_proba, ncol = 2, nrow = 1, labels=c("A","B"))

Sur nos deux graphes, la ligne noire pointillée correspond à la valeur de la statistique de test observée.

Sur le graphe de gauche, représentant la densité de la loi du \(\chi^2\) que suit notre statistique de test sous H0, la p-value correspond à l’aire coloriée en rouge : la probabilité que la statistique de test soit au moins aussi grande qu’observée.

Sur le graphe de gauche, représentant la probabilité que la statistique de test soit supérieure à une valeur donnée, la p-value se lit directement à l’intersection entre la ligne pointillée rouge et l’axe des ordonnées.

Puisqu’on dispose entre les mains toute la puissance de calcul d’un ordinateur moderne, on peut aussi se demander à quel point l’hypothèse asymptotique (\(n\) grand) est vérifiée pour \(n = 100\) comme on l’avait ici. On va pour cela simuler \(10^5\) datasets sous H0, et calculer pour chacun la valeur de la statistique de test. On pourra ensuite comparer visuellement la distribution empiriquement observée de nos statistiques de tests avec la distribution asymptotique du \(\chi^2\).

stats <- c()
pref <- rep(1/3, 3)

for (i in 1:1e5){
    # on simule notre dataset sous H0
    groups <- sample(size = n, x = c("Control", "Treatment1", "Treatment2"), prob = pref, replace = T)
    # on calcule les fréquences empiriques
    t <- table(groups)
    hatp <- t/sum(t)
    # puis on calcule la statistique de test qu'on enregistre
    D <- n * sum((hatp - pref) ^ 2 / pref)
    stats <- c(stats, D)
}

# on créé un dataframe avec nos statistiques de tests comme unique colonne
stats <- data.frame(D = stats)
# et on compare l'histogramme empirique à la densité asymptotique
plot_densite_vs_histogram <- ggplot() +
    geom_line(data=data_chisq, aes(x = x, y = dchisq)) +
    geom_histogram(data=stats, aes(x=D, y = ..density..), bins=50, alpha=0.5) +
    ylab("Densité de la statistique de test D sous H0") +
    xlab("d") +
    xlim(0, 12) +
    theme_bw()
plot_densite_vs_histogram

Pas trop mal, visuellement ! L’adéquation devrait être d’autant meilleure que \(n\) est grand dans vos simulations.

Nous savons à présent comparer une distribution empirique discrète à une distribution de référence fixée. Mais que se passe-t-il si on souhaite comparer notre distribution empirique à une famille paramétrée de distributions ?

Ajustement à une famille de loi paramétrée

Contexte formalisé

Comme précédemment, on dispose d’un échantillon \(X\) de taille \(n\), tiré dans une loi \(p\) inconnue prenant des valeurs dans un ensemble fini \(\mathcal{X}\) de \(d\) valeurs.

Cependant, plutôt que de comparer à une distribution de référence fixée \(p^\text{ref}\), on souhaite à présent savoir si notre échantillon a été tirée dans une loi de probabilité appartenant à une certaine famille de lois paramétrées par un ou plusieurs réels. On note cette famille \[\mathcal{P} = \{ p(\theta), \theta \in \Theta \subset \mathbb{R}^k \}.\]

On souhaite tester :

• H0 : \(p \in \mathcal{P}\),
• H1 : \(p \notin \mathcal{P}\).

Contexte pratique

Reprenons notre exemple précédent avec nos trois groupes : \(\mathcal{X} = \{ \text{Control, Treatment1, Treatment2} \}\). On pourrait se poser la question suivante :

• H0 : les deux groupes traitement ont la même probabilité d’être choisis
• H1 : l’un des deux groupes a une probabilité différente d’apparaître.

Une des façons de faire pourrait être de chercher à tester l’adéquation de nos observations avec la famille de lois paramétrées suivante : \[ \{ (\theta, \frac{1-\theta}{2}, \frac{1-\theta}{2}), \theta \in (0,1) \}. \]

La probabilité d’être choisi dans le groupe contrôle est un paramètre libre. Une fois fixé, les probabilités d’être dans chacun des deux groupes traitement sont égales et fixées également.

Ce que dit la théorie

Dans un premier temps, on va estimer le paramètre \(\theta\) par son estimateur de maximum de vraisemblance \(\hat{\theta}\).

On peut alors chercher à mesurer la “pseudo-distance du \(\chi^2\)” entre notre distribution empirique \(\hat{p}\) et une distribution de référence obtenue grâce à notre estimateur de maximum de vraisemblance : \(p(\hat{\theta})\). \[ D^2_n( \hat{p}, p(\hat{\theta}) ) = n \sum_{i=1}^d \frac{\left( \hat{p}_i - p_i(\hat{\theta}) \right)^2}{p_i(\hat{\theta})} \]

Il faut alors s’assurer que \(p\) vérifie un certain nombre de conditions techniques :

• c’est une fonction injective de \(\Theta\) dans l’ensemble des distributions de probabilités sur \(d\) éléments,
• c’est une fonction de classe \(\mathcal{C}^2\),
• aucune de ses composantes ne s’annule sur \(\Theta\),
• ses \(k\) dérivées partielles sont linéairement indépendantes en tout point de \(\Theta\).

Un théorème nous assure alors que, sous ces conditions,

• sous H0, \(D^2_n( \hat{p}, p(\hat{\theta}) )\) tend en loi vers une distribution du \(\chi^2\) à \(d-1-k\) degrés de liberté.
• sous H1, \(D^2_n( \hat{p}, p(\hat{\theta}) )\) tend vers l’infini presque sûrement.

Construction du test

Tout comme précédemment, on retiendra bien que ce théorème est un théorème asymptotique ! On ne pourra donc l’appliquer pour effectuer notre test sur de vraies données que pour de grands effectifs \(n\).

Si ces conditions d’effectifs suffisants sont respectées, le test fonctionne ainsi :

• l’estimateur de maximum de vraisemblance \(\hat{\theta}\) est calculé,
• la statistique de test \(D^2_n( \hat{p}, p(\hat{\theta}) )\) est calculée,
• on fixe l’erreur de premier ordre \(\alpha\) acceptable,
• on calcule \(q_{1-\alpha}\), le quantile de la loi du \(\chi^2\) à \(d-1-k\) degrés de liberté de niveau \(1-\alpha\),
• on rejette H0 si \(D^2_n( \hat{p}, p^\text{ref} ) > q_{1-\alpha}\).

Intuitivement, on mesure une pseudo-distance entre notre distribution empirique et notre distribution de référence, et on rejette l’hypothèse nulle d’égalité des deux distributions lorsque cette pseudo-distance devient trop grande.

Et avec R ? une construction ad-hoc

Etant donnée le nombre immense de possibilités de familles de lois, nous n’avons pas ici une unique fonction qui nous renverrait directement notre réponse sur un plateau d’argent comme précédemment. Mais à présent que nous avons bien compris les différentes étapes, il n’est pas trop difficile d’enchaîner celles-ci avec R.

Dans notre example pratique, nous devons tout d’abord trouver l’estimateur de maximum de vraisemblance de \(\theta\). On écrit donc la vraisemblance : \[ \mathcal{L}(\theta) = \theta^{n_1} \left( \frac{1 - \theta}{2} \right)^{n_2 + n_3} \] où \(n_1, n_2, n_3\) sont respectivement les effectifs observés dans les groupes contrôle, traitement1 et traitement2.

La log-vraisemblance est alors \[ \ln \mathcal{L}(\theta) = n_1 \ln \theta + (n_2 + n_3) \ln \frac{1 - \theta}{2} \]

On cherche à maximiser cette fonction, en cherchant comment annuler sa dérivée. Et on tombe sur un estimateur assez cohérent avec ce qu’on aurait pu imaginer intuitivement : \[ \hat{\theta} = \frac{n_1}{n} \]

# estimation de p chapeau et theta chapeau
t <- table(d$group)
hatp <- t/sum(t)
hat_theta <- hatp[1]

# dont on déduit p(theta chapeau)
p_hat_theta = c(hat_theta, (1-hat_theta)/2, (1-hat_theta)/2)

# qu'on utilise pour calculer la statistique de test
D <- n * sum((hatp - p_hat_theta) ^ 2 / p_hat_theta)

# on calcule et on affiche la p-value associée au test
pval <- 1 - pchisq(D, df=1)
pval

## [1] 0.151494

La p-value est grande : on ne peut pas rejeter H0 dans notre cas. Et pour cause, on avait bien simulé nos données telles que la probabilité d’être dans chacun des deux groupes traitement était égale.

Nous savons à présent comparer une distribution empirique discrète à une loi de référence fixée, ou à une famille de loi paramétrées. Nous allons maintenant nous intéresser à un autre cas de figure fréquent : tester l’indépendance de deux distributions discrètes.

Indépendance de deux distributions discrètes appariées

Contexte formalisé

On dispose cette fois de \(n\) couples d’observations de variables discrètes \((X_i, Y_i)\) mesurées sur chaque individu \(i\). Les \((X_i)\) sont à valeur dans un ensemble \(\mathcal{X}\) à \(d\) éléments, tandis que les \((Y_i)\) sont à valeur dans un ensemble \(\mathcal{Y}\) à \(e\) éléments.

On suppose que nos couples sont indépendants et identiquement distribués selon une loi \(\nu\) sur \(\mathcal{X} \times \mathcal{Y}\). Typiquement, une telle distribution serait paramétrée par \((d e - 1)\) valeurs.

On souhaite à présent savoir si les \(X\) et les \(Y\) prennent leurs valeurs indépendamment l’un de l’autre, c’est à dire si la loi \(\nu\) appartient à la famille des lois produit. On note cette famille : \[ \mathcal{F} = \{ p \otimes q, ~ p \in \mathcal{P}(\mathcal{X}), q \in \mathcal{P}(\mathcal{Y}) \} \] où \(\mathcal{P}(\mathcal{X})\) désigne l’ensemble des distributions sur l’ensemble discret \(\mathcal{X}\).

Typiquement, une distribution de \(\mathcal{F}\) n’est donc plus paramétrée que par \((d + e - 2)\) valeurs : les \((d-1)\) valeurs qui paramétrisent \(p\) et les \((e-1)\) valeurs qui paramétrisent \(q\).

Avec ces notations, on pourra tester :

• H0 : \(\nu \in \mathcal{F}\),
• H1 : \(\nu \notin \mathcal{F}\).

Contexte pratique

On continue de filer le même exemple que précédemment. Pour chaque individu \(i\), on observe le sexe \(X_i\) et le groupe \(Y_i\). On se demande si le choix du groupe est indépendant du sexe.

• H0 : le choix du groupe est indépendant du sexe (i.e. la loi jointe appartient à la famille de lois produit \(\mathcal{F}\)),
• H1 : le choix du groupe dépend du sexe.

Ce que dit la théorie

On peut estimer chacun des \(p_i\) par la fréquence empirique de l’observation \(x_i\), c’est à dire : \[ \hat{p}_i = \frac{1}{n} \sum_{k=1}^n \delta_{X_k = x_i} \] On peut également faire de même pour les fréquences empiriques des \(y_j\), en calculant : \[ \hat{q}_j = \frac{1}{n} \sum_{k=1}^n \delta_{Y_k = y_j} \] Et enfin, faire de même pour les fréquences empiriques des différents couples \((x_i, y_j)\) : \[ \hat{\nu}_{i,j} = \frac{1}{n} \sum_{k=1}^n \delta_{(X_k, Y_k) = (x_i, y_j)} \]

On souhaite ensuite mesurer la “pseudo-distance du \(\chi^2\)” entre la distribution empiriquement observée \(\hat{\nu}\) et celle donnée par le meilleur ajustement à une loi produit \(\hat{p} \otimes \hat{q}\), soit : \[ D^2_n( \hat{\nu}, \hat{p}\otimes\hat{q} ) = n \sum_{i=1}^d \sum_{j=1}^e \frac{\left( \hat{\nu}_{i,j} - \hat{p}_i \hat{q}_j \right)^2}{\hat{p}_i \hat{q}_j} \]

Un théorème nous assure alors que

• sous H0, \(D^2_n( \hat{\nu}, \hat{p}\otimes\hat{q} )\) tend en loi vers une distribution du \(\chi^2\) à \((d-1) (e-1)\) degrés de liberté.
• sous H1, \(D^2_n( \hat{\nu}, \hat{p}\otimes\hat{q} )\) tend vers l’infini presque sûrement.

Pour ceux qui souhaiteraient savoir d’où vient le nombre de degrés de liberté, notez que, comme précédemment, il correspond à la dimension de l’espace dans lequel vit la loi de probabilité \(\nu\) (soit \(de -1\)) moins la dimension de l’espace dans lequel vit la distribution à laquelle on souhaite se ramener (ici \(d + e - 2\)).

Construction du test

Comme toujours avec les tests du \(\chi^2\), notre test repose sur un théorème asymptotique : on ne peut donc l’appliquer pour effectuer notre test sur de vraies données que pour de grands effectifs \(n\) !

Si ces conditions d’effectifs suffisants sont respectées, le test fonctionne ainsi :

• la statistique de test \(D^2_n( \hat{\nu}, \hat{p}\otimes\hat{q} )\) est calculée,
• on fixe l’erreur de premier ordre \(\alpha\) acceptable,
• on calcule \(q_{1-\alpha}\), le quantile de niveau \(1-\alpha\) de la loi du \(\chi^2\) à \((d-1) (e-1)\) degrés de liberté,
• on rejette H0 si \(D^2_n( \hat{\nu}, \hat{p}\otimes\hat{q} ) > q_{1-\alpha}\).

Intuitivement, on mesure une pseudo-distance entre notre distribution empirique et notre meilleure estimation de loi produit, et on rejette l’hypothèse nulle d’égalité des deux distributions lorsque cette pseudo-distance devient trop grande.

Et avec R ? version courte

Dans notre example, on souhaite tester l’indépendance du sexe et du groupe chez nos individus. En une ligne, il nous suffit d’écrire :

chisq.test(table(d$sex, d$group))

## 
##  Pearson's Chi-squared test
## 
## data:  table(d$sex, d$group)
## X-squared = 2.4536, df = 2, p-value = 0.2932

Ou encore, pour un résultat équivalent avec une syntaxe différente,

chisq.test(d$sex, d$group)

## 
##  Pearson's Chi-squared test
## 
## data:  d$sex and d$group
## X-squared = 2.4536, df = 2, p-value = 0.2932

La p-value grande nous pousse ici à conserver H0 par défaut : on ne peut pas rejeter l’hypothèse nulle d’indépendance du sexe et du groupe.

Et avec R ? version longue

Pour se convaincre que tout fonctionne bien, il peut être utile de calculer et afficher certaines quantités.

# la table des fréquences empiriques est toujours intéressante à afficher
t <- table(d$sex, d$group)
t

##    
##     Control Treatment1 Treatment2
##   F      19         19         20
##   M      11         10         21

# on peut calculer les fréquences marginales
tot <- sum(t)
hatp <- rowSums(t)/tot
hatp

##    F    M 
## 0.58 0.42

hatq <- colSums(t)/tot
hatq

##    Control Treatment1 Treatment2 
##       0.30       0.29       0.41

# la loi produit obtenue grâge à hatp et hatq
hatpq <- hatp %o% hatq
hatpq

##   Control Treatment1 Treatment2
## F   0.174     0.1682     0.2378
## M   0.126     0.1218     0.1722

# l'estimation de hatnu
hatnu <- t/tot
hatnu

##    
##     Control Treatment1 Treatment2
##   F    0.19       0.19       0.20
##   M    0.11       0.10       0.21

# et enfin la statistique de test
D <- n * sum((hatnu - hatpq) ^ 2 / hatpq)
D

## [1] 2.45364

# qu'on peut utiliser pour calculer la p-value
pval <- 1 - pchisq(D, df=2)
pval

## [1] 0.2932235

Cette p-value pourrait être reportée comme précédemment sur le graphe de densité de la loi du \(\chi^2\) que suit notre statistique de test sous H0. Il se trouve cependant qu’ici c’est exactement la même loi, et ce n’est donc pas très intéressant à répéter.

Nous savons à présent utiliser notre test du \(\chi^2\) pour tester l’indépendance de deux distributions discrètes appariées. Voyons à présent comment on peut adapter ces résultats pour tester l’homogénéité de plusieurs échantillons discrets.

Homogénéité de plusieurs échantillons discrets

Contexte formalisé

On observe un échantillon de taille \(n_1\) \(X^{(1)} = (X^{(1)}_1, X^{(1)}_2, ..., X^{(1)}_{n_1})\) ainsi qu’un échantillon de taille \(n_2\), \(X^{(2)}\).

On suppose que les \((X_i^{(1)})\) proviennent de variables aléatoires indépendantes et identiquement distribuées selon une loi \(p^{(1)}\) et que les \((X_i^{(2)})\) proviennent de variables aléatoires indépendantes et identiquement distribuées selon une loi \(p^{(2)}\).

On souhaite tester :

• H0 : \(p^{(1)} = p^{(2)}\),
• H0 : \(p^{(1)} \neq p^{(2)}\).

Contexte pratique

Imaginons qu’un deuxième centre de recrutement ouvre, et qu’il recrute également des individus dans trois groupes comme précédemment : contrôle, traitement1 et traitement2.

n2 <- 80
d2 <- data.frame(sex = sample(size = n2, 
                             x = c("F", "M"), 
                             prob = c(0.5, 0.5), 
                             replace = T),
                group = sample(size = n2, 
                               x = c("Control", "Treatment1", "Treatment2"), 
                               prob = c(0.35, 0.35, 0.3), 
                               replace = T))

On peut chercher à tester :

• H0 : la distribution des trois groupes est identique dans le premier centre (dataset d) et dans le second centre (dataset d2).
• H1 : les deux distributions sont différentes.

Ce que dit la théorie

Pour répondre à la question, on va chercher à se ramener au cas précédent. On va former \(n_1 + n_2\) paires d’observations \((X_i, Y_i)\), avec :

• pour \(i \leq n_1\), \((X_i, Y_i) = (X_i^{(1)}, 1)\),
• pour \(n_1 < i \leq n_1+n_2\), \((X_i, Y_i) = (X_i^{(2)}, 2)\).

Notre test peut alors se reformuler en :

• H0 : les variables \(X\) et \(Y\) sont indépendantes, c’est à dire la loi jointe correspond à la loi produit des \(X\) et \(Y\).
• H1 : les deux variables ne sont pas indépendantes.

Et c’est tout. Et ça fonctionnerait encore de la même façon avec plus que deux échantillons. Finalement, tester l’homogénéité de plusieurs échantillons revient simplement à créer une variable correspondant au numéro de l’échantillon, et à tester l’indépendance entre le numéro de l’échantillon et la variable d’intérêt.

Et avec R ?

Nous avons déjà fait tout le boulot précédemment. Finalement, il ne s’agit ici que de manipuler nos données pour nous placer dans le même cas de figure.

# on rajoute à chacun des deux datasets une nouvelle colonne "centre"
d <- d %>% mutate(center = rep("Center1", n))
d2 <- d2 %>% mutate(center = rep("Center2", n2))

# on colle les lignes de d2 à la suite de d1
df <- rbind(d, d2)

# et on effectue notre test d'indépendance entre la variable de centre et la variable de groupe
table(df$group, df$center)

##             
##              Center1 Center2
##   Control         30      35
##   Treatment1      29      24
##   Treatment2      41      21

chisq.test(table(df$group, df$center))

## 
##  Pearson's Chi-squared test
## 
## data:  table(df$group, df$center)
## X-squared = 5.1493, df = 2, p-value = 0.07618

La p-value grande nous incite ici à conserver H0 par défaut. Notez que les simulations étaient pourtant réalisées à partir de loi sous-jacentes légèrement différentes. Trop légèrement différentes, visiblement, pour que notre test ait assez de puissance pour les discriminer. On obtiendra évidemment plus de puissance en augmentant les effectifs, ou en augmentant la différence entre les distributions de probabilité.

Conclusion

S’il y a bien deux choses à retenir à propos du test du \(\chi^2\), les voici :

1. il s’applique pour comparer des distributions discrètes : soit une distribution à une distribution de référence, soit une distribution à une famille de distributions, soit plusieurs distributions empiriques ensemble.
2. il repose sur des résultats asymptotiques : et requiert donc de grands effectifs !