Chap 3: Modèle de régression linéaire

U. Paris Ouest,
M1 - Cours de Modélisation Appliquée
Modèle de régression linéaire:
cas bivarié
Laurent Ferrara
Janvier 2014
U. Paris Ouest
L. Ferrara, 2013-14
Soit 2 variables continues X et Y. On observe les unités
expérimentales : (xi , yi), pour i = 1, …, n.
•
•
•
•
•
•
1. Existe-t-il un lien entre X et Y?
2. Comment le mesurer ?
3. Comment modéliser ce lien?
4. Comment estimer les paramètres de ce modèle?
5. Comment valider ce modèle ?
6. Comment tirer partie de ce modèle pour prévoir les
valeurs d’une variable d’après les valeurs de l’autre?
U. Paris Ouest
L. Ferrara, 2013-14
Exemple : données USA 1992 sur 50 états (state.x77)
40
45
50
Gra
55
60
65
Existe-t-il un lien entre :
les revenus d’un état et le nombre de ses « high-school graduates »?
3000
U. Paris Ouest
3500
4000
4500
5000
Inc
5500
6000
L. Ferrara, 2013-14
71
68
69
70
Life Exp
72
73
Causalité?
Existe-t-il un lien entre :
Le nombre de meurtres et l ’espérance de vie?
2
U. Paris Ouest
4
6
8
Murder
10
12
14
L. Ferrara, 2013-14
Quel type de lien?
• Mise en évidence un lien linéaire entre les 2 variables.
– Y est considérée comme la variable à expliquer, ou indépendante, ou
exogène
– X est considérée comme la variable explicative, ou dépendante, ou
endogène.
• Relation statistique entre les 2 variables (non-déterministe) :
la connaissance de X n’implique pas la connaissance parfaite
de Y : il existe une erreur aléatoire autour de la valeur
prédite
U. Paris Ouest
L. Ferrara, 2013-14
Comment mesurer un lien linéaire?
• Outil principal : Coefficient de corrélation linéaire
Cov( X , Y )
c( X , Y ) 
V ( X )V (Y )
Estimateur empirique :
n
 ( X ,Y ) 
 ( x  X )( y  Y )
i 1
i
n
 ( xi  X )
i 1
U. Paris Ouest
i
n
2
2
(
y

Y
)
 i
i 1
L. Ferrara, 2013-14
Comment mesurer un lien linéaire?
• Signification :
c( X , Y )  1   a, b t.q. : Y  aX  b
c( X , Y )  1 ?
• Test de Student
– H0 :
 ( X ,Y )  0
– H1 :
 ( X ,Y )  0
U. Paris Ouest
L. Ferrara, 2013-14
Comment mesurer un lien linéaire?
• Sous l’hypothèse nulle H0 :
 ( X ,Y )
(1   ( X , Y )) n  2
2
Donc, si
1 / 2
est tq : t* > tn 2
U. Paris Ouest
t* 
suit une loi de Student à (n-2) dl
 ( X ,Y )
(1   2 ( X , Y )) n  2
on rejette H0 au risque 
L. Ferrara, 2013-14
Attention au piège : dépendance non linéaire
le coeff de corrélation ne mesure que la dépendance linéaire.
> cor(x, y)
[1] 0.99
> cor(x, y2)
[1] 0.246
> cor(x, y3)
[1] 0.854
> cor(x, yexp)
[1] 0.898
• Effectuer une analyse graphique au préalable pour identifier
la forme de la dépendance.
• Un coeff de corrélation élevé ne signifie pas forcément une
dépendance linéaire.
U. Paris Ouest
L. Ferrara, 2013-14
3
y2
2
1
0
-1
1
0
y
-1
0
1
-1
6
2
3
yexp
4
5
4
2
0
0
-4
1
-2
y3
1
x
6
x
0
-1
U. Paris Ouest
0
x
1
-1
0
1
x
L. Ferrara, 2013-14
Attention au piège : Corrélation fallacieuse
Existence d’un coeff de corrélation non nul entre deux
variables qu’aucune théorie économique, physique … ne
relie.
2 cas :
– résultat purement aléatoire
– existence d’un troisième variable qui explique conjointement les 2
phénomènes (en général : le temps)
Exemple de Krugman :
lien désindustrialisation - délocalisation aux USA (Application à la France)
U. Paris Ouest
L. Ferrara, 2013-14
s91
ju
il91
no
vm 91
ar
s92
ju
il92
no
vm 92
ar
s93
ju
il93
no
vm 93
ar
s94
ju
il94
no
vm 94
ar
s95
ju
il95
no
vm 95
ar
s96
ju
il96
no
vm 96
ar
s97
ju
il97
no
vm 97
ar
s98
ju
il98
no
vm 98
ar
s99
ju
il99
no
vm 99
ar
s00
ju
il00
no
vm 00
ar
s01
ju
il01
no
vm 01
ar
s02
ju
il02
no
vm 02
ar
s03
m
ar
Evolution de l’emploi industriel France (Trimestriel 1991-2003)
empindus
4700,0
4600,0
4500,0
4400,0
4300,0
4200,0
4100,0
4000,0
U. Paris Ouest
L. Ferrara, 2013-14
s91
ju
il91
no
vm 91
ar
s92
ju
il92
no
vm 92
ar
s93
ju
il93
no
vm 93
ar
s94
ju
il94
no
vm 94
ar
s95
ju
il95
no
vm 95
ar
s96
ju
il96
no
vm 96
ar
s97
ju
il97
no
vm 97
ar
s98
ju
il98
no
vm 98
ar
s99
ju
il99
no
vm 99
ar
s00
ju
il00
no
vm 00
ar
s01
ju
il01
no
vm 01
ar
s02
ju
il02
no
vm 02
ar
s03
m
ar
Evolution des importations de biens en volume France 1991-2003
Imports
2,10
1,90
1,70
1,50
1,30
1,10
0,90
0,70
U. Paris Ouest
L. Ferrara, 2013-14
Corrélation = - 0,50, t de Student = 3,99
Conclusion statistique : on rejette l’hypothèse H0 de nullité de
la corrélation linéaire entre les 2 variables
Conclusion économique rapide : les pays à faibles coûts
salariaux détruisent les emplois dans l ’industrie Française
Or, Krugman a montré qu’en fait les destructions d’emplois
industriels étaient causées en partie par la baisse des
dépenses (en valeur) des ménages en produits manufacturés,
liée à la forte hausse de la productivité dans l’industrie par
comparaison avec celle dans les services
U. Paris Ouest
L. Ferrara, 2013-14
On remarque également que les coefficients de corrélation
entre chacune des variables et le temps sont de :
-0,75 pour l’emploi industriel
0,94 pour les imports
Exercice :
Proposer des exemples de corrélation fallacieuse
U. Paris Ouest
L. Ferrara, 2013-14
Attention au piège :
Un coeff de corrélation nul ne signifie pas que les variables
sont indépendantes (sauf dans le cas Gaussien)
En particulier, il peut exister une relation sur les moments
d’ordre supérieur du modèle
Exemple : lien linéaire entre les variances de X et Y
(cas des processus ARCH en séries chronologiques)
U. Paris Ouest
L. Ferrara, 2013-14
Autres outils de mesure de dépendance:
–
–
–
–
Concordance
Corrélation de rang (Tau de Kendall, coefficient de Spearman)
Corrélation conditionnelle
…
– L’expression générale de la dépendance ne peut se faire que par la
loi jointe.
 Si celle-ci n’est pas calculable: concept de copules
U. Paris Ouest
L. Ferrara, 2013-14
Comment modéliser un lien linéaire?
• Quel est le « meilleur » ajustement linéaire entre 2 v.a. ?
• Exemple : taux longs souverains / dette publique brute
U. Paris Ouest
L. Ferrara, 2013-14
Notation
yi
est la ième observation de la variable exogène
xi
est la ième observation de la variable endogène
yˆ i
Est la valeur ajustée (estimée) de la ième observation
Equation de la meilleure
droite d’ajustement:
U. Paris Ouest
yˆi  b0  b1 xi
L. Ferrara, 2013-14
Erreur de prévision
(ou erreur résiduelle)
En utilisant
yˆ i
pour prédire
yi ,
on fait une erreur de prévision:
ei  yi  yˆ i
La droite d’ajustement qui colle le mieux aux
données est celle pour laquelle les n erreurs de
prévisions sont les plus petites possibles au sens
d’un certain critère.
U. Paris Ouest
L. Ferrara, 2013-14
Critère des “Moindres Carrés”
yˆ i  b0  b1 xi
Equation de la droite :
Choisir les valeurs b0 et b1 qui minimise la somme
des carrés des erreurs.
i.e. : minimiser:
n
2
Q    yi  yˆ i 
U. Paris Ouest
i 1
L. Ferrara, 2013-14
La droite de régression
Par le calcul, minimiser (dériver, annuler et résoudre
pour b0 et b1):
2
n
Q    yi  b0  b1 xi 
i 1
et obtenir les estimateurs des moindres carrés
ordinaires (MCO) de b0 et b1:
 x  x y
n
bˆ1 
U. Paris Ouest
i 1
i
i
y
 x  x 
n
i 1
2

bˆ0  y  bˆ1 x
i
L. Ferrara, 2013-14
Remarques
En termes géométriques
• la droite de régression est celle qui minimise la distance
quadratique entre les points et les projections orthogonales
de ces points sur cette droite.
• la droite de régression est celle qui maximise la variance du
nuage de points projetés orthogonalement sur cette droite.
U. Paris Ouest
L. Ferrara, 2013-14
Formalisation
Hypothèses du modèle linéaire :
• H1 : E(Yi) fonction linéaire des xi (déterministes)
yi = b0 + b1 xi + i ,
pour i=1,…,n
• H2 : Les erreurs, i, sont indépendantes entre elles
• H3 : E(i) = 0, les erreurs sont d’espérance nulle
(en moyenne le modèle est bien spécifié)
U. Paris Ouest
L. Ferrara, 2013-14
• H4 : E(2i) = 2 , les erreurs sont de variance égale
pour toute valeur de X
(hypothèse d ’homoscédasticité)
• H5 : E(Xi i) = 0 , les erreurs,sont indépendantes des
valeurs de X
• H6 : Hypothèse de Normalité
Les erreurs, i, sont identiquement distribuées selon
la loi Normale.
U. Paris Ouest
L. Ferrara, 2013-14
Estimation des paramètres
Quels paramètres ?  b0 , b1 , 2 
bˆ0 , bˆ1 , ˆ 2
bˆ0 , bˆ1 estimés par MCO
ˆ 2
U. Paris Ouest
estimée par l’erreur quadratique moyenne
ou Mean Squared Error (MSE)
L. Ferrara, 2013-14
La MSE est définie par :
n
MSE  ˆ 2 

ˆ
Y

Y
 i i
i 1

2
n2
On pondère par le nombre de degrés de liberté du modèle
défini par :
degrés de liberté = nbre d’observations - nbre de paramètres
U. Paris Ouest
L. Ferrara, 2013-14
Loi asymptotique des paramètres
Les estimateurs MCO sont sans biais et convergents
• On montre que :
• On montre que :
E (bˆ1 )  b1
E (bˆ0 )  b0
V (bˆ1 ) 
ˆ 2
n
2
(
x

X
)
 i
i 1
Donc
U. Paris Ouest
V (bˆ1 )  0 si n  
L. Ferrara, 2013-14
Loi asymptotique des paramètres
• De même,




2
X
2 1

ˆ
V (b0 )  ˆ
 n
n

2
( xi  X ) 


i 1


V (bˆ0 )  0 si n  
U. Paris Ouest
L. Ferrara, 2013-14
Remarques
• Dans ce cadre, sous l ’hypothèse de normalité des erreurs,
estimateur MCO = estimateur EMV
• La variance estimée par le modèle est différente de la
variance empirique (valable pour tout échantillon qui suit le
modèle linéaire)
• La variance résiduelle mesure avec quelle amplitude les
valeurs de Y s ’écartent de la droite de régression.
– C ’est une mesure de la précision du modèle
– C ’est une mesure du risque associé au modèle
U. Paris Ouest
L. Ferrara, 2013-14
Exemple : 2 précisions différentes
U. Paris Ouest
L. Ferrara, 2013-14
Remarques
• Quel est le but du jeu de toute tentative de modélisation
d’une variable Y ?
 Minimiser la variance résiduelle
Y = partie déterministe + partie aléatoire
Y = f(X) + 
Par indépendance, V(Y) = V(f(X)) + V()
(Voir partie « Analyse de la Variance »)
U. Paris Ouest
L. Ferrara, 2013-14
Validation du modèle
On valide le modèle à l’aide des tests statistiques.
2 types de tests d’hypothèses sont développés :
1) Tests sur les paramètres du modèle
2) Tests sur les résidus du modèle
U. Paris Ouest
L. Ferrara, 2013-14
(1-) IC pour la pente bˆ1
Formule en mots:
Paramètre estimé ± (t-multiplier × standard error)
Formule en notations:


ˆ
b1  t1 ,n  2   
2


U. Paris Ouest
ˆ
 x  X 
2
i





L. Ferrara, 2013-14
Test sur la pente bˆ1
Null hypothesis H0: 1 =  (en général =0)
Alternative hypothesis H1: 1 ≠  (en général  0)
Test statistic
t* 
b1  

 MSE





2 
xi  x 


b1  
se b1 

P-value = Risque maximum d’accepter H1 à tort (à
comparer avec le risque de première espèce )
La P-value est déterminée par référence à une tdistribution avec n-2 degrés de liberté
U. Paris Ouest
L. Ferrara, 2013-14
(1-) IC pour la constante bˆ0
Formule en mots:
Paramètre estimé ± (t-multiplier × standard error)
Formule en notations:
bˆ0  t1
U. Paris Ouest
2
,n2
  ˆ
2
1
x

n  xi  X 2
L. Ferrara, 2013-14
Test sur la constante bˆ0
Null hypothesis H0: 0 =  (en général = 0)
Alternative hypothesis HA: 0 ≠  (en général  0)
Test statistic
b0  
t*
MSE
1

n
x

2
 x  x 
b0  
se b0 
2
i
P-value = Risque maximum d’accepter H1 à tort
(à comparer avec le risque de première espèce )
La P-value est déterminée par référence à une tdistribution avec n-2 degrés de liberté.
U. Paris Ouest
L. Ferrara, 2013-14
Test sur le terme d’erreur
Les intervalles et les tests précédents sont basés
sur la Normalité du terme d’erreur. Il importe
donc de tester les résidus.
– Test d’adéquation (Jarque-Bera, KS, …)
– Test graphiques (QQ-Plot)
Les résultats restent valides en cas d’écart à la loi
Normale si l’échantillon est grand. (résultats
asymptotiques)
U. Paris Ouest
L. Ferrara, 2013-14
Exemple : Poids / Taille
> w.fit <- lm(weight ~ 1 + height)
> summary(w.fit)
Call: lm(formula = weight ~ 1 + height)
Residuals:
Min
1Q Median
3Q Max
-13.2 -4.08 -0.0963 4.64 14.2
Coefficients:
Value Std. Error
(Intercept) -266.534
51.032
height
6.138
0.735
t value Pr(>|t|)
-5.223
0.001
8.347
0.000
Residual standard error: 8.64 on 8 degrees of freedom
Multiple R-Squared: 0.897
> resid(w.fit)
1
2
3
4
5
6
7
8
9
10
-5.27 -0.509 -13.2 5.04 3.45 0.0413 14.2 -0.234 6.87 -10.4
U. Paris Ouest
L. Ferrara, 2013-14
Mesure de la qualité du modèle
On mesure la qualité du modèle par l’analyse de la variance
On montre les 2 relations suivantes :
• la somme des résidus est nulle, i.e. :
n
e
i 1
i
0
• la moyenne de la variable et la moyenne de la
variable estimée sont égales, i.e. :
n
n
 y   yˆ
i 1
U. Paris Ouest
i
i 1
i
L. Ferrara, 2013-14
On en déduit l’équation de l’analyse de la variance:
2
ˆ
ˆ
(
y

y
)

(
y

y
)

e
 i
 i
i
2
2
i
i
i
Variance totale = Variance expliquée + Variance résiduelle
Objectif : Maximiser la variance expliquée
U. Paris Ouest
L. Ferrara, 2013-14
• R2 : mesure de la variance expliquée
R2  1
ˆ 2
n
2
(
Y

Y
)
 i
i 1
valeur entre 0 et 1
• Critères d’information : Akaike (1971)
U. Paris Ouest
L. Ferrara, 2013-14
Prévision
Que veut-on prévoir?
• La réponse «moyenne» de la population = E(Yh) pour
une valeur xh
– Ex : Quel est le poids moyen pour une taille donnée?
(Plus précis que le poids moyen de l’échantillon)
• La réponse Yh(new) à une nouvelle valeur donnée xh
– Ex : Quel est le poids estimé par le modèle d’un nouvel
individu choisi au hasard de taille donnée?
U. Paris Ouest
L. Ferrara, 2013-14
En fait les 2 prévisions sont égales :
Yˆh  b0  b1 xh
est le meilleur estimateur dans chaque cas.
Seuls les intervalles de confiance autour des réponses vont varier
U. Paris Ouest
L. Ferrara, 2013-14
Intervalle de confiance pour la
réponse moyenne de la population
E(Yh)
U. Paris Ouest
L. Ferrara, 2013-14
(1-) IC pour la réponse moyenne
E(Yh)
Formule en mots:
Sample estimate ± (t-multiplier × standard error)
Formule en notation:
2 
1

xh  X  
2 
yˆ h  t1 ,n 2   ˆ 

2
n

2


x

X
 i


U. Paris Ouest
L. Ferrara, 2013-14
Implications sur la précision
• Au plus les valeurs des xi sont étalées, au plus
l’intervalle de confiance est petit,
donc l’estimation de E(Yh) est plus précise.
• Suivant le même échantillon de xi, au plus la
valeur de xh est loin de la moyenne empirique, au
plus l’intervalle de confiance est grand,
donc l’estimation de E(Yh) est moins précise.
U. Paris Ouest
L. Ferrara, 2013-14
Remarques
• xh est une valeur correspondant au champ de l’étude
mais pas nécessairement une valeur de l’échantillon
• L’IC pour E(Yh) est correct même si le terme
d’erreur est seulement approché par une loi Normale
• Si le nombre d’observations est grand, l’IC pour
E(Yh) est correct même si le terme d’erreur s’écarte
fortement d’une loi Normale
U. Paris Ouest
L. Ferrara, 2013-14
Exemple : Estimation du poids moyen pour 2 tailles données
(60, proche de la moyenne, et 80, plus élevée que la moyenne)
> predict(w.fit, base2, type = "response", ci.fit = T, se.fit = T)
$fit:
1
2
102 224
$se.fit:
1
2
7.36 8.33
$residual.scale:
[1] 8.64
$df:
[1] 8
$ci.fit:
lower upper
1 84.7
119
2 205.3
244
attr(, "conf.level"):
[1] 0.95
U. Paris Ouest
L. Ferrara, 2013-14
Intervalle de Prévision pour la
réponse Yh(new) à une nouvelle
valeur xh(new)
U. Paris Ouest
L. Ferrara, 2013-14
Prévision de Yh(new)
si la moyenne E(Y) n’est pas connue
ie : si les paramètres sont estimés
 on rajoute une incertitude sur la moyenne
de Y
U. Paris Ouest
L. Ferrara, 2013-14
Propriété:
La prévision est non biaisée
Yˆh  bˆ0  bˆ1 xh
eh  Yh  Yˆh
 b0  b1 xh   h  (bˆ0  bˆ1 xh )
 (b  bˆ )  (b  bˆ ) x  
0
0
1
1
h
h
E (eh )  0
U. Paris Ouest
L. Ferrara, 2013-14
Variance de la prévision
Elle dépend de 2 composantes :
1. Variance due à l’estimation de E(Yh) par
yˆ h
2. Variance de Y inhérente à sa distribution
Estimation:
U. Paris Ouest




2 
2 
1
 1




x

x
x

x
  ˆ 2 1   n h

ˆ 2  ˆ 2   n h
n
2
2
n





x

x
x

x


i
i




i 1
i 1


L. Ferrara, 2013-14
(1-) IC pour la réponse Yh
Sample prediction ± (t-multiplier × standard error)
2 


xh  x  
1
2
yˆ h  t1 ,n 2   ˆ 1  
 n  x  x 2 
2
i


U. Paris Ouest
L. Ferrara, 2013-14
Regression Plot
Mortality = 389.189 - 5.97764 Latitude
S = 19.1150
R-Sq = 68.0 %
R-Sq(adj) = 67.3 %
Mortality
250
150
Regression
95% CI
95% PI
50
30
40
50
Latitude
U. Paris Ouest
L. Ferrara, 2013-14
Exemple d’application
Loi d’Okun (1962)
Relation entre marché du travail et production
• Différentes versions existent dans la littérature, ex:
• Voir exemple sous RATS
U. Paris Ouest
L. Ferrara, 2013-14
Exemple d’application
• Taux de croissance du PIB et Emploi privé US
U. Paris Ouest
L. Ferrara, 2013-14