LM et GLM

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)

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

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.
• une variable d’outcome liés à la prise de poids : diff_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.

On importe ces données dans R :

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 continues et binaires qu’il renferme.

Tests sur variables numériques

T test : adéquation à une norme

On observe des différences de poids \((X_i)\) chez nos individus entre le début et la fin de l’expérience. Est-ce que la différence de poids moyenne est différente de zéro ?

On souhaite réaliser le test d’hypothèse suivant :

• H0 : la moyenne des \((X_i)\) vaut zero.
• H1 : la moyenne des \((X_i)\) est différente de zéro.

Conditions d’application :

• les variables aléatoires sont indépendantes et de même loi.
• soit les variables aléatoires sont distribuées normalement, soit l’échantillon est “grand”.

On commence toujours par représenter nos données.

plot_diff_weight <- d %>% ggplot(aes(x=diff_weight)) +
    geom_histogram(alpha=0.7, position="identity", bins=20) +
    theme_bw() +
    ylab("effectif") + 
    xlab("différence de poids à 6 mois")
plot_diff_weight

Ici les conditions d’application semblent respectées.

Un résultat théorique certifie qu’une statistique de test calculée à partir de l’effectif, de la moyenne empirique et de la variance empirique, est distribuée suivant une loi de Student avec un certain nombre de “degrés de liberté”.

Il s’agit donc de calculer la valeur observée de cette statistique de test sur notre échantillon, et de la comparer aux quantiles de la loi théorique.

La fonction pour faire ça en R est la suivante :

t.test(x=d$diff_weight)

## 
##  One Sample t-test
## 
## data:  d$diff_weight
## t = 8.0269, df = 499, p-value = 7.219e-15
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.83606 3.02619
## sample estimates:
## mean of x 
##  2.431125

L’output de la fonction nous rappelle :

1. en première ligne l’échantillon sur lequel on a appliqué le test,
2. en seconde ligne on récupère la valeur observée de la statistique de test, ainsi que le nombre de “degrés de liberté” de la loi théorique.
3. en fin de seconde ligne, super importante : la p-value !!! Probabilité que la valeur de la statistique soit au moins aussi extrême qu’observée.
4. troisième ligne : un rappel de l’hypothèse alternative du test (modifiable avec l’argument alternative).
5. un intervalle de confiance à 95% de la moyenne (bornes modifiables via l’argument conf.level).
6. tout à la fin, l’estimation ponctuelle de la moyenne.

Que concluons nous ???

T test : comparaison de 2 moyennes

On observe des valeurs \((X_{i,j})\) de nos individus \(i\) dans deux groupes \(j=1\) et \(j=2\). Ici les deux groupes correspondent au control vs. treatment, et les observations correspondent à une différence de poids entre le début et la fin de l’expérience.

Est-ce que la différence de poids moyenne est différente dans les deux groupes ?

On souhaite réaliser le test d’hypothèse suivant :

• H0 : la moyenne des \((X_{i,1})\) est identique à celle des \((X_{i,2})\).
• H1 : les deux moyennes sont différentes.

Conditions d’application :

• les variables aléatoires sont indépendantes et de même loi.
• soit les variables aléatoires sont distribuées normalement, soit l’échantillon est “grand”.

On commence toujours par représenter nos données et quantités d’intérêt.

m <- d %>% group_by(group) %>% summarise(mean=mean(diff_weight)) %>% select(mean)

plot_diff_weight_by_group <- d %>% ggplot(aes(x=diff_weight, fill=group)) +
    geom_histogram(alpha=0.5, position="identity", bins=20) +
    theme_bw() +
    ylab("effectif") + 
    xlab("différence de poids à 6 mois") +
    scale_fill_manual(values=c(bleuclair, rose)) +
    geom_vline(xintercept = c(m$mean[1], m$mean[2]), color=c(bleuclair, rose))
plot_diff_weight_by_group

Ici les conditions d’application semblent respectées.

Un résultat théorique certifie qu’une statistique de test calculée à partir de l’effectif, des moyennes empiriques des deux groupes et de la variance empirique des deux groupes, est distribuée suivant une loi de Student avec un certain nombre de “degrés de liberté”.

Il s’agit donc de calculer la valeur observée de cette statistique de test sur notre échantillon, et de la comparer aux quantiles de la loi théorique.

La fonction pour faire ça en R est la suivante :

t.test(x=d$diff_weight[d$group == "control"], y=d$diff_weight[d$group == "treatment"])

## 
##  Welch Two Sample t-test
## 
## data:  d$diff_weight[d$group == "control"] and d$diff_weight[d$group == "treatment"]
## t = -8.9306, df = 457.61, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -6.133784 -3.921211
## sample estimates:
##   mean of x   mean of y 
## -0.08262375  4.94487391

Une syntaxe alternative pratique existe également :

t.test(diff_weight ~ group, data=d)

## 
##  Welch Two Sample t-test
## 
## data:  diff_weight by group
## t = -8.9306, df = 457.61, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group control and group treatment is not equal to 0
## 95 percent confidence interval:
##  -6.133784 -3.921211
## sample estimates:
##   mean in group control mean in group treatment 
##             -0.08262375              4.94487391

Ca se lit de façon très similaire à ce qu’on avait précédemment, les différences étant :

1. L’IC à 95% correspond à celui de la différence des deux moyennes.
2. Les estimations ponctuelles correspondent aux deux moyennes.

Que concluons-nous ???

T test : option “apparié”

Je rajoute une section à propos de cette option d’appariement, qui bien souvent embrouille plus qu’elle n’aide.

On possède des observations appariées, typiquement la même quantité mesurée “avant” \((X_i)\) et “après” \((Y_i)\) intervention sur chaque individu. Par exemple, le poids du patient avant intervention et après intervention. Y a t’il un changement de poids du patient avant et après intervention ?

Cela revient à tester :

• H0 : la moyenne des \((Y_i - X_i)\) vaut zero.
• H1 : la moyenne des \((Y_i - X_i)\) est différente de zéro.

C’est donc totalement similaire au premier test qu’on avait appliqué directement sur la différence de poids.

La façon de faire ça en R est la suivante (et on vérifie que ça donne bien le même résultat) :

t.test(x=d$weight_tf, y=d$weight_t0, paired=TRUE)

## 
##  Paired t-test
## 
## data:  d$weight_tf and d$weight_t0
## t = -8.0269, df = 499, p-value = 7.219e-15
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -3.02619 -1.83606
## sample estimates:
## mean difference 
##       -2.431125

ou

t.test(d$diff_weight)

## 
##  One Sample t-test
## 
## data:  d$diff_weight
## t = 8.0269, df = 499, p-value = 7.219e-15
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.83606 3.02619
## sample estimates:
## mean of x 
##  2.431125

Que concluons-nous ?

ANOVA : comparaison de moyennes

On souhaite plus précisément tester ici la différence d’au moins une moyenne parmi plusieurs groupes.

On dispose cette fois d’observations de la même quantité $(X_{i,j}) chez l’individu \(i\) dans le groupe \(j\), dans un nombre de groupes potentiellement supérieur à 2. Est-ce que la moyenne d’au moins un groupe diffère des autres groupes ?

Par exemple, ici, la différence de poids moyenne change-t-elle dans les trois classes d’âge ? (< 40, 40-60, > 60)

On teste formellement :

• H0 : les moyennes des \((X_{i,1})\), des \((X_{i,2})\) et des \((X_{i,3})\) sont identiques.
• H1 : au moins une des moyennes diffère.

Conditions d’application :

• les variables aléatoires sont indépendantes et suivent une loi normale.
• il y a égalité de variance entre les différents groupes.

Une fois n’est pas coutûme, commençons par représenter les données :

m <- d %>% group_by(age_group) %>% summarise(mean=mean(diff_weight)) %>% select(mean)

plot_diff_weight_by_age_group <- d %>% ggplot(aes(x=diff_weight, fill=as.factor(age_group))) +
    geom_histogram(alpha=0.5, position="identity", bins=20) +
    theme_bw() +
    ylab("effectif") + 
    xlab("différence de poids à 6 mois") +
    scale_fill_manual(values=c(bleufonce, bleuclair, rose)) +
    geom_vline(xintercept = c(m$mean[1], m$mean[2], m$mean[3]), color=c(bleufonce, bleuclair, rose))
plot_diff_weight_by_age_group

Les conditions d’application semblent respectées (au moins graphiquement).

Un résultat théorique certifie cette fois qu’une certaine statistique de test, calculée à partir des effectifs des différents groupes et des variances intra-groupes / inter-groupes, est distribuée suivant une loi de Fisher avec un certain nombre de “degrés de liberté”.

Il s’agit donc de calculer la valeur observée de cette statistique de test sur notre échantillon, et de la comparer aux quantiles de la loi théorique.

La fonction pour faire ça en R est la suivante :

res <- aov(diff_weight ~ age_group, data=d)
summary(res)

##              Df Sum Sq Mean Sq F value Pr(>F)
## age_group     1     12   11.99   0.261   0.61
## Residuals   498  22875   45.93

Les deux points les plus intéressants de l’output sont :

1. la “F value” qui est la valeur de la statistique de test,
2. et la p-value : la probabilité sous H0 d’obtenir une statistique de test au moins aussi extrême qu’observée.

Que concluons-nous ???

Finissons par un petit test, permettant de se convaincre qu’une ANOVA avec deux groupes correspond à un T test.

res <- aov(diff_weight ~ group, data=d)
summary(res)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## group         1   3159  3159.5   79.76 <2e-16 ***
## Residuals   498  19728    39.6                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

t.test(diff_weight ~ group, data=d)

## 
##  Welch Two Sample t-test
## 
## data:  diff_weight by group
## t = -8.9306, df = 457.61, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group control and group treatment is not equal to 0
## 95 percent confidence interval:
##  -6.133784 -3.921211
## sample estimates:
##   mean in group control mean in group treatment 
##             -0.08262375              4.94487391

Qu’est-ce qui nous permet de conclure que les deux tests sont équivalents avec deux groupes ?

Test de corrélation de Pearson

On dispose cette fois d’observations d’une variables \((Y_i)\) d’intérêt, par exemple encore une fois la différence de poids. Et d’une autre variable quantitative \((X_i)\) correspondant à quelque chose qui pourrait être corrélé à nos \((Y_i)\), par exemple ici le temps pendant lequel le patient a strictement observé son traitement (compliance).

Est-ce que le gain de poids dépend de la compliance, dans le groupe traité ?

Le coefficient de corrélation \(\rho\) correspond à la covariance des deux quantités normalisée par les écarts-types des deux quantités. Il est compris entre -1 et 1 et renseigne sur la dispersion des données en 2D. En cas de corrélation linéaire parfaite, il vaut 1 ou -1.

Formellement on peut tester, avec le test de Bravais-Pearson :

• H0 : \(\rho = 0\),
• H1 : \(\rho \neq 0\).

Conditions d’application : les \((Y_i, X_i)\) doivent être appariées et distribuées normalement.

Commençons par représenter les données :

plot_diff_weight_vs_compliance <- d %>% ggplot(aes(y=diff_weight, x=t_compliance, color = group)) +
    geom_point() +
    geom_smooth(method="lm", se=F, formula="y ~ x") +
    theme_bw() +
    xlab("t_compliance") + 
    ylab("différence de poids à 6 mois") +
    scale_color_manual(values=c(bleuclair, rose))
plot_diff_weight_vs_compliance

On peut à présent effectuer un test de corrélation en prenant tous les points ensemble :

cor.test(d$diff_weight, d$t_compliance)

## 
##  Pearson's product-moment correlation
## 
## data:  d$diff_weight and d$t_compliance
## t = 8.9827, df = 498, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2953913 0.4464803
## sample estimates:
##       cor 
## 0.3734095

Puis un test de corrélation en prenant uniquement le groupe control :

d_control <- d %>% filter(group == "control")
cor.test(d_control$diff_weight, d_control$t_compliance)

## 
##  Pearson's product-moment correlation
## 
## data:  d_control$diff_weight and d_control$t_compliance
## t = 0.77159, df = 248, p-value = 0.4411
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.07558871  0.17196019
## sample estimates:
##        cor 
## 0.04893724

Et enfin, le même test de corrélation en prenant uniquement le groupe treatment :

d_treatment <- d %>% filter(group == "treatment")
cor.test(d_treatment$diff_weight, d_treatment$t_compliance)

## 
##  Pearson's product-moment correlation
## 
## data:  d_treatment$diff_weight and d_treatment$t_compliance
## t = 13.982, df = 248, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5883172 0.7280229
## sample estimates:
##       cor 
## 0.6639239

Que concluons-nous ?

Le modèle linéaire

Le modèle linéaire a le bon goût d’unifier tous ces tests dans un unique framework.

Il consiste à modéliser une variable de réponse \(Y_i\) chez l’individu \(i\) comme une fonction linéaire de nos variables explicatrices \(X_{i,j}\), à laquelle on rajoute un bruit Gaussien centré. Une fonction linéaire n’est rien de plus qu’une somme pondérée de nos variables \(X_{i,j}\), soit :

\[ Y_i = a_1 X_{i,1} + a_2 X_{i,2} + ... + a_J X_{i,J} + \epsilon_i, ~~~\text{ avec } \epsilon_i \sim \mathcal{N}(0, \sigma^2) \]

On suppose que tous les \((Y_i)\) sont indépendantes et suivent la même loi. Le but du jeu, pour le statisticien, consiste ensuite à déterminer les valeurs des coefficients \(a_j\) grâce aux observations.

Bonne nouvelle : la syntaxe de R pour fitter un modèle linéaire est assez rapide et intuitive.

Un simple intercept

Premier cas simple : on aimerait modéliser la différence de poids comme étant uniquement une moyenne + un petit bruit Gaussien, et calculer cette moyenne comme on l’avait fait lors du premier test de Student.

res1 <- lm(formula="diff_weight ~ 1", data=d)
summary(res1)

## 
## Call:
## lm(formula = "diff_weight ~ 1", data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.9791  -4.6822  -0.1511   3.9584  25.4556 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.4311     0.3029   8.027 7.22e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.772 on 499 degrees of freedom

On retombe sur un test qu’on avait déjà réalisé : lequel ?

Un effet de variable discrète

Deuxième cas simple : on aimerait modéliser la différence de poids comme étant la somme d’un coefficient particulier chez les individus du groupe “control” et un autre coefficient chez les individus du groupe “treatment”. Soit exactement ce qu’on avait fait avec un test de Student comparant les deux moyennes.

Il se trouve que R, lorsqu’on lui donne à manger une variable à plusieurs modalités, va automatiquement la “scinder” en sous-variables indicatrices de chaque modalité. La syntaxe reste donc très simple :

# première version avec moyenne commune aux deux groupes
# et un petit ajout pour le groupe traitement
res2 <- lm(formula="diff_weight ~ group", data=d)
summary(res2)

## 
## Call:
## lm(formula = "diff_weight ~ group", data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.7634  -4.2507   0.0323   4.2829  22.9418 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -0.08262    0.39807  -0.208    0.836    
## grouptreatment  5.02750    0.56295   8.931   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.294 on 498 degrees of freedom
## Multiple R-squared:  0.138,  Adjusted R-squared:  0.1363 
## F-statistic: 79.76 on 1 and 498 DF,  p-value: < 2.2e-16

# deuxième version avec moyenne dans chaque groupe
res3 <- lm(formula="diff_weight ~ 0 + group", data=d)
summary(res3)

## 
## Call:
## lm(formula = "diff_weight ~ 0 + group", data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.7634  -4.2507   0.0323   4.2829  22.9418 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## groupcontrol   -0.08262    0.39807  -0.208    0.836    
## grouptreatment  4.94487    0.39807  12.422   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.294 on 498 degrees of freedom
## Multiple R-squared:  0.2366, Adjusted R-squared:  0.2335 
## F-statistic: 77.18 on 2 and 498 DF,  p-value: < 2.2e-16

On a exactement le même modèle sous-jacent, à un jeu de ré-écriture des paramètres près. Ici, les p-value indiquées correspondent à un test d’adéquation du paramètre à 0.

Si on souhaite comparer les deux moyennes et retomber sur le même résultat qu’au début du document, quelle version faudrait-il regarder ?

Troisième possibilité, notre différence de poids dépend d’une variable prenant 3 modalités ou plus (comme on avait proposé lors de la réalisation du test ANOVA ci-dessus). La même syntaxe fonctionne toujours aussi bien :

res4 <- lm(formula="diff_weight ~ 0 + age_group", data=d)
summary(res4)

## 
## Call:
## lm(formula = "diff_weight ~ 0 + age_group", data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.6381  -4.2684   0.1706   4.4074  25.7966 
## 
## Coefficients:
##           Estimate Std. Error t value Pr(>|t|)    
## age_group   1.0451     0.1469   7.112 3.98e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.857 on 499 degrees of freedom
## Multiple R-squared:  0.09204,    Adjusted R-squared:  0.09023 
## F-statistic: 50.59 on 1 and 499 DF,  p-value: 3.977e-12

L’output standard du modèle linéaire ne nous calcule des p-values que pour les tests d’adéquation de chaque paramètre à 0. Bien que l’ANOVA soit effectivement un test réalisé au sein du même framework de modèle linéaire, la p-value associée n’est récupérable, en pratique avec R, que via la fonction aov vue précédemment.

Un effet de variable continue

Dernier exemple similaire à ce qu’on avait fait précédemment : expliquer notre variable de réponse “différence de poids” par le “temps de compliance”.

res5 <- lm(formula="diff_weight ~ t_compliance", data=d)
summary(res5)

## 
## Call:
## lm(formula = "diff_weight ~ t_compliance", data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.6251  -4.2947   0.0167   4.1866  16.3475 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.04517    0.39382  -0.115    0.909    
## t_compliance  0.12592    0.01402   8.983   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.289 on 498 degrees of freedom
## Multiple R-squared:  0.1394, Adjusted R-squared:  0.1377 
## F-statistic: 80.69 on 1 and 498 DF,  p-value: < 2.2e-16

Quel test vu précédemment peut être remplacé par celui-ci ?

Vers plus de complexification

Et si on souhaite un effet additif du groupe et du temps de compliance :

res6 <- lm(formula="diff_weight ~ 0 + group + t_compliance", data=d)
summary(res6)

## 
## Call:
## lm(formula = "diff_weight ~ 0 + group + t_compliance", data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.7324  -3.8577   0.3067   3.9575  14.2405 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## groupcontrol   -2.38373    0.44163  -5.398 1.05e-07 ***
## grouptreatment  2.46635    0.45241   5.452 7.88e-08 ***
## t_compliance    0.12152    0.01295   9.382  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.807 on 497 degrees of freedom
## Multiple R-squared:  0.3515, Adjusted R-squared:  0.3476 
## F-statistic: 89.78 on 3 and 497 DF,  p-value: < 2.2e-16

Ce n’est pas exactement ce qu’on avait fait tout à l’heure sur le graphe : on avait alors deux régressions linéaires, une dans chaque groupe. R possède encore une fois une syntaxe pratique pour ça :

res7 <- lm(formula="diff_weight ~ 0 + group * t_compliance", data=d)
summary(res7)

## 
## Call:
## lm(formula = "diff_weight ~ 0 + group * t_compliance", data = d)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.9342  -3.7083   0.2186   3.8499  14.5558 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## groupcontrol                -0.31917    0.45725  -0.698    0.485    
## grouptreatment              -0.03519    0.48810  -0.072    0.943    
## t_compliance                 0.01249    0.01633   0.765    0.445    
## grouptreatment:t_compliance  0.23168    0.02380   9.734   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.326 on 496 degrees of freedom
## Multiple R-squared:  0.4555, Adjusted R-squared:  0.4511 
## F-statistic: 103.7 on 4 and 496 DF,  p-value: < 2.2e-16

On a introduit super facilement ce qu’on appelle des interactions entre nos variables. Le modèle linéaire, et la syntaxe de R, sont donc des outils de modélisation extrêmement puissants/versatiles. Ils demandent à bien réfléchir en amont à ce qu’on souhaite étudier pour éviter une trop grande multiplication de tests statistiques.

A vous de proposer : quel modèle linéaire vous semble le plus intéressant pour diff_weight ?

Intro GLM : la régression logistique

Jusqu’à présent, nous avons cherché à expliquer une variable de réponse numérique par des variables explicatives de type varié, numériques ou catégorielles. A présent, on souhaiterait expliquer des variables binaires, par exemple has_lost_weight.

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 un 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.

Que pouvons-nous en conclure ?

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

Et sur ce deuxième exemple, que pouvons-nous en conclure ?

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

Que pouvons-nous en conclure ?

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 continues. Les points principaux à retenir sont les suivants :

1. Pour comparer deux moyennes, ou une moyenne à une norme, le T test est la bonne option, à condition que les données soient en grand nombre, ou qu’elles soient distribuées normalement.
2. Pour comparer plus que deux moyennes, une ANOVA fait l’affaire, sous des conditions identiques.
3. Pour quantifier la corrélation entre deux variables continues distribuées normalement, un test de corrélation de Pearson peut être utilisé.
4. Le modèle linéaire unifie tous les tests précédents dans un unique framework, où une variable de réponse continue \(Y\) est modélisée comme une combinaison linéaire de variables explicatives \(X\) (discrètes ou continues).

Par ailleurs, nous avons ensuite élargi le framework du modèle linéaire (LM) au modèle linéaire généralisé (GLM) :

1. La régression logistique est un exemple de GLM qui permet d’analyser des variables de réponses binaires.
2. Comme pour un LM, le GLM nous permet d’intégrer l’influence de variables explicatives numériques ou catégorielles.
3. L’outcome est très similaire, on interprète : la p-value, et le signe des coefficients.

Notez pour finir qu’il existe d’autres types de régressions bien adaptées à des types de données différentes. Parmi les plus célèbres, vous rencontrerez très certainement 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.