10 Minutes to Learn APPELPY-Python's Applied Econometrics Library

** Learn in 10 minutes APPEL PY **

This post is a Japanese translation of this notebook.

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')

** Introducing the features of Appelpy. ** **

--Exploratory data analysis --Model estimation --Model diagnosis

Note: This notebook was run using the latest version of Appelpy v0.4.2 and its dependencies.

Model dataset settings

Now let's explore the features of Appelpy using one of the richest datasets in econometrics. The California Test Score Data Set (** Caschool **).

Load the dataset with Pandas and get the data with the Statsmodels method.

# 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()

Inspect the data

Now that we have the data, let's take a closer look at the distribution of some key variables before modeling.

Appelpy has a ** ʻeda` module ** (exploratory data analysis).

from appelpy.eda import statistical_moments

Get all the statistical moments at once.

1. Average
2. Distributed
3. Skewness (default is Fisher)
4. Kurtosis

statistical_moments(df)

Model estimation

At the core of Appelpy is to model the phenomena in the dataset.

The ** linear_model ** module contains classes for modeling data. Example: OLS (Ordinary Least Squares)

from appelpy.linear_model import OLS

** Here is a recipe for estimating a linear model. ** **

1. Create two lists, y_list and X_list, for the dependent variable (y) and the independent variable (X).

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

2. Fit the model in one line by passing the model's dataframe and two lists as input to the OLS object.

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

3. Get the result of the key model with the model_selection_stats attribute, including the root mean square error (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}

Get a summary of Statsmodels with the results_output attribute ...

model_nonrobust.results_output

... Now it's finally easy to get the results of a standardized model using the results_output_standardized attribute.

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

The output is similar to Stata's listcoef command.

Model diagnosis

Call the diagnostic_plot method of the model object to display many visual ** regression diagnostics: **.

--pp_plot: P-P plot -- qq_plot: Q-Q plot --rvp_plot: Plot of residuals and predicted values --rvf_plot: Plot of residuals and measured values

Appelpy's ** diagnostics ** module has a suite of classes and functions for regression diagnostics.

--If visual diagnostics are not enough, call other methods like heteroskedasticity_test to do a formal test.

--Are there any observations that pollute the model? Pass the model object to the BadApples class to see the points of influence, outliers, and high leverage.

drawing

** Here is a recipe for a 2x2 set of diagnostic plots. ** **

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()

Good grief. The specified model has a special problem with test scores above 1350.

There seems to be a non-constant variance in the RVF plot.

It's probably homoscedastic, but let's do a more specific test.

Homoscedasticity

from appelpy.diagnostics import heteroskedasticity_test

Call one of several ** homoscedasticity tests . Example: - Breusch-Pagan * (default for Stata's hettest command) - Breusch-Pagan studentized * (default for R's bptest command) -* White * (Stata's hettest, white command)

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

The Studentized Breusch-Pagan test is considered more robust and shows a significant difference from the initial hypothesis of homoscedasticity, so let's fit it to a model with the Heteroskedastic-Robust standard error (HC1 error). ..

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

model_hc1.results_output

What are the significant regressors (independent variables) now?

model_hc1.significant_regressors(0.05)

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

The str (student-teacher ratio) and ʻexpnstu` regressors are especially relevant for investing in students. You can do a ** Wald test ** to test if the two regressors are together and significantly different from 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

Not surprisingly, the standard error of some variables is higher for the HC1 error than for the non-robust error.

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

Anomalous observations: leverage, outliers, and influence

Let's inspect the observed results.

--High influence --High leverage --Outliers

from appelpy.diagnostics import BadApples

Set up an instance of BadApples using a model object.

bad_apples = BadApples(model_hc1).fit()

Let's plot the leverage value against the square of the normalized residuals.

The observations in the upper right quadrant have higher than average leverage and residual values (that is, the most influential observations in the model).

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

Leverage

Leverage values are stored in the measures_leverage attribute.

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

Let's take a look at the five most leveraged observations.

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

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

These highly leveraged points have at least one independent variable that is at least 4 standard deviations away from the mean.

Outliers

The BadApples object calculates a heuristic to determine if the scale value is "extreme".

Extreme indexes are stored in attributes, for example outliers are stored in ʻindices_outliers`.

bad_apples.indices_outliers

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

Outlier heuristics: If the absolute value of the standardized residuals or the absolute value of the studentized residuals is greater than 2, the observations are outliers.

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

Proportion of observations as outliers: 4.2857%

Due to its heuristics, less than about 5% of the observed values are expected to be considered outliers.

Influence

Also consider various influence measurements such as Cook's distance, dffits and dfbeta.

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

wrap up

--Load the dataset and preprocess it in Pandas. --Appelpy ** ʻeda**: Inspect the model dataset with useful functions such as statistical moments. --Appelpy ** modeling **.Linear_modelcontains classes such as OLS and WLS for model estimation. Use a short recipe for each model to access model estimates. --Appelpy **diagnostics**: Checks if model estimates are valid with diagnostic plots, statistical tests, and diagnostics such asBadApples` for impact analysis.

Let's check other features as well.

--Logistic regression (under discrete_model) --DummyEncoder and InteractionEncoder classes