1. régression linéaire simple

On travaille sur le format long de l’étude de santé mentale en prison qui contient 26 variables.

str(smp.l)
## 'data.frame':    799 obs. of  26 variables:
##  $ age         : int  31 49 50 47 23 34 24 52 42 45 ...
##  $ prof        : Factor w/ 8 levels "agriculteur",..: 3 NA 7 6 8 6 3 2 6 6 ...
##  $ duree       : int  4 NA 5 NA 4 NA NA 5 4 NA ...
##  $ discip      : int  0 0 0 0 1 0 0 0 1 0 ...
##  $ n.enfant    : int  2 7 2 0 1 3 5 2 1 2 ...
##  $ n.fratrie   : int  4 3 2 6 6 2 3 9 12 5 ...
##  $ ecole       : int  1 2 2 1 1 2 1 2 1 2 ...
##  $ separation  : int  0 1 0 1 1 0 1 0 1 0 ...
##  $ juge.enfant : int  0 0 0 0 NA 0 1 0 1 0 ...
##  $ place       : int  0 0 0 1 1 0 1 0 0 0 ...
##  $ abus        : int  0 0 0 0 0 0 0 0 1 1 ...
##  $ grav.cons   : int  1 2 2 1 2 1 5 1 5 5 ...
##  $ dep.cons    : int  0 0 0 0 1 0 1 0 1 0 ...
##  $ ago.cons    : int  1 0 0 0 0 0 0 0 0 0 ...
##  $ ptsd.cons   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ alc.cons    : int  0 0 0 0 0 0 0 0 1 1 ...
##  $ subst.cons  : int  0 0 0 0 0 0 1 0 1 0 ...
##  $ scz.cons    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ char        : int  1 1 1 1 1 1 1 1 4 1 ...
##  $ rs          : int  2 2 2 2 2 1 3 2 3 2 ...
##  $ ed          : int  1 2 3 2 2 2 3 2 3 2 ...
##  $ dr          : int  1 1 2 2 2 1 2 2 1 2 ...
##  $ suicide.s   : int  0 0 0 1 0 0 3 0 4 0 ...
##  $ suicide.hr  : int  0 0 0 0 0 0 1 0 1 0 ...
##  $ suicide.past: int  0 0 0 0 1 0 1 0 1 0 ...
##  $ dur.interv  : int  NA 70 NA 105 NA NA 105 84 78 60 ...

Une droite de régression est une droite qui passe au mieux par l’ensemble d’un nuage de points, c’est la droite qui minimise la somme des carrés des distances qui séparent les points et leur projection sur la droite.

plot(smp.l$age, smp.l$dur.interv, xlab = 'Âge', ylab = 'Durée de l\'entretien')

# Comme les variables âge et durée de l'entretien sont discrètes, plusieurs valeurs sont superposées, donc on utilise jitter()
plot(jitter(smp.l$age), jitter(smp.l$dur.interv), xlab = 'Âge', ylab = 'Durée de l\'entretien')

# On peut augmenter l'intensité du bruit dans jitter
plot(jitter(smp.l$age), jitter(smp.l$dur.interv, factor = 4), xlab = 'Âge', ylab = 'Durée de l\'entretien')

# On ajoute la droite de régression logistique abline()
plot(jitter(smp.l$age), jitter(smp.l$dur.interv, factor = 4), xlab = 'Âge', ylab = 'Durée de l\'entretien')
abline(lm(smp.l$dur.interv ~ smp.l$age), lwd = 2) # lm() permet de calculer les coordonnées de la droite de régression, abline() de tracer la ligne

On a donc obtenu une droite \(durée = a + b*âge\) et la question est \(b ≠ 0\)?. Pour y répondre il faut passer au modèle de régression linéaire \(durée = a+B*âge + bruit\) avec un bruit qui va suivre une loi normale.
C’est encore lm()(linear modèle) qui va permettre d’effectuer la régression logistique.

mod1 <- lm(dur.interv ~ age, data = smp.l) # on peut ajouter dans le lm une option subset =
summary(mod1)
## 
## Call:
## lm(formula = dur.interv ~ age, data = smp.l)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -62.470 -14.402  -1.712  12.341  60.055 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 57.04091    2.22028  25.691   <2e-16 ***
## age          0.12625    0.05375   2.349   0.0191 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.57 on 745 degrees of freedom
##   (52 observations deleted due to missingness)
## Multiple R-squared:  0.00735,    Adjusted R-squared:  0.006018 
## F-statistic: 5.516 on 1 and 745 DF,  p-value: 0.0191

On obtient ainsi le résultat de notre équation \(durée = 57.04091 + 0.12625 * âge + bruit\) et au bout de la ligne âge le p de savoir si \(b≠0\) : 0.0191, au risque de 5% la pente de la droite de régression est effectivement différente de zéro.
Quand l’âge augmente, la durée de l’entretien avec le clinicien augmente également.

On peut également tester la nullité de la corrélation (Pearson : corrélation de deux variables dont au moins une est quantitative et au moins une dont la distribution est normale) entre durée de l’entretien et âge du patient (cf. D-Tests statistiques) et on obtient les mêmes résultats.

cor.test(smp.l$dur.interv, smp.l$age)
## 
##  Pearson's product-moment correlation
## 
## data:  smp.l$dur.interv and smp.l$age
## t = 2.3487, df = 745, p-value = 0.0191
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.01408787 0.15650345
## sample estimates:
##        cor 
## 0.08573358

Finalement entre test \(b=0\) et test de nullité de corrélation on obtient le même p, à quoi sert de faire une corrélation linéaire simple ? :
- d’abord on obtient de b (si on prend 10 ans, c’est 10*b soit 1,2 minutes en plus) - on peut l’utiliser avec une variable qualitative

avec une variable qualitative (par exemple binaire)

plot(smp.l$dep.cons, jitter(smp.l$dur.interv))
abline(lm(smp.l$dur.interv ~ smp.l$dep.cons), lwd = 2)

On obtient la droite \(durée = a + b*dep\) et on se demande si \(b≠0\) ?
Pourrait-on utiliser un test t de student qui comparerait la durée moyenne d’interview selon que les détenus sont déprimés ou non ?

mod2 <- lm(dur.interv ~ dep.cons, data = smp.l)
summary(mod2)
## 
## Call:
## lm(formula = dur.interv ~ dep.cons, data = smp.l)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -62.538 -13.923   1.077  12.077  61.077 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  58.9234     0.9041  65.171  < 2e-16 ***
## dep.cons      7.6143     1.4481   5.258  1.9e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.33 on 747 degrees of freedom
##   (50 observations deleted due to missingness)
## Multiple R-squared:  0.03569,    Adjusted R-squared:  0.0344 
## F-statistic: 27.65 on 1 and 747 DF,  p-value: 1.9e-07
t.test(smp.l$dur.interv ~ smp.l$dep.cons, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  smp.l$dur.interv by smp.l$dep.cons
## t = -5.2583, df = 747, p-value = 1.9e-07
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -10.457001  -4.771515
## sample estimates:
## mean in group 0 mean in group 1 
##        58.92341        66.53767

On obtient exactement les mêmes p et les coefficients \(a = 58.9234\) et \(b = 66.53767 - 58.92341 = 7.6143\) (la différence des deux moyennes).
On peut dire que la régression linéaire est une généralisation du test de nullité d’un coefficient de corrélation (Pearson) et également une généralisation du test t de Student.

2. régression linéaire multiple et analyse de variance

La durée de l’entretien est-elle liée à plusieurs variables (qui sont toutes plus au moins liées les unes des autres, par exemple la dépression est plus ou moins liée à l’âge, l’abus de substance en partie liée à la dépression…) ?
- âge
- existence d’une dépression
- existence d’un abus de substance
- existence d’un trouble schizophrénique

La formule pourrait s’écrire \(durée = a + b*âge + c*dep + d*subst + e*scz +bruit\) et chercher si \(b≠0\) qui revient à dire “est-ce que la durée de l’entretien est liée à l’âge, toute choses égales par ailleurs (à dépression, abus de substance et schizophrénie constantes) ?”, ou encore, “de combien augmente la durée de l’entretien quand l’âge augmente d’un an à dépression, abus de substance et schizophrénie constantes ?”.
On peut aussi dire que b caractérise la durée de l’entretien après ajustement sur la dépression, l’abus de substance et la schizophrénie.

mod3 <- lm(dur.interv ~ age + dep.cons + subst.cons + scz.cons, data = smp.l)
summary(mod3)
## 
## Call:
## lm(formula = dur.interv ~ age + dep.cons + subst.cons + scz.cons, 
##     data = smp.l)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -63.654 -14.522  -1.193  11.482  62.482 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 48.90105    2.62213  18.649  < 2e-16 ***
## age          0.22096    0.05708   3.871 0.000118 ***
## dep.cons     7.38932    1.44783   5.104 4.24e-07 ***
## subst.cons   5.25157    1.74318   3.013 0.002678 ** 
## scz.cons     2.27256    2.52323   0.901 0.368062    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.1 on 742 degrees of freedom
##   (52 observations deleted due to missingness)
## Multiple R-squared:  0.05833,    Adjusted R-squared:  0.05325 
## F-statistic: 11.49 on 4 and 742 DF,  p-value: 4.692e-09

Un détenu déprimé (par rapport à un non déprimé), toutes choses égales par ailleurs (en annulant l’effet de l’âge, de l’abus de substance et de la schizophrénie) a une durée d’interview qui durera 7.38932 minutes de plus.
Un détenu déprimé avec un abus de substance aura une durée d’interview de 7.38932 + 5.25157 = 12,63549 min de plus.
Au risque alpha de 5%, on ne peut pas dire que la présence d’une schizophrénie augmente la durée de l’interview (\(p = 0.368062\)).

avec variable explicative multicatégorielle

entre les différentes catégories de la variable

Attention, dans ce modèle, la variable à expliquer est quantitative, les variables explicatives peuvent être quantitiatives ou binaires.
Il existe des situations où l’on pourrait souhaiter mettre des variables explicatives à plus de deux classes (par exemple la profession). Pour intégrer la variable profession, il faudra la recoder en sept variables binaires :
- f1 : artisan = 1 vs 0
- f2 : cadre = 1 vs 0
- f3 : prof inter = 1 vs 0
- f4 : ouvrier = 1 vs 0
- f5 : employé = 1 vs 0
- f6 : autre = 1 vs 0
- f7 : sans prof = 1 vs 0
- la variable agriculteur est codée implicitement (c’est tous les f = 0)
En fait, R va faire lui même la conversion en autant de variables que le nombre de catégories moins une. Attention donc, lors du codage de la profession à ne pas le coder en chiffre mais bien en mots pour que R reconnaisse le côté catégoriel.

mod4 <- lm(dur.interv ~ age + dep.cons + subst.cons + scz.cons + prof, data = smp.l)
summary(mod4)
## 
## Call:
## lm(formula = dur.interv ~ age + dep.cons + subst.cons + scz.cons + 
##     prof, data = smp.l)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -63.280 -14.164  -1.337  10.959  63.184 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             62.79202   10.20779   6.151 1.26e-09 ***
## age                      0.21289    0.05884   3.618 0.000317 ***
## dep.cons                 7.36792    1.45840   5.052 5.53e-07 ***
## subst.cons               5.34589    1.76902   3.022 0.002599 ** 
## scz.cons                 2.50439    2.54734   0.983 0.325863    
## profartisan            -11.48515    9.82936  -1.168 0.243005    
## profautre              -10.28748   10.33482  -0.995 0.319862    
## profcadre              -19.29636   10.38568  -1.858 0.063574 .  
## profemploye            -13.55809    9.76340  -1.389 0.165358    
## profouvrier            -14.01270    9.72111  -1.441 0.149880    
## profprof.intermediaire -13.01926    9.96911  -1.306 0.191977    
## profsans emploi        -14.27866    9.71782  -1.469 0.142174    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.11 on 731 degrees of freedom
##   (56 observations deleted due to missingness)
## Multiple R-squared:  0.06595,    Adjusted R-squared:  0.05189 
## F-statistic: 4.692 on 11 and 731 DF,  p-value: 5.825e-07

On interprète les résultats de la façon suivante, par rapport aux agriculteurs (qui est la profession de référence - car la première par ordre alphabétique), les artisants passent 11.48515 min de moins en interview. Il peut être très utile de changer la modalité de référence.

smp.l$prof <- relevel(smp.l$prof, ref = 'ouvrier') # attention à ce que la variable smp.l$prof soit bien au format as.factor
mod5 <- lm(dur.interv ~ age + dep.cons + subst.cons + scz.cons + prof, data = smp.l)
summary(mod5)
## 
## Call:
## lm(formula = dur.interv ~ age + dep.cons + subst.cons + scz.cons + 
##     prof, data = smp.l)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -63.280 -14.164  -1.337  10.959  63.184 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            48.77932    2.83938  17.180  < 2e-16 ***
## age                     0.21289    0.05884   3.618 0.000317 ***
## dep.cons                7.36792    1.45840   5.052 5.53e-07 ***
## subst.cons              5.34589    1.76902   3.022 0.002599 ** 
## scz.cons                2.50439    2.54734   0.983 0.325863    
## profagriculteur        14.01270    9.72111   1.441 0.149880    
## profartisan             2.52755    2.48989   1.015 0.310381    
## profautre               3.72522    3.99637   0.932 0.351567    
## profcadre              -5.28366    4.25567  -1.242 0.214798    
## profemploye             0.45460    2.12659   0.214 0.830785    
## profprof.intermediaire  0.99344    2.95809   0.336 0.737089    
## profsans emploi        -0.26596    1.87727  -0.142 0.887375    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.11 on 731 degrees of freedom
##   (56 observations deleted due to missingness)
## Multiple R-squared:  0.06595,    Adjusted R-squared:  0.05189 
## F-statistic: 4.692 on 11 and 731 DF,  p-value: 5.825e-07

mesurer l’effet global d’une variable multicatégorielle

On n’a toujours pas la réponse à la question de savoir si la profession (quelle quelle soit) a un effet sur la durée de l’interview. Il faut ajouter drop1(), on voit qu’il n’y a pas d’effet significatif de la profession en général sur la durée de l’interview.

drop1(mod5, .~., test = 'F')
## Single term deletions
## 
## Model:
## dur.interv ~ age + dep.cons + subst.cons + scz.cons + prof
##            Df Sum of Sq    RSS    AIC F value    Pr(>F)    
## <none>                  266846 4395.6                      
## age         1    4778.4 271624 4406.8 13.0899 0.0003173 ***
## dep.cons    1    9317.1 276163 4419.1 25.5233 5.527e-07 ***
## subst.cons  1    3333.6 270180 4402.8  9.1322 0.0025992 ** 
## scz.cons    1     352.8 267199 4394.6  0.9666 0.3258633    
## prof        7    2295.5 269142 4388.0  0.8983 0.5071556    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mesurer l’intéraction entre deux variables explicatives

mod3 <- lm(dur.interv ~ age + dep.cons + subst.cons + scz.cons, data = smp.l)
summary(mod3)
## 
## Call:
## lm(formula = dur.interv ~ age + dep.cons + subst.cons + scz.cons, 
##     data = smp.l)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -63.654 -14.522  -1.193  11.482  62.482 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 48.90105    2.62213  18.649  < 2e-16 ***
## age          0.22096    0.05708   3.871 0.000118 ***
## dep.cons     7.38932    1.44783   5.104 4.24e-07 ***
## subst.cons   5.25157    1.74318   3.013 0.002678 ** 
## scz.cons     2.27256    2.52323   0.901 0.368062    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.1 on 742 degrees of freedom
##   (52 observations deleted due to missingness)
## Multiple R-squared:  0.05833,    Adjusted R-squared:  0.05325 
## F-statistic: 11.49 on 4 and 742 DF,  p-value: 4.692e-09

On a dit précédement qu’un détenu déprimé avec un abus de substance aura une durée d’interview de 7.38932 + 5.25157 = 12,63549 min de plus.
Ceci est vrai s’il n’y a pas de synergie, de potentialisation entre la dépression et l’abus de substance sur la durée de l’entretien.

mod6 <- lm(dur.interv ~ age + dep.cons * subst.cons + scz.cons, data = smp.l)
summary(mod6)
## 
## Call:
## lm(formula = dur.interv ~ age + dep.cons * subst.cons + scz.cons, 
##     data = smp.l)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -62.032 -14.251  -1.163  11.472  62.313 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         49.51693    2.65788  18.630  < 2e-16 ***
## age                  0.21728    0.05711   3.805 0.000154 ***
## dep.cons             6.15780    1.69775   3.627 0.000306 ***
## subst.cons           3.17244    2.29849   1.380 0.167931    
## scz.cons             1.97233    2.53094   0.779 0.436059    
## dep.cons:subst.cons  4.49688    3.24296   1.387 0.165963    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.08 on 741 degrees of freedom
##   (52 observations deleted due to missingness)
## Multiple R-squared:  0.06077,    Adjusted R-squared:  0.05443 
## F-statistic: 9.588 on 5 and 741 DF,  p-value: 7.024e-09

Cela ajoute la ligne pour dépression et abus de substance avec un \(p = 0.165963\), cela veut dire que statistiquement il n’y a pas d’interaction entre dépression et abus de substance sur la durée de l’interview et que donc quand un détenu a une dépression et un abus de substance, la durée est bien augmentée de 7.38932 + 5.25157 = 12,63549 min.
Un piège très classique, lorsqu’on met un terme d’interaction dans un modèle de régression c’est qu’on ne peut plus interpréter les effets principaux, c’est à dire que les lignes dépression et abus de substance ne sont plus interprétables.

Analyse de variance

C’est un cas particulier de régression linéaire multiple où les variables explicatives sont toutes catégorielles. Le terme est d’origine historique (de simplification de calcul).

On peut cependant garder ce nom d’ANOVA dans le cas particulier de l’analyse de variance à un facteur qui revient à une généralisation du test t (comparaison de la moyenne d’une variable aléatoire quantitative à plus de deux groupes). C’est une méthode qui permet donc de comparer les moyennes lorsqu’il y a plus de deux groupes. Par exemple, si on souhaite comparer la durée moyenne d’interview en fonction de la profession, on ne peut pas réaliser un test t car il y a plus de deux professions et il faut réaliser une ANOVA à un facteur.

mod7 <- lm(dur.interv ~ prof, data = smp.l)
summary(mod7)
## 
## Call:
## lm(formula = dur.interv ~ prof, data = smp.l)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -61.731 -13.826  -1.731  12.947  58.912 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             61.7315     1.3359  46.211   <2e-16 ***
## profagriculteur         17.0185     9.9071   1.718   0.0863 .  
## profartisan              2.0941     2.5033   0.837   0.4031    
## profautre                2.4993     4.0755   0.613   0.5399    
## profcadre               -4.7750     4.3063  -1.109   0.2679    
## profemploye              0.3220     2.1742   0.148   0.8823    
## profprof.intermediaire   1.3440     3.0096   0.447   0.6553    
## profsans emploi         -0.6432     1.9168  -0.336   0.7373    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 19.63 on 735 degrees of freedom
##   (56 observations deleted due to missingness)
## Multiple R-squared:  0.008295,   Adjusted R-squared:  -0.001149 
## F-statistic: 0.8783 on 7 and 735 DF,  p-value: 0.5231

Nous obtenons la durée de l’interview par rapport à la profession de référence qui est l’ouvrier (modifiée plus haut).
Pour obtenir l’effet global, il faut utiliser drop1()

drop1(mod7, .~., test = 'F') # F pour Fisher
## Single term deletions
## 
## Model:
## dur.interv ~ prof
##        Df Sum of Sq    RSS    AIC F value Pr(>F)
## <none>              283316 4432.1               
## prof    7    2369.9 285686 4424.3  0.8783 0.5231

Avec un \(p = 0.5231\) il n’y a pas de différence moyenne de durée d’interview en fonction des professions.

3. conditions de validités du modèle de régression

Elles ne sont pas très faciles à vérifier en pratique :
- la variable à expliquer doit être quantitative
- normalité du bruit
- la variance du bruit ne doit dépendre ni des valeurs de la variable à expliquer (= “indépendante”), ni des valeurs des variables explicatives (= “dépendantes”).
- le bruit doit être un “vrai” bruit (pas de structure de corrélation évidente)

Si les conditions de validité ne sont pas vérifiées on réalisera un test de Kruskal-Wallis (qui est une ANOVA non paramétrique).

En pratique, il faut au moins vérifier la normalité du bruit, en faisant un histogramme des résidus du modèle (le bruit)

mod3 <- lm(dur.interv ~ age + dep.cons + subst.cons + scz.cons, data = smp.l)
hist(resid(mod3)) # histogramme des résidus de ce modèle

NB : les résidus sont obtenus en effectuant la soustraction de la valeur prédite de y à la valeur observée de y pour un x particulier. Il est important de noter que la valeur prédite provient de notre droite de régression. La valeur observée provient de notre ensemble de données. \(Résiduel=y_{observé}−y_{prédit}\)
cf. https://www.greelane.com/fr/science-technologie-math%C3%A9matiques/math/what-are-residuals-3126253/

Lire la suite

A-Introduction et représentations graphiques
B-Dispersion et intervalles de confiance
C-Coefficient de corrélation
D-Tests statistiques
E-Régression linéaire
F-Régression logistique
G-Données de survie
H-Statistique exploratoire multidimensionnelle
I-Multiplicité des tests