Go to Qiita Advent Calendar 2024 Top
LoginSignup
4
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@lystahi(lystahi)

10分で学ぶAPPELPY - Pythonの応用計量経済学ライブラリ

Last updated at Posted at 2020-05-04

10分で学ぶAPPELPY

この投稿はこちらのノートブックの和訳です。

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

Appelpyの機能をご紹介します。

  • 探索的データ分析
  • モデル推定
  • モデル診断

注意: このノートブックは、Appelpy v0.4.2とその依存関係の最新バージョンを使用して実行されました。

モデルデータセットの設定

それでは、計量経済学の中で最も豊富なデータセットの1つを使って、Appelpyの機能を探ってみましょう。The California Test Score Data Set (Caschool)です。

Pandasでデータセットを読み込み、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()
distcod county district grspan enrltot teachers calwpct mealpct computer testscr compstu expnstu str avginc elpct readscr mathscr scrtot
0 75119 Alameda Sunol Glen Unified KK-08 195 10.900000 0.510200 2.040800 67 690.799988 0.343590 6384.911133 17.889910 22.690001 0.000000 691.599976 690.000000 1381.599976
1 61499 Butte Manzanita Elementary KK-08 240 11.150000 15.416700 47.916698 101 661.200012 0.420833 5099.380859 21.524664 9.824000 4.583333 660.500000 661.900024 1322.400024
2 61549 Butte Thermalito Union Elementary KK-08 1550 82.900002 55.032299 76.322601 169 643.599976 0.109032 5501.954590 18.697226 8.978000 30.000002 636.299988 650.900024 1287.200012
3 61457 Butte Golden Feather Union Elementary KK-08 243 14.000000 36.475399 77.049202 85 647.700012 0.349794 7101.831055 17.357143 8.978000 0.000000 651.900024 643.500000 1295.400024
4 61523 Butte Palermo Union Elementary KK-08 1335 71.500000 33.108601 78.427002 171 640.849976 0.128090 5235.987793 18.671329 9.080333 13.857677 641.799988 639.900024 1281.700012

データを検査する

データを手に入れたので、モデリングを行う前に、いくつかの主要な変数の分布について詳しく調べてみましょう。

Appelpyには、edaモジュール (探索的データ分析)があります。

from appelpy.eda import statistical_moments

すべての統計的モーメントを一度に取得します。

  1. 平均
  2. 分散
  3. 歪度（デフォルトはフィッシャー）
  4. 尖度
statistical_moments(df)
mean var skew kurtosis
distcod 67472.8 1.20201e+07 -0.0307844 -1.09626
enrltot 2628.79 1.53124e+07 2.87017 10.2004
teachers 129.067 35311.2 2.93254 11.219
calwpct 13.246 131.213 1.68306 4.58959
mealpct 44.7052 735.678 0.183954 -0.999802
computer 303.383 194782 2.87169 10.7729
testscr 654.157 363.03 0.0916151 -0.254288
compstu 0.135927 0.00421926 0.922369 1.43113
expnstu 5312.41 401876 1.0679 1.87571
str 19.6404 3.57895 -0.0253655 0.609597
avginc 15.3166 52.2135 2.21516 6.53212
elpct 15.7682 334.375 1.4268 1.4354
readscr 654.97 404.331 -0.0583962 -0.358049
mathscr 653.343 351.72 0.255082 -0.159811
scrtot 1308.31 1452.12 0.0916149 -0.254288

モデルの推定

Appelpyの中核となるのは、データセットの中の現象をモデル化することです。

linear_model モジュールは、データをモデル化するためのクラスを含んでいます。例：OLS（Ordinary Least Squares）

from appelpy.linear_model import OLS

線形モデルを推定するためのレシピを紹介します。

  1. 従属変数(y)と独立変数(X)のために、y_list と X_list の2つのリストを作成します。
y_list = ['scrtot']
X_list = ['avginc', 'str', 'enrltot', 'expnstu']
  1. モデルのデータフレームと2つのリストをOLSオブジェクトの入力として渡すことで、1行でモデルをフィットします。
model_nonrobust = OLS(df, y_list, X_list).fit()
  1. 二乗平均平方根誤差(Root MSE)を含む、model_selection_stats属性でキーモデルの結果を取得します。
model_nonrobust.model_selection_stats

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

Statsmodels のサマリーを results_output 属性で取得します...

model_nonrobust.results_output
OLS Regression Results
Dep. Variable: scrtot R-squared: 0.544
Model: OLS Adj. R-squared: 0.540
Method: Least Squares F-statistic: 123.7
Date: Mon, 04 May 2020 Prob (F-statistic): 2.05e-69
Time: 18:32:26 Log-Likelihood: -1959.6
No. Observations: 420 AIC: 3929.
Df Residuals: 415 BIC: 3949.
Df Model: 4
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
const 1312.9203 27.833 47.171 0.000 1258.208 1367.632
avginc 3.8640 0.185 20.876 0.000 3.500 4.228
str -1.3954 0.893 -1.563 0.119 -3.150 0.359
enrltot -0.0016 0.000 -4.722 0.000 -0.002 -0.001
expnstu -0.0061 0.003 -2.316 0.021 -0.011 -0.001
Omnibus: 7.321 Durbin-Watson: 0.763
Prob(Omnibus): 0.026 Jarque-Bera (JB): 7.438
Skew: -0.326 Prob(JB): 0.0243
Kurtosis: 2.969 Cond. No. 1.39e+05


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.39e+05. This might indicate that there are
strong multicollinearity or other numerical problems.


...これでようやく results_output_standardized 属性を使って標準化されたモデルの結果を簡単に取得できるようになりました。


model_nonrobust.results_output_standardized
# This is a Styler object; add .data to access the underlying Pandas dataframe
Unstandardized and Standardized Estimates
coef t P>|t| coef_stdX coef_stdXy stdev_X
scrtot
avginc +3.8640 +20.876 0.000 +27.9209 +0.7327 7.2259
str -1.3954 -1.563 0.119 -2.6399 -0.0693 1.8918
enrltot -0.0016 -4.722 0.000 -6.3035 -0.1654 3913.1050
expnstu -0.0061 -2.316 0.021 -3.8364 -0.1007 633.9371

出力はStataのlistcoefコマンドと似ています。

モデルの診断

モデルオブジェクトの diagnostic_plot メソッドを呼び出して、多くの視覚的な 回帰診断を表示します:。

  • pp_plot: P-P プロット
  • qq_plot: Q-Q プロット
  • rvp_plot: 残差と予測値のプロット
  • rvf_plot: 残差と実測値のプロット

Appelpyの diagnostics モジュールは、回帰診断のためのクラスと関数のスイートを持っています。

  • 視覚的な診断だけでは十分でない場合は、heteroskedasticity_testのような他のメソッドを呼び出して正式なテストを行います。

  • モデルを汚している観測値はありますか？ モデルオブジェクトを BadApples クラスに渡して、影響力のある点、外れ値、高レバレッジの点を調べます。

描画

診断プロットの2x2セットのレシピがこちらにあります。

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

output_30_0.png

やれやれ。指定されたモデルには1350以上のテストスコアに 特別な問題があります。

RVFプロットに一定でない分散があるようです。

おそらく等分散性もあると思いますが、もっと具体的なテストをしてみましょう。

等分散性

from appelpy.diagnostics import heteroskedasticity_test

いくつかの 等分散性検定 のいずれかを呼び出します。例：

  • Breusch-Pagan (Stataの hettest コマンドのデフォルト)
  • Bresuch-Pagan studentized (Rの bptest コマンドのデフォルト)
  • White (Stataの hettest, white コマンド)
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

Studentized Breusch-Pagan検定は、よりロバストであると見なされ、等分散性の初期仮説との有意差を示しているので、Heteroskedastic-Robust標準誤差（ HC1誤差）をもつモデルに適合させましょう。

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

model_hc1.results_output
OLS Regression Results
Dep. Variable: scrtot R-squared: 0.544
Model: OLS Adj. R-squared: 0.540
Method: Least Squares F-statistic: 78.60
Date: Mon, 04 May 2020 Prob (F-statistic): 1.36e-49
Time: 18:32:27 Log-Likelihood: -1959.6
No. Observations: 420 AIC: 3929.
Df Residuals: 415 BIC: 3949.
Df Model: 4
Covariance Type: HC1
coef std err z P>|z| [0.025 0.975]
const 1312.9203 29.364 44.711 0.000 1255.367 1370.474
avginc 3.8640 0.240 16.119 0.000 3.394 4.334
str -1.3954 0.932 -1.498 0.134 -3.221 0.430
enrltot -0.0016 0.000 -5.794 0.000 -0.002 -0.001
expnstu -0.0061 0.003 -2.166 0.030 -0.012 -0.001
Omnibus: 7.321 Durbin-Watson: 0.763
Prob(Omnibus): 0.026 Jarque-Bera (JB): 7.438
Skew: -0.326 Prob(JB): 0.0243
Kurtosis: 2.969 Cond. No. 1.39e+05


Warnings:
[1] Standard Errors are heteroscedasticity robust (HC1)
[2] The condition number is large, 1.39e+05. This might indicate that there are
strong multicollinearity or other numerical problems.

今、有意な回帰子（独立変数）は何でしょうか？

model_hc1.significant_regressors(0.05)

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

str (student-teacher ratio) と expnstu の回帰子は、特に生徒への投資に関連しています。 2つの回帰子が一緒になって0から有意に異なるかどうかを検定するために、Wald検定 を行うことができます。

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

一部の変数の標準誤差が、非ロバスト誤差と比較して、HC1誤差の方が高いのは当然のことです。

pd.concat([pd.Series(model_nonrobust.results.bse, name='std_error_nonrobust'),
           pd.Series(model_hc1.results.bse, name='std_error_hc1')], axis='columns')
std_error_nonrobust std_error_hc1
const 27.833430 29.364467
avginc 0.185093 0.239715
str 0.892703 0.931541
enrltot 0.000341 0.000278
expnstu 0.002613 0.002794

異常な観測値：レバレッジ、外れ値、そして影響力

観察した結果を検査してみましょう。

  • 高い影響力
  • 高いレバレッジ
  • 外れ値
from appelpy.diagnostics import BadApples

モデルオブジェクトを用いて BadApples のインスタンスを設定する。

bad_apples = BadApples(model_hc1).fit()

正規化された残差の2乗に対してレバレッジ値をプロットしてみましょう。

右上の象限にある観測値は、平均よりも高いレバレッジ値と残差値を持っています（つまり、モデルで最も影響力のある観測値です）。

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

output_53_0.png

レバレッジ

レバレッジの値は measures_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

レバレッジが最も高い5つの観測値を見てみましょう。

bad_apples.show_extreme_observations().loc[[48,34,404,413,159]]
scrtot avginc str enrltot expnstu
48 1261.099976 12.109128 19.017494 27176 5864.366211
34 1252.200012 11.722225 21.194069 25151 5117.039551
404 1387.900024 55.327999 16.262285 1059 6460.657227
413 1398.200012 50.676998 15.407042 687 7217.263184
159 1294.500000 14.298300 20.291372 21338 5123.474121 
model_hc1.X_standardized.loc[[48,34,404,413,159]]
avginc str enrltot expnstu
48 -0.443884 -0.329278 6.273077 0.870684
34 -0.497428 0.821246 5.755585 -0.308182
404 5.537230 -1.785664 -0.401163 1.811299
413 4.893572 -2.237740 -0.496228 3.004803
159 -0.140922 0.344087 4.781167 -0.298032

これらの高レバレッジな点は、少なくとも1つの独立変数が平均値から4標準偏差以上離れた値を持っています。

外れ値

BadApplesオブジェクトは、尺度の値が「極端」であるかどうかを判断するためのヒューリスティックを計算します。

極端なインデックスは属性に格納され、例えば外れ値は indices_outliers に格納されます。

bad_apples.indices_outliers

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

外れ値のヒューリスティック：標準化された残差の絶対値またはスチューデント化された残差の絶対値が2より大きい場合、観測値は外れ値です。

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

Proportion of observations as outliers: 4.2857%

そのヒューリスティックにより、観測値の約5％以下が外れ値と見なされると予想されます。

影響力

また、クックの距離、dffits、dfbetaなどの様々な影響力測定も考慮してください。

bad_apples.measures_influence.sort_values('cooks_d', ascending=False).head()
dfbeta_const dfbeta_avginc dfbeta_str dfbeta_enrltot dfbeta_expnstu cooks_d dffits_internal dffits
404 0.010742 -0.811324 0.068295 0.061925 0.039761 0.152768 -0.873981 -0.882753
403 0.024082 -0.529623 0.016131 0.031978 0.026944 0.063718 -0.564440 -0.567338
413 0.074297 -0.373064 0.013021 0.033955 -0.099557 0.044116 -0.469661 -0.470876
38 -0.277356 -0.193382 0.217223 -0.136237 0.315817 0.030955 -0.393412 -0.397701
415 -0.157629 0.099707 0.045684 -0.026836 0.249169 0.025414 0.356468 0.357936

要約

  • データセットをロードし、Pandasで前処理を行います。
  • Appelpy eda: 統計モーメントなどの便利な関数でモデルデータセットを検査します。
  • Appelpy modelling. linear_modelには、モデル推定のためのOLSやWLSなどのクラスが含まれています。 各モデルの短いレシピを使用して、モデルの推定値にアクセスします。
  • Appelpy diagnostics: モデル推定値が有効かどうかを、診断プロット、統計的検定、影響力分析のためのBadApplesなどの診断でチェックします。
その他の機能もチェックしてみましょう。
  • ロジスティック回帰 (discrete_model の下)
  • DummyEncoder と InteractionEncoder クラス
4
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
4
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?