【統計学】分散不均一だと何が問題なのか|不偏性とガウスマルコフ性について
こんにちは、青の統計学です。
今回は、分散均一と分散不均一について解説いたします。
推定量期待値の分散に関わる問題で、検定方法についても触れようと思います。
各種検定のチートシートは以下をクリック！
【最短合格】統計検定準一級のチートシート|難易度や出題範囲について
分散均一と分散不均一(Dispersion homogeneity and dispersion heterogeneity)
分散均一とは、誤差項の分散が均一という意味ですが、主に最小二乗法(OLS)によって出力された推定量が不偏性を持つための条件として登場します。
母集団モデルを単回帰または重回帰モデルでかける 誤差項の平均は0(平均独立) 説明変数間に完全共線性がない 誤差項は互いに無相関 誤差項の分散は均一
上から4つの仮定が満たされると、下のOLS推定量は不偏推定量となります。
$$β = (X’X)^{-1}X’y$$
不偏推定量とは
不偏推定量とは、真のパラメータを期待値として推定する性質を持つ推定量です。
$$E[\overline{X}]=μ$$
【n-1で割る理由】不偏分散と不偏性についてわかりやすく解説
しかし、5つ目の仮定が満たされる場合、OLS推定量は最小分散不偏推定量となります。
不偏性を満たし、かつ分散が最小ということですね。
「最小」というのは、すべての推定量の中で最も分散が小さい(=効率性がある)という意味です。
最良線形不偏推定量とは少し違うので注意してください。
【MSEを最小化】ガウス・マルコフの定理と最良線形不偏推定量について
分散均一は重回帰分析の場合、以下のように表せます。
$$V[u_i | X_1,X_2,..x_k]=V[u_i]=s_i$$
このように、説明変数と独立に誤差項が決定されています。
誤差項の分散って何に使うの？
そもそも説明変数ではない誤差項にどんな役割があるのでしょうか。
実は回帰係数の分散を求める際に、誤差項の分散を利用します。
$$V[\hat{β_1}|x_1,x_2,…x_n]=\frac{s^2}{\sum_{i=1}^n(x_i-\overline{x})^2}$$
ただし、誤差項の分散sなんて直接求めることはできません。
分散均一の場合、以下のように誤差項の分散は表すことができます。
$$\hat{s}=\frac{1}{n-2}\sum_{i=1}^{n}\hat{u}_{i}^2$$
分散を直接求めることはできないので残差(予測値と実際の観測データの差)を使っているということです。
分散不均一だと何が問題なのか？
仮定4つが満たされていると推定量の不偏性は担保されていますが、5つ目の分散不均一性がないとどう困るのでしょうか？
効率性の低下
先ほど簡単に触れましたが、分散不均一だと「効率性」が保証されません。
分散不均一は、OLS推定量の効率性を低下させます。
$$V[\hat{θ}]≦V[\hat(θ)*]$$
上のような、どの不偏推定量の分散よりも分散が小さい場合、「効率性がある」と評価されます。
効率性とは、推定量の精度のことで、効率性が低下すると推定量の分散が大きくなります。
結果として、推定結果が不確かになります。
以下のコンテンツでは、推定量の分散の下限を設定できる「クラメールラオの下限」について言及しています。
【統計検定】フィッシャー情報量とクラメール・ラオの不等式について解説|python
標準誤差の過小評価
分散不均一がある場合、OLS推定量の標準誤差が過小評価される可能性があります。
標準誤差は推定量の信頼性を示す指標であり、過小評価されると、信頼区間や仮説検定の結果が不正確になります。
また、先ほど紹介した「最小二乗推定量の分散の求め方」を分散不均一な場合には適用できません。
$$V[\hat{β_1}|x_1,x_2,…x_n]=\frac{s^2}{\sum_{i=1}^n(x_i-\overline{x})^2}$$
そもそも不均一だとsがもとまりませんね。
この時に使うのが残差で、残差の二乗 \(u_{i}^2\) を計算し、これを誤差項の分散の推定値とします。
つまり、分散不均一な場合の最小二乗推定量の分散の求め方は以下のようになります。
$$Var(\hat{β}) = \frac{\sum_{i=1}^n(u_{i}^2 * x_{i}^2)} {\sum_{i=1}^{n}(x_i – \overline{x})^2}^2$$
標準誤差については以下のコンテンツで詳細にまとめています。
分散不均一に頑健な標準誤差を求めよう
分散不均一に対して頑健な標準誤差の一つとして、ホワイトの標準誤差というものがあります。
手順を示します。Xは説明変数を格納した行列です。
1：OLS（Ordinary Least Squares）によってパラメータを推定し、残差 \(\hat{ε}\) を計算します。
いうまでもないかもしれませんが、ホワイトの標準誤差は線形回帰で使用可能です。
2：残差 \(\hat{ε}\) の二乗と説明変数 X の外積の行列を計算します。
各観測値について、この行列の要素は以下のようになります。
これらの行列の和を計算し、それを S とします。
$$S = \sum \hat{ε}_i^2 * (X_i X^T_i)]$$
この外積和を使って分散不均一に頑健な共分散行列を計算します。
$$V(\hat{β}) = (X^TX)^{-1} * S * (X^TX)^{-1}$$
最後にホワイトの頑健な標準誤差は、この共分散行列 \(V(\hat{β})\) の対角要素の平方根として計算されます
対角要素なので、ただの分散です。
$$\pmatrix{ V(X_1)&COV(X_1,X_2)&…&COV(X_1,X_n) \\ COV(X_2,X_1)&V(X_2)&…&COV(X_n,X_1) \\ COV(X_n,X_1)&COV(X_n,X_2)&…&V(X_n)}$$
例えば、共分散行列\(V(\hat{β})\) のある対角要素を\(v_{ii}\)とすると
$$SE(\hat{β}_i) = \sqrt{v_{ii}}$$
では、共分散行列を使うことで頑健になるのでしょうか？
対角行列つまり分散の平方根（つまり標準偏差）は、値のばらつきの度合いを表す統計的な尺度であり、元のデータの単位と同じ単位で解釈できるのが特徴です。
よって共分散行列の対角要素の平方根を取ることで、各パラメータの推定値のばらつきの度合い（つまり不確実性）を直感的に解釈可能な形で表現することができます。
頑健な標準誤差を使うことで分散不均一が回避されるということではありません。
標準誤差に、分散不均一という情報を織り込むことで標準誤差が過小に評価されることを防いでいるというわけです。
CODE|ブルーシュ=ペーガン検定とホワイト検定
これまで分散不均一の問題点や特徴を解説してきましたが、最後に分散不均一かどうかを判定することが重要でしょう。
ここでは、ブルーシュ-ペーガン検定とホワイト検定について解説します。
ブルーシュ=ペイガン検定
回帰分析の残差の二乗を従属変数とする新たな回帰モデルを設定します。
$$Y = Xβ + ε$$
$$ε^2 = Zγ + u$$
元の回帰式の誤差項の二乗が目的変数になっていますね。
Zは通常Xと同じ説明変数を格納した行列ですが、説明変数が誤差項のうちどのくらいを説明できるかを考えます。
次に、\(ε^2 = Zγ + u\)の決定係数を使って、LM統計量(Lagrange Multiplier Statistics)を\(nR^2\)を計算し、カイ二乗分布を用いてp値を計算します。
帰無仮説が棄却されると分散不均一と言えます。
→棄却できないことが望ましい。
ホワイト検定
誤差の分散が説明変数の非線形関数として表現されるようになっております。
具体的には、回帰分析の残差の二乗を従属変数とする新たな回帰モデルを設定しますが、このとき説明変数(Z*)としては、元の説明変数、その二乗項、および交互作用項を含むことが一般的です。
$$ε^2 = Z*δ + v$$
ここがブルーシュ ペーガンの検定との違いですね。表現力が増しました。
次に、この新たな回帰モデルのF統計量（または\(χ^2\)統計量）を計算し、帰無仮説（誤差の分散が一定である）を検証します。
帰無仮説が棄却されると、分散不均一の存在が示唆されます。
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_white
np.random.seed(42) # 乱数のシードを設定して再現性を保証します
# データを生成
n = 100
x1 = np.random.randn(n)
x2 = np.random.randn(n)
error_variance = np.exp(x1)
error = np.random.randn(n) * np.sqrt(error_variance)
y = 3 + 2 * x1 + 0.5 * x2 + error
data = pd.DataFrame({"y": y, "x1": x1, "x2": x2})
# 元の回帰モデルを推定
X = sm.add_constant(data[["x1", "x2"]])
model = sm.OLS(y, X).fit()
# ホワイト検定を実行
white_test_stat, white_test_p_value, _, _ = het_white(model.resid, model.model.exog)
print("White test statistic:", white_test_stat)
print("White test p-value:", white_test_p_value)
ホワイト検定を実装するコードを紹介します。
White test statistic: 21.12266604682107 White test p-value: 0.0007679500542173898
p_valueは相当小さくなりましたね。
通常の有意水準(??)であれば、有意なので分散不均一だと言えそうです。
つまり、設定した回帰式から生まれた推定パラメータの分散は、まだまだ小さくできる余地があるということです。
有意水準とp値について復習したい方はこちらをご覧ください。
