Statistiques au LAPEC – Analyse de données de comptage

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 les données :

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 des données de comptage.

Le GLM

Principe

Le GLM : Generalized Linear Model permet de généraliser une grande part de ce qu’on a vu jusqu’à présent.

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

Régression linéaire

Dans le cas d’une variable de réponse \(Y\) numérique, on suppose que la loi de \(Y\) sachant \(X\) est une loi Normale de paramètre \(\mathbb{E}(Y|X)\).

On modélise directement cette moyenne comme une combinaison linéaire de variables explicatives, i.e. \(g\) est la fonction identité.

Avec R, on pourra faire un modèle linéaire en passant par les fonctions lm ou glm, en choisissant la bonne famille (“gaussian”).

res_lm_1 <- lm("diff_weight ~ group*t_compliance", data=d)
summary(res_lm_1)

## 
## Call:
## lm(formula = "diff_weight ~ group*t_compliance", data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2949 -1.5933  0.1597  1.6771  7.8364 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 -0.960433   0.233093  -4.120 4.43e-05 ***
## grouptreatment               0.216735   0.328371   0.660  0.50954    
## t_compliance                -0.004425   0.007897  -0.560  0.57549    
## grouptreatment:t_compliance -0.032257   0.011248  -2.868  0.00431 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.522 on 496 degrees of freedom
## Multiple R-squared:  0.04863,    Adjusted R-squared:  0.04288 
## F-statistic: 8.451 on 3 and 496 DF,  p-value: 1.74e-05

res_lm_2 <- glm("diff_weight ~ group*t_compliance", family=gaussian, data=d)
summary(res_lm_2)

## 
## Call:
## glm(formula = "diff_weight ~ group*t_compliance", family = gaussian, 
##     data = d)
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 -0.960433   0.233093  -4.120 4.43e-05 ***
## grouptreatment               0.216735   0.328371   0.660  0.50954    
## t_compliance                -0.004425   0.007897  -0.560  0.57549    
## grouptreatment:t_compliance -0.032257   0.011248  -2.868  0.00431 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 6.362231)
## 
##     Null deviance: 3317.0  on 499  degrees of freedom
## Residual deviance: 3155.7  on 496  degrees of freedom
## AIC: 2350.1
## 
## Number of Fisher Scoring iterations: 2

Ces deux régressions identiques correspondent au graphique suivant :

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

Régression logistique

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 transformation suivante de notre combinaison linéaire de variables explicatives :

\[ 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. C’est elle qu’on choisit comme fonction de lien.

Avec R, la syntaxe de base pour faire une régression logistique illustre bien le fait que ce soit un GLM :

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

## 
## Call:
## glm(formula = "is_depressed ~ group*weight_t0", family = binomial, 
##     data = d)
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -15.21067    1.91745  -7.933 2.14e-15 ***
## grouptreatment            -0.16950    2.74736  -0.062    0.951    
## weight_t0                  0.21993    0.02769   7.943 1.97e-15 ***
## grouptreatment:weight_t0  -0.01093    0.03896  -0.281    0.779    
## ---
## 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: 427.93  on 496  degrees of freedom
## AIC: 435.93
## 
## Number of Fisher Scoring iterations: 5

Cette régression logistique peut s’illustrer avec le graphe suivant :

p_logistic <- d %>% ggplot(aes(x=weight_t0, y=as.numeric(is_depressed), col=group)) + 
    geom_point() +
    theme_bw() +
    xlab("Poids initial (kg)") +
    ylab("Indicatrice de dépression") +
    geom_smooth(formula = "y~x",
                method = "glm",
                method.args = list(family = binomial),
                se = FALSE,
                fullrange = T) +
    scale_color_manual(values=c(bleuclair, rose))
p_logistic

Régression de Poisson

Principe

La régression de Poisson est encore un exemple classique de GLM, utilisé pour modéliser une variable de réponse typiquement discrète positive, typiquement issue d’un processus de comptage d’événements.

On suppose que notre \(Y\) suit une loi de Poisson, de paramètre \(\lambda(X)\), et on cherche à modéliser la dépendance de \(\lambda\) aux variables explicatives \(X\). Comme \(\lambda\) doit être strictement positif, on a l’idée d’appliquer la fonction exponentielle à notre combinaison linéaire de variables explicatives :

\[ \lambda(X) = \exp \left( a_0 + a_1 x_1 + a_2 x_2 + ... + a_n x_n \right) \]

Cas univarié

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

On souhaiterait dans un premier temps expliquer n_meals_glycemia_out_of_target à partir de l’IMC bmi. Première chose à faire : un graphe !

p <- d %>% ggplot(aes(x=bmi, y=n_meals_glycemia_out_of_target)) + 
    geom_point() +
    theme_bw() +
    xlab("IMC (kg/m^2)") +
    ylab("Nombre de repas avec glycémie hors cible")
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 = poisson),
                se = FALSE,
                color = bleuclair) +
    geom_smooth(formula = "y~x", method = "lm", color = rose, se = F) 

res_poisson <- glm("n_meals_glycemia_out_of_target ~ bmi", family=poisson, data=d)
summary(res_poisson)

## 
## Call:
## glm(formula = "n_meals_glycemia_out_of_target ~ bmi", family = poisson, 
##     data = d)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.543143   0.048252   -52.7   <2e-16 ***
## bmi          0.204633   0.001578   129.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 14325.7  on 499  degrees of freedom
## Residual deviance:  2173.2  on 498  degrees of freedom
## AIC: 4186.9
## 
## Number of Fisher Scoring iterations: 5

Que concluons-nous ?

Cas multivarié

Regardons ce qu’on obtient en proposant un effet joint du bras de traitement et de l’IMC :

res_poisson_2 <- glm("n_meals_glycemia_out_of_target ~ bmi*group", family=poisson, data=d)
summary(res_poisson_2)

## 
## Call:
## glm(formula = "n_meals_glycemia_out_of_target ~ bmi*group", family = poisson, 
##     data = d)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -1.812655   0.052625 -34.445  < 2e-16 ***
## bmi                 0.188905   0.001670 113.102  < 2e-16 ***
## grouptreatment     -1.298835   0.148460  -8.749  < 2e-16 ***
## bmi:grouptreatment  0.018963   0.005252   3.611 0.000305 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 14325.7  on 499  degrees of freedom
## Residual deviance:  1094.9  on 496  degrees of freedom
## AIC: 3112.5
## 
## Number of Fisher Scoring iterations: 4

Ce modèle s’illustre avec le graphique suivant :

p_poisson <- d %>% ggplot(aes(x=bmi, y=n_meals_glycemia_out_of_target, col=group)) + 
    geom_point() +
    theme_bw() +
    xlab("IMC (kg/m^2)") +
    ylab("Nombre de repas avec glycémie hors cible") +
    geom_smooth(formula = "y~x",
                method = "glm",
                method.args = list(family = poisson),
                se = FALSE,
                fullrange = T) +
    scale_color_manual(values=c(bleuclair, rose))
p_poisson

Que concluons-nous ?

Quel autre modèle auriez-vous envie d’envisager ?

Conclusion

Les take-home messages sont les suivants :

• le modèle linéaire (LM) est un GLM avec fonction de lien identité, qui sert à modéliser des données numériques.
• la régression logistique est un GLM avec fonction de lien logit, qui sert à modéliser des données binaires.
• la régression de Poisson est un GLM avec fonction de lien log, qui sert à modéliser des données de comptage.
• tous ces modèles partagent une théorie commune et s’interprètent de façon similaire.