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

**10分で学ぶAPPELPY**

この投稿は[こちらのノートブック](https://github.com/mfarragher/appelpy-examples/blob/master/00_ten-minutes-to-appelpy.ipynb)の和訳です。

```python
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](https://rdrr.io/cran/Ecdat/man/Caschool.html) (**Caschool**)です。

Pandasでデータセットを読み込み、Statsmodelsメソッドでデータを取得します。


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


```python
df.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>distcod</th>
      <th>county</th>
      <th>district</th>
      <th>grspan</th>
      <th>enrltot</th>
      <th>teachers</th>
      <th>calwpct</th>
      <th>mealpct</th>
      <th>computer</th>
      <th>testscr</th>
      <th>compstu</th>
      <th>expnstu</th>
      <th>str</th>
      <th>avginc</th>
      <th>elpct</th>
      <th>readscr</th>
      <th>mathscr</th>
      <th>scrtot</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>0</td>
      <td>75119</td>
      <td>Alameda</td>
      <td>Sunol Glen Unified</td>
      <td>KK-08</td>
      <td>195</td>
      <td>10.900000</td>
      <td>0.510200</td>
      <td>2.040800</td>
      <td>67</td>
      <td>690.799988</td>
      <td>0.343590</td>
      <td>6384.911133</td>
      <td>17.889910</td>
      <td>22.690001</td>
      <td>0.000000</td>
      <td>691.599976</td>
      <td>690.000000</td>
      <td>1381.599976</td>
    </tr>
    <tr>
      <td>1</td>
      <td>61499</td>
      <td>Butte</td>
      <td>Manzanita Elementary</td>
      <td>KK-08</td>
      <td>240</td>
      <td>11.150000</td>
      <td>15.416700</td>
      <td>47.916698</td>
      <td>101</td>
      <td>661.200012</td>
      <td>0.420833</td>
      <td>5099.380859</td>
      <td>21.524664</td>
      <td>9.824000</td>
      <td>4.583333</td>
      <td>660.500000</td>
      <td>661.900024</td>
      <td>1322.400024</td>
    </tr>
    <tr>
      <td>2</td>
      <td>61549</td>
      <td>Butte</td>
      <td>Thermalito Union Elementary</td>
      <td>KK-08</td>
      <td>1550</td>
      <td>82.900002</td>
      <td>55.032299</td>
      <td>76.322601</td>
      <td>169</td>
      <td>643.599976</td>
      <td>0.109032</td>
      <td>5501.954590</td>
      <td>18.697226</td>
      <td>8.978000</td>
      <td>30.000002</td>
      <td>636.299988</td>
      <td>650.900024</td>
      <td>1287.200012</td>
    </tr>
    <tr>
      <td>3</td>
      <td>61457</td>
      <td>Butte</td>
      <td>Golden Feather Union Elementary</td>
      <td>KK-08</td>
      <td>243</td>
      <td>14.000000</td>
      <td>36.475399</td>
      <td>77.049202</td>
      <td>85</td>
      <td>647.700012</td>
      <td>0.349794</td>
      <td>7101.831055</td>
      <td>17.357143</td>
      <td>8.978000</td>
      <td>0.000000</td>
      <td>651.900024</td>
      <td>643.500000</td>
      <td>1295.400024</td>
    </tr>
    <tr>
      <td>4</td>
      <td>61523</td>
      <td>Butte</td>
      <td>Palermo Union Elementary</td>
      <td>KK-08</td>
      <td>1335</td>
      <td>71.500000</td>
      <td>33.108601</td>
      <td>78.427002</td>
      <td>171</td>
      <td>640.849976</td>
      <td>0.128090</td>
      <td>5235.987793</td>
      <td>18.671329</td>
      <td>9.080333</td>
      <td>13.857677</td>
      <td>641.799988</td>
      <td>639.900024</td>
      <td>1281.700012</td>
    </tr>
  </tbody>
</table>
</div>



# データを検査する

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

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


```python
from appelpy.eda import statistical_moments
```

すべての統計的モーメントを一度に取得します。
1. 平均
2. 分散
3. 歪度（デフォルトはフィッシャー）
4. 尖度


```python
statistical_moments(df)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>mean</th>
      <th>var</th>
      <th>skew</th>
      <th>kurtosis</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>distcod</td>
      <td>67472.8</td>
      <td>1.20201e+07</td>
      <td>-0.0307844</td>
      <td>-1.09626</td>
    </tr>
    <tr>
      <td>enrltot</td>
      <td>2628.79</td>
      <td>1.53124e+07</td>
      <td>2.87017</td>
      <td>10.2004</td>
    </tr>
    <tr>
      <td>teachers</td>
      <td>129.067</td>
      <td>35311.2</td>
      <td>2.93254</td>
      <td>11.219</td>
    </tr>
    <tr>
      <td>calwpct</td>
      <td>13.246</td>
      <td>131.213</td>
      <td>1.68306</td>
      <td>4.58959</td>
    </tr>
    <tr>
      <td>mealpct</td>
      <td>44.7052</td>
      <td>735.678</td>
      <td>0.183954</td>
      <td>-0.999802</td>
    </tr>
    <tr>
      <td>computer</td>
      <td>303.383</td>
      <td>194782</td>
      <td>2.87169</td>
      <td>10.7729</td>
    </tr>
    <tr>
      <td>testscr</td>
      <td>654.157</td>
      <td>363.03</td>
      <td>0.0916151</td>
      <td>-0.254288</td>
    </tr>
    <tr>
      <td>compstu</td>
      <td>0.135927</td>
      <td>0.00421926</td>
      <td>0.922369</td>
      <td>1.43113</td>
    </tr>
    <tr>
      <td>expnstu</td>
      <td>5312.41</td>
      <td>401876</td>
      <td>1.0679</td>
      <td>1.87571</td>
    </tr>
    <tr>
      <td>str</td>
      <td>19.6404</td>
      <td>3.57895</td>
      <td>-0.0253655</td>
      <td>0.609597</td>
    </tr>
    <tr>
      <td>avginc</td>
      <td>15.3166</td>
      <td>52.2135</td>
      <td>2.21516</td>
      <td>6.53212</td>
    </tr>
    <tr>
      <td>elpct</td>
      <td>15.7682</td>
      <td>334.375</td>
      <td>1.4268</td>
      <td>1.4354</td>
    </tr>
    <tr>
      <td>readscr</td>
      <td>654.97</td>
      <td>404.331</td>
      <td>-0.0583962</td>
      <td>-0.358049</td>
    </tr>
    <tr>
      <td>mathscr</td>
      <td>653.343</td>
      <td>351.72</td>
      <td>0.255082</td>
      <td>-0.159811</td>
    </tr>
    <tr>
      <td>scrtot</td>
      <td>1308.31</td>
      <td>1452.12</td>
      <td>0.0916149</td>
      <td>-0.254288</td>
    </tr>
  </tbody>
</table>
</div>



# モデルの推定

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

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


```python
from appelpy.linear_model import OLS
```

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

1. 従属変数(y)と独立変数(X)のために、`y_list` と `X_list` の2つのリストを作成します。


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

2. モデルのデータフレームと2つのリストをOLSオブジェクトの入力として渡すことで、1行でモデルをフィットします。


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

3. 二乗平均平方根誤差(Root MSE)を含む、`model_selection_stats`属性でキーモデルの結果を取得します。


```python
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` 属性で取得します...


```python
model_nonrobust.results_output
```



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

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

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


```python

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



<div>
<style  type="text/css" >
</style><table id="T_2aa4e7e4_8dea_11ea_8105_acde48001122" ><caption>Unstandardized and Standardized Estimates</caption><thead>    <tr>        <th class="blank level0" ></th>        <th class="col_heading level0 col0" >coef</th>        <th class="col_heading level0 col1" >t</th>        <th class="col_heading level0 col2" >P>|t|</th>        <th class="col_heading level0 col3" >coef_stdX</th>        <th class="col_heading level0 col4" >coef_stdXy</th>        <th class="col_heading level0 col5" >stdev_X</th>    </tr>    <tr>        <th class="index_name level0" >scrtot</th>        <th class="blank" ></th>        <th class="blank" ></th>        <th class="blank" ></th>        <th class="blank" ></th>        <th class="blank" ></th>        <th class="blank" ></th>    </tr></thead><tbody>
                <tr>
                        <th id="T_2aa4e7e4_8dea_11ea_8105_acde48001122level0_row0" class="row_heading level0 row0" >avginc</th>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row0_col0" class="data row0 col0" >+3.8640</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row0_col1" class="data row0 col1" >+20.876</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row0_col2" class="data row0 col2" >0.000</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row0_col3" class="data row0 col3" >+27.9209</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row0_col4" class="data row0 col4" >+0.7327</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row0_col5" class="data row0 col5" >7.2259</td>
            </tr>
            <tr>
                        <th id="T_2aa4e7e4_8dea_11ea_8105_acde48001122level0_row1" class="row_heading level0 row1" >str</th>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row1_col0" class="data row1 col0" >-1.3954</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row1_col1" class="data row1 col1" >-1.563</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row1_col2" class="data row1 col2" >0.119</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row1_col3" class="data row1 col3" >-2.6399</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row1_col4" class="data row1 col4" >-0.0693</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row1_col5" class="data row1 col5" >1.8918</td>
            </tr>
            <tr>
                        <th id="T_2aa4e7e4_8dea_11ea_8105_acde48001122level0_row2" class="row_heading level0 row2" >enrltot</th>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row2_col0" class="data row2 col0" >-0.0016</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row2_col1" class="data row2 col1" >-4.722</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row2_col2" class="data row2 col2" >0.000</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row2_col3" class="data row2 col3" >-6.3035</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row2_col4" class="data row2 col4" >-0.1654</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row2_col5" class="data row2 col5" >3913.1050</td>
            </tr>
            <tr>
                        <th id="T_2aa4e7e4_8dea_11ea_8105_acde48001122level0_row3" class="row_heading level0 row3" >expnstu</th>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row3_col0" class="data row3 col0" >-0.0061</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row3_col1" class="data row3 col1" >-2.316</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row3_col2" class="data row3 col2" >0.021</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row3_col3" class="data row3 col3" >-3.8364</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row3_col4" class="data row3 col4" >-0.1007</td>
                        <td id="T_2aa4e7e4_8dea_11ea_8105_acde48001122row3_col5" class="data row3 col5" >633.9371</td>
            </tr>
    </tbody></table>
</div>

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

# モデルの診断

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

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


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

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

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


## 描画

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


```python
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449916/aba24d07-33d1-864c-02ba-4fe6753f9607.png)


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

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

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

## 等分散性


```python
from appelpy.diagnostics import heteroskedasticity_test
```

いくつかの **等分散性検定** のいずれかを呼び出します。例：
- *Breusch-Pagan* (Stataの `hettest` コマンドのデフォルト)
- *Bresuch-Pagan studentized* (Rの `bptest` コマンドのデフォルト)
- *White* (Stataの `hettest, white` コマンド)


```python
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



```python
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`誤差）をもつモデルに適合させましょう。


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


```python
model_hc1.results_output
```



<div>
<table class="simpletable">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>scrtot</td>      <th>  R-squared:         </th> <td>   0.544</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.540</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   78.60</td>
</tr>
<tr>
  <th>Date:</th>             <td>Mon, 04 May 2020</td> <th>  Prob (F-statistic):</th> <td>1.36e-49</td>
</tr>
<tr>
  <th>Time:</th>                 <td>18:32:27</td>     <th>  Log-Likelihood:    </th> <td> -1959.6</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>   420</td>      <th>  AIC:               </th> <td>   3929.</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   415</td>      <th>  BIC:               </th> <td>   3949.</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     4</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>         <td>HC1</td>       <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable">
<tr>
     <td></td>        <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>const</th>   <td> 1312.9203</td> <td>   29.364</td> <td>   44.711</td> <td> 0.000</td> <td> 1255.367</td> <td> 1370.474</td>
</tr>
<tr>
  <th>avginc</th>  <td>    3.8640</td> <td>    0.240</td> <td>   16.119</td> <td> 0.000</td> <td>    3.394</td> <td>    4.334</td>
</tr>
<tr>
  <th>str</th>     <td>   -1.3954</td> <td>    0.932</td> <td>   -1.498</td> <td> 0.134</td> <td>   -3.221</td> <td>    0.430</td>
</tr>
<tr>
  <th>enrltot</th> <td>   -0.0016</td> <td>    0.000</td> <td>   -5.794</td> <td> 0.000</td> <td>   -0.002</td> <td>   -0.001</td>
</tr>
<tr>
  <th>expnstu</th> <td>   -0.0061</td> <td>    0.003</td> <td>   -2.166</td> <td> 0.030</td> <td>   -0.012</td> <td>   -0.001</td>
</tr>
</table>
<table class="simpletable">
<tr>
  <th>Omnibus:</th>       <td> 7.321</td> <th>  Durbin-Watson:     </th> <td>   0.763</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.026</td> <th>  Jarque-Bera (JB):  </th> <td>   7.438</td>
</tr>
<tr>
  <th>Skew:</th>          <td>-0.326</td> <th>  Prob(JB):          </th> <td>  0.0243</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 2.969</td> <th>  Cond. No.          </th> <td>1.39e+05</td>
</tr>
</table><br/><br/>Warnings:<br/>[1] Standard Errors are heteroscedasticity robust (HC1)<br/>[2] The condition number is large, 1.39e+05. This might indicate that there are<br/>strong multicollinearity or other numerical problems.<br></br>
</div>



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


```python
model_hc1.significant_regressors(0.05)
```




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



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


```python
from appelpy.diagnostics import wald_test
```


```python
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誤差の方が高いのは当然のことです。


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




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>std_error_nonrobust</th>
      <th>std_error_hc1</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>const</td>
      <td>27.833430</td>
      <td>29.364467</td>
    </tr>
    <tr>
      <td>avginc</td>
      <td>0.185093</td>
      <td>0.239715</td>
    </tr>
    <tr>
      <td>str</td>
      <td>0.892703</td>
      <td>0.931541</td>
    </tr>
    <tr>
      <td>enrltot</td>
      <td>0.000341</td>
      <td>0.000278</td>
    </tr>
    <tr>
      <td>expnstu</td>
      <td>0.002613</td>
      <td>0.002794</td>
    </tr>
  </tbody>
</table>
</div>



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

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

- 高い影響力
- 高いレバレッジ
- 外れ値


```python
from appelpy.diagnostics import BadApples
```

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


```python
bad_apples = BadApples(model_hc1).fit()
```

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

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


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



![output_53_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/449916/e10393d2-ea10-093c-10c9-31d61a334b7a.png)


### レバレッジ

レバレッジの値は `measures_leverage` 属性に格納されます。


```python
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つの観測値を見てみましょう。


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




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>scrtot</th>
      <th>avginc</th>
      <th>str</th>
      <th>enrltot</th>
      <th>expnstu</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>48</td>
      <td>1261.099976</td>
      <td>12.109128</td>
      <td>19.017494</td>
      <td>27176</td>
      <td>5864.366211</td>
    </tr>
    <tr>
      <td>34</td>
      <td>1252.200012</td>
      <td>11.722225</td>
      <td>21.194069</td>
      <td>25151</td>
      <td>5117.039551</td>
    </tr>
    <tr>
      <td>404</td>
      <td>1387.900024</td>
      <td>55.327999</td>
      <td>16.262285</td>
      <td>1059</td>
      <td>6460.657227</td>
    </tr>
    <tr>
      <td>413</td>
      <td>1398.200012</td>
      <td>50.676998</td>
      <td>15.407042</td>
      <td>687</td>
      <td>7217.263184</td>
    </tr>
    <tr>
      <td>159</td>
      <td>1294.500000</td>
      <td>14.298300</td>
      <td>20.291372</td>
      <td>21338</td>
      <td>5123.474121</td>
    </tr>
  </tbody>
</table>
</div>




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




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>avginc</th>
      <th>str</th>
      <th>enrltot</th>
      <th>expnstu</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>48</td>
      <td>-0.443884</td>
      <td>-0.329278</td>
      <td>6.273077</td>
      <td>0.870684</td>
    </tr>
    <tr>
      <td>34</td>
      <td>-0.497428</td>
      <td>0.821246</td>
      <td>5.755585</td>
      <td>-0.308182</td>
    </tr>
    <tr>
      <td>404</td>
      <td>5.537230</td>
      <td>-1.785664</td>
      <td>-0.401163</td>
      <td>1.811299</td>
    </tr>
    <tr>
      <td>413</td>
      <td>4.893572</td>
      <td>-2.237740</td>
      <td>-0.496228</td>
      <td>3.004803</td>
    </tr>
    <tr>
      <td>159</td>
      <td>-0.140922</td>
      <td>0.344087</td>
      <td>4.781167</td>
      <td>-0.298032</td>
    </tr>
  </tbody>
</table>
</div>



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

### 外れ値

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

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


```python
bad_apples.indices_outliers
```




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



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


```python
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などの様々な影響力測定も考慮してください。


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




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>dfbeta_const</th>
      <th>dfbeta_avginc</th>
      <th>dfbeta_str</th>
      <th>dfbeta_enrltot</th>
      <th>dfbeta_expnstu</th>
      <th>cooks_d</th>
      <th>dffits_internal</th>
      <th>dffits</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>404</td>
      <td>0.010742</td>
      <td>-0.811324</td>
      <td>0.068295</td>
      <td>0.061925</td>
      <td>0.039761</td>
      <td>0.152768</td>
      <td>-0.873981</td>
      <td>-0.882753</td>
    </tr>
    <tr>
      <td>403</td>
      <td>0.024082</td>
      <td>-0.529623</td>
      <td>0.016131</td>
      <td>0.031978</td>
      <td>0.026944</td>
      <td>0.063718</td>
      <td>-0.564440</td>
      <td>-0.567338</td>
    </tr>
    <tr>
      <td>413</td>
      <td>0.074297</td>
      <td>-0.373064</td>
      <td>0.013021</td>
      <td>0.033955</td>
      <td>-0.099557</td>
      <td>0.044116</td>
      <td>-0.469661</td>
      <td>-0.470876</td>
    </tr>
    <tr>
      <td>38</td>
      <td>-0.277356</td>
      <td>-0.193382</td>
      <td>0.217223</td>
      <td>-0.136237</td>
      <td>0.315817</td>
      <td>0.030955</td>
      <td>-0.393412</td>
      <td>-0.397701</td>
    </tr>
    <tr>
      <td>415</td>
      <td>-0.157629</td>
      <td>0.099707</td>
      <td>0.045684</td>
      <td>-0.026836</td>
      <td>0.249169</td>
      <td>0.025414</td>
      <td>0.356468</td>
      <td>0.357936</td>
    </tr>
  </tbody>
</table>
</div>



# 要約

- データセットをロードし、Pandasで前処理を行います。
- Appelpy **`eda`**: 統計モーメントなどの便利な関数でモデルデータセットを検査します。
- Appelpy **modelling**. `linear_model`には、モデル推定のためのOLSやWLSなどのクラスが含まれています。 各モデルの短いレシピを使用して、モデルの推定値にアクセスします。
- Appelpy **`diagnostics`**: モデル推定値が有効かどうかを、診断プロット、統計的検定、影響力分析のための`BadApples`などの診断でチェックします。

###### その他の機能もチェックしてみましょう。
- ロジスティック回帰 (`discrete_model` の下)
- DummyEncoder と InteractionEncoder クラス