Le modèle linéaire
Ana Fermin
Exercice
Nous considérons les données salaires.tex.
salary : salaire brut actuel en $ par an
salbegin : salaire de départ en $ par an
minority : appartenance à une minorité (Non, Oui)
jobtime : nombre de mois depuis l’entrée dans l’entreprise
prevexp : nombre de mois de travail avant l’entrée dans l’entreprise
educ : nombre d’années d’étude
sex : H,F
Lecture et codage des données
Salaire <- read.table('salaires.txt',header=TRUE)
summary(Salaire)## salary salbegin jobtime prevexp
## Min. : 15750 Min. : 9000 Min. :63.00 Min. : 0.00
## 1st Qu.: 24000 1st Qu.:12488 1st Qu.:72.00 1st Qu.: 19.25
## Median : 28875 Median :15000 Median :81.00 Median : 55.00
## Mean : 34420 Mean :17016 Mean :81.11 Mean : 95.86
## 3rd Qu.: 36938 3rd Qu.:17490 3rd Qu.:90.00 3rd Qu.:138.75
## Max. :135000 Max. :79980 Max. :98.00 Max. :476.00
## educ minority sex
## Min. : 8.00 Min. :0.0000 Min. :0.0000
## 1st Qu.:12.00 1st Qu.:0.0000 1st Qu.:0.0000
## Median :12.00 Median :0.0000 Median :0.0000
## Mean :13.49 Mean :0.2194 Mean :0.4557
## 3rd Qu.:15.00 3rd Qu.:0.0000 3rd Qu.:1.0000
## Max. :21.00 Max. :1.0000 Max. :1.0000Devant ces données,que faut-il faire ? Quelles sont les bonnes habitudes à adquerir ? Comment re-coder et transformer des variables dans votre tableau des données ?
Manipuler les données à l’aide des fonctions mutate et select du package dplyr (dans tidyverse).
library(tidyverse) # Organisation des données et plotsSalaire <-Salaire %>%
mutate(minority = factor(minority,labels=c("Non","Oui")),
sex = factor(sex,labels=c("H","F")))
summary(Salaire)## salary salbegin jobtime prevexp
## Min. : 15750 Min. : 9000 Min. :63.00 Min. : 0.00
## 1st Qu.: 24000 1st Qu.:12488 1st Qu.:72.00 1st Qu.: 19.25
## Median : 28875 Median :15000 Median :81.00 Median : 55.00
## Mean : 34420 Mean :17016 Mean :81.11 Mean : 95.86
## 3rd Qu.: 36938 3rd Qu.:17490 3rd Qu.:90.00 3rd Qu.:138.75
## Max. :135000 Max. :79980 Max. :98.00 Max. :476.00
## educ minority sex
## Min. : 8.00 Non:370 H:258
## 1st Qu.:12.00 Oui:104 F:216
## Median :12.00
## Mean :13.49
## 3rd Qu.:15.00
## Max. :21.00Régression linéaire simple
- On souhaite expliquer la variable salary en fonction de salbegin à l’aide de la régression linéaire simple.
ggplot(data = Salaire, aes(x=salbegin,y=salary)) + geom_point()
- Ajuster le modèle de règression linèaire simple en exècutant la commande suivante (vous pouvez consulter l’aide avec la commande help(lm)).
fit.lm1 <- lm(salary~salbegin, data=Salaire)- Analyser les résultats de la régression ajustée obtenues à l’aide de la commande ci-dessous
summary(fit.lm1)##
## Call:
## lm(formula = salary ~ salbegin, data = Salaire)
##
## Residuals:
## Min 1Q Median 3Q Max
## -35424 -4031 -1154 2584 49293
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.928e+03 8.887e+02 2.17 0.0305 *
## salbegin 1.909e+00 4.741e-02 40.28 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8115 on 472 degrees of freedom
## Multiple R-squared: 0.7746, Adjusted R-squared: 0.7741
## F-statistic: 1622 on 1 and 472 DF, p-value: < 2.2e-16La variable salbegin est-elle significative pour expliquer salary ?
Relever et interpréter la valeur observée du coefficient \({\rm R}^2\).
Donner l’équation du modèle ajusté.
ggplot(data = Salaire, aes(x=salbegin,y=salary)) +
geom_point() + geom_smooth(method = "lm", fullrange = TRUE,se = FALSE) 
Régression linéaire multiple
- On souhaite expliquer la variable salary en fonction de toutes les autres variables à l’aide de la régression linéaire multiple.
Quelles sont les variables significatives au seuil 5% ?
Tester la nullité de chaque coefficient au risque \(5\%\) en n’oubliant pas de donner les hypothèses nulle et alternative du test et la loi de la statistique du test. Que peut-on conclure ?
Tester la significativité globale du modèle à un niveau de risque de \(5\%\) en n’oubliant pas de donner les hypothèses nulle et alternative du test et la loi de la statistique du test. Que peut-on conclure ?
Relever et interpréter la valeur observée du coefficient \({\rm R}^2\) (R-squared).
Relever et interpréter la valeur observée du coefficient \({\rm R}^2\)-ajusté (Adjusted R-squared).
Analiser les résidus. Tracer quelques graphiques pour la validation du résidus (installer le package ggfortify et utliser la fonction autoplot)
Pensez vous que le modèle ajusté est pertinent ? Justifier votre réponse.
lm(salary ~ salbegin+jobtime+prevexp+educ+minority+sex, data=Salaire)##
## Call:
## lm(formula = salary ~ salbegin + jobtime + prevexp + educ + minority +
## sex, data = Salaire)
##
## Coefficients:
## (Intercept) salbegin jobtime prevexp educ
## -12059.761 1.707 156.511 -18.764 589.259
## minorityOui sexF
## -1377.447 -2418.777lm(salary ~ ., data=Salaire)##
## Call:
## lm(formula = salary ~ ., data = Salaire)
##
## Coefficients:
## (Intercept) salbegin jobtime prevexp educ
## -12059.761 1.707 156.511 -18.764 589.259
## minorityOui sexF
## -1377.447 -2418.777fit.lm2 <- lm(salary ~ ., data=Salaire)
summary(fit.lm2)##
## Call:
## lm(formula = salary ~ ., data = Salaire)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28359 -3917 -908 2550 45568
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.206e+04 3.482e+03 -3.463 0.000582 ***
## salbegin 1.706e+00 6.124e-02 27.868 < 2e-16 ***
## jobtime 1.565e+02 3.405e+01 4.597 5.53e-06 ***
## prevexp -1.876e+01 3.601e+00 -5.210 2.83e-07 ***
## educ 5.893e+02 1.664e+02 3.542 0.000437 ***
## minorityOui -1.377e+03 8.513e+02 -1.618 0.106317
## sexF -2.419e+03 7.990e+02 -3.027 0.002605 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7398 on 467 degrees of freedom
## Multiple R-squared: 0.8147, Adjusted R-squared: 0.8123
## F-statistic: 342.2 on 6 and 467 DF, p-value: < 2.2e-16fit.lm3 <- lm(salary ~ . - minority, data=Salaire)
summary(fit.lm3)##
## Call:
## lm(formula = salary ~ . - minority, data = Salaire)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28446 -4048 -1014 2598 45788
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.255e+04 3.475e+03 -3.612 0.000337 ***
## salbegin 1.723e+00 6.051e-02 28.472 < 2e-16 ***
## jobtime 1.545e+02 3.408e+01 4.534 7.37e-06 ***
## prevexp -1.944e+01 3.583e+00 -5.424 9.36e-08 ***
## educ 5.930e+02 1.666e+02 3.559 0.000410 ***
## sexF -2.233e+03 7.921e+02 -2.819 0.005021 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7410 on 468 degrees of freedom
## Multiple R-squared: 0.8137, Adjusted R-squared: 0.8117
## F-statistic: 408.7 on 5 and 468 DF, p-value: < 2.2e-16library("ggfortify") #autoplot (graphiques des residus)
autoplot(fit.lm3,2:3)
Que faut-il faire ? Y a–t-il des transformations qui s’imposent ?
- Appliquer une transformation logarithmique aux variables salary et salbegin. Relever le log(salaire) actuel en dessous duquel se situent 50% des employés de l’échantillon et le log(salaire) actuel au-dessus duquel se situent 25% des employés de l’échantillon.
Devant ces données,que faut-il faire ? Quelles sont les bonnes habitudes à adquerir ? Comment creer et eliminer des variables dans un autre tableau des données ?
- Manipuler les données à l’aide des fonctions mutate et select.
Salaire1 <- Salaire %>% mutate(logsalary = log(salary),
logsalbegin = log(salbegin))
Salaire1 <- Salaire1 %>% select(-salary, -salbegin)
summary(Salaire1)## jobtime prevexp educ minority sex
## Min. :63.00 Min. : 0.00 Min. : 8.00 Non:370 H:258
## 1st Qu.:72.00 1st Qu.: 19.25 1st Qu.:12.00 Oui:104 F:216
## Median :81.00 Median : 55.00 Median :12.00
## Mean :81.11 Mean : 95.86 Mean :13.49
## 3rd Qu.:90.00 3rd Qu.:138.75 3rd Qu.:15.00
## Max. :98.00 Max. :476.00 Max. :21.00
## logsalary logsalbegin
## Min. : 9.665 Min. : 9.105
## 1st Qu.:10.086 1st Qu.: 9.432
## Median :10.271 Median : 9.616
## Mean :10.357 Mean : 9.669
## 3rd Qu.:10.517 3rd Qu.: 9.769
## Max. :11.813 Max. :11.290- On souhaite expliquer la variable log(salary) en fonction de toutes les autres variables (log(salbegin), jobtime, prevexp, educ et sex) à l’aide de la régression linéaire.
Tester la significativité globale du modèle à un niveau de risque de 5%. Que peut-on conclure ?
Relever et interpréter la valeur observée du coefficient \({\rm R}^2\).
Pensez vous que le modèle ajusté est pertinent ? Justifier votre réponse.
fit.lm4 <- lm(logsalary ~ .- minority , data=Salaire1)
summary(fit.lm4)##
## Call:
## lm(formula = logsalary ~ . - minority, data = Salaire1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.44031 -0.11068 -0.00731 0.09581 0.86024
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.116e+00 3.125e-01 3.571 0.000392 ***
## jobtime 4.517e-03 7.579e-04 5.960 4.97e-09 ***
## prevexp -5.527e-04 7.932e-05 -6.968 1.10e-11 ***
## educ 1.071e-02 3.912e-03 2.737 0.006431 **
## sexF -4.995e-02 1.844e-02 -2.708 0.007019 **
## logsalbegin 9.107e-01 3.382e-02 26.924 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1639 on 468 degrees of freedom
## Multiple R-squared: 0.8317, Adjusted R-squared: 0.8299
## F-statistic: 462.6 on 5 and 468 DF, p-value: < 2.2e-16autoplot(fit.lm4,2:3)
- Quel modèle preferez vous ? Justifier votre réponse.
Selection de modèles
library(leaps) # Selection de modeles
selection <- regsubsets(logsalary ~ ., data=Salaire1)
plot(selection,scale="r2") # qui donne le modèle complet
plot(selection,scale="adjr2")
plot(selection,scale="bic")
library(MASS) # Selection de modeles
reg <- lm(logsalary~.,data=Salaire1)
selection <- stepAIC(reg,direction="backward")## Start: AIC=-1709
## logsalary ~ jobtime + prevexp + educ + minority + sex + logsalbegin
##
## Df Sum of Sq RSS AIC
## <none> 12.506 -1709.0
## - minority 1 0.0615 12.567 -1708.7
## - educ 1 0.2087 12.714 -1703.2
## - sex 1 0.2317 12.738 -1702.3
## - jobtime 1 0.9637 13.470 -1675.8
## - prevexp 1 1.2124 13.718 -1667.1
## - logsalbegin 1 18.1355 30.641 -1286.2summary(selection)##
## Call:
## lm(formula = logsalary ~ jobtime + prevexp + educ + minority +
## sex + logsalbegin, data = Salaire1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.44407 -0.10824 -0.01063 0.09517 0.85214
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.226e+00 3.204e-01 3.827 0.000147 ***
## jobtime 4.541e-03 7.570e-04 5.999 3.98e-09 ***
## prevexp -5.373e-04 7.986e-05 -6.729 5.03e-11 ***
## educ 1.091e-02 3.909e-03 2.792 0.005459 **
## minorityOui -2.881e-02 1.901e-02 -1.516 0.130276
## sexF -5.509e-02 1.873e-02 -2.941 0.003430 **
## logsalbegin 8.996e-01 3.457e-02 26.024 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1636 on 467 degrees of freedom
## Multiple R-squared: 0.8325, Adjusted R-squared: 0.8304
## F-statistic: 386.9 on 6 and 467 DF, p-value: < 2.2e-16