(chap:14-hetero)=
import matplotlib.pyplot as plt
import lmdiag
import numpy as np
import pandas as pd
import wooldridge
from seaborn import residplot
from statsmodels.formula.api import ols
from statsmodels.stats.api import het_breuschpagan, het_white
from statsmodels.stats.outliers_influence import reset_ramsey
# 警告メッセージを非表示
import warnings
warnings.filterwarnings("ignore")
<仮定5(均一分散; Homogeneity)が満たされない場合>
仮定5の均一分散の下では説明変数は誤差項の分散に影響を与えない。即ち,
$$\text{Var}\left(u|x\right)=\sigma^2$$この仮定に関連して次の点について留意する必要がある。
仮説を検証するということを目的とすると,検定が無効というのは致命的な問題である。特に,不均一分散(Heteroskedasticity)の問題は,横断面データを使うと頻繁に出てくる問題である。ではどのように対応すればよいのか。
不均一分散頑健的推定(Heteroskedasticity-Robust Inference)
この手法を使うと,OLS推定量の標準誤差が調整され未知の不均一分散であっても,$t$検定,$F$検 定が有効になるというものである。
(理由)均一分散であっても不均一分散であっても,$n\rightarrow\infty$の場合,不均一分散頑健的推定の$t$($F$)値は$t$($F$)分布に従う。言い換えると,標本の大きさが十分に大きければ,$t$($F$)値の分布は$t$($F$)分布で近似できるということである。
更なる利点は,通常のOLS推定の後に標準誤差の調整が施され,計算はstatsmodels
を使うと簡単におこなうことが可能である。
(注意)大標本でのみ有効。
不均一分散頑健的推定では,OLS推定の共分散行列(covariance matrix)と呼ばれる箇所を調整し,OLS推定量の標準誤差を修正する。その調整方法が複数提案されていおり,statsmodels
では以下の種類に対応している。
HC0
: White (1980)の不均一分散頑健共分散行列推定HC1
: MacKinnon and White (1985)の不均一分散頑健共分散行列推定v1HC2
: MacKinnon and White (1985)の不均一分散頑健共分散行列推定v2HC3
: MacKinnon and White (1985)の不均一分散頑健共分散行列推定v3ここでHC
はH
eteroskedasticity-C
onsistent Covariance Matrix EstimatorsのH
とC
。
不均一分散頑健共分散行列推定を使い計算した$t$値や$F$値を頑健的$t$値,頑健的$F$値と呼ぶことにする。
{note}
`HC0`などは不均一分散に対応する推定方法である。一方で,時系列分析では残差に不均一分散と自己相関の両方が存在する場合がある。この両方の問題に頑健的な推定を不均一分散・自己相関頑健推定(Heteroskedasticity-Autocorrelation Robust Inference)と呼ぶ。以下で説明する引数`cov_type`に`HAC`を次のように指定することにより,不均一分散・自己相関頑健標準誤差を計算することができる。
```
cov_type='HAC'
```
ちなみに,`HC`はHeteroskedasticity Consistentの略であり,`HAC`はHeteroskedasticity-Autocorrelation Consistent`の略である。
OLS推定量の不均一分散頑健標準偏差が簡単に計算できるのであれば,通常の標準偏差を使う必要はないのではないか,という疑問が生じる。この問を関して以下の点を考える必要がある。
従って,標本の大きさが「大標本」と判断できる場合(例えば,$n\geq 1000$)以外は通常の標準偏差と不均一分散頑健標準偏差の両方を表示することを勧める。
wooldridge
パッケージのデータセットgpa3
を使い説明する。この例では大学のGPAと高校の成績や性別,人種などがどのような関係にあるかを探る。
gpa3 = wooldridge.data('gpa3').query('spring == 1') # 春学期だけを抽出
wooldridge.data('gpa3', description=True)
gpa2
に一部の変数の説明が続いている。
wooldridge.data('gpa2', description=True)
被説明変数:
cumgpa
:累積GPA説明変数
sat
:SATの成績hsperc
:高校の成績の%点(上位から)tothrs
:データ抽出時から学期までの時間?(gpa3
の定義)female
:女性ダミー変数(女性=1
)black
:人種ダミー変数(黒人=1
)white
:人種ダミー変数(白人=1
)form_ols = 'cumgpa ~ sat + hsperc + tothrs + female + black + white'
mod_ols = ols(form_ols, data=gpa3)
res_ols = mod_ols.fit()
print(res_ols.summary().tables[1])
上のOLSの結果を使い頑健$t$値を計算するために,res_ols
のメソッド.get_robustcov_results()
を使う。
cov_type
は頑健性の計算法の指定(デフォルトはCH1
)。use_t
は$t$検定を指定(デフォルトはNone
で「自動」に決められる)res_robust = res_ols.get_robustcov_results(cov_type='HC3', use_t=True)
print(res_robust.summary().tables[1])
coef
の値は同じであり,必ずそうなる。不均一分散頑健推定は,表の中で標準誤差,$t$値,$p$値,信頼区間に影響を与える。std err
を比べると,帰無仮説$\hat{\beta}_j=0$の棄却判断を覆すほど大きく変わる変数はない。これは不均一分散がそれほど大きな問題ではないことを示唆している。この点を確かめるために,res_ols
の誤差項を図示してみる。誤差項を図示する方法として2つを紹介する。一つ目は,res_ols
の属性.resid
を使う。res_ols.resid
はPandas
のSeries
(シリーズ)なので,そのメソッドplot()
を使い図示する。style
はマーカーを指定するオプション。
res_ols.resid.plot(style='o')
pass
2つ目の方法はplt.plot()
を使う。オップション'o'
はマーカの指定である。
plt.plot(res_ols.resid, 'o')
pass
3つ目はscatter()
を使う。.index
はres_ols.resd
の属性でインデックスを示す。
plt.scatter(res_ols.resid.index, res_ols.resid)
pass
OLS推定をする際,fit()
の関数に.get_robustcov_results()
で使った同じオプションを追加すると頑健$t$値などを直接出力できる。
res_HC3 = ols(form_ols, data=gpa3).fit(cov_type='HC3', use_t=True)
print(res_HC3.summary().tables[1])
同じデータgpa3
を使い,黒人ダミーと白人ダミーの係数は両方とも0
という仮説を検定する。
hypotheses = 'black = white = 0'
まず通常の$F$検定を考える。
f_test_ols = res_ols.f_test(hypotheses)
f_test_ols.summary()
返り値(左から)
F statistic
:$F$統計量F p-value
:$F$の$p$値df_denom
:分母の自由度df_num
:分子の自由度次に頑健$F$検定の方法を説明する。上の不均一分散頑健推定の方法2で使ったf_test_HC3
を使う。
f_test_HC3 = res_HC3.f_test(hypotheses)
f_test_HC3.summary()
$t$検定の場合と同じように,大きく変わる結果につながってはない。
均一分散の場合 $t$($F$)値は厳密に$t$($F$)分散に従う。それが故に,均一分散が好まれる理由である。ここでは均一分散の検定について考える。帰無仮説と対立仮説は以下となる。
$\text{H}_0$:誤差項は均一分散
$\text{H}_A$:帰無仮説は成立しない
2つの検定方法を考える。
データhprice1
を使って,住宅価格の決定要因を検討する。ここで考える均一分散の検定にBreusch-Pagan検定と呼ばれるもので,statsmodels
のサブパッケージstats
の関数het_breuschpagan
を使う。
hprice1 = wooldridge.data('hprice1')
wooldridge.data('hprice1', description=True)
以下で使う変数について。
被説明変数
price
:住宅価格(単位:1000ドル)説明変数
lotsize
:土地面積(単位:平方フィート)sqrft
:家の面積(単位:平方フィート)bdrms
:寝室の数まず変数の変換をしない場合を考える。
form_h = 'price ~ lotsize + sqrft + bdrms'
res_h = ols(form_h, data=hprice1).fit()
print(res_h.summary().tables[1])
この結果に対してBreusch-Pagan検定をおこなう。het_breuschpagan()
の引数について:
res_h
の属性.resid
で誤差項の値res_h
の属性.model
の属性exog
を使い定数項を含む説明変数の値het_breuschpagan(res_h.resid, res_h.model.exog)
返り値(上から)
LM statistic
:$LM$統計量LM p-value
:$LM$の$p$値F statistics
:$F$統計量F p-value
:$F$の$p$値$LM$検定とはLagrange Multiplier検定のことで,大標本の場合に仮定1〜4(GM仮説)のもとで成立する。一般にBreusch-Pagan検定は$LM$統計量を使ったものを指すが,$F$統計量としても計算できる。
5%の有意水準で帰無仮説($\text{H}_0$:誤差項は均一分散)を棄却でき,不均一分散の可能性が高い。対処法として変数の変換が挙げられる。
bdrms
以外を対数変換する。
form_h_log = 'np.log(price) ~ np.log(lotsize) + np.log(sqrft) + bdrms'
res_h_log = ols(form_h_log, data=hprice1).fit()
print(res_h_log.summary().tables[1])
het_breuschpagan(res_h_log.resid, res_h_log.model.exog)
5%の有意水準でお帰無仮説を棄却できない。即ち,対立仮説($\text{H}_A$:帰無仮説は成立しない)を採択し,均一分散の可能性が高い。
この検定はOLS推定量の標準誤差と検定統計量を無効にする不均一分散を主な対象としており,Breusch-Pagan検定よりもより一般的な式に基づいて検定をおこなう。statsmodels
のサブパッケージstats
の関数het_white
を使う。
hprice
のデータを使った上の例を使う。
het_white()
の引数:
res_h
の属性.resid
で誤差項の値res_h
の属性.model
の属性exog
を使い定数項を含む説明変数の値het_white(res_h.resid, res_h.model.exog)
返り値(上から)
LM statistic
:$LM$統計量LM p-value
:$LM$の$p$値F statistics
:$F$統計量F p-value
:$F$の$p$値一般にWhite検定は$LM$統計量を使ったものを指すが,$F$統計量としても計算できる。
5%の有意水準で帰無仮説($\text{H}_0$:誤差項は均一分散)を棄却でき,不均一分散の可能性が高い。対処法として変数の変換が挙げられる。
het_white(res_h_log.resid, res_h_log.model.exog)
5%の有意水準で帰無仮説を棄却できない。即ち,対立仮説($\text{H}_A$:帰無仮説は成立しない)を採択し,均一分散の可能性が高い。
仮定4〜6は残差に関するものであり,残差をプロットし不均一分散や非線形性を確認することは回帰分析の重要なステップである。
残差を図示する方法としてlmdiag
以外に以下を紹介する。
matplotlib
を直接使うseaborn
というパッケージの中にある関数residplot
を使う上で計算したres_h
とres_h_log
を利用し
.fittedvalues
).resid
)となる図を作成する。
lmdiag
¶対数変換前
plt.figure(figsize=(8,7))
lmdiag.plot(res_h)
pass
対数変換後
plt.figure(figsize=(8,7))
lmdiag.plot(res_h_log)
pass
対数変換により残差の変化がより均一的になり,Residuals vs. LeverageのCook's Distanceを見ても外れ値がなくなっている。
Matplotlib
のplot()
¶対数変換前
plt.scatter(res_h.fittedvalues, res_h.resid)
pass
対数変換後
plt.scatter(res_h_log.fittedvalues, res_h_log.resid)
pass
対数変換により残差の変化がより均一的になったのが確認できる。
seaborn
のresidplot()
¶seaborn
はmatplotlib
を利用し様々な図を描ける。seaborn
についてはこのサイトを参照。
通常import seaborn as sns
でインポートすることが多いようであるが,ここではresidplot
のみをインポートしている。
residplot()
は散布図を作成する関数である。
x=
:横軸の変数を指定y=
:縦軸の変数を指定lowerss=True
(デフォルトはFalse
)にすると,散布図にベスト・フィットする曲線を表示する。対数変換前
residplot(x=res_h.fittedvalues, y=res_h.resid, lowess=True)
pass
対数変換後
residplot(x=res_h_log.fittedvalues, y=res_h_log.resid, lowess=True)
pass
仮定1で回帰式は線形となっているが,その仮定が正しければ,上の残差の散布図は概ね平行にそして0を中心に上下等間隔に散らばることになる。そのパターンに比較的に近いのは対数変換後の図である。不均一分散の場合は,そのパターンが大きく崩れている場合であり,その原因のその1つに回帰式の特定化の間違いがある。例えば,説明変数の2乗が含まれるべきなのに欠落している場合が挙げられる。極端な場合,残差の散布図は$U$字または逆$U$字のような形になりえる。
一方,線形性の検定もある。以下ではその1つである RESET (Regression Specification Error Test) 検定を考える。
この検定の考え方はそれほど難しくはない。次の回帰式を推定するとしよう。
$$y=\beta_0+\beta_1x_1+\beta_2x_2+u$$この線形式が正しければ,$x_1^2$や$x_2^3$等を追加しても統計的に有意ではないはずである。さらに,この線形回帰式の予測値$\hat{y}$は$x_1$や$x_2$の線形になっているため、$\hat{y}^2$や$\hat{y}^3$は$x_1$や$x_2$の非線形となっている。従って,次式を推計し,もし非線形の効果がなければ$\delta_1$も$\delta_2$も有意ではないはずである。
$$y=\beta_0+\beta_1x_1+\beta_2x_2+\delta_1\hat{y}^2+\delta_2\hat{y}^3$$この考えに基づいて以下の仮説を検定する。
$\text{H}_0:\;\delta_1=\delta_2=0$(線形回帰式が正しい)
$\text{H}_A$: $\text{H}_0$は成立しない
<コメント>
statsmodels
のサブパッケージ.stats.outliers_influence
のRESET検定用の関数reset_ramsey
を使う。(ramsey
はこの検定を考案した学者名前)
reset_ramsey()
の使い方:
degree
(デフォルトは5)は$y$の何乗までを含めるかを指定する。対数変換前
reset_ramsey(res_h,degree=3).summary()
返り値
F
: $F$統計量p
: $p$値df_denom
: 分母の自由度df_num
: 分子の自由度(2乗以上の制約数)対数変換後
reset_ramsey(res_h_log,degree=3)
5%の有意水準のもとで,res_h
の帰無仮説を棄却できるが,res_h_log
では棄却できない。後者の推計式がより適しているようである。