Statistiques au LAPEC – Analyse de données binaires

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)
# pour faire des graphes multiples
library(cowplot)
# pour afficher joliment les tableaux dans un document rmarkdown
library(knitr)

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

On importe ensuite nos données pour travailler dessus :

d <- read.csv("etude_fictive_diabetes.csv")
head(d)

##   X id   group sex age height weight_t0 diabetes steps_t0 t_compliance
## 1 1  1 control   F  72   1.61     80.25       II     2642           32
## 2 2  2 control   M  71   1.81     79.66        I     1675           48
## 3 3  3 control   F  32   1.72     67.25       II     2698            8
## 4 4  4 control   M  76   1.74     73.59        I     2519            2
## 5 5  5 control   F  34   1.53     83.98        I     2764           12
## 6 6  6 control   F  70   1.71     69.14       II     2606           24
##     age_group      bmi bmi_group diff_weight diff_steps is_depressed
## 1  plus de 60 30.95945     obese  -3.6168091  -26.66633            1
## 2  plus de 60 24.31550    normal  -2.6491665  656.62903            1
## 3 moins de 40 22.73188    normal  -1.5288284   17.31293            0
## 4  plus de 60 24.30638    normal  -3.0515988 -589.42996            1
## 5 moins de 40 35.87509     obese   0.8933788  761.91017            1
## 6  plus de 60 23.64488    normal  -0.3074968  216.47307            1
##   feels_sleepy n_meals_glycemia_out_of_target weight_tf steps_tf
## 1      morning                             53  76.63319 2615.334
## 2      morning                             19  77.01083 2331.629
## 3      morning                              4  65.72117 2715.313
## 4      morning                             17  70.53840 1929.570
## 5    afternoon                            188  84.87338 3525.910
## 6       always                              9  68.83250 2822.473

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

Comparaison de proportions

Visuellement

Nous avons deux colonnes binaires qui se prêtent bien à ce premier exercice : sex et is_depressed. 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 en dépression.

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", alpha=0.7) +
    theme_bw() + 
    scale_fill_manual(values=c(bleuclair, rose))

p2 <- d %>% 
    mutate(is_depressed = as.factor(ifelse(is_depressed==1, "patient dépressif", "patient non dépressif"))) %>%
    ggplot(aes(x = group, fill = is_depressed)) +
    geom_bar(position = "dodge", alpha=0.7) +
    theme_bw() + 
    ylab("nombre") +
    xlab("") +
    scale_fill_manual(values=c(bleuclair, rose))

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

# Les deux lignes commentées suivantes permettent d'enregistrer les figures
#ggsave("../../Figures/illu_diabete_is_depressed_barplot.pdf", plot=p2, width=8, height=5)
#ggsave("../../Figures/illu_diabete_sex_barplot.pdf", plot=p1, width=8, height=5)

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	250
M	250

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

## 
##  Exact binomial test
## 
## data:  t
## number of successes = 250, number of trials = 500, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.4552856 0.5447144
## sample estimates:
## probability of success 
##                    0.5

Que concluons-nous ?

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 et largement connu !

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 = 0, df = 1, p-value = 1

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 without continuity correction
## 
## data:  t, null probability 0.5
## X-squared = 0, df = 1, p-value = 1
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4563413 0.5436587
## sample estimates:
##   p 
## 0.5

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 = 0, df = 1, p-value = 1
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4563413 0.5436587
## sample estimates:
##   p 
## 0.5

Que concluons-nous ?

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	127	123
M	123	127

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 
## 250 250

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.

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.

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.7885
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.7393262 1.5372832
## sample estimates:
## odds ratio 
##   1.065954

Que concluons-nous ?

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.25      0.25
## M    0.25      0.25

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

t2 / sum(t2)

##    
##     control treatment
##   F   0.254     0.246
##   M   0.246     0.254

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 exemple, 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.072, df = 1, p-value = 0.7884

Là encore, une fonction alternative existe :

prop.test(t2)

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  t2
## X-squared = 0.072, df = 1, p-value = 0.7884
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.07564103  0.10764103
## sample estimates:
## prop 1 prop 2 
##  0.508  0.492

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

Que concluons-nous ?

Vérifier que les tests de \(\chi^2\) et de Fisher semblent cohérents.

Sur l’outcome de dépression

Il ne reste qu’à ré-appliquer les même fonctions que précedemment avec une variable différente :

t3 <- table(d$is_depressed, d$group)
fisher.test(t3)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  t3
## p-value = 0.0001962
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.3413613 0.7260068
## sample estimates:
## odds ratio 
##  0.4986961

chisq.test(t3)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  t3
## X-squared = 13.801, df = 1, p-value = 0.0002032

Que concluons-nous ?

Pour 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 :

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

Var1	Freq
afternoon	174
always	190
morning	136

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

## 
##  Chi-squared test for given probabilities
## 
## data:  t4
## X-squared = 9.232, df = 2, p-value = 0.009892

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(t4)

## 
##  Chi-squared test for given probabilities
## 
## data:  t4
## X-squared = 9.232, df = 2, p-value = 0.009892

Que concluons-nous ?

Entre groupes : version Fisher

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

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

	control	treatment
afternoon	103	71
always	74	116
morning	73	63

fisher.test(t5)

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

Que concluons-nous ?

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

chisq.test(t5)

## 
##  Pearson's Chi-squared test
## 
## data:  t5
## X-squared = 15.905, df = 2, p-value = 0.0003519

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

Que concluons-nous ?

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. Nous n’étions également pas en mesure de regarder l’impact d’une variable continue sur un outcome binaire.

Le cadre approprié pour faire ça est celui de la régression logistique.

Cadre formel

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

On suppose que la loi de \(Y\) sachant \(X\) est une loi de Bernoulli, de paramètre \(p(X)\).

On souhaite modéliser \(p(X)\) en prenant en compte n’importe quelle somme de variables explicatives, mais toujours de façon à avoir \(p(X) \in (0,1)\).

Onn 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 is_depressed à partir du poids initial weight_t0. Première chose à faire : un graphe !

p <- d %>% ggplot(aes(x=weight_t0, y=as.numeric(is_depressed))) + 
    geom_point() +
    theme_bw() +
    xlab("Poids initial (kg)") +
    ylab("Indicatrice de dépression")
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("is_depressed ~ weight_t0", family=binomial, data=d)
summary(res)

## 
## Call:
## glm(formula = "is_depressed ~ weight_t0", family = binomial, 
##     data = d)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -14.82559    1.32129  -11.22   <2e-16 ***
## weight_t0     0.20773    0.01873   11.09   <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: 680.29  on 499  degrees of freedom
## Residual deviance: 443.24  on 498  degrees of freedom
## AIC: 447.24
## 
## 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’être dépressif à 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.230876

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.262221   0.265448   0.988    0.323
## age         -0.005438   0.005182  -1.049    0.294
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 693.15  on 499  degrees of freedom
## Residual deviance: 692.04  on 498  degrees of freedom
## AIC: 696.04
## 
## Number of Fisher Scoring iterations: 3

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

Exemple avec une covariable discrète

Tout l’intérêt d’une régression logistique est de pouvoir utiliser tout ce qu’on souhaite comme variable explicative. Qu’obtient-on si on cherche à expliquer is_depressed à partir de group ?

Un dessin tout d’abord :

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

Et la régression logistique, ensuite :

res <- glm("is_depressed ~ group", family=binomial, data=d)
summary(res)

## 
## Call:
## glm(formula = "is_depressed ~ group", family = binomial, data = d)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      0.0160     0.1265   0.126 0.899344    
## grouptreatment  -0.6972     0.1842  -3.785 0.000154 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 680.29  on 499  degrees of freedom
## Residual deviance: 665.73  on 498  degrees of freedom
## AIC: 669.73
## 
## 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 is_depressed à partir du group et de weight_t0 ?

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

## 
## Call:
## glm(formula = "is_depressed ~ group + weight_t0", family = binomial, 
##     data = d)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -14.83797    1.35251  -10.97  < 2e-16 ***
## grouptreatment  -0.93792    0.24555   -3.82 0.000134 ***
## weight_t0        0.21453    0.01945   11.03  < 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: 680.29  on 499  degrees of freedom
## Residual deviance: 428.01  on 497  degrees of freedom
## AIC: 434.01
## 
## Number of Fisher Scoring iterations: 5

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 par rapport à 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. la régression logistique permet d’étendre le framework de la régression linéaire à des données binaires.