Session 3 : Statistiques pratiques sur données discrètes

Pour les utilisateurs de R, on commence par charger un certain nombre de modules intéressants.

# le package moderne standard pour manipuler ses données
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# pour faire des graphes multiples
library(cowplot)

## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp

# pour afficher joliment les tableaux dans un document rmarkdown
library(knitr)

# quelques couleurs manuelles
bleufonce <- "#3d5468"
bleuclair <- "#5b7c98"
rose <- "#ff5555"

Préparation des données

Simulation

Pour se focaliser sur les stats et uniquement les stats, on travaille sur un jeu de données idéalisé. Et puisqu’on bosse dans les essais cliniques, on va simuler le RCT (Randomized Controlled Trial) ayant les caractéristiques suivantes :

• l’étude possède deux bras: control et treatment.
• on enregistre des variables démographiques à l’inclusion, type sex, age, age_group, height, weight_t0.
• des variables d’outcomes liés à la prise de poids : diff_weight, has_lost_weight.
• une variable d’outcome à trois modalités : feels_sleepy.
• une variable de temps de suivi du traitement : t_compliance.
• une variable d’outcome de nombre de chutes pendant le suivi : n_falls.

Pour ceux que ça intéresse, voici la loi dans laquelle les données sont simulées :

# data that do not depend on the group
d <- data.frame(id = seq(1,500),
                group = as.factor(c(rep("control",250), rep("treatment",250))),
                sex = ifelse(rbinom(n=500, size=1, prob=0.4) == 0, "F", "M"),
                age = round(runif(n=500, min=19, max=77)),
                height = round(rnorm(n=500, mean=1.7, sd=0.1), digits=2),
                weight_t0 = round(rnorm(n=500, mean=70, sd=10)),
                t_compliance = round(rexp(n=500, rate=0.05))  )

# possibilities for the "feels_sleepy" variable
modalites <- c("morning", "afternoon", "always")

# data that do depend on the group
d <- d %>% 
    mutate(diff_weight = c(rnorm(n=250, mean=0.5*(d$weight_t0[1:250]-70), sd=2.5), 
                           rnorm(n=250, mean=0.5*(d$weight_t0[251:500]-70)+0.25*d$t_compliance[251:500], sd=2.5)),
           feels_sleepy = c( sample(modalites, 250, prob = c(0.3, 0.4, 0.3), replace = TRUE),
                                    sample(modalites, 250, prob = c(0.2, 0.3, 0.5), replace = TRUE) ),
           n_falls = c(rpois(n=250, lambda=6), rpois(n=250, lambda=5.5)) )

d <- d %>% mutate(weight_tf = weight_t0 - diff_weight,
                  age_group = as.factor(ifelse(age < 40, 1, ifelse(age < 60, 2, 3))),
                  has_lost_weight = (diff_weight > 0) )

Ecriture/Lecture du dataset

Pour être sûr de travailler sur le fichier tel qu’il sera lu par tout le monde, on enregistre le dataset simulé puis on écrase le dataset avec l’import.

write.csv(d, "etude_fictive.csv")
d <- read.csv("etude_fictive.csv")
head(d)

##   X id   group sex age height weight_t0 t_compliance diff_weight feels_sleepy
## 1 1  1 control   M  38   1.83        62           35   -5.202373      morning
## 2 2  2 control   F  73   1.69        66           13   -2.903321       always
## 3 3  3 control   F  46   1.54        79           21    2.600586    afternoon
## 4 4  4 control   F  36   1.71        72           55    1.120590    afternoon
## 5 5  5 control   F  30   1.66        64            3   -2.330265    afternoon
## 6 6  6 control   F  75   1.75        63           12   -2.758114      morning
##   n_falls weight_tf age_group has_lost_weight
## 1       8  67.20237         1           FALSE
## 2       6  68.90332         3           FALSE
## 3       9  76.39941         2            TRUE
## 4       6  70.87941         1            TRUE
## 5       8  66.33026         1           FALSE
## 6      11  65.75811         3           FALSE

Nous sommes prêts à travailler sur ce jeu de données, et, pour ce qui nous intéresse aujourd’hui, à analyser les données discrètes qu’il renferme.

Comparer des proportions

Visuellement

Nous avons deux colonnes binaires qui se prêtent bien à ce premier exercice : sex et has_lost_weight. On souhaiterait savoir si

1. la proportion de femmes dans l’ensemble de l’étude est différente de 50%,
2. les deux groupes contiennent la même proportion de femmes,
3. les deux groupes contiennent la même proportion de patients ayant perdu du poids.

Comme toujours en stats, il est intéressant de représenter les données pour se faire une idée intuitive !

p1 <- d %>% 
    ggplot(aes(x = group, fill = sex)) +
    geom_bar(position = "dodge") +
    theme_bw() + 
    scale_fill_manual(values=c(bleuclair, rose))

p2 <- d %>% 
    ggplot(aes(x = group, fill = has_lost_weight)) +
    geom_bar(position = "dodge") +
    theme_bw() + 
    scale_fill_manual(values=c(bleuclair, rose))

plot_grid(p1, p2, nrow=1, ncol=2, label_size = 12)

A présent qu’on voit un peu mieux à quoi ressemblent nos données, on aimerait savoir si les différences de hauteur de barres sont dues à des fluctuations d’échantillonnage ou non.

Avec un attendu théorique binomial

On souhaite tout d’abord comparer notre proportion observée avec une certaine proportion fixée, et répondre à la question suivante :

La proportion de femmes dans l’ensemble de l’étude est-elle différente de 50% ?

La première solution – appelons la “solution exacte” – consiste à s’appuyer sur la distribution attendue sous H0 du nombre de succès dans notre suite d’essais, avec une probabilité de succès fixée.

Attendu théorique

Le nombre de femmes (succès) dans un échantillon de 500 individus (essais) avec une proportion de 50% de femmes (probabilité de succès) est distribué selon une loi binomiale de paramètre \(\mathcal{B}(n=500, p=0.5)\).

On va tout simplement comparer le nombre de femmes observées dans notre échantillon à cette distribution théorique sous H0, de façon à rejeter H0 si la valeur observée est trop extrême, ou à conserver H0 dans le cas inverse.

Et avec R ?

La syntaxe est la suivante :

t <- table(d$sex)
kable(t)

Var1	Freq
F	305
M	195

res <- binom.test(t, p=0.5, conf.level=0.95)
res

## 
##  Exact binomial test
## 
## data:  t
## number of successes = 305, number of trials = 500, p-value = 9.909e-07
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.5656996 0.6529917
## sample estimates:
## probability of success 
##                   0.61

Très petite p-value donc : on peut rejeter H0.

Compréhension visuelle du test

Pour que ça semble moins tombé du ciel, il nous est possible de représenter facilement cette distribution théorique, sous H0 : profitons-en !

# La p-value est la probabilité d'avoir une valeur de statistique de test supérieure à la valeur observée, sous H0
pval <- res$p.value
nF <- t[1]

# On représente la densité de la loi binomiale qui nous intéresse
x <- seq(150, 350, 1)
data_binom <- data.frame(x = x,
                         density = dbinom(x, size=500, prob=0.5),
                         cdf = pbinom(x, size=500, prob=0.5))

plot_density <- data_binom %>% ggplot(aes(x = x, y = density)) +
    geom_line() +
    geom_vline(xintercept = nF, lty="dashed") +
    geom_area(data=data_binom[x>nF,], fill="red") +
    ylab("Densité sous H0") +
    xlab("nombre de femmes") +
    theme_bw()

plot_tail <- data_binom %>% ggplot() +
    geom_line(aes(x = x, y = 1-cdf)) +
    geom_hline(yintercept = pval, lty="dashed", color="red") +
    geom_vline(xintercept = nF, lty="dashed") +
    ylab("Queue de distribution") +
    xlab("nombre de femmes") +
    theme_bw()

plot_grid(plot_density, plot_tail, ncol = 1, nrow = 2)

Sur le premier graphe, la valeur observée est extrême par rapport à la densité attendue sous H0. Sur le second graphe, la probabilité de voir une valeur aussi extrême ou plus extrême qu’observée se lit directement en ordonnée de la courbe pointillée rouge : cette valeur est notre p-value !

Avec un attendu théorique : version \(\chi^2\)

La seconde solution pour répondre à la même question consiste à faire un test du \(\chi^2\), qui n’est pas un test exact mais repose sur un résultat théorique asymptotique.

• Inconvénient : non-exact et repose sur un théorème compliqué,
• Avantage : très général, on en reparlera !

Ce que dit la théorie

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}\).

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.

Utilisation de la théorie pour construire un 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 ?

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

pobs <- t
pth <- c(0.5, 0.5)
res <- chisq.test(x=pobs, p=pth)
res

## 
##  Chi-squared test for given probabilities
## 
## data:  pobs
## X-squared = 24.2, df = 1, p-value = 8.683e-07

Note pour les utilisateurs de R : une fonction alternative existe pour faire la même chose :

prop.test(t, p=0.5, conf.level=0.95)

## 
##  1-sample proportions test with continuity correction
## 
## data:  t, null probability 0.5
## X-squared = 23.762, df = 1, p-value = 1.09e-06
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.5655522 0.6527314
## sample estimates:
##    p 
## 0.61

On remarque une très légère différence de p-value, qui est due à un paramètre par défaut de “correction de continuité” dans cette fonction. On retrouvera un résultat tout à fait similaire à un bon vieux \(\chi^2\) avec la commande suivante :

prop.test(t, p=0.5, conf.level=0.95, correct=F)

## 
##  1-sample proportions test without continuity correction
## 
## data:  t, null probability 0.5
## X-squared = 24.2, df = 1, p-value = 8.683e-07
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.5665640 0.6517587
## sample estimates:
##    p 
## 0.61

Compréhension visuelle du test

Dobs <- res$statistic
# 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,
                         density = dchisq(x, df=res$parameter),
                         cdf = pchisq(x, df=res$parameter))

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

plot_tail <- data_chisq %>% ggplot() +
    geom_line(aes(x = x, y = 1-cdf)) +
    geom_hline(yintercept = res$p.value, lty="dashed", color="red") +
    geom_vline(xintercept = Dobs, lty="dashed") +
    ylab("Queue de distribution sous H0") +
    xlab("statistique de test") +
    theme_bw()

plot_grid(plot_density, plot_tail, ncol = 1, nrow = 2)

Entre groupes : version Fisher

A présent, on souhaite comparer nos proportions de femmes dans les deux bras du RCT, on construit donc la table de contingence suivante :

t2 <- table(d$sex, d$group)
kable(t2)

	control	treatment
F	157	148
M	93	102

On considère ensuite qu’on connaît très exactement le nombre d’hommes et de femmes, ainsi que le nombre de patients dans les groupes traitement et contrôle :

rowSums(t2)

##   F   M 
## 305 195

colSums(t2)

##   control treatment 
##       250       250

La distribution des tables de contingence qu’on attend sous H0 est celle qui correspond à l’indépendance des deux marginales, c’est à dire celle où on obtient le contenu de la première colonne en piochant sans remise 250 fois dans le pool connu de femmes et d’hommes.

En R, on applique ce test avec la commande suivante :

fisher.test(t2, conf.level=0.95)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  t2
## p-value = 0.4633
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.7990228 1.6945366
## sample estimates:
## odds ratio 
##   1.163114

Note : ça ne correspond pas totalement à l’expérience réelle, dans laquelle techniquement seul le nombre de patients dans les deux groupes est connu. Il existe des alternatives à ça, où on ne considère qu’une seule des marginales comme étant fixée. Mais comme en pratique personne n’utilise ces alternatives, on n’en parlera pas plus.

Entre groupes : version \(\chi^2\)

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 est la loi produit des deux marginales),
• H1 : le choix du groupe dépend du sexe.

L’idée intuitive est similaire à ce qui précède : on considère encore comme hypothèse nulle l’indépendance des marginales.

La distribution théorique attendue est donc :

p_sex <- rowSums(t2)/sum(t2)
p_group <- colSums(t2)/sum(t2)
pth <- p_sex %o% p_group
pth

##   control treatment
## F   0.305     0.305
## M   0.195     0.195

Que l’on compare à la distribution observée :

t2 / sum(t2)

##    
##     control treatment
##   F   0.314     0.296
##   M   0.186     0.204

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}\).

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\)).

Utilisation de la théorie pour construire un 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 ?

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 :

res <- chisq.test(t2)
res

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  t2
## X-squared = 0.53804, df = 1, p-value = 0.4632

Là encore, une fonction alternative existe :

prop.test(t2)

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  t2
## X-squared = 0.53804, df = 1, p-value = 0.4632
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.05615244  0.13181448
## sample estimates:
##    prop 1    prop 2 
## 0.5147541 0.4769231

Cette fois-ci, le comportement des deux fonctions est le même par défaut, avec cette “correction de continuité” dans les deux cas.

Notons enfin la très grande proximité des résultats de chisq.test et fisher.test : ouf !

Comparer des données catégorielles

La bonne nouvelle, c’est que tout est très similaire à ce qu’on vient de voir.

D’ailleurs, tout ce que je pouvais écrire de façon générale précédemment a été écrit de façon générale, avec plus de deux catégories notamment.

Il ne reste donc qu’à s’entraîner ! On s’intéresse à présent aux données d’outcome feels_sleepy.

d %>% ggplot(aes(x = group, fill = feels_sleepy)) +
    geom_bar(position = "dodge") +
    theme_bw() + 
    scale_fill_manual(values=c(bleuclair, bleufonce, rose))

Avec une loi donnée

On peut utiliser un test du \(\chi^2\), par exemple pour tester la différence de la distribution observée par rapport à une distribution homogène sur les trois catégories :

t3 <- table(d$feels_sleepy)
kable(t3)

Var1	Freq
afternoon	169
always	210
morning	121

chisq.test(t3, p=c(1/3, 1/3, 1/3))

## 
##  Chi-squared test for given probabilities
## 
## data:  t3
## X-squared = 23.812, df = 2, p-value = 6.75e-06

Pour les utilisateurs de R, notons que, par défaut, le comportement de la fonction chisq.test à qui on ne donne pas d’argument p est justement de comparer à une distribution uniforme :

chisq.test(t3)

## 
##  Chi-squared test for given probabilities
## 
## data:  t3
## X-squared = 23.812, df = 2, p-value = 6.75e-06

Entre groupes : version Fisher

Le même principe que dans le cas 2x2 est appliqué sur notre table de contingence 3x2.

t4 <- table(d$feels_sleepy, d$group)
kable(t4)

	control	treatment
afternoon	96	73
always	82	128
morning	72	49

fisher.test(t4)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  t4
## p-value = 0.0001462
## alternative hypothesis: two.sided

Cette fois, la p-value est petite ! On peut conclure à une association entre l’outcome feels_sleepy et le bras de traitement.

Entre groupes : version \(\chi^2\)

chisq.test(t4)

## 
##  Pearson's Chi-squared test
## 
## data:  t4
## X-squared = 17.578, df = 2, p-value = 0.0001524

Encore un résultat bien similaire à celui du test de Fisher : ouf !

Intro GLM : la régression logistique

Jusqu’à présent, nous avons cherché à rejeter une hypothèse nulle d’homogénéité d’un certain nombre de distributions discrètes, définies selon la valeur d’une covariable. Mais nous n’avons pas quantifié l’impact de cette covariable sur notre distribution d’intérêt.

Le cadre approprié pour faire ça est celui de la régression logistique, qui se trouve être un exemple de GLM : Generalized Linear Model.

Cadre formel

Dans le GLM, on s’intéresse à une variable de réponse \(Y\), qu’on souhaite expliquer avec des variables \(X_1, X_2, ..., X_n\) discrètes ou continues.

Il se trouve simplement que la réponse \(Y\) ne réagit pas forcément linéairement aux modifications de \(X\). On introduit donc ce qu’on appelle une fonction de lien \(g\), et on suppose que la variable de réponse est distribuée suivant une loi dont l’espérance est une transformation d’une combinaison linéaire de nos variables explicatives \((X_i)\).

\[ \mathbb{E}( Y | X ) = g^{-1}( a_0 + a_1 x_1 + a_2 x_2 + ... + a_n x_n ) \]

Dans le cas d’une variable de réponse \(Y\) binaire, la loi de \(Y\) sachant \(X\) est une loi de Bernoulli, de paramètre \(p(X)\). Ce paramètre \(p(X)\) se trouve être la probabilité d’obtenir un succès (codé par un 1), et est donc également l’espérance de \(Y | X\).

Que prendre, alors, comme fonction de lien \(g\) ? Il nous faut une fonction \(g\) de \([0,1]\) dans \(\mathbb{R}\). Pour d’autres raisons qu’on ne détaillera pas ici, on choisit la fonction suivante :

\[ p(X) = \frac{1}{1 + e^{-(a_0 + a_1 x_1 + a_2 x_2 + ... + a_n x_n)}} \\ \ln \left( \frac{p}{1-p} \right) = a_0 + a_1 x_1 + a_2 x_2 + ... + a_n x_n \]

Sur la première ligne, la fonction de droite \(u \mapsto \frac{1}{1+e^{-u}}\) est appelée fonction logistique. Sur la seconde ligne, la fonction de gauche \(u \mapsto \ln \frac{u}{1-u}\) est la réciproque de la fonction logistique, appelée fonction logit.

Exemple avec une covariable continue

On commence avec une covariable continue, de façon à se faire une idée visuelle de la forme de la fonction logistique.

On souhaiterait dans un premier temps expliquer has_lost_weight à partir du poids initial weight_t0. Première chose à faire : un graphe !

p <- d %>% ggplot(aes(x=weight_t0, y=as.numeric(has_lost_weight))) + 
    geom_point() +
    theme_bw() +
    xlab("Poids initial (kg)") +
    ylab("Indicatrice de perte de poids")
p

L’illustration suivante superpose à ce premier graphique le fit d’un modèle logistique, en comparaison avec ce qui pourrait être une mauvaise utilisation d’un simple modèle linéaire.

p + geom_smooth(formula = "y~x",
                method = "glm",
                method.args = list(family = binomial),
                se = FALSE,
                color = bleuclair) +
    geom_smooth(formula = "y~x", method = "lm", color = rose, se = F) 

L’avantage de la régression logistique saute aux yeux : on prédit bien une variable entre 0 et 1, et la valeur de la courbe bleue peut directement s’interpréter comme la probabilité, pour un patient ayant le poids en abscisse, d’avoir perdu du poids à l’issue de l’essai.

On peut récupérer les valeurs de notre fit avec un appel à :

res <- glm("has_lost_weight ~ weight_t0", family=binomial, data=d)
summary(res)

## 
## Call:
## glm(formula = "has_lost_weight ~ weight_t0", family = binomial, 
##     data = d)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -15.12487    1.39353  -10.85   <2e-16 ***
## weight_t0     0.22973    0.02058   11.16   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 658.96  on 499  degrees of freedom
## Residual deviance: 407.41  on 498  degrees of freedom
## AIC: 411.41
## 
## Number of Fisher Scoring iterations: 5

Les coefficients renvoyés correspondent à ceux de la combinaison linéaire de variables explicatives, c’est à dire aux \((a_i)\) des formules pré-citées. Dans notre cas, le coefficient en face de weight_t0 est positif, donc un patient plus lourd initialement est plus susceptible d’avoir perdu du poids à l’issue de l’essai.

Pour se donner une idée quantitative de l’impact d’un kg supplémentaire au début de l’essai, on peut calculer rapidement l’odds ratio associé à ce kg supplémentaire en prenant l’exponentielle de notre paramètre. En effet,

\[ \begin{align*} \frac{ \frac{ p(X_1 = x+1) }{ 1 - p(X_1 = x + 1) } }{ \frac{ p(X_1=x) }{ 1 - p(X_1 = x) } } &= \exp \ln \frac{ \frac{ p(X_1 = x+1) }{ 1 - p(X_1 = x + 1) } }{ \frac{ p(X_1=x) }{ 1 - p(X_1 = x) } } \\ &= \exp \left( \ln \frac{ p(X_1 = x+1) }{ 1 - p(X_1 = x + 1) } - \ln \frac{ p(X_1=x) }{ 1 - p(X_1 = x) } \right) \\ &= \exp \left( a_0 + a_1(x + 1) - a_0 - a_1x \right) \\ &= \exp a_1 \end{align*} \]

Dans notre cas, ça donne :

exp(coef(res)[2])

## weight_t0 
##  1.258264

L’output de la régression logistique nous donne aussi la p-value associée au test de non-nullité de nos coefficients de régression. Ici, les p-values sont petites, on peut donc rejeter l’hypothèse nulle de nullité de nos deux coefficients.

Sur un exemple différent, en cherchant à expliquer sex à partir de age:

dbis <- d %>% mutate(sex = ifelse(sex == "F", 1, 0))
res <- glm("sex ~ age", family=binomial, data=dbis)
summary(res)

## 
## Call:
## glm(formula = "sex ~ age", family = binomial, data = dbis)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.823832   0.276540   2.979  0.00289 **
## age         -0.007917   0.005456  -1.451  0.14679   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 668.75  on 499  degrees of freedom
## Residual deviance: 666.64  on 498  degrees of freedom
## AIC: 670.64
## 
## Number of Fisher Scoring iterations: 4

Cette fois les p-values sont grandes, et on ne peut pas conclure à une différence significative de nos paramètres avec zéro.

Exemple avec une covariable discrète

Tout l’intérêt du GLM est de pouvoir utiliser tout ce qu’on souhaite comme variable explicative. Qu’obtient-on si on cherche à expliquer has_lost_weight à partir de group ?

Un dessin tout d’abord :

d %>% 
    ggplot(aes(x = group, fill = has_lost_weight)) +
    geom_bar(position = "dodge") +
    theme_bw() + 
    scale_fill_manual(values=c(bleuclair, rose))

Et la régression logistique, ensuite :

res <- glm("has_lost_weight ~ group", family=binomial, data=dbis)
summary(res)

## 
## Call:
## glm(formula = "has_lost_weight ~ group", family = binomial, data = dbis)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      0.0160     0.1265   0.126    0.899    
## grouptreatment   1.1149     0.1941   5.743  9.3e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 658.96  on 499  degrees of freedom
## Residual deviance: 624.38  on 498  degrees of freedom
## AIC: 628.38
## 
## Number of Fisher Scoring iterations: 4

La p-value petite nous pousse à conclure à un effet significativement différent de 0 du bras de traitement sur l’outcome.

Exemple multivarié

Comme dans le modèle linéaire, on peut proposer des effets additifs de différentes covariables. Par exemple, que se passe-t-il si on souhaite expliquer has_lost_weight à partir du group et de weight_t0 ?

res <- glm("has_lost_weight ~ group + weight_t0", family=binomial, data=dbis)
summary(res)

## 
## Call:
## glm(formula = "has_lost_weight ~ group + weight_t0", family = binomial, 
##     data = dbis)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -20.49007    1.90096 -10.779  < 2e-16 ***
## grouptreatment   2.48993    0.33260   7.486 7.08e-14 ***
## weight_t0        0.29197    0.02685  10.872  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 658.96  on 499  degrees of freedom
## Residual deviance: 332.59  on 497  degrees of freedom
## AIC: 338.59
## 
## Number of Fisher Scoring iterations: 6

Le modèle avec une influence des deux variables est bien supporté par nos données ! On obtient donc une estimation de l’effet du bras de traitement tout en ayant pris en compte l’influence d’une autre variable explicative.

Conclusion

Nous avons vu aujourd’hui un certain nombre de méthodes permettant d’analyser des variables discrètes. Les points principaux à retenir sont les suivants :

1. pour des données binaires, un test exact de différence à une proportion fixée existe, mais un test du \(\chi^2\) fait aussi l’affaire.
2. pour des données bivariées discrètes, un test de Fisher d’indépendance entre les deux variables existe, mais un test du \(\chi^2\) fait aussi l’affaire.
3. certains préfèreront toujours effectuer un test de Fisher, considéré comme “exact” sous ses hypothèses. D’autres préfèrent le test du \(\chi^2\), plus général et “robuste”, mais reposant sur un résultat théorique asymptotique.
4. le GLM permet d’étendre le framework de la régression linéaire à des données non continues.
5. la régression logistique est un exemple de GLM, qui permet d’expliquer des données binaires par des données discrètes ou continues.

Notez pour finir qu’il existe d’autres types de régressions bien adaptées à des données discrètes différentes. Parmi les plus célèbres, vous rencontrerez peut-être un jour la régression softmax, pour des données catégorielles, ou la régression de Poisson, pour des données de comptage d’événements. Elles fonctionnent toutes sur le même principe, avec une fonction de lien différente.