-
Notifications
You must be signed in to change notification settings - Fork 45
/
3_regression.qmd
596 lines (438 loc) · 18.3 KB
/
3_regression.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
---
title: "Régression : une introduction"
date: 2020-07-15T13:00:00Z
draft: false
weight: 40
slug: regression
tags:
- scikit
- statsmodels
- Machine Learning
- US elections
- Regression
- Exercice
- Modelisation
categories:
- Modélisation
- Exercice
type: book
description: |
La régression linéaire est la première modélisation statistique
qu'on découvre dans un cursus quantitatif. Il s'agit en effet d'une
méthode très intuitive et très riche. Le _Machine Learning_ permet de
l'appréhender d'une autre manière que l'économétrie. Avec `scikit` et
`statsmodels`, on dispose de tous les outils pour satisfaire à la fois
data scientists et économistes.
image: featured_regression.png
echo: false
---
::: {.cell .markdown}
```{python}
#| echo: false
#| output: 'asis'
#| include: true
#| eval: true
import sys
sys.path.insert(1, '../../') #insert the utils module
from utils import print_badges
#print_badges(__file__)
print_badges("content/modelisation/3_regression.qmd")
```
:::
Le précédent chapitre visait à proposer un premier modèle pour comprendre
les comtés où le parti Républicain l'emporte. La variable d'intérêt étant
bimodale (victoire ou défaite), on était dans le cadre d'un modèle de
classification.
Maintenant, sur les mêmes données, on va proposer un modèle de régression
pour expliquer le score du parti Républicain. La variable est donc continue.
Nous ignorerons le fait que ses bornes se trouvent entre 0 et 100 et donc
qu'il faudrait, pour être rigoureux, transformer l'échelle afin d'avoir
des données dans cet intervalle.
{{< include _import_data_ml.qmd >}}
Ce chapitre va utiliser plusieurs _packages_
de modélisation, les principaux étant `Scikit` et `Statsmodels`.
Voici une suggestion d'import pour tous ces _packages_.
```{python}
#| echo: true
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
import sklearn.metrics
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
```
## Principe général
Le principe général de la régression consiste à trouver une loi $h_\theta(X)$
telle que
$$
h_\theta(X) = \mathbb{E}_\theta(Y|X)
$$
Cette formalisation est extrêmement généraliste et ne se restreint d'ailleurs
par à la régression linéaire.
En économétrie, la régression offre une alternative aux méthodes de maximum
de vraisemblance et aux méthodes des moments. La régression est un ensemble
très vaste de méthodes, selon la famille de modèles
(paramétriques, non paramétriques, etc.) et la structure de modèles.
### La régression linéaire
C'est la manière la plus simple de représenter la loi $h_\theta(X)$ comme
combinaison linéaire de variables $X$ et de paramètres $\theta$. Dans ce
cas,
$$
\mathbb{E}_\theta(Y|X) = X\beta
$$
Cette relation est encore, sous cette formulation, théorique. Il convient
de l'estimer à partir des données observées $y$. La méthode des moindres
carrés consiste à minimiser l'erreur quadratique entre la prédiction et
les données observées (ce qui explique qu'on puisse voir la régression comme
un problème de *Machine Learning*). En toute généralité, la méthode des
moindres carrés consiste à trouver l'ensemble de paramètres $\theta$
tel que
$$
\theta = \arg \min_{\theta \in \Theta} \mathbb{E}\bigg[ \left( y - h_\theta(X) \right)^2 \bigg]
$$
Ce qui, dans le cadre de la régression linéaire, s'exprime de la manière suivante :
$$
\beta = \arg\min \mathbb{E}\bigg[ \left( y - X\beta \right)^2 \bigg]
$$
Lorsqu'on amène le modèle théorique ($\mathbb{E}_\theta(Y|X) = X\beta$) aux données,
on formalise le modèle de la manière suivante :
$$
Y = X\beta + \epsilon
$$
Avec une certaine distribution du bruit $\epsilon$ qui dépend
des hypothèses faites. Par exemple, avec des
$\epsilon \sim \mathcal{N}(0,\sigma^2)$ i.i.d., l'estimateur $\beta$ obtenu
est équivalent à celui du Maximum de Vraisemblance dont la théorie asymptotique
nous assure l'absence de biais, la variance minimale (borne de Cramer-Rao).
::: {.cell .markdown}
```{=html}
<div class="alert alert-success" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-pencil"></i> Exercice 1a : Régression linéaire avec scikit</h3>
```
Cet exercice vise à illustrer la manière d'effectuer une régression linéaire avec `scikit`.
Dans ce domaine,
`statsmodels` est nettement plus complet, ce que montrera l'exercice suivant.
L'intérêt principal de faire
des régressions avec `scikit` est de pouvoir comparer les résultats d'une régression linéaire
avec d'autres modèles de régression. Cependant, le chapitre sur les
pipelines montrera qu'on peut très bien insérer, avec quelques efforts
de programmation orientée objet, une régression `statsmodels` dans
un pipeline `scikit`.
L'objectif est d'expliquer le score des Républicains en fonction de quelques
variables. Contrairement au chapitre précédent, où on se focalisait sur
un résultat binaire (victoire/défaite des Républicains), cette
fois on va chercher à modéliser directement le score des Républicains.
1. A partir de quelques variables, par exemple, *'Unemployment_rate_2019', 'Median_Household_Income_2019', 'Percent of adults with less than a high school diploma, 2015-19', "Percent of adults with a bachelor's degree or higher, 2015-19"*, expliquer la variable `per_gop` à l'aide d'un échantillon d'entraînement `X_train` constitué au préalable.
⚠️ Utiliser la variable `Median_Household_Income_2019`
en `log` sinon son échelle risque d'écraser tout effet.
2. Afficher les valeurs des coefficients, constante comprise
3. Evaluer la pertinence du modèle avec le $R^2$ et la qualité du fit avec le MSE.
4. Représenter un nuage de points des valeurs observées
et des erreurs de prédiction. Observez-vous
un problème de spécification ?
```{=html}
</div>
```
:::
```{python}
#| output: false
# 1. Régression linéaire de per_gop sur différentes variables explicatives
xvars = ['Unemployment_rate_2019', 'Median_Household_Income_2019', 'Percent of adults with less than a high school diploma, 2015-19', "Percent of adults with a bachelor's degree or higher, 2015-19"]
df2 = votes[["per_gop"] + xvars].copy()
df2['log_income'] = np.log(df2["Median_Household_Income_2019"])
df2 = df2.dropna().astype(np.float64)
X_train, X_test, y_train, y_test = train_test_split(
df2.drop(["Median_Household_Income_2019","per_gop"], axis = 1),
100*df2[['per_gop']].values.ravel(), test_size=0.2, random_state=0
)
ols = LinearRegression().fit(X_train, y_train)
y_pred = ols.predict(X_test)
#y_pred[:10]
```
```{python}
#| output: false
#2. Afficher les valeurs des coefficients
print(ols.intercept_, ols.coef_)
```
```{python}
#| output: false
# 3. Evaluer la pertinence du modèle
rmse = sklearn.metrics.mean_squared_error(y_test, y_pred, squared = False)
rsq = sklearn.metrics.r2_score(y_test, y_pred)
print('Mean squared error: %.2f'
% rmse)
# The coefficient of determination: 1 is perfect prediction
print('Coefficient of determination: %.2f'
% rsq)
```
```{python}
#| output: false
#4. Nuage de points des valeurs observées
tempdf = pd.DataFrame({"prediction": y_pred, "observed": y_test,
"error": y_test - y_pred})
fig = plt.figure()
g = sns.scatterplot(data = tempdf, x = "observed", y = "error")
g.axhline(0, color = "red")
# La répartition des erreurs n'est clairement pas
# aléatoire en fonction de $X$.
# Le modèle souffre
# donc d'un problème de spécification.
```
Voici le nuage de points de nos erreurs:
```{python}
g.figure.get_figure()
```
```{python}
#| output: false
g.figure.get_figure().savefig("featured_regression.png")
```
Clairement, le modèle présente un problème de spécification.
```{python}
# packages utiles
import statsmodels.api as sm
import statsmodels.formula.api as smf
```
::: {.cell .markdown}
```{=html}
<div class="alert alert-success" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-pencil"></i> Exercice 1b : Régression linéaire avec statsmodels</h3>
```
Cet exercice vise à illustrer la manière d'effectuer une régression linéaire avec `statsmodels` qui offre des fonctionnalités plus proches de celles de `R`, et moins orientées *Machine Learning*.
L'objectif est toujours d'expliquer le score des Républicains en fonction de quelques
variables.
1. A partir de quelques variables, par exemple, *'Unemployment_rate_2019', 'Median_Household_Income_2019', 'Percent of adults with less than a high school diploma, 2015-19', "Percent of adults with a bachelor's degree or higher, 2015-19"*, expliquer la variable `per_gop`.
⚠️ utiliser la variable `Median_Household_Income_2019`
en `log` sinon son échelle risque d'écraser tout effet.
2. Afficher un tableau de régression.
3. Evaluer la pertinence du modèle avec le R^2.
4. Utiliser l'API `formula` pour régresser le score des républicains en fonction de la variable `Unemployment_rate_2019`, de `Unemployment_rate_2019` au carré et du log de
`Median_Household_Income_2019`.
```{=html}
</div>
```
:::
```{python}
#| output: false
# 1. Régression linéaire de per_gop sur différentes variables explicatives
xvars = ['Unemployment_rate_2019', 'Median_Household_Income_2019', 'Percent of adults with less than a high school diploma, 2015-19', "Percent of adults with a bachelor's degree or higher, 2015-19"]
xvars = ['Unemployment_rate_2019', 'Median_Household_Income_2019', 'Percent of adults with less than a high school diploma, 2015-19', "Percent of adults with a bachelor's degree or higher, 2015-19"]
df2 = votes[["per_gop"] + xvars].copy()
df2['log_income'] = np.log(df2["Median_Household_Income_2019"])
df2 = df2.dropna().astype(np.float64)
X = sm.add_constant(df2.drop(["Median_Household_Income_2019","per_gop"], axis = 1))
results = sm.OLS(df2[['per_gop']], X).fit()
```
```{python}
#| output: false
# 2. Afficher le tableau de régression
print(results.summary())
# html_snippet = results.summary().as_html()
```
```{python}
#| output: false
# 3. Calcul du R^2
print("R2: ", results.rsquared)
```
```{python}
#| output: false
# 4. Nouvelle régression avec l'API formula
results = smf.ols('per_gop ~ Unemployment_rate_2019 + I(Unemployment_rate_2019**2) + np.log(Median_Household_Income_2019)', data=df2).fit()
print(results.summary())
```
::: {.cell .markdown}
```{=html}
<div class="alert alert-warning" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-lightbulb"></i> Hint</h3>
```
Pour sortir une belle table pour un rapport sous $\LaTeX$, il est possible d'utiliser
la méthode [`Summary.as_latex`](https://www.statsmodels.org/devel/generated/statsmodels.iolib.summary.Summary.as_latex.html#statsmodels.iolib.summary.Summary.as_latex). Pour un rapport HTML, on utilisera [`Summary.as_html`](https://www.statsmodels.org/devel/generated/statsmodels.iolib.summary.Summary.as_latex.html#statsmodels.iolib.summary.Summary.as_latex)
```{=html}
</div>
```
:::
::: {.cell .markdown}
```{=html}
<div class="alert alert-info" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-comment"></i> Note</h3>
```
Les utilisateurs de `R` retrouveront des éléments très familiers avec `statsmodels`,
notamment la possibilité d'utiliser une formule pour définir une régression.
La philosophie de `statsmodels` est similaire à celle qui a influé sur la construction
des packages `stats` et `MASS` de `R`: offrir une librairie généraliste, proposant
une large gamme de modèles. Néanmoins, `statsmodels` bénéficie de sa jeunesse
par rapport aux packages `R`. Depuis les années 1990, les packages `R` visant
à proposer des fonctionalités manquantes dans `stats` et `MASS` se sont
multipliés alors que `statsmodels`, enfant des années 2010, n'a eu qu'à
proposer un cadre général (les *generalized estimating equations*) pour
englober ces modèles.
```{=html}
</div>
```
:::
### La régression logistique
Ce modèle s'applique à une distribution binaire.
Dans ce cas, $\mathbb{E}_{\theta} (Y|X) = \mathbb{P}_{\theta} (Y = 1|X)$.
La régression logistique peut être vue comme un modèle linéaire en probabilité :
$$
\text{logit}\bigg(\mathbb{E}_{\theta}(Y|X)\bigg) = \text{logit}\bigg(\mathbb{P}_{\theta}(Y = 1|X)\bigg) = X\beta
$$
La fonction $\text{logit}$ est $]0,1[ \to \mathbb{R}: p \mapsto \log(\frac{p}{1-p})$.
Elle permet ainsi de transformer une probabilité dans $\mathbb{R}$.
Sa fonction réciproque est la sigmoïde ($\frac{1}{1 + e^{-x}}$),
objet central du *Deep Learning*.
Il convient de noter que les probabilités ne sont pas observées, c'est l'*outcome*
binaire (0/1) qui l'est. Cela amène à voir la régression logistique de deux
manières différentes :
* En économétrie, on s'intéresse au modèle latent qui détermine le choix de
l'outcome. Par exemple, si on observe les choix de participer ou non au marché
du travail, on va modéliser les facteurs déterminant ce choix ;
* En *Machine Learning*, le modèle latent n'est nécessaire que pour classifier
dans la bonne catégorie les observations
L'estimation des paramètres $\beta$ peut se faire par maximum de vraisemblance
ou par régression, les deux solutions sont équivalentes sous certaines
hypothèses.
::: {.cell .markdown}
```{=html}
<div class="alert alert-info" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-comment"></i> Note</h3>
```
Par défaut, `scikit` applique une régularisation pour pénaliser les modèles
peu parcimonieux (comportement différent
de celui de `statsmodels`). Ce comportement par défaut est à garder à l'esprit
si l'objectif n'est pas de faire de la prédiction.
```{=html}
</div>
```
:::
```{python}
# packages utiles
from sklearn.linear_model import LogisticRegression
import sklearn.metrics
```
::: {.cell .markdown}
```{=html}
<div class="alert alert-success" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-pencil"></i> Exercice 2a : Régression logistique avec scikit</h3>
```
Avec `scikit`, en utilisant échantillons d'apprentissage et d'estimation :
1. Evaluer l'effet des variables déjà utilisées sur la probabilité des Républicains
de gagner. Affichez la valeur des coefficients.
2. Déduire une matrice de confusion et
une mesure de qualité du modèle.
3. Supprimer la régularisation grâce au paramètre `penalty`. Quel effet sur les paramètres estimés ?
```{=html}
</div>
```
:::
```{python}
#| output: false
#1. Modèle logit avec les mêmes variables que précédemment
xvars = ['Unemployment_rate_2019', 'Median_Household_Income_2019', 'Percent of adults with less than a high school diploma, 2015-19', "Percent of adults with a bachelor's degree or higher, 2015-19"]
df2 = votes[["per_gop"] + xvars].copy()
df2['log_income'] = np.log(df2["Median_Household_Income_2019"])
df2 = df2.dropna().astype(np.float64)
df2['y'] = (df2['per_gop']>0.5).astype(int)
X_train, X_test, y_train, y_test = train_test_split(
df2.drop(["Median_Household_Income_2019","y"], axis = 1),
1*df2[['y']].values.ravel(), test_size=0.2, random_state=0
)
clf = LogisticRegression().fit(X_train, y_train)
y_pred = clf.predict(X_test)
print(clf.intercept_, clf.coef_)
```
```{python}
#| output: false
from sklearn.metrics import ConfusionMatrixDisplay
# 2. Matrice de confusion
ConfusionMatrixDisplay.from_predictions(y_test, y_pred)
sc_accuracy = sklearn.metrics.accuracy_score(y_pred, y_test)
sc_f1 = sklearn.metrics.f1_score(y_pred, y_test)
sc_recall = sklearn.metrics.recall_score(y_pred, y_test)
sc_precision = sklearn.metrics.precision_score(y_pred, y_test)
print(sc_accuracy)
print(sc_f1)
print(sc_recall)
print(sc_precision)
```
```{python}
#| output: false
#3. Supprimer la régularisation
clf2 = LogisticRegression(penalty='none').fit(X_train, y_train)
y_pred2 = clf.predict(X_test)
print(clf2.intercept_, clf2.coef_)
# Les coefficients sont complètement différents
```
```{python}
# packages utiles
from scipy import stats
```
::: {.cell .markdown}
```{=html}
<div class="alert alert-success" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-pencil"></i> Exercice 2b : Régression logistique avec statmodels</h3>
```
En utilisant échantillons d'apprentissage et d'estimation :
1. Evaluer l'effet des variables déjà utilisées sur la probabilité des Républicains
de gagner.
2. Faire un test de ratio de vraisemblance concernant l'inclusion de la variable de (log)-revenu.
```{python}
#| output: false
#1. Modèle logit avec les mêmes variables que précédemment
xvars = [
'Unemployment_rate_2019', 'Median_Household_Income_2019',
'Percent of adults with less than a high school diploma, 2015-19',
"Percent of adults with a bachelor's degree or higher, 2015-19"]
df2 = votes[["per_gop"] + xvars]
df2['log_income'] = np.log(df2["Median_Household_Income_2019"])
df2 = df2.dropna().astype(np.float64)
df2['y'] = (df2['per_gop']>0.5).astype(int)
mylogit = smf.logit(
formula = 'y ~ Unemployment_rate_2019 + I(Unemployment_rate_2019**2) + np.log(Median_Household_Income_2019)',
data=df2[df2['Median_Household_Income_2019']>0]).fit()
print(mylogit.summary())
```
```{python}
#| output: false
#2. Faire un test de ratio de vraisemblance
logit_h0 = smf.logit(
formula = 'y ~ Unemployment_rate_2019 + I(Unemployment_rate_2019**2)',
data=df2[df2['Median_Household_Income_2019']>0]).fit()
# print(logit_h0.summary())
lr = -2*(mylogit.llf - logit_h0.llf)
lrdf = (logit_h0.df_resid - mylogit.df_resid)
lr_pvalue = stats.chi2.sf(lr, df=lrdf)
#print(lr_pvalue)
# La pvalue du test de maximum de ratio
# de vraisemblance étant proche de 1,
# cela signifie que la variable log revenu ajoute,
# presque à coup sûr, de l'information au modèle.
```
```{=html}
</div>
```
:::
::: {.cell .markdown}
```{=html}
<div class="alert alert-warning" role="alert">
<h3 class="alert-heading"><i class="fa-solid fa-lightbulb"></i> Hint</h3>
```
La statistique du test est :
$$
LR = -2\log\bigg(\frac{\mathcal{L}_{\theta}}{\mathcal{L}_{\theta_0}}\bigg) = -2(\mathcal{l}_{\theta} - \mathcal{l}_{\theta_0})
$$
```{=html}
</div>
```
:::
## Pour aller plus loin
Ce chapitre n'évoque les enjeux de la régression
que de manière très introductive. Pour compléter ceci,
il est recommandé d'explorer les champs suivants:
- Les modèles linéaires généralisés pour découvrir la régression
avec des hypothèses plus générales que celles que nous avons posées
jusqu'à présent ;
- Les autres modèles de régression de _machine learning_ comme les forêts
aléatoires ;
- Les tests d'hypothèses pour aller plus loin sur ces questions que notre
test de ratio de vraisemblance.