Apprenez APPELPY en 10 minutes - Bibliothèque d'économie métrique appliquée Python

** Apprenez en 10 minutes APPEL PY **

Cet article est une traduction japonaise de ce cahier.

import pandas as pd
import numpy as np
import statsmodels.api as sm  # for datasets
import matplotlib.pyplot as plt
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')

** Présentation des fonctionnalités d'Appelpy. ** **

--L'analyse exploratoire des données --Estimation du modèle --Diagnostic de modèle

Remarque: ce notebook a été exécuté avec la dernière version d'Appelpy v0.4.2 et ses dépendances.

Paramètres du jeu de données du modèle

Explorons maintenant les fonctionnalités d'Appelpy en utilisant l'un des ensembles de données les plus riches en économie métrique. The California Test Score Data Set (** Caschool **).

Chargez l'ensemble de données avec Pandas et récupérez les données avec la méthode Statsmodels.

# Load data
df = sm.datasets.get_rdataset('Caschool', 'Ecdat').data
# Add together the two test scores (total score)
df['scrtot'] = df['readscr'] + df['mathscr']

df.head()

Inspectez les données

Maintenant que nous avons les données, examinons de plus près la distribution de certaines variables clés avant la modélisation.

Appelpy a un module ** ʻeda` ** (analyse exploratoire des données).

from appelpy.eda import statistical_moments

Obtenez tous les moments statistiques à la fois.

1. Moyenne
2. Distribué
3. Distorsion (la valeur par défaut est Fisher)
4. Netteté

statistical_moments(df)

Estimation du modèle

Au cœur d'Appelpy se trouve la modélisation de ce qui se passe dans l'ensemble de données.

Le module ** linear_model ** contient des classes pour la modélisation des données. Exemple: OLS (moindres carrés ordinaires)

from appelpy.linear_model import OLS

** Voici une recette pour estimer un modèle linéaire. ** **

1. Créez deux listes, y_list et X_list, pour la variable dépendante (y) et la variable indépendante (X).

y_list = ['scrtot']
X_list = ['avginc', 'str', 'enrltot', 'expnstu']

2. Ajustez le modèle sur une ligne en passant le bloc de données du modèle et deux listes en entrée de l'objet OLS.

model_nonrobust = OLS(df, y_list, X_list).fit()

3. Obtenez le résultat du modèle de clé avec l'attribut model_selection_stats, y compris l'erreur de racine carrée moyenne carrée (Root MSE).

model_nonrobust.model_selection_stats

{'root_mse': 25.85923055850994,
 'r_squared': 0.5438972040191434,
 'r_squared_adj': 0.5395010324916171,
 'aic': 3929.1191674301117,
 'bic': 3949.320440986499}

Obtenez un résumé des Statsmodels avec l'attribut results_output ...

model_nonrobust.results_output

... Maintenant, il est enfin facile d'obtenir les résultats d'un modèle standardisé en utilisant l'attribut results_output_standardized.

model_nonrobust.results_output_standardized
# This is a Styler object; add .data to access the underlying Pandas dataframe

La sortie est similaire à la commande listcoef de Stata.

Diagnostic de modèle

Appelez la méthode diagnostic_plot de l'objet modèle pour afficher de nombreux diagnostics de régression visuels **: **.

--pp_plot: tracé P-P -- qq_plot: tracé Q-Q --rvp_plot: tracé des valeurs résiduelles et prédites --rvf_plot: tracé des valeurs résiduelles et mesurées

Le module ** diagnostics ** d'Appelpy a une suite de classes et de fonctions pour le diagnostic de régression.

• Si les diagnostics visuels ne suffisent pas, appelez d'autres méthodes comme heteroskedasticity_test pour faire un test formel.

• Y a-t-il des observations qui polluent le modèle? Passez l'objet modèle à la classe BadApples pour voir les points d'influence, les valeurs aberrantes et l'effet de levier élevé.

dessin

** Voici une recette pour un ensemble 2x2 de tracés de diagnostic. ** **

fig, axarr = plt.subplots(2, 2, figsize=(10, 10))
model_nonrobust.diagnostic_plot('pp_plot', ax=axarr[0][0])
model_nonrobust.diagnostic_plot('qq_plot', ax=axarr[1][0])
model_nonrobust.diagnostic_plot('rvf_plot', ax=axarr[1][1])
axarr[0, 1].axis('off')
plt.tight_layout()

Bon chagrin. Le modèle spécifié a un problème particulier avec des scores de test supérieurs à 1350.

Il semble y avoir une variance non constante dans le graphique RVF.

Il est probablement réparti uniformément, mais faisons un test plus spécifique.

Dispersion égale

from appelpy.diagnostics import heteroskedasticity_test

Appelez l'un des nombreux ** tests d'isodispersité **. Exemple:

• • Breusch-Pagan * (commande par défaut de Stata hettest)
• • Bresuch-Pagan studentized * (par défaut pour la commande bptest de R)
• • White * (la commande `` hettest, white` de Stata '')

bp_stats = heteroskedasticity_test('breusch_pagan', model_nonrobust)
print('Breusch-Pagan test :: {}'.format(bp_stats['distribution'] + '({})'.format(bp_stats['nu'])))
print('Test statistic: {:.4f}'.format(bp_stats['test_stat']))
print('Test p-value: {:.4f}'.format(bp_stats['p_value']))

Breusch-Pagan test :: chi2(1)
Test statistic: 2.5579
Test p-value: 0.1097

bps_stats = heteroskedasticity_test('breusch_pagan_studentized', model_nonrobust)
print('Breusch-Pagan test (studentized) :: {}'.format(bps_stats['distribution'] + '({})'.format(bps_stats['nu'])))
print('Test statistic: {:.4f}'.format(bps_stats['test_stat']))
print('Test p-value: {:.4f}'.format(bps_stats['p_value']))

Breusch-Pagan test (studentized) :: chi2(1)
Test statistic: 11.3877
Test p-value: 0.0225

Le test Studentized Breusch-Pagan est considéré comme plus robuste et montre des différences significatives par rapport à l'hypothèse d'hypodispersité initiale, donc adaptons-le à un modèle avec l'erreur standard Heteroskedastic-Robust (erreur HC1). ..

model_hc1 = OLS(df, y_list, X_list, cov_type='HC1').fit()

model_hc1.results_output

Quels sont les retributeurs importants (variables indépendantes) actuellement?

model_hc1.significant_regressors(0.05)

['avginc', 'enrltot', 'expnstu']

Les retours de «str» (ratio élèves-enseignant) et «expnstu» sont particulièrement pertinents pour investir dans les étudiants. Vous pouvez faire un ** test Wald ** pour tester si les deux rétributeurs sont ensemble et significativement différents de 0.

from appelpy.diagnostics import wald_test

wald_stats = wald_test(model_hc1, ['str', 'expnstu'])
print('Wald test :: {}'.format(wald_stats['distribution'] + '({})'.format(wald_stats['nu'])))
print('Test statistic: {:.4f}'.format(wald_stats['test_stat']))
print('Test p-value: {:.4f}'.format(wald_stats['p_value']))

Wald test :: chi2(2)
Test statistic: 4.7352
Test p-value: 0.0937

Sans surprise, l'erreur standard de certaines variables est plus élevée avec l'erreur HC1 qu'avec l'erreur non robuste.

pd.concat([pd.Series(model_nonrobust.results.bse, name='std_error_nonrobust'),
           pd.Series(model_hc1.results.bse, name='std_error_hc1')], axis='columns')

Observations anormales: effet de levier, valeurs aberrantes et influence

Inspectons les résultats observés.

--Haute influence --Effet de levier élevé --Plus

from appelpy.diagnostics import BadApples

Configurez une instance de BadApples à l'aide d'un objet modèle.

bad_apples = BadApples(model_hc1).fit()

Tracons la valeur de levier par rapport au carré du résidu normalisé.

Les observations dans le quadrant supérieur droit ont un effet de levier et des valeurs résiduelles supérieurs à la moyenne (c'est-à-dire les observations les plus influentes du modèle).

fig, ax = plt.subplots(figsize=(11,7))
bad_apples.plot_leverage_vs_residuals_squared(annotate=True)
plt.show()

Influence

La valeur de l'effet de levier est stockée dans l'attribut mesures_leverage.

bad_apples.measures_leverage.sort_values(ascending=False).head()

48     0.113247
34     0.085633
404    0.076762
413    0.065741
159    0.061353
Name: leverage, dtype: float64

Jetons un coup d'œil aux cinq observations les plus exploitées.

bad_apples.show_extreme_observations().loc[[48,34,404,413,159]]

model_hc1.X_standardized.loc[[48,34,404,413,159]]

Ces points à fort effet de levier ont des valeurs où au moins une variable indépendante est à au moins 4 écarts-types de la moyenne.

Outlier

L'objet BadApples calcule une heuristique pour déterminer si la valeur d'échelle est "extrême".

Les index extrêmes sont stockés dans des attributs, par exemple les valeurs aberrantes sont stockées dans ʻindices_outliers`.

bad_apples.indices_outliers

[5, 6, 7, 8, 9, 10, 13, 14, 37, 38, 47, 134, 370, 392, 403, 404, 415, 417]

Heuristique des valeurs aberrantes: si la valeur absolue des résidus standardisés ou la valeur absolue des résidus étudiés est supérieure à 2, la valeur observée est une valeur aberrante.

print(f"Proportion of observations as outliers: {len(bad_apples.indices_outliers) / len(df):.4%}")

Proportion of observations as outliers: 4.2857%

En raison de son heuristique, moins d'environ 5% des valeurs observées devraient être considérées comme des valeurs aberrantes.

Influence

Considérez également diverses mesures d'influence telles que la distance de cuisson, les dffits et le dfbeta.

bad_apples.measures_influence.sort_values('cooks_d', ascending=False).head()

emballer

--Chargez l'ensemble de données et prétraitez-le avec Pandas. --Appelpy ** ʻeda**: Inspectez l'ensemble de données du modèle avec des fonctions utiles telles que les moments statistiques. --Appelpy ** modélisation **.Linear_modelcontient des classes telles que OLS et WLS pour l'estimation du modèle. Utilisez une courte recette pour chaque modèle pour accéder aux estimations du modèle. --Appelpy **diagnostics **: vérifie si les estimations du modèle sont valides avec des graphiques de diagnostic, des tests statistiques et des diagnostics tels que BadApples` pour l'analyse d'impact.

Vérifions également d'autres fonctionnalités.

--Régression logistique (sous discrete_model)

• Classes DummyEncoder et InteractionEncoder