La régression linéaire simple

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 quantitative Y à l’aide d’une variable quantitative X en utilisant un modèle linéaire :

Y=a + b \times X + \textrm{bruit}

En pratique, on devra vérifier que le bruit suit une loi normale.

Un exemple

Étudions si le score de relation varie en fonction de l’âge des personnes.

Calculons les paramètres avec la commande lm(Y~X) :

mod.sr.a<-lm(score.relation~age,data=satis) 
summary(mod.sr.a)


Call:
lm(formula = score.relation ~ age, data = satis)

Residuals:
Min 1Q Median 3Q Max
-22.047 -2.447 1.339 3.620 5.673

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 33.76634 0.84496 39.962 <2e-16 *** age 0.02668 0.01485 1.796 0.0734 .
— –
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 4.609 on 347 degrees of freedom
(185 observations deleted due to missingness)
Multiple R-squared: 0.00921, Adjusted R-squared: 0.006355
F-statistic: 3.226 on 1 and 347 DF, p-value: 0.07336

On y voit les coefficients avec leur p-value respectifs :
– a = 33.76634 avec p=<2e-16 – b = 0.02668 avec p = 0.0734 Mais on y trouve aussi la valeur du p-value (0.07336), que l’on avait dans la test de corrélation linéaire :

cor.test(satis$score.relation,satis$age,use="complete.obs")


Pearson's product-moment correlation


data: satis$score.relation and satis$age
t = 1.796, df = 347, p-value = 0.07336
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.009102243 0.198945290
sample estimates:
cor
0.09596955

Cette valeur étant trop élevée (supérieure à 0.05), on peut pas affirmer que les deux variables sont liés.

Test de normalité du bruit

Pour juger de la validité d’un modèle de régression, on peut tracer un histogramme :

hist(resid(mod.sr.a),col="brown")

histo_resid_01.png

On voit que le bruit ne suit pas une loi normale. Le modèle n’est pas valide.

Représentation graphique

plot(satis$age,satis$score.relation,las=1,xlab="âge",ylab="Score de relation")
abline(mod.sr.a,lwd=2,col="red")

reglin_01.png

Cas particulier : la variable X est binaire

Cherchons à expliquer le score de relation par le fait de recommander l’établissement.

Préparons la variable binaire :

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

Calculons les coefficients :

mod.sr.r<-lm(score.relation~recommander.b,data=satis) summary(mod.sr.r)


Call:
lm(formula = score.relation ~ recommander.b, data = satis)

Residuals:
Min 1Q Median 3Q Max
-19.1091 -1.6826 0.8909 2.9975 7.8909

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 32.1091 0.3918 81.959 <2e-16 *** recommander.b 4.5735 0.4763 9.602 <2e-16 ***
— –
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Residual standard error: 4.109 on 338 degrees of freedom
(194 observations deleted due to missingness)
Multiple R-squared: 0.2143, Adjusted R-squared: 0.212
F-statistic: 92.19 on 1 and 338 DF, p-value: < 2.2e-16

On voit que les valeurs de p sont faibles.

On lit aussi b=4.5735.

Si la variable recommande augmente de 1 (passe de 0 à 1), le score de relation augmente de 4.5735.

Comparons avec le test de Student :

t.test(score.relation~recommander.b,data=satis,var.equal=TRUE)


Two Sample t-test


data: score.relation by recommander.b
t = -9.6016, df = 338, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-5.510459 -3.636576
sample estimates:
mean in group 0 mean in group 1
32.10909 36.68261

On y retrouve la p-value très faible.

On voit aussi l’écart des moyennes entre les deux sous-groupes :

36.68261-32.10909


[1] 4.57352

Une autre façon de trouver ces valeurs est de réaliser deux prédictions :

predict(mod.sr.r,data=satis,data.frame(recommander.b=0))


1
32.10909

predict(mod.sr.r,data=satis,data.frame(recommander.b=1))


1
36.68261

36.68261-32.10909


[1] 4.57352