「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第４章「変数の工夫」②回帰の仮定（多重共線性・不均一分散）

第４章「変数の工夫」

書籍の著者　山澤成康 先生

この記事は、書籍「回帰分析から学ぶ計量経済学」第４章「変数の工夫」の Python写経活動 を取り扱います。

今回は「回帰の仮定」を吟味します！
具体的には、多重共線性と不均一分散に取り組みます！
では書籍を開いて回帰分析の旅に出発です🚀

はじめに

書籍「回帰分析から学ぶ計量経済学」のご紹介

このシリーズは書籍「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」（オーム社、「テキスト」と呼びます）の Python 写経です。

テキストは、2023年11月に発売され、副題「Excelで読み解く経済のしくみ」のとおり、主に Excel を用いて、計量経済学を平易に学べる素晴らしい書籍です。
テキストの「はじめに」に著者の先生が執筆の動機を書かれています。

私にとって計量経済学は高嶺の花ですが、このテキストでさまざまな回帰分析のアプローチを知ることができました。
また、書籍の Excel 処理を Python に置き換える「寄り道写経」の実践を通じて、回帰分析のお気持ちに少し近づけた感じがいたします。

回帰分析に慣れ親しむのに丁度良いレベル感と内容ですので、これはぜひともブログにしたい！と思って現在に至ります。
計量経済学の色を薄め、データ分析の色を濃いめに書いてまいります！

引用表記

この記事は、出典に記載の書籍に掲載された文章と配布データを引用し、適宜、掲載文章・配布データを改変して書いています。
【出典】
「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」
第１版第１刷、著者 山澤成康、オーム社

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

第４章 変数の工夫

この記事は第４章の以下の節を取り扱います。

記事に用いるデータは、オーム社の書籍紹介サイトからダウンロードできる Excel ファイル内のデータをもとにしてCSVファイルを作成し、data フォルダに格納しています。

第４章で用いるライブラリをインポートします。

### インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計処理
import scipy.stats as stats
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.api import het_white     # ホワイトの検定
from statsmodels.stats.outliers_influence import variance_inflation_factor # VIF
import lmdiag  # 残差プロット

# 描画
import matplotlib.pyplot as plt
from matplotlib.ticker import PercentFormatter  # 軸目盛りを％表示
import seaborn as sns
import japanize_matplotlib

# ワーニング非表示
import warnings
warnings.simplefilter('ignore')

4.5 多重共線性の数値例

■ 多重共線性の概要 p.130～
テキストの多重共線性の説明をサマリーします。

多重共線性の有無は VIF で確認できます。
１つの説明変数を目的変数に、残りの説明変数を説明変数にした回帰分析を実行して決定係数$${R^2}$$を求め、 VIF を算出します。
テキストによると説明変数が３つある場合、$${X_1}$$の VIF を求める際の$${X_1}$$を目的変数にする回帰は次の式のようになります。

VIF は以下の式で算出します。
$${R^2}$$は上記の回帰の決定係数です。

多重共線性の判定に関しては、テキストによると、VIF が 10 よりも大きい場合に多重共線性がある、と考えるようです。

多重共線性の解決策は説明変数の相関を減らすことであり、VIF が最も高い変数を除外することが１つの解決策になります。

■ 多重共線性の数値例 p.132～
テキストの仮想データを読み込みます。
説明変数$${X_1}$$と$${X_2}$$の相関が高く、目的変数$${Y}$$は以下の数式（新のモデル）で算出されています。

### データの読み込み

# CSVファイルの読み込み
df5 = pd.read_csv('./data/04_05_corr.csv')
# データフレームの表示
print('df5.shape:', df5.shape)
display(df5.head())

【実行結果】

◆ 説明変数間の相関の確認
相関係数を確認します。

### 相関係数の算出
display(df5.corr().style.map(
    lambda x: 'background-color: lightpink' if (abs(x)>=0.9) & (x!=1)  else ''))

【実行結果】
$${X_1}$$と$${X_2}$$の相関係数は$${0.968}$$です！

相関係数を可視化しましょう。
seaborn の pairplot で散布図行列を描画します。

### 散布図行列の描画
sns.pairplot(data=df5);

【実行結果】
$${X_1}$$と$${X_2}$$、$${Y}$$と$${X_3}$$には強い正の相関がみられます。

◆ 回帰分析
多重共線性がある状態で回帰分析を実行します。

### 回帰分析の実行
result = smf.ols(formula='Y ~ X1 + X2 + X3', data=df5).fit()
display(result.summary())

【実行結果】
$${X_1, X_2}$$の標準誤差（std err）は$${X_3}$$と比べてかなり大きな値になっており、かつ、係数の推定値$${0.3725, 0.5997}$$は真値$${0.5, 0.5}$$とは乖離しています。
多重共線性の影響を受けている感じです。

◆ VIFの計算
VIF を計算します。
statsmodels の variance_inflation_factor を利用します。

### VIFの計算 by statsmodelsのvariance_inflation_factor
# 注意事項：必要に応じて説明変数に定数項を加える

# 説明変数X1～X3と定数項constのデータフレームを作成 
X = sm.add_constant(df5.iloc[:, 1:-1].copy())
# VIFを格納するデータフレームの初期化
vif = pd.DataFrame()
# VIFの算出 ※説明変数ごとにforループでVIF計算を繰り返し処理
vif['VIF'] = [variance_inflation_factor(X, i) for i in range(1, X.shape[1])]
# データフレームのインデックスに説明変数名を設定
vif.index = X.columns[1:]
# 結果の表示
display(vif)

【実行結果】
$${X_1, X_2}$$の VIF は 10 を超えています。
多重共線性あり、です！

以下のコードは$${1 / (1 - R^2)}$$で VIF を計算しています。
statsmodels の回帰分析の結果 result から属性 rsquared で決定係数を取得できます。

### VIF計算 回帰分析⇒決定係数取得⇒VIF計算

## 設定と準備
exp_vars = df5.columns[1:-1]  # 説明変数の列名
vif2 = pd.DataFrame()         # VIFを格納するデータフレーム

## VIF計算
# 説明変数ごとに回帰分析・VIF計算を繰り返し処理
for target in exp_vars:
    # 回帰分析の実行 ※formulaには 'target ~ exp1 + ･･･ ' を設定している
    result = smf.ols(
        formula=target + ' ~ ' + ' + '.join(exp_vars[exp_vars != target]),
        data=df5[exp_vars]).fit()
    # VIFを算出してデータフレームに追加  ※result.rsquaredで決定係数を取得
    vif2.loc[target, 'VIF'] = 1 / (1 - result.rsquared)

## 結果の表示
display(vif2)

【実行結果】

◆ $${X_2}$$を除外して回帰分析
説明変数を$${X_1}$$と$${X_3}$$だけにして回帰分析を実行します。

### X2を除いた回帰分析の実行
result = smf.ols(formula='Y ~ X1 + X3', data=df5).fit()
display(result.summary())

【実行結果】
なんと自由度調整済み決定係数が$${1}$$！

$${X_1}$$の係数の標準誤差が小さくなりました。
多重共線性がない状態だからですね！

テキストによると「$${X_1}$$の係数は$${0.96}$$」であり「$${X_1}$$と$${X_2}$$の係数の和に近い数値になって」います。
係数の真値は$${X_1}$$・$${X_2}$$ともに$${0.5}$$であり、和は$${1}$$。
確かに近い数値になっています。

4.6 都道府県分析の注意点（不均一分散）

■ 不均一の概要 p.134～
テキストの不均一分散の説明をサマリーします。

経済分析では、都道府県別の平均値を使って分析する場合に、不均一分散に注意する必要があるそうです。
各都道府県のサンプル数を揃えていないと、分散が異なる場合があるようです。

◆ ◆ ◆

不均一分散データの分析を進める際には、ぜひ、著者の山澤先生のWeb記事をご覧ください！
テキストよりも詳細な説明を確認できます！
この記事の分析もこのWeb記事を参考にさせていただきました。
ありがとうございます！

■ 不均一分散データ p.134～
データを読み込みます。
データは以下の数式に従って作成されています。

### データの読み込み

# CSVファイルの読み込み
df6 = pd.read_csv('./data/04_06_heteroscedastic.csv')
# データフレームの表示
print('df6.shape:', df6.shape)
display(df6.head())

【実行結果】

$${X, Y}$$の散布図を描画します。

### XとYの散布図の描画

# 描画領域の設定
plt.figure(figsize=(6, 4))
# 散布図の描画
plt.scatter(df6.X, df6.Y, s=70, ec='blue', alpha=0.3)
# 修飾
plt.xlabel('X')
plt.ylabel('Y')
plt.grid(lw=0.5)

【実行結果】
$${X}$$が大きくなるにつれて、$${Y}$$の分散が大きくなっている感じです。

続いて$${X}$$と誤差$${e}$$の散布図を描画します。

### Xとeの散布図の描画

# 描画領域の設定
plt.figure(figsize=(6, 4))
# 散布図の描画
plt.scatter(df6.X, df6.e, s=70, ec='blue', alpha=0.3)
# 縦軸0の水平線の描画
plt.axhline(0, color='tab:red', ls='--')
# 修飾
plt.xlabel('X')
plt.ylabel('e')
plt.grid(lw=0.5)

【実行結果】
$${X}$$が大きくなるにつれて、誤差の分散が大きくなっています。
誤差の分散が不均一です！

■ 通常の回帰分析の実行
不均一分散データで通常の回帰分析を実行します。

### 通常の回帰分析の実行
formula = 'Y ~ X'
result = smf.ols(formula=formula, data=df6).fit()
display(result.summary())

【実行結果】
$${X}$$の係数$${0.5087}$$、標準誤差$${0.028}$$を覚えておきましょう。

ひとまず残差を可視化します。

### 残差の描画

# 描画領域の設定
plt.figure(figsize=(6, 4))
# 散布図の描画
plt.scatter(df6.X, result.resid, s=70, ec='blue', alpha=0.3)
# 縦軸0の水平線の描画
plt.axhline(0, color='tab:red', ls='--')
# 修飾
plt.xlabel('X')
plt.ylabel('残差')
plt.grid(lw=0.5)

【実行結果】
$${X}$$の値が大きくなるにつれて、残差の分散は大きくなっています。
先に見た誤差の散布図の特徴をうまく捉えている感じです。

■ 閑話休題：残差診断プロット
lmdiag ライブラリで statsmodels の回帰分析結果から簡単に残差診断プロットを描画できます。

### 残差診断
plt.figure(figsize=(10, 8))
lmdiag.plot(result);

【実行結果】
左上が残差プロットです。
この他に、正規Q-Qプロット、標準化残差、レバレッジ（外れ値）を可視化できます。

■ 不均一分散の検定：ホワイトの検定 p.136
ホワイトの検定は、帰無仮説「誤差項は均一分散」とし、誤差項の二乗が説明変数および説明変数の二乗と相関があるかを検定するものだそうです。
手順は・・・

① 通常の回帰分析を実施

② 誤差の二乗$${v_i}$$を目的変数、説明変数と説明変数の二乗を説明変数にして回帰分析を実施

③ 上記②に係る$${\beta_1=\beta_2 =0}$$の$${F}$$検定で、不均一分散を判定

先程の通常の回帰の結果を利用して②の回帰を実行し、ホワイトの検定を実践してみます。

### ホワイトの検定 書籍の計算方法
# 残差の二乗を定数項・X^2・Xで回帰分析する
#  ⇒F値のp値0.007は有意⇒不均一分散があると判断する

result_resid = smf.ols(formula='I(result.resid**2) ~ I(X**2) + X', data=df6).fit()
display(result_resid.summary().tables[0])

【実行結果】
$${F}$$値（F-statistic）$${=5.490}$$、$${p}$$値$${=0.0072}$$であり、有意水準$${5\%}$$で帰無仮説「誤差項は均一分散」は棄却され、対立仮説「誤差項は不均一分散」を採択します。
誤差項は不均一分散でした！

◆ 参考：statsmodels の het_white 関数
statsmodels の het_white 関数を利用してホワイトの検定を実行してみましょう。
こちらの関数も通常の回帰分析の結果（残差）を利用します。

### ホワイトの検定 statsmodelsのhet_white関数を利用
#  F値のp値0.007は有意⇒不均一分散があると判断する

# ホワイトの検定の実行：引数 残差の値、説明変数の値
LM_stats, LM_pvalue, F_stats, F_pvalue = het_white(result.resid, result.model.exog)
# 結果表示
print('LM値: ', LM_stats)
print('p値:  ', LM_pvalue, '\n', '-'*30)
print('F値:  ', F_stats)
print('p値:  ', F_pvalue)

【実行結果】
先に実施したホワイトの検定の$${F}$$値・$${p}$$値と同じになりました。

■ 不均一分散の対応１：ロバスト標準誤差（ホワイトの推定量）
テキストは不均一分散の対応策の１つとして「標準誤差を計算しなおす方法」に言及しています。
不均一分散の問題は係数の標準誤差に生じるので標準誤差の計算を工夫する考えだそうです。

テキストによると・・・

もともと、不均一分散の場合、通常の回帰分析では標準誤差が大きめに、$${t}$$値は小さめに、ということでした。
ホワイトの推定量も分散大きめ・$${t}$$値小さめという点に、少々モヤモヤが残ります。

で、こちらのサイトより、不均一分散のロバスト標準誤差の言及を引用させていただきます。

大標本の場合、ロバスト標準誤差の意義はありそうです。

では不均一分散データ（小標本）のロバスト標準誤差を算出しましょう。
statsmodels の ols で 引数 cov_type を指定します。

### 不均一分散頑健的推定
# 参考サイト：不均一分散頑健的推定 https://py4etrics.github.io/14_Hetero.html
# 　cov_type: HC0, HC1, HC2, HC3がある
# 係数の推定値の標準誤差が通常の回帰の結果よりも大きくなっている

result_robust = smf.ols(formula=formula, data=df6).fit(cov_type='HC1', use_t=True)
display(result_robust.summary())

【実行結果】
通常の回帰の結果、$${X}$$の係数$${0.0587}$$、標準誤差$${0.028}$$と比べて、ロバスト標準誤差の方は$${X}$$の係数$${0.5087}$$（同じ値）、標準誤差$${0.029}$$（大きな値）となりました。

■ 不均一分散の対応２：加重最小二乗法
不均一分散の原因が分かっている場合に加重最小二乗法が有効、とテキストは紹介しています。

誤差が説明変数の大きさに比例して大きくなることが分かっている場合、目的変数・説明変数ともに説明変数で割ると均一分散になるようです。
テキストの式をお借りします。

上の２番目の式を利用して、加重最小二乗法を実行します。
statsmodels の ols を利用します。

### 加重最小二乗法の実行（書籍の式：Y/X = a 1/X + b + u/X を使用）
# 定数項⇒Xの係数b、1/X項⇒切片a になります
# 係数の推定値の標準誤差が小さくなってt値が大きくなることを確認すること

result_weighted2 = smf.ols(formula='I(Y/X) ~ I(1/X)', data=df6).fit()
display(result_weighted2.summary())

【実行結果】

statsmodels の 加重最小二乗法 wls を利用してみましょう。
wls では重みが誤差分散の逆数に比例することが必要らしく、引数 weights に説明変数$${X}$$の二乗の逆数$${1/X^2}$$を設定しました。

### 加重最小二乗法の実行 （不均一分散の原因が分かっている）
# WLS では重みが誤差分散の逆数に比例することが必要
# 係数の推定値の標準誤差が小さくなってt値が大きくなることを確認すること

result_weighted = smf.wls(formula=formula, data=df6, weights=1/df6.X**2).fit()
display(result_weighted.summary())

【実行結果】
こちらは切片、傾きともに正しい位置に表示されています。

先の通常の回帰で加重最小二乗法を実行したときの推定値と同じ結果を得られました。
以下は先の通常の回帰で加重最小二乗法を実行したときの推定値です。

不均一分散の原因が分かっているケースが多いかは不明ですが、分かっている場合には加重最小二乗法を使いたいです。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。