POL9590A - Méthodes quantitatives

1.1 Aperçu des données

data(swiss)
head(swiss)

##              Fertility Agriculture Examination Education Catholic
## Courtelary        80.2        17.0          15        12     9.96
## Delemont          83.1        45.1           6         9    84.84
## Franches-Mnt      92.5        39.7           5         5    93.40
## Moutier           85.8        36.5          12         7    33.77
## Neuveville        76.9        43.5          17        15     5.16
## Porrentruy        76.1        35.3           9         7    90.57
##              Infant.Mortality
## Courtelary               22.2
## Delemont                 22.2
## Franches-Mnt             20.2
## Moutier                  20.3
## Neuveville               20.6
## Porrentruy               26.6

Le jeu de données contient 47 observations et 6 variables:

• Fertility: taux standardisé de fécondité (variable dépendante)
• Agriculture: % d’hommes travaillant dans l’agriculture
• Education: % de recrues militaires ayant dépassé l’école primaire
• Catholic: % de catholiques (vs. protestants)

1.2 Estimation du modèle

mod <- lm(Fertility ~ Agriculture + Education + Catholic, data = swiss)
summary(mod)

## 
## Call:
## lm(formula = Fertility ~ Agriculture + Education + Catholic, 
##     data = swiss)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.178  -6.548   1.379   5.822  14.840 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 86.22502    4.73472  18.211  < 2e-16 ***
## Agriculture -0.20304    0.07115  -2.854  0.00662 ** 
## Education   -1.07215    0.15580  -6.881 1.91e-08 ***
## Catholic     0.14520    0.03015   4.817 1.84e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.728 on 43 degrees of freedom
## Multiple R-squared:  0.6423, Adjusted R-squared:  0.6173 
## F-statistic: 25.73 on 3 and 43 DF,  p-value: 1.089e-09

2.1 La matrice chapeau et les leviers

Le levier \(h_{ii}\) mesure l’influence potentielle d’une observation sur les coefficients — basé uniquement sur sa position dans l’espace des prédicteurs \(\mathbf{X}\), indépendamment de \(y_i\).

\[\mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}', \quad h_{ii} = \mathbf{x}_i'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i\]

Seuil empirique: \(h_{ii} > 2(k+1)/n\)

n <- nrow(swiss)
k <- 3  # Agriculture, Education, Catholic
seuil_levier <- 2 * (k + 1) / n
round(seuil_levier, 3)

## [1] 0.17

hat_vals <- hatvalues(mod)
hat_vals

##   Courtelary     Delemont Franches-Mnt      Moutier   Neuveville   Porrentruy 
##   0.09508174   0.06608974   0.12663318   0.05522591   0.04101996   0.12561333 
##        Broye        Glane      Gruyere       Sarine      Veveyse        Aigle 
##   0.05813563   0.06178212   0.07532280   0.06677811   0.06530904   0.06289975 
##      Aubonne     Avenches     Cossonay    Echallens     Grandson     Lausanne 
##   0.07406457   0.06460185   0.07639677   0.05968830   0.06133548   0.09395597 
##    La Vallee       Lavaux       Morges       Moudon        Nyone         Orbe 
##   0.07675957   0.09978790   0.05484826   0.05879560   0.03284225   0.04848798 
##         Oron      Payerne Paysd'enhaut        Rolle        Vevey      Yverdon 
##   0.08570120   0.04914880   0.06662873   0.05398449   0.04684066   0.04094843 
##      Conthey    Entremont       Herens     Martigwy      Monthey   St Maurice 
##   0.08950952   0.09179174   0.09923514   0.07504439   0.07416766   0.07686721 
##       Sierre         Sion       Boudry La Chauxdfnd     Le Locle    Neuchatel 
##   0.08639200   0.06679670   0.03904427   0.15527423   0.09035896   0.12652423 
##   Val de Ruz ValdeTravers V. De Geneve  Rive Droite  Rive Gauche 
##   0.05757509   0.12781384   0.45439561   0.13589989   0.10860138

leviers_eleves <- which(hat_vals > seuil_levier)

round(sort(hat_vals[leviers_eleves], decreasing = TRUE), 3)

## V. De Geneve 
##        0.454

Interprétation: Avec \(n = 47\) et \(k+1 = 4\), le seuil est \(2 \times 4/47 \approx 0.17\). Une seule province dépasse ce seuil: V. De Geneve (\(h_{ii} = 0.454\)). Genève est radicalement atypique dans l’espace des prédicteurs — très peu d’agriculture, éducation élevée, quasi-aucun catholique. Son profil est si éloigné de la moyenne des prédicteurs qu’elle détient à elle seule un pouvoir potentiel considérable sur les coefficients. Cela dit, le levier ne mesure que le potentiel d’influence — la distance de Cook déterminera si ce potentiel est réellement exercé.

2.2 Distance de Cook et graphique d’influence

La distance de Cook mesure l’influence réelle en combinant levier et résidu:

\[D_i = \frac{r_i^2}{k+1}\left(\frac{h_{ii}}{1-h_{ii}}\right)\]

Seuil: \(D_i > 4/n\)

influencePlot(mod,
              main = "Graphique d'influence — swiss",
              sub  = "Taille des cercles proportionnelle à la distance de Cook")

##                  StudRes        Hat        CookD
## Neuveville    2.03100984 0.04101996 0.0411227119
## La Chauxdfnd -1.30151919 0.15527423 0.0766076148
## Neuchatel     1.94289778 0.12652423 0.1284116005
## V. De Geneve -0.05288177 0.45439561 0.0005960716
## Rive Gauche  -2.16798412 0.10860138 0.1318152674

cook_vals <- cooks.distance(mod)
seuil_cook <- 4 / n
round(seuil_cook, 3)

## [1] 0.085

obs_influentes <- which(cook_vals > seuil_cook)

print(round(sort(cook_vals[obs_influentes], decreasing = TRUE), 4))

## Rive Gauche   Neuchatel 
##      0.1318      0.1284

Interprétation: Avec un seuil de \(4/47 \approx 0.085\), deux provinces dépassent le seuil:

• Rive Gauche (\(D_i = 0.132\)): résidu studentisé de \(-2.17\) (le modèle sur-prédit fortement sa fécondité) combiné à un levier de 0.109. C’est l’observation la plus influente du jeu de données.
• Neuchatel (\(D_i = 0.128\)): résidu de \(+1.94\) (sous-prédit) et levier de 0.127. Les deux provinces exercent une traction mesurable sur les coefficients.

Le cas de V. De Geneve illustre parfaitement que levier \(\neq\) influence: malgré son levier extrême (\(h_{ii} = 0.454\), le plus élevé de loin), son \(D_i\) est quasi-nul (\(0.0006\)) car son résidu studentisé est de seulement \(-0.053\) — le modèle prédit parfaitement sa fécondité. Un fort levier sans résidu ne déplace pas les coefficients.

2.3 Résidus studentisés et test de Bonferroni

Les résidus studentisés externes (jackknife) permettent de comparer les résidus sur une échelle commune, en réestimant \(\hat{\sigma}\) sans l’observation \(i\) à chaque fois:

\[r_i^* = \frac{e_i}{\hat{\sigma}_{(-i)}\sqrt{1-h_{ii}}}\]

Le test de Bonferroni teste formellement si le plus grand \(|r_i^*|\) constitue une valeur aberrante, en corrigeant pour les \(n\) comparaisons multiples:

outlierTest(mod)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##              rstudent unadjusted p-value Bonferroni p
## Rive Gauche -2.167984           0.035875           NA

rstud <- rstudent(mod)
print(round(sort(abs(rstud), decreasing = TRUE)[1:5], 3))

## Rive Gauche  Neuveville       Glane   Neuchatel      Sierre 
##       2.168       2.031       1.990       1.943       1.646

Interprétation: Rive Gauche a le plus grand résidu en valeur absolue (\(r^* = -2.168\)) et une \(p\)-valeur brute de \(0.036\) — ce qui est inférieur à 0.05 et pourrait, à première vue, signaler une valeur aberrante. Cependant, il s’agit d’une \(p\)-valeur non-corrigée, calculée comme si on ne testait qu’une seule observation. Or, on teste simultanément les \(n = 47\) résidus. La correction de Bonferroni multiplie cette \(p\)-valeur par \(n\): \(0.036 \times 47 = 1.69\). Puisqu’une probabilité ne peut dépasser 1, R affiche NA — c’est sa façon d’indiquer que la \(p\)-valeur ajustée est techniquement \(> 1\), donc certainement non significative. Le message No Studentized residuals with Bonferroni p < 0.05 confirme qu’aucune valeur aberrante n’est formellement identifiée. Les provinces de Neuveville (\(r^* = 2.031\)), Glane (\(r^* = 1.990\)) et Neuchatel (\(r^* = 1.943\)) ont aussi des résidus relativement élevés. Ce sont des provinces difficiles à prédire par le modèle, non des erreurs de mesure.

3.1 Graphiques CPR (Component + Residual Plots)

Les CPR plots tracent, pour chaque prédicteur \(x_j\), la quantité \(\hat{\beta}_j x_{ij} + e_i\) — c’est-à-dire \(y_i\) moins la contribution de toutes les autres variables. Si la courbe loess (rouge) suit la droite de régression (bleue), la relation est linéaire.

crPlots(mod,
        col.lines = c("blue", "#de1d5b"),
        main = "CPR plots — swiss")

Interprétation:

• Agriculture: la courbe loess présente une légère inflexion aux valeurs intermédiaires, mais reste proche de la droite. Pas de non-linéarité flagrante.
• Education: la courbe loess suit de très près la droite de régression. Relation linéaire bien spécifiée.
• Catholic: déviation modeste.

3.2 Graphiques AVP (Added Variable Plots)

Les AVP montrent l’effet net et partiel de chaque prédicteur après contrôle pour toutes les autres variables. La pente de la droite est exactement \(\hat{\beta}_j\).

avPlots(mod,
        col.lines = "#de1d5b",
        main = "AVP — swiss")

Interprétation:

• Education: pente fortement négative (\(\hat{\beta}_2 = -1.072\)), nuage de points bien aligné. L’effet de l’éducation est la relation la plus robuste du modèle.
• Catholic: pente positive bien définie (\(\hat{\beta}_3 = 0.145\)), avec les provinces très catholiques (>90%) qui ancrent la droite à droite.
• Agriculture: pente négative modeste (\(\hat{\beta}_1 = -0.203\)) — l’effet est réel mais moins marqué que pour les deux autres prédicteurs.

3.3 Test RESET de Ramsey

Le test RESET ajoute \(\hat{y}^2\) et \(\hat{y}^3\) comme régresseurs supplémentaires et teste si leurs coefficients sont conjointement nuls:

\[H_0: \text{la spécification linéaire est correcte (les termes non-linéaires n'améliorent pas le modèle)}\]

resettest(mod, power = 2:3, type = "fitted")

## 
##  RESET test
## 
## data:  mod
## RESET = 0.055965, df1 = 2, df2 = 41, p-value = 0.9456

Interprétation: \(F(2, 41) = 0.056\), \(p = 0.946\) — on ne rejette clairement pas \(H_0\). Malgré les légères déviations visuelles dans les CPR plots (notamment pour Agriculture et Catholic aux extrêmes), elles ne sont pas suffisamment importantes pour invalider la spécification linéaire. Les termes \(\hat{y}^2\) et \(\hat{y}^3\) n’ajoutent rien au modèle.

# Votre code ici

mod2 <- lm(Fertility ~ Agriculture + I(Agriculture^2) + Education + Catholic,
           data = swiss)
summary(mod2)

## 
## Call:
## lm(formula = Fertility ~ Agriculture + I(Agriculture^2) + Education + 
##     Catholic, data = swiss)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.266  -5.624   1.450   4.988  13.873 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      81.168342   6.811567  11.916 4.66e-15 ***
## Agriculture       0.024091   0.231306   0.104    0.918    
## I(Agriculture^2) -0.002344   0.002272  -1.032    0.308    
## Education        -1.028327   0.161372  -6.372 1.16e-07 ***
## Catholic          0.151721   0.030779   4.929 1.34e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.722 on 42 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6179 
## F-statistic: 19.59 on 4 and 42 DF,  p-value: 3.658e-09

resettest(mod2, power = 2:3, type = "fitted")

## 
##  RESET test
## 
## data:  mod2
## RESET = 0.30219, df1 = 2, df2 = 40, p-value = 0.7409

4.1 Graphique Scale-Location

Le graphique Scale-Location trace \(\sqrt{|r_i|}\) en fonction de \(\hat{y}_i\). Une courbe loess horizontale indique une variance constante des résidus (homoscédasticité).

plot(mod, which = 3,
     col = "#de1d5b",
     main = "Scale-Location — swiss")

Interprétation: La courbe loess est relativement plate, sans tendance systématique à la hausse ou à la baisse. Il n’y a pas d’évidence visuelle d’hétéroscédasticité. La dispersion des résidus semble constante sur l’ensemble de la plage des valeurs ajustées.

4.2 Test de Breusch-Pagan

Le test de Breusch-Pagan régresse \(e_i^2\) sur les prédicteurs et teste si les variances résiduelles dépendent de \(\mathbf{x}_i\):

\[H_0: \text{Var}(\varepsilon_i \mid \mathbf{x}_i) = \sigma^2 \quad \forall i \quad \text{(homoscédasticité)}\]

Sous \(H_0\), \(BP \sim \chi^2(k)\) avec \(k = 3\) degrés de liberté.

bptest(mod)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod
## BP = 2.7503, df = 3, p-value = 0.4317

Interprétation: \(\chi^2(3) = 2.750\), \(p = 0.432\) — on ne rejette pas \(H_0\). Il n’y a aucune évidence statistique d’hétéroscédasticité. Les erreurs-types OLS classiques sont valides pour l’inférence dans ce modèle.

5.1 Graphique Q-Q

plot(mod, which = 2,
     col = "#de1d5b",
     main = "Q-Q plot des résidus studentisés — swiss")

Interprétation: Les points suivent bien la diagonale sur la majeure partie de la distribution. Les légères déviations aux extrêmes (Rive Gauche en queue inférieure et Neuveville en queue supérieure) correspondent aux observations à résidus élevés identifiées en section 2, mais restent dans des limites acceptables. Il n’y a pas de déviation en S (asymétrie) ni de valeurs extrêmes très lourdes.

5.2 Test de Shapiro-Wilk

\[H_0: \text{les résidus sont normalement distribués}\]

shapiro.test(residuals(mod))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod)
## W = 0.97153, p-value = 0.3027

Interprétation: \(W = 0.972\), \(p = 0.303\). On ne rejette pas \(H_0\). Les résidus sont tout à fait compatibles avec une distribution normale. L’hypothèse est satisfaite, ce qui valide les tests \(t\) et les intervalles de confiance malgré le petit \(n = 47\).

# Votre code ici

swiss_sans <- swiss[rownames(swiss) != "Rive Gauche", ]
mod_sans   <- lm(Fertility ~ Agriculture + Education + Catholic,
                 data = swiss_sans)

round(coef(mod), 3)

## (Intercept) Agriculture   Education    Catholic 
##      86.225      -0.203      -1.072       0.145

round(coef(mod_sans), 3)

## (Intercept) Agriculture   Education    Catholic 
##      85.544      -0.204      -0.997       0.152

shapiro.test(residuals(mod_sans))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_sans)
## W = 0.9611, p-value = 0.1265