La régression logistique

Un fichier de données

Travaillons sur une enquête de satisfactions dans un hôpital, récupéré lors d’un cours de FUN (France Unité Numérique). Les fichiers sont disponibles ici.

Par la suite, ces données seront stockées dans la variable satis .

Dans ces données comme souvent, les observations sont en lignes et les variables en colonnes.

Présentation

On cherche à expliquer une variable binaire Y à l’aide d’autres variables.

Pour que déboucher sur une valeur numérique, on utilise le modèle suivant :

\ln \left( \frac{P(Y=1)}{P(Y=0} \right)=\ln \left( \frac{P(Y=1)}{1-P(Y=1} \right) =a+b \times X_1+c \times X_2 + \cdots =K 

On aura alors pour les prédictions :

p(Y=1)=\frac{e^K}{1+e^K} 

La procédure est voisine de la celle de la régression multiple.

Un exemple

On cherche à expliquer la variable binaire recommander.b en fonction de plusieurs variables :
– sexe.cat
– age
– score.relation
– profession.cat

Les variables catégorielles

Il faudra faire très attention à bien convertir les variables catégorielles pour qu’elles ne soient pas interprétées comme des variables quantitatives.

satis$sexe.cat<-factor(satis$sexe,labels=c("H","F")) 
satis$profession.cat<-factor(satis$profession) 
satis$service.cat<-factor(satis$service)

Dans le logiciel, la première catégorie servira de référence. Par exemple, pour le sexe codé (“H”,”F”), “H” sera la référence. On cherchera donc à savoir si le fait d’être une femme modifie le score de relation.

Pour changer cette référence, on peut le faire de la sorte :

satis$sexe.cat <- relevel(satis$sexe.cat,ref="F")

La variable binaire

satis$recommander.b<-ifelse(satis$recommander==2,1,0)

La régression logistique

mod.r.m<-glm(recommander.b~sexe.cat+age+score.relation+profession.cat,data=satis,family=binomial("logit"))

Calculons les valeurs de p par catégorie :

drop1(mod.r.m,.~.,test="Chisq")


Single term deletions


Model:
recommander.b ~ sexe.cat + age + score.relation + profession.cat
Df Deviance AIC LRT Pr(>Chi)
332.48 354.48
sexe.cat 1 336.74 356.74 4.259 0.03905 *
age 1 333.00 353.00 0.516 0.47243
score.relation 1 409.14 429.14 76.665 < 2e-16 ***
profession.cat 7 341.19 349.19 8.711 0.27408


— -
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

La valeur de p pour l’âge est faible (inférieure à 0.05).

Regardons l’ensemble des résultats :

summary(mod.r.m)


Call:
glm(formula = recommander.b ~ sexe.cat + age + score.relation +
profession.cat, family = binomial("logit"), data = satis)

Deviance Residuals:
Min 1Q Median 3Q Max
-2.2358 -0.8142 0.5003 0.6937 2.0298

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -23.745655 882.744298 -0.027 0.9785
sexe.catF 0.595364 0.293048 2.032 0.0422 *
age 0.006309 0.008796 0.717 0.4732
score.relation 0.256624 0.034247 7.493 6.71e-14 ***
profession.cat2 15.641720 882.743601 0.018 0.9859
profession.cat3 15.338964 882.743436 0.017 0.9861
profession.cat4 14.630198 882.743434 0.017 0.9868
profession.cat5 14.881076 882.743450 0.017 0.9866
profession.cat6 15.129875 882.743504 0.017 0.9863
profession.cat7 15.301040 882.743609 0.017 0.9862
profession.cat8 14.363204 882.743485 0.016 0.9870

— –
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 419.70 on 334 degrees of freedom
Residual deviance: 332.48 on 324 degrees of freedom
(199 observations deleted due to missingness)
AIC: 354.48


Number of Fisher Scoring iterations: 13

Interprétation

Servons-nous de deux prédictions pour mieux comprendre ce qui se passe :

predict(mod.r.m,data=satis,data.frame(sexe.cat="H",age=30,score.relation=40,profession.cat="5"))


1
1.589661

predict(mod.r.m,data=satis,data.frame(sexe.cat="F",age=30,score.relation=40,profession.cat="5"))


1
2.185025

2.185025-1.589661


[1] 0.595364

Les autres paramètres étant inchangés, on a une augmentation de la valeur prédite de 0.595364.

Mais interprétons les prédictions.

Pour un homme :

exp(1.589661)/(1+exp(1.589661))


[1] 0.8305684

Pour une femme :

exp(2.185025)/(1+exp(2.185025))


[1] 0.8988967

La probabilité de recommander est de 0.8305684 pour homme et de 0.8988967 pour une femme (avec les autres paramètres choisis fixés).

Cas particulier : une seule variable explicative binaire

Cherchons à expliquer la variable recommander.b à l’aide de la variable binaire sexe.

mod.r.s<-glm(recommander.b~sexe,data=satis,family=binomial("logit")) 
summary(mod.r.s)


Call:
glm(formula = recommander.b ~ sexe, family = binomial("logit"),
data = satis)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.4964 -1.4602 0.8890 0.9188 0.9188

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.64401 0.14454 4.456 8.37e-06 ***
sexe 0.08039 0.21085 0.381 0.703

— –
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 516.95 on 404 degrees of freedom
Residual deviance: 516.81 on 403 degrees of freedom
(129 observations deleted due to missingness)
AIC: 520.81


Number of Fisher Scoring iterations: 4

Comparons avec l’odds-ratio.

twoby2(1-satis$sexe,1-satis$recommander.b)


2 by 2 table analysis:


Outcome : 0
Comparing : 0 vs. 1

0 1 P(0) 95% conf. interval
0 130 63 0.6736 0.6043 0.7360
1 139 73 0.6557 0.5892 0.7165

95% conf. interval
Relative Risk: 1.0273 0.8945 1.1798
Sample Odds Ratio: 1.0837 0.7169 1.6383
Conditional MLE Odds Ratio: 1.0835 0.7021 1.6745
Probability difference: 0.0179 -0.0740 0.1088

Exact P-value: 0.7523
Asymptotic P-value: 0.703



exp(0.08039)


[1] 1.08371

Calculons le p-value avec le test du Chi2 :

chisq.test(satis$recommander.b,satis$sexe,correct=FALSE)


Pearson's Chi-squared test


data: satis$recommander.b and satis$sexe
X-squared = 0.1454, df = 1, p-value = 0.703

On retrouve la même p-value.