Retour à la page d’accueil

Cette page web accompagne les quatre premières séances du cours de Méthodes quantitatives pour la sociologie 2 (année universitaire 2015-2016). Elle reprend les codes R utilisés pendant les séances de cours ainsi que des exercices supplémentaires avec solution. L’ensemble des données nécessaires sont accessibles à l’adresse http://teaching.slmc.fr/mqs2/.

Télécharger le support complet du cours en version imprimable (en date du 27 avril 2016)

Télécharger le document de travail de l’Insee sur la régression logistique (de Cédric Afsa, paru en mars 2016)

Cette page a été conçue avec R Markdown sous RStudio et compilée le 2016-05-14 à 15:45:09. Elle peut être enregistrée sous la forme d’un fichier .html autonome pour être consultée hors connexion.

Séance 1 : Introduction aux méthodes de régression

L’objectif de cette première séance est d’introduire les spécificités des méthodes de régression à partir d’un exemple fictif suivi tout du long : l’analyse de la relation entre poids et activité physique régulière. Le cours procède en trois temps :

  1. Lien entre le poids et la taille : régression linéaire simple
  2. Impact de l’activité physique : régression linéaire multiple
  3. L’efficacité des régimes : endogénéité et causalité inverse

Simulation des données utilisées dans le cours

L’échantillon de données utilisé a été entièrement simulé à l’aide des fonctions rbinom() et rnorm() de R (taper ?rbinom ou ?rnorm dans R pour en savoir plus). Ces fonctions font appel à un générateur de nombres pseudo-aléatoires qu’il est possible d’initialiser avec la fonction set.seed() (lien pour en savoir plus sur la génération de nombres aléatoires sous R)

# Initialisation du générateur de nombres pseudo-aléatoires de R
set.seed(1)

# Taille de l'échantillon simulé
n <- 100

# Simulation des variables
femme <- rbinom(n, 1, 0.5)
sport <- femme*rbinom(n, 1, 0.65) + (1-femme)*rbinom(n,1,0.45)
taille <- (1-femme)*rnorm(n,177,8) + femme*rnorm(n,165,8)
imc <- sport*rnorm(n,22,1) + (1-sport)*rnorm(n,26,1)
poids <- imc*(taille/100)**2
regime <- (imc > 23 + rnorm(n,0,2.5))*1
l.poids <- regime*(poids + rnorm(n,2,2)) + (1-regime)*(poids + rnorm(n,0,2))

# Arrondi pour les variables numériques
taille <- round(taille,0)
poids <- round(poids,0)
l.poids <- round(l.poids,0)

# Affichage des 10 premiers individus
data.frame(poids,taille,sport)[1:10,]
##    poids taille sport
## 1     88    181     0
## 2     79    177     0
## 3     61    164     1
## 4     65    163     0
## 5     67    165     0
## 6     69    174     1
## 7     44    146     1
## 8     63    170     1
## 9     71    168     0
## 10    80    192     1

Lien entre le poids et la taille : régression linéaire simple

Le poids des individus étant fortement lié à leur taille, on va chercher à contrôler l’effet de la taille dans l’analyse de la relation entre poids et activité physique régulière.

On commence donc par analyser la relation entre le poids et la taille : comme il s’agit de deux variables de nature quantitative, on peut représenter leur association par un nuage de points et la synthétiser à l’aide du coefficient de corrélation linéaire de Pearson.

# Construction du nuage de points
plot(taille,poids,xlab="Taille en cm",ylab="Poids en kg")

# Ajout d'un titre général
title(main="Relation entre le poids et la taille des individus")

# Calcul du coefficient de corrélation avec la fonction cor()
cor(taille,poids)
## [1] 0.8645734

On peut également estimer le modèle de régression linéaire simple suivant: \[\text{(m1)} \quad poids_i = \beta_0 + \beta_1 \times taille_i + \epsilon_i\]

Pour ce faire, on utilise la fonction lm(): on l’utilise pour effectuer l’estimation puis on stocke ses résultats dans un objet (m1 dans l’exemple) ; la fonction summary permet ensuite d’afficher une synthèse de ces résultats.

# Estimation du modèle m1 et affectation à l'objet m1
m1 <- lm(poids ~ taille)

# Résumé du modèle m1
summary(m1)
## 
## Call:
## lm(formula = poids ~ taille)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.204  -5.020   1.726   4.492  16.837 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -107.16662   10.46301  -10.24   <2e-16 ***
## taille         1.04095    0.06112   17.03   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.483 on 98 degrees of freedom
## Multiple R-squared:  0.7475, Adjusted R-squared:  0.7449 
## F-statistic: 290.1 on 1 and 98 DF,  p-value: < 2.2e-16

En moyenne dans l’échantillon, une taille supérieure de 1 cm est associée à un poids supérieur de 1,04 kg. Les valeurs prédites par le modèle sont stockées dans l’objet m1 et accessibles via l’entrée fitted.value :

# Récupération et affichage des valeurs prédites par m1
predict1 <- round(m1$fitted.values,2)
data.frame(poids,taille,sport,predict1)[1:10,]
##    poids taille sport predict1
## 1     88    181     0    81.24
## 2     79    177     0    77.08
## 3     61    164     1    63.55
## 4     65    163     0    62.51
## 5     67    165     0    64.59
## 6     69    174     1    73.96
## 7     44    146     1    44.81
## 8     63    170     1    69.79
## 9     71    168     0    67.71
## 10    80    192     1    92.70

On peut également récupérer la valeur des coefficients de m1 pour tracer la droite correspondant à ces valeurs prédites :

# Construction du nuage de points
plot(taille,poids,xlab="Taille en cm",ylab="Poids en kg")

# Récupération des valeurs des coefficients avec la fonction (coef())
m1beta0 <- coef(m1)[1]
m1beta1 <- coef(m1)[2]

# Ajout d'une droite d'équation y = m1beta0 + m1beta1 * x
abline(a = m1beta0, b = m1beta1)

# Ajout d'un titre général
title(main="Représentation graphique de la prédiction du modèle 1")

Impact de l’activité physique : régression linéaire multiple

Quand on compare le poids moyen des individus selon qu’ils pratiquent ou non une activité physique régulière, la relation apparaît assez nettement (- 18 kg pour ceux qui pratiquent une activité physique régulière).

data.frame(
  sport=c(0,1)
  , poids=f(as.numeric(by(poids,sport,mean)))
  , taille=f(as.numeric(by(taille,sport,mean)))
)
##   sport poids taille
## 1     0 78,76 174,24
## 2     1 60,84 166,76

Cependant, les personnes pratiquant une activité physique régulière étant en moyenne plus petites que les autres (- 8 cm), l’intensité de cette relation serait certainement moins forte si on contrôlait par la taille.

Dans un premier temps, on peut revenir au graphique et y différencier les individus selon qu’ils pratiquent ou non une activité physique régulière.

# Construction du nuage de points
plot(taille,poids,pch=sport,xlab="Taille en cm",ylab="Poids en kg")

# Ajout de la légende
legend("topleft",legend=c("Pas d'activité physique régulière","Activité physique régulière"),pch=c(0,1))

# Ajout d'un titre général
title(main="Poids, taille et activité physique des individus")

Cette représentation semble indiquer l’idée que l’activité physique régulière est associée, à taille égale, à un poids moindre. Pour approfondir cette analyse, on met en oeuvre un modèle de régression multiple : \[\text{(m2)} \quad poids_i = \beta_0 + \beta_1 \times taille_i + \beta_2 \times sport_i + \epsilon_i\] où \(sport_i\) est une variable indicatrice.

# Estimation du modèle m2
m2 <- lm(poids ~ taille + sport)

# Résumé du modèle m2
summary(m2)
## 
## Call:
## lm(formula = poids ~ taille + sport)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.1134 -2.4723 -0.5312  2.4404 13.1832 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -69.48486    6.24333  -11.13   <2e-16 ***
## taille        0.85085    0.03573   23.82   <2e-16 ***
## sport       -11.55417    0.76182  -15.17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.549 on 97 degrees of freedom
## Multiple R-squared:  0.9251, Adjusted R-squared:  0.9236 
## F-statistic:   599 on 2 and 97 DF,  p-value: < 2.2e-16

L’estimation du coefficient \(\beta_2\) semble indiquer que la pratique d’une activité physique régulière est associée à un poids moyen d’environ 11,6 kg plus faible. Cette valeur est moindre que les 18 kg qu’indiquaient les statistiques univariées : ceci est dû à l’effet de structure de la taille sur l’activité physique régulière (le fait que dans l’échantillon les personnes qui pratiquent une activité physique régulière mesurent en moyenne 8 cm de moins que les autres).

À nouveau on peut calculer les valeurs prédites par le modèle et les représenter sur le graphique.

predict2 <- round(m2$fitted.values,2)
data.frame(poids,taille,sport,predict1,predict2)[1:10,]
##    poids taille sport predict1 predict2
## 1     88    181     0    81.24    84.52
## 2     79    177     0    77.08    81.12
## 3     61    164     1    63.55    58.50
## 4     65    163     0    62.51    69.20
## 5     67    165     0    64.59    70.90
## 6     69    174     1    73.96    67.01
## 7     44    146     1    44.81    43.18
## 8     63    170     1    69.79    63.60
## 9     71    168     0    67.71    73.46
## 10    80    192     1    92.70    82.32

Le modèle étant de meilleure qualité que le précédent, les valeurs prédites sont en général plus proches des valeurs effectives.

# Construction du nuage de points
plot(taille,poids,pch=sport,xlab="Taille en cm",ylab="Poids en kg")

# Ajout des droites représentant la prédiction du modèle m2
abline(a = coef(m2)[1], b = coef(m2)[2])
abline(a = coef(m2)[1] + coef(m2)[3], b = coef(m2)[2],lty=2)

# Ajout d'une légende mixte
legend("topleft",legend=c("Pas d'activité physique régulière","Activité physique régulière"),pch=c(0,1),lty=c(1,2))

# Ajout d'un titre général
title(main="Représentation graphique de la prédiction du modèle 2")

Au-delà de la valeur des coefficients, il est possible d’interpréter la significativité des coefficients à l’aide de leur p-valeur ou d’un intervalle de confiance (cf. le cours).

L’efficacité des régimes : endogénéité et causalité inverse

On cherche enfin à mesurer la relation entre poids et régimes amaigrissants, en contrôlant par la taille et par l’activité physique régulière. Le modèle de régression multiple que l’on estime est donc le suivant : \[\text{(m3)} \quad poids_i = \beta_0 + \beta_1 \times taille_i + \beta_2 \times sport_i + \beta_3 \times regime_i + \epsilon_i\]

m3 <- lm(poids ~ taille + sport + regime)
summary(m3)
## 
## Call:
## lm(formula = poids ~ taille + sport + regime)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3520 -2.5752 -0.3603  2.3478 12.8767 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -70.15622    6.09121 -11.518   <2e-16 ***
## taille        0.84487    0.03491  24.204   <2e-16 ***
## sport       -10.56913    0.84272 -12.542   <2e-16 ***
## regime        2.04703    0.82824   2.472   0.0152 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.459 on 96 degrees of freedom
## Multiple R-squared:  0.9296, Adjusted R-squared:  0.9274 
## F-statistic: 422.4 on 3 and 96 DF,  p-value: < 2.2e-16

Le coefficient associé au fait d’avoir suivi un régime amaigrissant durant les 12 derniers mois est estimé à la valeur 2,05 : autrement dit, en moyenne les personnes qui ont suivi un régime amaigrissant durant les 12 derniers moids pèsent 2,05 kg de plus que les autres.

Interrogées, les personnes déclarant avoir suivi un régime amaigrissant au cours des 12 derniers mois déclarent pourtant avoir perdu du poids, ce qui est confirmé par la comparaison des poids moyens au moment de l’enquête et un an avant (variable l.poids) :

data.frame(
  regime=c("Pas de régime amaigrissant","Régime amaigrissant")
  , l.poids=as.numeric(by(l.poids,regime,mean))
  , poids=as.numeric(by(poids,regime,mean))
  , delta=as.numeric(by(poids,regime,mean)-by(l.poids,regime,mean))
)
##                       regime  l.poids    poids      delta
## 1 Pas de régime amaigrissant 63.07692 63.43590  0.3589744
## 2        Régime amaigrissant 76.80328 75.34426 -1.4590164

Ce phénomène illustre le mécanisme de causalité inverse (cf. cours) : même si individuellement les personnes ayant suivi un régime amaigrissant ont perdu du poids, comme par ailleurs les personnes en surpoids ont une probabilité plus élevée de suivre de tels régimes, le coefficient estimé par la régression est globalement positif.

C’est ce type de situation qui justifie l’expression “corrélation n’est pas causalité” : le modèle de régression ne révèle que des associations statistiques, et non des relations de cause à effet. Pour mesurer l’efficacite des régimes en s’affranchissant de ce type d’endogénéité, on peut par exemple utiliser des données de panel (observations répétées dans le temps).

 

Séance 2 : Régression linéaire

En parallèle des éléments théoriques introduits par le cours, un exemple est suivi tout du long : l’analyse des déterminants du salaire à partir des données de l’enquête Emploi en continu (EEC) du quatrième trimestre 2012 (extraction au 1/20ème). Les fichiers dont sont extraites ces données ainsi que leur documentation sont librement téléchargeables sur le site de l’Insee.

Les différentes variables explicatives sont ajoutées progressivement au modèle : âge, sexe, niveau de diplôme atteint, temps partiel, secteur d’emploi (public ou privé). Elles permettent de revenir sur certains concepts fondamentaux vus en cours :

  1. Introduction de la variable d’âge et choix de spécification
  2. Introduction de la variable de sexe et problèmes de colinéarité
  3. Introduction de la variable de diplôme et choix de la modalité de référence
  4. Comparaison des modèles : R2, R2 ajusté et test de significativité globale des coefficients de Fisher
  5. Introduction des autres variables : utilisation de variables croisées et test d’hypothèses complexes

Note : Ce plan diffère de celui suivi pendant le cours. Les modèles de régression estimés avec SAS et qui figurent dans le support de cours ont été adaptés en conséquence.

Chargement et examen des données

Les données utilisées sont stockées au format RDS dans le fichier http://teaching.slmc.fr/mqs2/ee12t4.rds. Le format RDS est un format propre à R de stockage des données sous une forme compressée (taper ?readRDS dans R pour en savoir plus).

  # Lecture des données et assignation à l'objet e
  e <- readRDS(gzcon(url("http://teaching.slmc.fr/mqs2/ee12t4.rds")))
  
  # Passage en minuscule des noms de variables et affichage de la structure de l'objet e
  colnames(e) <- tolower(colnames(e))
  str(e)
## 'data.frame':    1783 obs. of  19 variables:
##  $ ident    : chr  "G0C5B116" "G0C5B1I6" "G0D5GIX6" "G0D5OGWD" ...
##  $ noi      : int  1 1 2 1 2 1 1 2 4 1 ...
##  $ extri1613: num  1416 1524 1594 1988 2522 ...
##  $ sexe     : int  1 1 2 1 2 1 1 2 1 1 ...
##  $ age      : int  64 43 50 24 30 30 51 51 16 82 ...
##  $ cse      : int  37 23 46 56 NA 54 48 68 NA NA ...
##  $ acteu    : int  1 1 1 1 3 1 1 1 3 3 ...
##  $ stc      : int  2 1 2 2 NA 2 2 2 NA NA ...
##  $ tam1d    : int  NA 1 NA NA NA NA NA NA NA NA ...
##  $ aidref   : int  NA NA NA NA 4 NA NA NA 5 NA ...
##  $ salred   : int  3350 NA 2508 1500 NA 1517 3163 950 NA NA ...
##  $ tpp      : int  1 1 1 1 NA 1 1 1 NA NA ...
##  $ ddipl    : int  4 1 3 4 5 1 5 7 6 7 ...
##  $ nbagenf  : int  0 1 0 0 0 0 1 1 1 0 ...
##  $ duhab    : int  7 7 7 6 NA 6 6 6 NA NA ...
##  $ pub3fp   : int  4 NA 4 4 NA 4 4 4 NA NA ...
##  $ naia     : int  1948 1969 1962 1988 1982 1982 1961 1961 1995 1930 ...
##  $ fordat   : int  1970 1992 1982 2007 2000 2006 1978 1977 NA NA ...
##  $ ancentr  : int  108 72 50 30 NA 5 372 132 NA NA ...

Passer le nom des variables des tables importées en minuscule n’est absolument pas obligatoire, mais cela permet de faire moins d’erreurs par la suite.

La variable dépendante de l’étude, salred, est le salaire mensuel en euros redressé : certains enquêtés préfèrent en effet déclarer leur salaire en tranches, voire refusent de le déclarer. Des techniques plus ou moins complexes sont utilisées pour imputer une valeur à l’ensemble des salariés.

# Distribution de salred
summary(e$salred)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      31    1250    1600    1851    2168   11420    1004

# On supprime les 1 004 observations pour lesquelles salred est manquant (inactifs, chômeurs, indépendants, etc.)
e <- e[!is.na(e$salred),]
nrow(e)
## [1] 779

Introduction de la variable d’âge et choix de spécification

On commence par estimer un modèle de régression linéaire simple de la forme : \[\text{(m1)} \quad salaire_i = \beta_0 + \beta_1 \times age_i + \epsilon_i\] en utilisant la fonction lm():

# Estimation du modèle m1
m1 <- lm(salred ~ age, data = e)

# Résumé du modèle m1
summary(m1)
## 
## Call:
## lm(formula = salred ~ age, data = e)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2123.9  -603.6  -166.3   315.2  9864.3 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  904.887    160.700   5.631 2.51e-08 ***
## age           22.544      3.697   6.099 1.68e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1168 on 777 degrees of freedom
## Multiple R-squared:  0.04568,    Adjusted R-squared:  0.04445 
## F-statistic: 37.19 on 1 and 777 DF,  p-value: 1.684e-09

Âge et salaire semblent ainsi positivement corrélés : une augmentation de l’âge de 1 an dans l’échantillon est associée en moyenne à une augmentation de salaire de 22,5 euros par mois. Sa p-valeur étant inférieure à 0,01, ce coefficient est significatif au seuil de 1 % (cf. cours). Graphiquement, il est possible de représenter la prédiction du modèle par une droite :

# Construction du nuage de point
plot(e$age,e$salred,xlab="Âge en années",ylab="Salaire en euros")

# Ajout d'une droite d'équation y =  coef(m1)[1] +  coef(m1)[2] * x
abline(a = coef(m1)[1], b = coef(m1)[2])

# Ajout d'un titre général
title(main="Représentation graphique de la prédiction du modèle m1")

Changement de spécification

Certaines valeurs extrêmes de salaire polarisent très fortement le nuage de points. Pour limiter cet effet, on peut modifier la spécification du modèle en régressant le salaire en logarithme : \[\text{(m2)} \quad \ln(salaire_i) = \beta_0 + \beta_1 \times age_i + \epsilon_i\]

# Création de la variable lnsalred
e$lnsalred <- log(e$salred)

# Estimation du modèle m2
m2 <- lm(lnsalred ~ age, data = e)

# Affichage des coefficients du modèle m2
summary(m2)$coefficients
##               Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept) 6.84842934 0.086897843 78.810119 0.000000e+00
## age         0.01186784 0.001998874  5.937264 4.367229e-09

Le coefficient de l’âge s’interprète désormais comme l’augmentation de salaire moyenne en pourcentage associé à 1 an supplémentaire dans l’échantillon, ici + 1,19 %.

Introduire une variable explicative et son carré

Il est également fréquent d’introduire certaines variables explicatives avec leur carré (tout particulièrement dans le cas de l’âge), ce qui permet de capter certaines relations non-linéaires. On estime alors le modèle : \[\text{(m3)} \quad \ln(salaire_i) = \beta_0 + \beta_1 \times age_i + \beta_2 \times age_i^2 + \epsilon_i\]

# Création de la variable d'âge au carré age2
e$age2 <- e$age^2

# Estimation du modèle m3
m3 <- lm(lnsalred ~ age + age2, data = e)

# Affichage des coefficients du modèle m3
summary(m3)$coefficients
##                  Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept)  5.4874355803 0.2605441145 21.061445 3.102072e-78
## age          0.0839229496 0.0131807866  6.367067 3.296151e-10
## age2        -0.0008801731 0.0001592134 -5.528259 4.419753e-08

\(\beta_1\) et \(\beta_2\) doivent alors être interprétés conjointement (cf. cours). Graphiquement, la prédiction du modèle m3 peut être représentée par une courbe :

# Construction du nuage de points
plot(e$age,e$lnsalred,xlab="Âge en années",ylab="Logarithme du salaire")

# Ajout de la courbe représentant la prédiction du modèle m3
curve(coef(m3)[1] + coef(m3)[2]*x + coef(m3)[3]*x^2, add=TRUE)

# Ajout d'un titre général
title(main="Représentation graphique de la prédiction du modèle m3")

Introduction de la variable de sexe et problèmes de colinéarité

De nombreux travaux d’économie et de sociologie du travail montrent l’existence d’un phénomène de discrimination salariale : à caractéristiques égales, être une femme est associé à un salaire en moyenne plus faible. La méthode de la régression linéaire multiple permet précisément de mettre en évidence le lien entre sexe et salaire en contrôlant par l’ensemble des autres variables disponibles dans les données.

Pour introduire la variable sexe dans la régression, on commence par la dichotomiser : on crée la variable homme valant 1 si l’individu est un homme et 0 sinon et la variable femme valant 1 si l’individu est une femme et 0 sinon. Ces variables permettent tout d’abord d’indiquer le sexe des individus sur le nuage de points du logarithme du salaire en fonction de l’âge :

# Dichotomisation de la variable codant le sexe
e$homme <- (e$sexe == 1) * 1
e$femme <- (e$sexe == 2) * 1

# Construction du nuage de points en distinguant les hommes et les femmes
plot(e$age,e$lnsalred,pch = e$femme,xlab="Âge en années",ylab="Logarithme du salaire")

# Ajout d'une légende
legend("topleft",legend=c("Hommes","Femmes"),pch=c(0,1))

# Ajout d'un titre général
title(main="Logarithme du salaire, âge et sexe")

Il semble qu’il y ait davantage de femmes dans les valeurs extrêmes faibles et d’hommes dans les valeurs extrêmes élevées. Pour confirmer cette impression, on estime le modèle m4 : \[\text{(m4)} \quad \ln(salaire_i) = \beta_0 + \beta_1 age_i + \beta_2 age_i^2 + \beta_3 femme_i + \beta_4 homme_i + \epsilon_i\]

# Estimation du modèle m4
m4 <- lm(lnsalred ~ age + age2 + femme + homme, data = e)

# Affichage des coefficients du modèle m4
summary(m4)
## 
## Call:
## lm(formula = lnsalred ~ age + age2 + femme + homme, data = e)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2079 -0.2628  0.0385  0.3483  1.9731 
## 
## Coefficients: (1 not defined because of singularities)
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.7038824  0.2501340  22.803  < 2e-16 ***
## age          0.0828141  0.0125918   6.577 8.83e-11 ***
## age2        -0.0008742  0.0001521  -5.748 1.30e-08 ***
## femme       -0.3686416  0.0424589  -8.682  < 2e-16 ***
## homme               NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.592 on 775 degrees of freedom
## Multiple R-squared:  0.1612, Adjusted R-squared:  0.158 
## F-statistic: 49.66 on 3 and 775 DF,  p-value: < 2.2e-16

La ligne correspondant à la variable homme ne comporte que des NA (pour Not Assigned). R a en effet automatiquement exclu cette variable du modèle en raison du phénomène de colinéarité (cf. cours). De manière générale, pour introduire une variable qualitative il convient de la dichotomiser et d’introduire toutes les variables indicatrices sauf une (cf. infra avec le diplôme sur le choix de cette modalité). Le modèle m4 estimé est donc en fait : \[\text{(m4)} \quad \ln(salaire_i) = \beta_0 + \beta_1 age_i + \beta_2 age_i^2 + \beta_3 femme_i + \epsilon_i\]

# Ré-estimation du modèle m4 (sans la colinéarité)
m4 <- lm(lnsalred ~ age + age2 + femme, data = e)

# Affichage des coefficients du modèle m4
summary(m4)$coefficients
##                  Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept)  5.7038824061 0.2501339896 22.803308 1.801019e-88
## age          0.0828140895 0.0125917871  6.576834 8.834572e-11
## age2        -0.0008741721 0.0001520925 -5.747633 1.301026e-08
## femme       -0.3686415575 0.0424589295 -8.682309 2.273092e-17

Dans le cas des variables dichotomiques comme le sexe, l’interprétation du coefficient est simple : celui-ci correspond à l’écart moyen associé au fait d’avoir la caractéristique codée par la variable dichotomique. Dans le cas du modèle m4, le fait d’être une femme est associé, à âge égal, à un salaire en moyenne de 36,9 % inférieur. La prédiction de ce modèle peut être représentée par deux courbes :

# Construction du nuage de points
plot(e$age,e$lnsalred,pch = e$femme,xlab="Âge en années",ylab="Logarithme du salaire")

# Ajout des courbes représentant la prédiction du modèle m4
curve(coef(m4)[1] + coef(m4)[2]*x + coef(m4)[3]*x^2, add=TRUE)
curve(coef(m4)[1] + coef(m4)[4] + coef(m4)[2]*x + coef(m4)[3]*x^2, lty=2, add=TRUE)

# Ajout d'une légende mixte
legend("topleft",legend=c("Hommes","Femmes"),pch=c(0,1),lty=c(1,2))

# Ajout d'un titre général
title(main="Représentation graphique de la prédiction du modèle m4")

Introduction de la variable de diplôme et choix de la modalité de référence

Pour approfondir la question de la discrimination salariale, on peut chercher à neutraliser l’effet du diplôme sur le salaire. Pour ce faire, il suffit de rajouter les modalités de la variable ddipl dans le modèle de régression (sauf la dernière) : \[\text{(m5)} \quad \ln(salaire_i) = \beta_0 + \beta_1 age_i + \beta_2 age_i^2 + \beta_3 femme_i + \beta_4 ddipl1_i + \beta_5 ddipl3_i + ... + \beta_8 ddipl6_i +\epsilon_i\]

# Dichotomisation de la variable de niveau de diplôme
e$ddipl1 <- (e$ddipl == 1) * 1 # Diplôme supérieur à baccalauréat + 2 ans
e$ddipl3 <- (e$ddipl == 3) * 1 # Baccalauréat + 2 ans
e$ddipl4 <- (e$ddipl == 4) * 1 # Baccalauréat ou brevet professionnel ou autre diplôme de ce niveau
e$ddipl5 <- (e$ddipl == 5) * 1 # CAP, BEP ou autre diplôme de ce niveau
e$ddipl6 <- (e$ddipl == 6) * 1 # Brevet des collèges
e$ddipl7 <- (e$ddipl == 7) * 1 # Aucun diplôme ou certificat d'études primaires

# Estimation du modèle m5
m5 <- lm(lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + ddipl4 + ddipl5 + ddipl6, data = e)

# Résumé du modèle m5
summary(m5)
## 
## Call:
## lm(formula = lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + 
##     ddipl4 + ddipl5 + ddipl6, data = e)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2307 -0.2146  0.0671  0.2927  1.9488 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.4264080  0.2330520  23.284  < 2e-16 ***
## age          0.0691972  0.0116498   5.940 4.32e-09 ***
## age2        -0.0006555  0.0001412  -4.641 4.07e-06 ***
## femme       -0.3961442  0.0391407 -10.121  < 2e-16 ***
## ddipl1       0.7714419  0.0699169  11.034  < 2e-16 ***
## ddipl3       0.6990360  0.0730079   9.575  < 2e-16 ***
## ddipl4       0.5127076  0.0667608   7.680 4.84e-14 ***
## ddipl5       0.3052744  0.0644052   4.740 2.55e-06 ***
## ddipl6       0.3301545  0.0994454   3.320 0.000943 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.54 on 770 degrees of freedom
## Multiple R-squared:  0.3067, Adjusted R-squared:  0.2995 
## F-statistic: 42.58 on 8 and 770 DF,  p-value: < 2.2e-16

Avant toute chose, on constate que même en intégrant le niveau de diplôme le coefficient associé au fait d’être une femme reste significatif et du même ordre de grandeur que dans le modèle précédent (en moyenne un salaire inférieur de 39,6 %)

Dans le cas d’une variable catégorielle à plusieurs modalités, les coefficients s’interprètent relativement à la modalité de différence. Ainsi dans le modèle m6, le fait d’avoir le brevet des collèges (ddipl6) par rapport à aucun diplôme ou le certificat d’étude primaire (ddipl7, la référence) est associé, à sexe et âge égaux, à un salaire en moyenne 33,0 % supérieur.

Le choix de la modalité de référence est donc crucial pourl’interprétation des résultats, dans la mesure où il affecte la valeur des coefficients ainsi que les tests de significativité (cf. cours). Les modèles m6b à m6f sont des variantes du modèle m6 où la modalité de référence pour le diplôme est modifiée :

# Estimation des modèles m5b à m5f
m5b <- lm(lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + ddipl4 + ddipl5 + ddipl7, data = e)
m5c <- lm(lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + ddipl4 + ddipl6 + ddipl7, data = e)
m5d <- lm(lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + ddipl5 + ddipl6 + ddipl7, data = e)
m5e <- lm(lnsalred ~ age + age2 + femme + ddipl1 + ddipl4 + ddipl5 + ddipl6 + ddipl7, data = e)
m5f <- lm(lnsalred ~ age + age2 + femme + ddipl3 + ddipl4 + ddipl5 + ddipl6 + ddipl7, data = e)

# Affichage des coefficients des modèles m5b, m5d et m5f
summary(m5b)$coefficients
##                 Estimate   Std. Error    t value     Pr(>|t|)
## (Intercept)  5.756562484 0.2392113212  24.064758 7.444386e-96
## age          0.069197201 0.0116497884   5.939782 4.319187e-09
## age2        -0.000655489 0.0001412384  -4.641012 4.074485e-06
## femme       -0.396144237 0.0391406976 -10.121032 1.086937e-22
## ddipl1       0.441287341 0.0979900521   4.503389 7.723714e-06
## ddipl3       0.368881430 0.1001984990   3.681507 2.480444e-04
## ddipl4       0.182553033 0.0955541404   1.910467 5.644427e-02
## ddipl5      -0.024880171 0.0942519921  -0.263975 7.918698e-01
## ddipl7      -0.330154527 0.0994454102  -3.319957 9.426885e-04
summary(m5d)$coefficients
##                 Estimate   Std. Error    t value      Pr(>|t|)
## (Intercept)  5.939115517 0.2305332000  25.762517 4.861661e-106
## age          0.069197201 0.0116497884   5.939782  4.319187e-09
## age2        -0.000655489 0.0001412384  -4.641012  4.074485e-06
## femme       -0.396144237 0.0391406976 -10.121032  1.086937e-22
## ddipl1       0.258734308 0.0621396268   4.163757  3.484282e-05
## ddipl3       0.186328398 0.0652342668   2.856296  4.401572e-03
## ddipl5      -0.207433204 0.0568835020  -3.646632  2.835791e-04
## ddipl6      -0.182553033 0.0955541404  -1.910467  5.644427e-02
## ddipl7      -0.512707559 0.0667607629  -7.679774  4.842889e-14
summary(m5f)$coefficients
##                 Estimate   Std. Error    t value      Pr(>|t|)
## (Intercept)  6.197849825 0.2360574201  26.255687 5.212520e-109
## age          0.069197201 0.0116497884   5.939782  4.319187e-09
## age2        -0.000655489 0.0001412384  -4.641012  4.074485e-06
## femme       -0.396144237 0.0391406976 -10.121032  1.086937e-22
## ddipl3      -0.072405911 0.0681599288  -1.062294  2.884351e-01
## ddipl4      -0.258734308 0.0621396268  -4.163757  3.484282e-05
## ddipl5      -0.466167512 0.0598961190  -7.782934  2.282870e-14
## ddipl6      -0.441287341 0.0979900521  -4.503389  7.723714e-06
## ddipl7      -0.771441868 0.0699169027 -11.033696  2.209062e-26

La valeur (et notamment le signe) des paramètres estimés relatifs au diplôme change du tout au tout selon la modalité de référence choisie, de même pour les tests de significativité : en règle générale, quand la variable est ordonnée on préfère choisir comme modalité de référence une modalité plutôt centrale (cf. cours).

Comparaison des modèles : R2, R2 ajusté et test de significativité globale des coefficients de Fisher

Pouvoir prédictif du modèle : R2 et R2 ajusté

Un modèle de régression est d’autant meilleur qu’il parvient à minimiser la somme des carrés des résidus \(\sum_{i=1}^n \epsilon_i\). Le R2 est un indicateur construit sur ce principe, qui permet de rendre compte de la qualité de la prédiction produite par le modèle.

Cependant, quand on ajoute des variables explicatives mécaniquement le R2 augmente : le R2 ajusté permet de prendre en compte ce phénomène en privilégiant, à qualité de prédiction égale, les modèles les plus parcimonieux en variables explicatives.

La fonction summary() affiche par défaut la valeur du R2, du R2 ajusté. Pour comparer ces statistiques dans les 10 premiers modèles, on les centralise dans un tableau en utilisant la fonction sapply() (?sapply pour plus de détails).

liste <- list(m1,m2,m3,m4,m5,m5b,m5c,m5d,m5e,m5f)
data.frame(
  modele=c("m1","m2","m3","m4","m5","m5b","m5c","m5d","m5e","m5f")
  , R2=sapply(liste,function(x){summary(x)$r.square})
  , R2adj=sapply(liste,function(x){summary(x)$adj.r.square})
)
##    modele         R2      R2adj
## 1      m1 0.04568262 0.04445442
## 2      m2 0.04339927 0.04216812
## 3      m3 0.07964610 0.07727405
## 4      m4 0.16123123 0.15798439
## 5      m5 0.30670561 0.29950255
## 6     m5b 0.30670561 0.29950255
## 7     m5c 0.30670561 0.29950255
## 8     m5d 0.30670561 0.29950255
## 9     m5e 0.30670561 0.29950255
## 10    m5f 0.30670561 0.29950255

Par construction, le R2 ajusté est toujours inférieur au R2 (mais d’assez peu ici car il y a peu de variables), et au fur et à mesure qu’on ajoute des variables explicatives, les deux statistiques augmentent. Avec les variables d’âge, de sexe et de niveau de diplôme, le modèle parvient à expliquer 30 % de la variance du salaire en logarithme, ce qui est acceptable étant donnée sa parcimonie. Au passage, on observe qu’un changement dans la modalité de référence n’induit aucun changement dans la qualité de prédiction globale du modèle, d’où les valeurs identiques pour les modèles m5 à m5f.

Pouvoir explicatif du modèle : test de significativité globale de Fisher

R2 et R2 ajusté constituent de bons indicateurs du pouvoir prédictif d’un modèle (sur les différentes finalités de la modélisation économétrique, cf. la conclusion de la séance 1). Pour juger de leur pouvoir explicatif, il est possible de tester conjointement la significativité de toutes les variables introduites dans le modèle. Si l’on peut rejeter à un seuil raisonnable l’hypothèse qu’aucune des variables n’est significative, alors cela signifie que le modèle présente un certain pouvoir explicatif.

Dès le modèle m1, la p-valeur du test de Fisher est inférieure à 0,001 : avec la seule variable d’âge l’association avec la variable de salaire est suffisante pour qu’on puisse estimer avec un risque très faible de se tromper (très inférieur à 0, %) que le modèle apporte des éléments d’explication de la ditribution des salaires.

Introduction des autres variables : utilisation de variables croisées et test d’hypothèses complexes

Les deux autres variables que l’on souhaite intégrer dans le modèle sont le fait de travailler à temps partiel et le secteur d’emploi (public ou privé).

Influence du travail à temps partiel

On estime tout d’abord le modèle m6 (en reprenant toutes les variables du modèle m5d) : \[ \text{(m6)} \quad \ln(salaire_i) = \beta_0 + \beta_1 age_i ... + \beta_8 ddipl7_i + \beta_9 tpp_i +\epsilon_i\]

# Recodage de tpp sous forme d'indicatrice
e$tpp <- (e$tpp == 2) * 1

# Estimation du modèle m6
m6 <- lm(lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + ddipl5 + ddipl6 + ddipl7 + tpp, data = e)

# Résumé du modèle m6
summary(m6)
## 
## Call:
## lm(formula = lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + 
##     ddipl5 + ddipl6 + ddipl7 + tpp, data = e)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.71034 -0.21731  0.01168  0.24895  1.91424 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.0577416  0.2002346  30.253  < 2e-16 ***
## age          0.0641887  0.0101166   6.345 3.80e-10 ***
## age2        -0.0005827  0.0001227  -4.750 2.43e-06 ***
## femme       -0.2300553  0.0355410  -6.473 1.71e-10 ***
## ddipl1       0.2717526  0.0539415   5.038 5.87e-07 ***
## ddipl3       0.1348533  0.0567137   2.378   0.0177 *  
## ddipl5      -0.2059405  0.0493732  -4.171 3.38e-05 ***
## ddipl6      -0.1351951  0.0829914  -1.629   0.1037    
## ddipl7      -0.4712866  0.0580047  -8.125 1.78e-15 ***
## tpp         -0.7340425  0.0461421 -15.908  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4687 on 769 degrees of freedom
## Multiple R-squared:  0.4784, Adjusted R-squared:  0.4723 
## F-statistic: 78.36 on 9 and 769 DF,  p-value: < 2.2e-16

Le R2 de la régression augmente très fortement (à 0,4784), ce qui indique que la variable de temps partiel est déterminante pour comprendre les écarts de salaire dans notre échantillon (ce qui est assez logique). L’amplitude de l’écart de salaire entre hommes et femmes est plus faible que précédemment (-23,0 %) : cela est dû au fait que les femmes travaillent bien plus souvent à temps partiel que les hommes (29,0 % contre 7,1 %), et que le temps partiel est associé à un salaire en moyenne beaucoup plus faible (1 006 euros contre 2 035 euros pour un temps complet). Il s’agit d’un effet de structure.

# Temps partiel selon le sexe
data.frame(
  sexe=c("Homme","Femme")
  ,tpp=as.numeric(by(e$tpp,e$femme,mean))
)
##    sexe        tpp
## 1 Homme 0.07070707
## 2 Femme 0.28981723
# Salaire et temps partiel
data.frame(
  tpp=c("Temps complet","Temps partiel")
  ,salred=as.numeric(by(e$salred,e$tpp,mean))
)
##             tpp   salred
## 1 Temps complet 2034.777
## 2 Temps partiel 1005.727

Influence du secteur d’activité

On peut formuler l’hypothèse que les écarts de salaire entre hommes et femmes sont plus faibles dans le secteur public que dans le secteur privé, dans la mesure où l’évolution salariale y est fortement déterminée par des grilles communes.

On peut commencer par estimer le modèle : \[ \text{(m7)} \quad \ln(salaire_i) = \beta_0 + \beta_1 age_i ... + \beta_9 tpp_i + \beta_{10} public_i +\epsilon_i\]

# Création de la variable public à partir de pub3fp
e$public <- (e$pub3fp %in% c('1','2','3')) * 1

# Estimation du modèle m7
m7 <- lm(lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + ddipl5 + ddipl6 + ddipl7 + tpp + public, data = e)

# Résumé du modèle m7
summary(m7)
## 
## Call:
## lm(formula = lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + 
##     ddipl5 + ddipl6 + ddipl7 + tpp + public, data = e)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.69361 -0.21361  0.01218  0.24854  1.92101 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.0492156  0.2003470  30.194  < 2e-16 ***
## age          0.0643376  0.0101158   6.360 3.46e-10 ***
## age2        -0.0005858  0.0001227  -4.775 2.16e-06 ***
## femme       -0.2368066  0.0360443  -6.570 9.28e-11 ***
## ddipl1       0.2670390  0.0540971   4.936 9.77e-07 ***
## ddipl3       0.1345860  0.0567050   2.373   0.0179 *  
## ddipl5      -0.2020448  0.0494879  -4.083 4.92e-05 ***
## ddipl6      -0.1284620  0.0831960  -1.544   0.1230    
## ddipl7      -0.4632515  0.0584387  -7.927 7.89e-15 ***
## tpp         -0.7372579  0.0462241 -15.950  < 2e-16 ***
## public       0.0476773  0.0426322   1.118   0.2638    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4686 on 768 degrees of freedom
## Multiple R-squared:  0.4792, Adjusted R-squared:  0.4724 
## F-statistic: 70.67 on 10 and 768 DF,  p-value: < 2.2e-16

Le coefficient associé à la variable public est faible et non significatif. Cependant, l’interprétation de ce seul coefficient ne permet pas d’apporter une réponse à notre hypothèse : il rend en effet compte de l’association moyenne entre secteur et salaire, sans qu’il soit possible que cette association diffère entre hommes et femmes. En d’autres termes, pour répondre à la question posée on souhaiterait croiser les variables sexe et secteur en estimant le modèle : \[ \text{(m8)} \quad \ln(salaire_i) = \beta_0 + \beta_1 age_i ... + \beta_8 tpp_i + \beta_9 femme\_prive_i + \beta_{10} homme\_public_i + \beta_{11} femme\_public _i +\epsilon_i\]

# Création des variables croisées
e$homme_prive <- e$homme*(1-e$public)
e$femme_prive <- e$femme*(1-e$public)
e$homme_public <- e$homme*e$public
e$femme_public <- e$femme*e$public

# Estimation du modèle m8
m8 <- lm(lnsalred ~ age + age2 + ddipl1 + ddipl3 + ddipl5 + ddipl6 + ddipl7 + tpp + femme_prive + homme_public + femme_public, data = e)

# Résumé du modèle m8
summary(m8)
## 
## Call:
## lm(formula = lnsalred ~ age + age2 + ddipl1 + ddipl3 + ddipl5 + 
##     ddipl6 + ddipl7 + tpp + femme_prive + homme_public + femme_public, 
##     data = e)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.65817 -0.21630  0.02006  0.25250  1.90068 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.0549018  0.1992901  30.382  < 2e-16 ***
## age           0.0648655  0.0100635   6.446 2.03e-10 ***
## age2         -0.0005866  0.0001220  -4.807 1.85e-06 ***
## ddipl1        0.2712860  0.0538275   5.040 5.81e-07 ***
## ddipl3        0.1342370  0.0564034   2.380  0.01756 *  
## ddipl5       -0.2097652  0.0492902  -4.256 2.34e-05 ***
## ddipl6       -0.1356673  0.0827874  -1.639  0.10168    
## ddipl7       -0.4736262  0.0582279  -8.134 1.67e-15 ***
## tpp          -0.7373260  0.0459782 -16.036  < 2e-16 ***
## femme_prive  -0.2881022  0.0396268  -7.270 8.85e-13 ***
## homme_public -0.1213313  0.0699339  -1.735  0.08315 .  
## femme_public -0.1463200  0.0526571  -2.779  0.00559 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4661 on 767 degrees of freedom
## Multiple R-squared:  0.4854, Adjusted R-squared:  0.478 
## F-statistic: 65.77 on 11 and 767 DF,  p-value: < 2.2e-16

Pour les salariés du privé, à âge et diplôme égaux et en contrôlant par le fait d’être à temps partiel ou non, le fait d’être une femme est associé à un salaire moyen inférieur de 28,8 %. Dans le secteur public, cet écart est de l’ordre de 2,5 %.

Test d’hypothèses complexes avec des variables croisées qualitatives

Quand des variables croisées interviennent dans un modèle, il est pertinent de mener d’autres tests que le seul test de significativité :

  • \(\beta_9\) = \(\beta_{11}\) : le salaire des femmes est-il significativement différent selon le secteur ?
  • \(\beta_{10}\) = \(\beta_{11}\) : le salaire des hommes et des femmes est-il significativement différent dans le secteur public ?

Pour mener à bien ces tests d’hypothèse, il suffit de modifier la modalité de référence (modèles m8b et m8c).

# Estimation du modèle m8b avec femme_prive comme référence (pour tester l'écart avec femme_public)
m8b <- lm(lnsalred ~ age + age2 + ddipl1 + ddipl3 + ddipl5 + ddipl6 + ddipl7 + tpp + homme_prive + homme_public + femme_public, data = e)

# Affichage des coefficients pertinents du modèle m8b
summary(m8b)$coefficients[10:12,]
##               Estimate Std. Error  t value     Pr(>|t|)
## homme_prive  0.2881022 0.03962682 7.270384 8.851702e-13
## homme_public 0.1667709 0.07136585 2.336845 1.970377e-02
## femme_public 0.1417821 0.05250716 2.700244 7.081603e-03

Toutes les autres variables du modèle contrôlées par ailleurs, en moyenne le salaire des femmes dans le secteur public est supérieur de 14,2 % à celui des femmes dans le secteur privé (femme_prive est la référence dans le modèle m8b). La p-valeur du test de significativité de femme_public étant inférieure à 0.01 (égale à 0,007), cet écart est significatif au seuil de 1 %.

# Estimation du modèle m8c avec homme_public comme référence (pour tester l'écart avec femme_public)
m8c <- lm(lnsalred ~ age + age2 + ddipl1 + ddipl3 + ddipl5 + ddipl6 + ddipl7 + tpp + homme_prive + femme_prive  + femme_public, data = e)

# Affichage des coefficients pertinents du modèle m8c
summary(m8c)$coefficients[10:12,]
##                 Estimate Std. Error    t value   Pr(>|t|)
## homme_prive   0.12133127 0.06993388  1.7349428 0.08315264
## femme_prive  -0.16677090 0.07136585 -2.3368447 0.01970377
## femme_public -0.02498877 0.07837733 -0.3188266 0.74994484

Dans le modèle 8c où homme_public est la modalité de référence, le coefficient associé à femme_public n’est pas significativement différent de 0 (p-valeur de 0,750) : toutes les variables du modèle étant contrôlées par ailleurs, on ne peut pas distinguer d’écart de salaire significatif entre hommes et femmes dans le secteur public.

Croisement de l’âge et du sexe

Croiser des variables explicatives dans un modèle permet souvent d’apporter des éléments de réponse à des hypothèses de travail très concrètes. Au-delà des seules variables qualitatives, il est également possible de croiser une variable qualitative et une variable quantitative. Afin de pouvoir représenter graphiquement l’impact du croisement de ces variables sur la modélisation, on repart du modèle m4 : \[\text{(m4)} \quad \ln(salaire_i) = \beta_0 + \beta_1 age_i + \beta_2 age_i^2 + \beta_3 femme_i + \epsilon_i\]

Pour croiser sexe et âge, il suffit d’autoriser les coefficients relatifs à l’âge à varier en fonction du sexe : \[\text{(m9)} \quad \ln(salaire_i) = \beta_0 + \beta_1 femme_i + homme_i \times (\beta_2 age_i + \beta_3 age_i^2) + femme_i \times (\beta_4 age_i + \beta_5 age_i^2) + \epsilon_i\]

En pratique, on définit la variable age_homme qui vaut l’âge pour les hommes et 0 pour les femmes et la variable age_femme qui faut l’âge pour les femmes et 0 pour les hommes, ainsi que leur carré respectif (age2_homme et age2_femme).

# Construction des variables de croisement
e$age_homme <- e$age*e$homme
e$age_femme <- e$age*e$femme
e$age2_homme <- e$age_homme^2
e$age2_femme <- e$age_femme^2

# Estimation du modèle m9
m9 <- lm(lnsalred ~ femme + age_homme + age2_homme + age_femme + age2_femme, data = e)

# Résumé du modèle m9
summary(m9)
## 
## Call:
## lm(formula = lnsalred ~ femme + age_homme + age2_homme + age_femme + 
##     age2_femme, data = e)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3030 -0.2592  0.0454  0.3369  2.0352 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.5364945  0.3454229  16.028  < 2e-16 ***
## femme       -0.0368156  0.4967968  -0.074    0.941    
## age_homme    0.0858810  0.0174282   4.928 1.02e-06 ***
## age2_homme  -0.0008548  0.0002105  -4.061 5.38e-05 ***
## age_femme    0.0797460  0.0181087   4.404 1.21e-05 ***
## age2_femme  -0.0008941  0.0002188  -4.087 4.83e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5904 on 773 degrees of freedom
## Multiple R-squared:  0.168,  Adjusted R-squared:  0.1626 
## F-statistic: 31.21 on 5 and 773 DF,  p-value: < 2.2e-16

On voit que les coefficients de age_homme et age_femme diffèrent sensiblement, avec age_homme > age_femme : cela semble indiquer qu’en moyenne dans l’échantillon, le gain de salaire associé à 1 an d’âge supplémentaire est supérieur pour les hommes.

# Construction du nuage de points
plot(e$age,e$lnsalred,pch = e$femme,xlab="Âge en années",ylab="Logarithme du salaire")

# Ajout des courbes représentant la prédiction du modèle m9
curve(coef(m9)[1] + coef(m9)[3]*x + coef(m9)[4]*x^2, add=TRUE)
curve(coef(m9)[1] + coef(m9)[2] + coef(m9)[5]*x + coef(m9)[6]*x^2, lty=2, add=TRUE)

# Ajout d'une légende mixte
legend("topleft",legend=c("Hommes","Femmes"),pch=c(0,1),lty=c(1,2))

# Ajout d'un titre général
title(main="Représentation graphique de la prédiction du modèle m9")

Le croisement de l’âge et du sexe permet ainsi de mettre en évidence un écart de salaire entre hommes et femmes qui se creuse avec l’âge.

Test d’hypothèses complexes avec des variables croisées qualitatives et quantitatives

Il est également possible de mener à bien des tests d’hypothèsees complexes quand les variables croisées sont qualitatives et quantitatives. Dans l’exemple, cela revient à tester si les coefficients associés à l’âge sont significativement différents pour les hommes et les femmes. Une analyse de la variance permet de construire le test approprié : si le logarithme du salaire est bien mieux expliqué par le modèle m9 que par le modèle m4, alors ces coefficients sont significativement différents les uns des autres.

# Analyse de la variance des modèles m4 et m9
anova(m4,m9)
## Analysis of Variance Table
## 
## Model 1: lnsalred ~ age + age2 + femme
## Model 2: lnsalred ~ femme + age_homme + age2_homme + age_femme + age2_femme
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1    775 271.62                              
## 2    773 269.43  2    2.1871 3.1374 0.04395 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La p-valeur du test de Fisher est de 0,044 < 0,050 : on peut donc rejeter l’hypothèse nulle d’égalité des coefficients \(\beta_2\) et \(\beta_4\) d’une part, \(\beta_3\) et \(\beta_5\) d’autre part (en reprenant les notations du modèle m9).

 

Exercice : Adhésion à un programme de retraite d’entreprise

Dans un article de 1995, l’économiste L.E. Papke propose une analyse empirique de l’adhésion des salariés américains à un certain type de programmes de retraite d’entreprise (programmes 401(k)) en fonction de la générosité de l’employeur (au sens de la part que représente sa contribution au programme). Le fichier http://teaching.slmc.fr/mqs2/401k.dta comporte une partie de l’échantillon d’entreprises utilisé dans cette étude ; il est décrit par le fichier http://teaching.slmc.fr/mqs2/401k.des. En particulier:

  • prate est le pourcentage de salariés de l’entreprise adhérant à un programme de retraite 401(k)
  • mrate est la part que représente l’abondement de l’employeur. Si mrate = 0.50, alors pour 1 dollar épargné l’employeur abonde à hauteur de 50 cents.

L’objectif de cet exercice, repris de l’ouvrage de J.M. Wooldridge cité en bibliographie (cf. introduction du cours), est d’étudier la relation entre prateet mrate.

 

Q1 Utilisez la fonction library() pour charger le package foreign. Utilisez sa fonction read.dta() pour charger les données du fichier http://teaching.slmc.fr/mqs2/401k.dta dans R. Affichez la structure des données importées (nombre de variables, d’observations, etc.)

 

Q2 Calculez le pourcentage moyen de salariés couverts par un programme de retraite 401(k) ainsi que la moyenne de l’abondement des employeurs.

 

Q3 Utilisez la fonction lm() pour estimer le modèle \[\text{(m1)} \quad prate = \beta_0 + \beta_1mrate + \epsilon\] Interprétez la valeur des estimations des coefficients \(\beta_0\) et \(\beta_1\).

 

Q4 Représentez la relation entre prate et mrate par un nuage de points ainsi que la prédiction du modèle m1 (pensez à vous inspirer des graphiques produits lors de la séance 1).

 

Q5 À partir de l’analyse du graphique et des statistiques d’ajustement du modèle (R2, R2 ajusté), quelle opinion formez-vous sur la qualité de cette modélisation ?

 

Q6 Introduisez le carré de l’abondement de l’employeur dans la régression pour estimer le modèle : \[\text{(m2)} \quad prate = \beta_0 + \beta_1 mrate + \beta_2 mrate^2 + \epsilon\] Interprétez les coefficients et représentez la prédiction du modèle sur le nuage de points. Que pensez-vous de cette nouvelle modélisation ?

 

Séance 3 : Régression sur données qualitatives

# Lecture des données et assignation à l'objet e
e <- readRDS(gzcon(url("http://teaching.slmc.fr/mqs2/ee12t4_new.rds")))

# Passage en minuscule des noms de variables et recodage de la variable d'âge
colnames(e) <- tolower(colnames(e))
e$age <- as.numeric(e$age)

# Création d'une variable de stabilité du contrat de travail
# et restriction de l'échantillon aux individus pour lesquels
# elle est pertinente
e$stable[e$chpub == "1" & e$titc %in% c("1","2")] <- 1
e$stable[e$contra == "1"] <- 1
e$stable[e$contra %in% as.character(2:5)] <- 0
e <- e[!is.na(e$stable),]

# Définition de la fonction logit()
logit <- function(x) exp(x)/(1 + exp(x))

Nature de la variable dépendante et choix du modèle

# Représentation de la relation entre âge et stabilité du contrat
plot(e$age,e$stable,xlab = "Âge en années", ylab = 
"Contrat de travail stable ou non")
title(main = "Relation entre âge et stabilité du contrat")

# Installation, chargement et utilisation du package data.table pour agréger les informations efficacement
# install.packages("data.table")
library(data.table)
ag <-data.table(e)[,j=list(count=.N),by=list(stable,age)]

# Représentation de la relation entre âge et travail
plot(e$age,e$stable,type="n",xlab = "Âge en années", ylab = 
"Contrat de travail stable ou non")
symbols(x=ag$age, y=ag$stable, circles=ag$count, inches=1/4, ann=F, bg="cornflowerblue", fg=NULL, add = T)
title(main = "Relation entre âge et stabilité du contrat")

# Représentation de la relation entre âge et travail (troisième essai)
plot(e$age,e$stable,type="n",xlab = "Âge en années", ylab = 
"Contrat de travail stable ou non")
symbols(x=ag$age, y=ag$stable, circles=sqrt(ag$count), inches=1/4, ann=F, bg="cornflowerblue", fg=NULL, add = T)

# Sauvegarde du graph pour la suite (sans son titre)
g1 <- recordPlot()

# Ajout du titre
title(main = "Relation entre âge et stabilité du contrat")

Modéliser des données qualitatives

# Génération de données aléatoires
set.seed(1)
x <- rnorm(100,0,1)
x2 <- x^2
y <- 2*x + rnorm(100,0,1)
z <- ( exp(y)/(1+exp(y)) > 0.5 ) * 1
plot(x,z,pch=19,col="cornflowerblue")

# Estimation du modèle sim1 (linéaire) et représentation de sa prédiction
plot(x,z,pch=19,col="cornflowerblue")
sim1 <- glm(z ~ x,family=gaussian)
abline(a = coef(sim1)[1], b = coef(sim1)[2])
title(main = "Modèle linéaire appliqué à des données simulées")

# Estimation du modèle sim2 (linéaire général logit) et représentation de sa prédiction
plot(x,z,pch=19,col="cornflowerblue")
sim2 <- glm(z ~ x,family=binomial(link="logit"))
curve(logit(coef(sim2)[1] + coef(sim2)[2]*x), add=TRUE)
title(main = "Modèle logistique appliqué à des données simulées")

# Estimation du modèle sim3 (linéaire général probit) et représentation de sa prédiction
plot(x,z,pch=19,col="cornflowerblue")
sim3 <- glm(z ~ x,family=binomial(link="probit"))
curve(pnorm(coef(sim3)[1] + coef(sim3)[2]*x), add=TRUE)
title(main = "Modèle probit appliqué à des données simulées")

# Représentation des 3 modélisation sur le même graph
plot(x,z,pch=19,col="cornflowerblue")
abline(a = coef(sim1)[1], b = coef(sim1)[2],lty=1)
curve(logit(coef(sim2)[1] + coef(sim2)[2]*x), add=TRUE,lty=2)
curve(pnorm(coef(sim3)[1] + coef(sim3)[2]*x), add=TRUE,lty=3)
legend("right",c("Modèle linéaire","Modèle logit","Modèle probit"),lty=1:3)
title(main = "Comparaison des modèles sur données simulées")

# Estimation du modèle 1 (linéaire) et représentation de sa prédiction
m1 <- glm(stable~age,data=e,family=gaussian)
g1
abline(a = coef(m1)[1],b = coef(m1)[2],lty=1)
title(main = "Modèle linéaire simple appliqué aux données de l'EEC")

# Estimation du modèle 2 (linéaire général logit) et représentation de sa prédiction
m2 <- glm(stable ~ age,data=e,family=binomial(link="logit"))
g1
curve(logit(coef(m2)[1] + coef(m2)[2]*x), add=TRUE)
title(main = "Modèle logit simple appliqué aux données de l'EEC")

# Estimation du modèle 3 (linéaire général probit) et représentation de sa prédiction
m3 <- glm(stable ~ age,data=e,family=binomial(link="probit"))
g1
curve(pnorm(coef(m2)[1] + coef(m2)[2]*x), add=TRUE)
title(main = "Modèle probit simple appliqué aux données de l'EEC")

# Ajout d'une légende et d'un titre
g1
abline(a = coef(m1)[1],b = coef(m1)[2],lty=1)
curve(logit(coef(m2)[1] + coef(m2)[2]*x), add=TRUE,lty=2)
curve(pnorm(coef(m2)[1] + coef(m2)[2]*x), add=TRUE,lty=3)
legend("right",c("Modèle linéaire","Modèle logit","Modèle probit"),lty=1:3)
title(main = "Comparaison des modèles sur les données de l'EEC")

Le modèle linéaire général et son estimation

# Fonction logistique et fonction Phi
plot(x=1,type="n",xlim=c(-6,6),ylim=c(0,1),xlab="",ylab="")
curve(logit(x), add=TRUE,lty=2)
curve(pnorm(x), add=TRUE,lty=3)
legend("right",c("Fonction logit","Fonction Phi"),lty=2:3)
title(main = "Fonction logistique et fonction Phi")

# Option graphique pour avoir les deux graphiques côte-à-côte
# par(mfrow=c(1,2))

# Représentation de la prédiction d'une modélisation vraisemblable de z
col1 <- rep("cornflowerblue",100)
col1[c(59,6)] <- "red"
col1[c(84,95)] <- "green"
plot(x,z,pch=19,col=col1)
curve(logit(coef(sim2)[1] + coef(sim2)[2]*x), add=TRUE)
text(x = 0.6, y = 0.5, "logit(b0A + b1A*x)")
title(main="Modélisation vraisemblable de z")

# Représentation de la prédiction d'une modélisation peu vraisemblable de z
col2 <- rep("cornflowerblue",100)
col2[(x < -1 & z == 0) | (x > 1 & z == 1)] <- "red"
plot(x,z,pch=19,col=col2)
curve(logit(-1.3*x), add=TRUE)
text(x = 0.8, y = 0.5, "logit(b0B + b1B*x)")
title(main="Modélisation peu vraisemblable de z")

# Réinitialisation de l'option graphique pour avoir les deux graphiques côte-à-côte et ajout du titre
# par(mfrow=c(1,1))
# title(main="Modélisation vraisemblable (à gauche) et peu vraisemblable \n (à droite) de z")

Interprétation des coefficients d’un modèle logistique

# Construction des variables explicatives supplémentaires
e <- within(e[!is.na(e$ddipl),],{

  ddipl1 <- (ddipl == "1") * 1
  ddipl3 <- (ddipl == "3") * 1
  ddipl4 <- (ddipl == "4") * 1
  ddipl5 <- (ddipl == "5") * 1
  ddipl6 <- (ddipl == "6") * 1
  ddipl7 <- (ddipl == "7") * 1

  femme <- (sexe == "2") * 1

  age1530 <- (age >= 15 & age <= 30) * 1
  age3150 <- (age >= 31 & age <= 50) * 1
  age5165 <- (age >= 51 & age <= 65) * 1
  
  agri <- (nafg4n == "ES") * 1
  indus <- (nafg4n == "ET") * 1
  cons <- (nafg4n == "EU") * 1
  tert <- (nafg4n == "EV") * 1
  
})

# Estimation d'un modèle un peu plus complet
m4 <- glm(stable ~ age1530 + age5165 + femme + ddipl1 + ddipl3 + ddipl5 + ddipl6 + ddipl7 + agri + cons + tert,data=e,family= binomial(link="logit"))

# Affichage des informations sur le modèle m4 
(t1a <- cbind(
  Var=rownames(summary(m4)$coefficients)
  ,data.frame(summary(m4)$coefficients,row.names=NULL)
))
##            Var    Estimate Std..Error     z.value     Pr...z..
## 1  (Intercept)  2.67652698  0.2551488  10.4900633 9.596409e-26
## 2      age1530 -1.62082391  0.1522592 -10.6451633 1.836619e-26
## 3      age5165  0.40179432  0.1942094   2.0688712 3.855818e-02
## 4        femme -0.42142108  0.1428129  -2.9508605 3.168900e-03
## 5       ddipl1  0.66964571  0.2178210   3.0742943 2.110012e-03
## 6       ddipl3  0.72026357  0.2384301   3.0208581 2.520595e-03
## 7       ddipl5  0.35229052  0.1989732   1.7705428 7.663677e-02
## 8       ddipl6  0.07985644  0.2737797   0.2916814 7.705302e-01
## 9       ddipl7 -0.27566101  0.2138604  -1.2889761 1.974064e-01
## 10        agri -0.65268362  0.5409504  -1.2065499 2.276055e-01
## 11        cons -0.37960462  0.3232267  -1.1744221 2.402260e-01
## 12        tert -0.55170790  0.2159040  -2.5553395 1.060843e-02
# Calcul des odds-ratio
t1b <- t1a[,c("Var","Estimate")]
t1b$OR <- exp(t1b$Estimate)
t1b$low <- exp(t1b$Estimate - qnorm(0.975) * t1a$Std..Error)
t1b$high <- exp(t1b$Estimate + qnorm(0.975) * t1a$Std..Error)
t1b
##            Var    Estimate         OR       low       high
## 1  (Intercept)  2.67652698 14.5345268 8.8149094 23.9653591
## 2      age1530 -1.62082391  0.1977357 0.1467175  0.2664946
## 3      age5165  0.40179432  1.4945039 1.0213761  2.1867968
## 4        femme -0.42142108  0.6561138 0.4959256  0.8680440
## 5       ddipl1  0.66964571  1.9535451 1.2747174  2.9938702
## 6       ddipl3  0.72026357  2.0549748 1.2878177  3.2791296
## 7       ddipl5  0.35229052  1.4223217 0.9630117  2.1007004
## 8       ddipl6  0.07985644  1.0831316 0.6333438  1.8523493
## 9       ddipl7 -0.27566101  0.7590702 0.4991644  1.1543042
## 10        agri -0.65268362  0.5206467 0.1803372  1.5031451
## 11        cons -0.37960462  0.6841318 0.3630853  1.2890534
## 12        tert -0.55170790  0.5759653 0.3772407  0.8793749

Evaluation de la pertinence du modèle

# Centralisation des statistiques construites à partir de la vraisemblance
(t2 <- data.frame(
  modele=c(paste0("sim",1:3),paste0("m",1:4))
  , source=c(rep("Simulations",3),rep("EEC",4))
  , description=c("Lineaire simple","Logit simple","Probit simple","Lineaire simple","Logit simple","Probit simple","Logit multiple")
  , logLik=sapply(list(sim1,sim2,sim3,m1,m2,m3,m4),logLik)
  , aic=sapply(list(sim1,sim2,sim3,m1,m2,m3,m4),function(x){x$aic})
))
##   modele      source     description     logLik        aic
## 1   sim1 Simulations Lineaire simple  -39.68448   85.36896
## 2   sim2 Simulations    Logit simple  -35.26943   74.53886
## 3   sim3 Simulations   Probit simple  -35.20471   74.40942
## 4     m1         EEC Lineaire simple -660.80047 1327.60093
## 5     m2         EEC    Logit simple -762.66522 1529.33044
## 6     m3         EEC   Probit simple -768.51554 1541.03109
## 7     m4         EEC  Logit multiple -744.08332 1512.16664
# Installation et chargement du package pROC (attention à la casse !)
# install.packages("pROC")
library(pROC)
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
# Courbe ROC sur les données simuleés
plot(roc(z ~ predict(sim2)),xlab="Spécificité : Part des négatifs classés négatifs",ylab="Sensibilité : Part des positifs classés positifs",main="Courbe ROC sur les données simulées")

## 
## Call:
## roc.formula(formula = z ~ predict(sim2))
## 
## Data: predict(sim2) in 41 controls (z 0) < 59 cases (z 1).
## Area under the curve: 0.9177
# Courbe ROC sur les données de l'EEC
plot(roc(e$stable ~ predict(m4)),xlab="Spécificité : Part des négatifs classés négatifs",ylab="Sensibilité : Part des positifs classés positifs",main="Courbe ROC sur les données de l'EEC")

## 
## Call:
## roc.formula(formula = e$stable ~ predict(m4))
## 
## Data: predict(m4) in 288 controls (e$stable 0) < 1818 cases (e$stable 1).
## Area under the curve: 0.7426

 

Séance 4 : Aspects transversaux et approfondissements

Cette séance aborde plusieurs questions communes aux modèles de régression sur données quantitatives et qualitatives :

  1. Hétéroscédasticité et erreurs standards robustes
  2. Causalité et endogénéité
  3. Utiliser des pondérations dans une régression

Hétéroscédasticité et erreurs standards robustes

# Génération de données aléatoires
n <- 500
set.seed(1)
age <- round(runif(n)*60+20,0)
prop <- F
prop[age < 30] <- (runif(sum(age<30)) < 0.10)
prop[age >= 30 & age < 50] <- (runif(sum(age >= 30 & age < 50)) < 0.50)
prop[age >= 50]<- (runif(sum(age >= 50)) < 0.70)
pat <- pmax(100*age + rnorm(n,10000,age*1000) + prop * rnorm(n,100000,age*6000),0)

# Représentation du patrimoine brut en fonction de l'âge
plot(age,pat,xlab = "Âge en années", ylab = 
"Patrimoine brut")
title(main = "Patrimoine brut en fonction de l'âge")

# Estimation du modèle linéaire simple m1
m1 <- lm(pat ~ age)
summary(m1)
## 
## Call:
## lm(formula = pat ~ age)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -170916 -105808  -61397    5536 1087919 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2894.5    26209.1   0.110    0.912    
## age           2100.3      498.5   4.214 2.98e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 189400 on 498 degrees of freedom
## Multiple R-squared:  0.03442,    Adjusted R-squared:  0.03248 
## F-statistic: 17.75 on 1 and 498 DF,  p-value: 2.985e-05
# Représentation de la prédiction du modèle m1
plot(age,pat,xlab = "Âge en années", ylab = 
"Patrimoine brut")
abline(coef(m1)[1],coef(m1)[2])
title(main = "Prédiction du modèle m1")

# Représentation des résidus du modèle m1
plot(age,m1$residuals,xlab = "Âge en années", ylab = 
"Résidus du modèle m1")
abline(0,0)
title(main = "Résidus du modèle m1")

sd(m1$residuals[age < 30])
## [1] 63013.76
sd(m1$residuals[age > 50])
## [1] 233291.5
# Estimation du modèle linéaire multiple m2
age2 <- age^2
m2 <- lm(pat ~ age + age2)
summary(m2)
## 
## Call:
## lm(formula = pat ~ age + age2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -147589 -114678  -62547    3534 1075524 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -94956.50   78284.70  -1.213   0.2257  
## age           6447.70    3315.34   1.945   0.0524 .
## age2           -42.85      32.31  -1.326   0.1853  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 189300 on 497 degrees of freedom
## Multiple R-squared:  0.03783,    Adjusted R-squared:  0.03396 
## F-statistic:  9.77 on 2 and 497 DF,  p-value: 6.889e-05
# Représentation de la prédiction du modèle m2
plot(age,pat,xlab = "Âge en années", ylab = 
"Patrimoine brut")
curve(coef(m2)[1] + coef(m2)[2]*x + coef(m2)[3] *x^2,add = T)
title(main = "Prédiction du modèle m2")

# Représentation des résidus du modèle m1
plot(age,m2$residuals,xlab = "Âge en années", ylab = 
"Résidus du modèle m2")
abline(0,0)
title(main = "Résidus du modèle m2")

sd(m2$residuals[age < 30])
## [1] 62283.16
sd(m2$residuals[age > 50])
## [1] 232932.4
# Installation et chargement du package sandwich
# install.packages("sandwich")
library(sandwich)

# Fonction pour calculée automatiquement les erreurs standards robustes
# Inspiré de : https://thetarzan.wordpress.com/2011/05/28/heteroskedasticity-robust-and-clustered-standard-errors-in-r/
robust <- function(m){
  require(sandwich)
  n <- nrow(m$model)
  p <- ncol(m$model)
  std <- (n / (n - p) )*sqrt(diag(sandwich(m)))
  t <- m$coefficients/std
  pval <- 2*pnorm(-abs(t))
  r <- matrix(
    c(m$coefficients, std, t, pval)
    , nrow=p
    , dimnames=dimnames(summary(m)$coefficients)
  )
  dimnames(r)[[2]][2] <- "Robust std. err."
  return(r)
}

# Comparaison des erreurs standards
(t1 <- summary(m2)$coefficients)
##                Estimate  Std. Error   t value   Pr(>|t|)
## (Intercept) -94956.4960 78284.70409 -1.212964 0.22571994
## age           6447.6972  3315.33867  1.944808 0.05236177
## age2           -42.8501    32.30654 -1.326360 0.18532971
(t2 <- robust(m2))
##                Estimate Robust std. err.   t value   Pr(>|t|)
## (Intercept) -94956.4960      63720.53507 -1.490202 0.13617102
## age           6447.6972       3090.62782  2.086210 0.03695965
## age2           -42.8501         32.81246 -1.305909 0.19158337

Causalité et endogénéité

# Initialisation du générateur de nombres aléatoires de R
set.seed(1)

# Simulation des variables explicatives
n <- 1000
cons <- 1 * (runif(n) < 0.2)
eff <- round(rchisq(n,14)^2,0)
hr <- eff*rnorm(n,1478,50)
rp <- 1 * (runif(n) + (eff > 49) * 0.1 > 0.6)
pctsy <- pmin(round(rchisq(n,4.5)^1.5,0),100)

# Génération des variables de prévention et d'accident du travail 
# l'une et l'autre corrélée au risque inhérent au métier 
risque <- rnorm(n)

# Utilisation de la fonction logit pour générer l'indicatrice de prévention
logit <- function(x) exp(x)/(1 + exp(x))
prev <- 1 * (0.50 < logit(-1 + 0.2*risque + 0.6*rp + 0.03*pctsy + 0.2*(eff > 250)) )
summary(prev)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.274   1.000   1.000
cor(prev,pctsy)
## [1] 0.534491
cor(prev,rp)
## [1] 0.3805093

# Utilisation d'une loi de Poisson pour générer le nombre d'accidents du travail par établissement
nbAT <- rpois(n,eff*pmax(0.25 + 0.15*risque - 0.05*prev + 0.3*cons - 0.0001*eff,0))
txAT <- nbAT * 1000000/hr
summary(nbAT)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   18.00   41.00   60.11   79.00  618.00
summary(nbAT/eff)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.1459  0.2566  0.2819  0.3930  1.0740
summary(txAT)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0   100.6   174.5   190.9   266.9   757.6
cor(nbAT,pctsy)
## [1] -0.01258487
cor(nbAT,rp)
## [1] -0.002652238
# Représentation de la relation entre prévention et txAT
plot(prev,txAT,xlab="Existence d'un dispositif de prévention",ylab="Nombre d'accidents du travail par million d'heures travaillées")
title(main = "Prévention des risques et accidents du travail")

# Estimation du modèle m3
m3 <- lm(txAT ~ prev + cons + eff)
summary(m3)
## 
## Call:
## lm(formula = txAT ~ prev + cons + eff)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -299.26  -71.97   -4.26   62.39  419.95 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 161.78882    5.53076  29.253  < 2e-16 ***
## prev         38.00129    7.20440   5.275 1.63e-07 ***
## cons        187.26146    7.89109  23.731  < 2e-16 ***
## eff          -0.08413    0.01870  -4.499 7.64e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 99.81 on 996 degrees of freedom
## Multiple R-squared:  0.3763, Adjusted R-squared:  0.3744 
## F-statistic: 200.3 on 3 and 996 DF,  p-value: < 2.2e-16
# Régression de première étape
m4 <- lm(prev ~ rp + pctsy + cons + eff)
summary(m4)
## 
## Call:
## lm(formula = prev ~ rp + pctsy + cons + eff)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06562 -0.24927 -0.03601  0.18195  0.85225 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.612e-01  2.274e-02 -11.487  < 2e-16 ***
## rp           3.618e-01  2.034e-02  17.786  < 2e-16 ***
## pctsy        2.246e-02  9.183e-04  24.459  < 2e-16 ***
## cons        -8.374e-03  2.536e-02  -0.330    0.741    
## eff          4.622e-04  5.906e-05   7.827 1.28e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3208 on 995 degrees of freedom
## Multiple R-squared:  0.4854, Adjusted R-squared:  0.4833 
## F-statistic: 234.6 on 4 and 995 DF,  p-value: < 2.2e-16
# Récupération de la prédiction du modèle 4
prevhat <- m4$fitted.values
# Régression de seconde étape
m5 <- lm(txAT ~ prevhat + cons + eff)
summary(m5)
## 
## Call:
## lm(formula = txAT ~ prevhat + cons + eff)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -309.27  -69.41    0.56   67.28  403.77 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 172.08063    5.74262  29.966  < 2e-16 ***
## prevhat     -23.69201   10.66651  -2.221  0.02656 *  
## cons        186.85375    7.98097  23.412  < 2e-16 ***
## eff          -0.05412    0.01929  -2.806  0.00512 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 100.9 on 996 degrees of freedom
## Multiple R-squared:  0.3621, Adjusted R-squared:  0.3601 
## F-statistic: 188.4 on 3 and 996 DF,  p-value: < 2.2e-16
# Installation et chargement du package AER
#install.packages("AER") 
library("AER")

# Estimation du modèle m6
m6 <- ivreg(txAT ~ prev + cons + eff | rp + pctsy + cons + eff)
summary(m6)
## 
## Call:
## ivreg(formula = txAT ~ prev + cons + eff | rp + pctsy + cons + 
##     eff)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -315.299  -73.529   -2.021   71.044  406.016 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 172.08063    5.88337  29.249  < 2e-16 ***
## prev        -23.69201   10.92795  -2.168  0.03039 *  
## cons        186.85375    8.17659  22.852  < 2e-16 ***
## eff          -0.05412    0.01976  -2.739  0.00628 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 103.4 on 996 degrees of freedom
## Multiple R-Squared: 0.3304,  Adjusted R-squared: 0.3284 
## Wald test: 179.5 on 3 and 996 DF,  p-value: < 2.2e-16

Utiliser des pondérations dans une régression

# Lecture des données et assignation à l'objet e
e <- readRDS(gzcon(url("http://teaching.slmc.fr/mqs2/ee12t4.rds")))
  
# Passage en minuscule des noms de variables et restriction aux salariés
colnames(e) <- tolower(colnames(e))
e <- e[!is.na(e$salred),]

# Création des variables utilisées dans les modèles
e$lnsalred <- log(e$salred)
e$age2 <- e$age^2
e$femme <- (e$sexe == 2) * 1
e$ddipl1 <- (e$ddipl == 1) * 1 # Diplôme supérieur à baccalauréat + 2 ans
e$ddipl3 <- (e$ddipl == 3) * 1 # Baccalauréat + 2 ans
e$ddipl4 <- (e$ddipl == 4) * 1 # Baccalauréat ou brevet professionnel ou autre diplôme de ce niveau
e$ddipl5 <- (e$ddipl == 5) * 1 # CAP, BEP ou autre diplôme de ce niveau
e$ddipl6 <- (e$ddipl == 6) * 1 # Brevet des collèges
e$ddipl7 <- (e$ddipl == 7) * 1 # Aucun diplôme ou certificat d'études primaires
# Estimation du modèle de référence
m7 <- lm(lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + ddipl5 + ddipl6 + ddipl7, data = e)
summary(m7)
## 
## Call:
## lm(formula = lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + 
##     ddipl5 + ddipl6 + ddipl7, data = e)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2307 -0.2146  0.0671  0.2927  1.9488 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.9391155  0.2305332  25.763  < 2e-16 ***
## age          0.0691972  0.0116498   5.940 4.32e-09 ***
## age2        -0.0006555  0.0001412  -4.641 4.07e-06 ***
## femme       -0.3961442  0.0391407 -10.121  < 2e-16 ***
## ddipl1       0.2587343  0.0621396   4.164 3.48e-05 ***
## ddipl3       0.1863284  0.0652343   2.856 0.004402 ** 
## ddipl5      -0.2074332  0.0568835  -3.647 0.000284 ***
## ddipl6      -0.1825530  0.0955541  -1.910 0.056444 .  
## ddipl7      -0.5127076  0.0667608  -7.680 4.84e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.54 on 770 degrees of freedom
## Multiple R-squared:  0.3067, Adjusted R-squared:  0.2995 
## F-statistic: 42.58 on 8 and 770 DF,  p-value: < 2.2e-16
# Estimation du modèle pondéré
m7w <- lm(lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + ddipl5 + ddipl6 + ddipl7, data = e, weights = extri1613)
summary(m7w)
## 
## Call:
## lm(formula = lnsalred ~ age + age2 + femme + ddipl1 + ddipl3 + 
##     ddipl5 + ddipl6 + ddipl7, data = e, weights = extri1613)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -135.344   -8.141    2.004   10.491   90.410 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.9313057  0.2393566  24.780  < 2e-16 ***
## age          0.0666245  0.0121210   5.497 5.27e-08 ***
## age2        -0.0006081  0.0001470  -4.136 3.93e-05 ***
## femme       -0.3767886  0.0398547  -9.454  < 2e-16 ***
## ddipl1       0.3073208  0.0604662   5.083 4.68e-07 ***
## ddipl3       0.2283322  0.0661951   3.449 0.000592 ***
## ddipl5      -0.1805393  0.0578196  -3.122 0.001861 ** 
## ddipl6      -0.1388831  0.0933496  -1.488 0.137220    
## ddipl7      -0.4939485  0.0689624  -7.163 1.85e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.17 on 770 degrees of freedom
## Multiple R-squared:  0.3069, Adjusted R-squared:  0.2997 
## F-statistic: 42.61 on 8 and 770 DF,  p-value: < 2.2e-16