線形回帰#
本稿では、基本的な回帰モデルの一つである線形回帰について、その設定や基本的性質・仮説検定の方法をscikit-learnやstatsmodelsによる実行例を交えつつ解説します。
モデルの概要#
線形回帰では、目的変数\(y \in \mathbb{R}\)が説明変数\(\mathbf{x} = (x_1, x_2, \dots, x_P) \in \mathbb{R}^P\)の線形結合に定数項\(b_0 \in \mathbb{R}\)と平均0、分散\(\sigma^2\)の正規分布に従うノイズ \(e \sim \mathcal{N}(0, \sigma^2)\) が付加されたものである、すなわち
というモデルに従うと仮定します。線形回帰では、観測された\(N\)個のデータ \(D = \left\{ (\mathbf{x}_i, y_i) \right\}_{i=1}^N\) に最もよく適合するような係数 \(\mathbf{b} = (b_1, b_2, \dots, b_P, b_{P+1})^T \in \mathbb{R}^{P+1}\) を最小二乗法により推定します。
また、本稿では線形回帰において一般的な以下の仮定もおくこととします：
ノイズ \(e\) はサンプルごとに独立
説明変数 \(\mathbf{x}\) とノイズ \(e\) は無相関
最小二乗法による係数の推定#
データ \(D = \left\{ (\mathbf{x}_i, y_i) \right\}_{i=1}^N\) の \(i\) 番目のサンプルの説明変数が \(\mathbf{x}_i=(x_{i, 1}, x_{i, 2}, \dots, x_{i, P}) \in \mathbb{R}^P\) 、目的変数が \(y_i \in \mathbb{R}\) であるとき、最小二乗法ではこのデータ \(D\) に対する二乗和誤差
が最小になるような \(\mathbf{b} \in \mathbb{R}^{P+1}\) を係数の推定値とします。
注釈
ノイズ \(e\) が平均0、分散 \(\sigma^2\) の正規分布に従う場合、最小二乗法と最尤推定は等価になります。このことは、対数尤度 \(\ell (\mathbf{b}; D)\) が
とも表されるため対数尤度 \(\ell (\mathbf{b}; D)\) の最大化と二乗和誤差 \(E(\mathbf{b})\) の最小化が等価になることから確認できます。
最小二乗法による係数の推定値 \(\hat{\mathbf{b}}\) は、行列を用いて解析的に書くことができます。 実際、
とおくと二乗和誤差 \(E(\mathbf{b})\) は
と書くことができますので、 \(E(\mathbf{b})\) の \(\mathbf{b}\) についての勾配 \(\nabla_\mathbf{b}E\) は
となります。ゆえに最小二乗法による係数の推定値 \(\hat{\mathbf{b}}\) を求めるには \(\nabla_\mathbf{b}E = \mathbf{0}\) 、すなわち
を解けばよいことになります。この方程式(2)は正規方程式と呼ばれており、その解は \(\mathbf{X}^T \mathbf{X}\) が正則であれば
と表すことができます。最小二乗法では、この \(\hat{\mathbf{b}}\) を係数の推定値とします。この \(\hat{\mathbf{b}}\) は比較的緩い仮定の下で不偏性や一致性などの望ましい性質を有していますが、それらの性質の説明は一旦末尾の付録に譲ることにし、先にPythonを用いた実行例の紹介を行うことにします。
注釈
線形回帰モデル(1)は、目的変数を任意の基底関数 \(\phi_j: \mathbb{R}^P \to \mathbb{R} ~ (j=1, 2, \dots, M)\) の線形結合
で表現する線形基底関数モデルに拡張することが可能です。
これは、基底関数を \(\phi_p(\mathbf{x})=x_p \quad (p=1, 2, \dots, P), \quad \phi_{P+1}(\mathbf{x})=1\) とすれば線形回帰モデル(1)と一致します。
線形基底関数モデルの係数 \(\mathbf{w} = (w_1, w_2, \dots, w_M)^T \in \mathbb{R}^M\) の推定は、 \(\mathbf{X}\) を
に置き換えたのち、線形回帰モデルと同様の方法で係数推定を行うことが可能です。すなわち、
となります。
実行例#
本節ではscikit-learnおよびstatsmodelsを用いた線形回帰の実行例を紹介します。ここではまずscikit-learnを用いた係数 \(\mathbf{b}\) の推定および予測の実行例を紹介し、その後statsmodelsを用いたt検定・F検定の実行例を紹介します。なお、その他の検定や理論的背景については末尾の付録にていくつか補足していますので、そちらもご参照ください。
準備#
# 必要なライブラリのインポート
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import sklearn
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.linear_model import LinearRegression
import statsmodels.api as sm
print(f"NumPy version: {np.__version__}")
print(f"pandas version: {pd.__version__}")
print(f"scikit-learn version: {sklearn.__version__}")
print(f"matplotlib version: {matplotlib.__version__}")
print(f"seaborn version: {sns.__version__}")
print(f"statsmodels version: {sm.__version__}")
NumPy version: 1.26.4
pandas version: 2.2.3
scikit-learn version: 1.5.2
matplotlib version: 3.9.2
seaborn version: 0.13.2
statsmodels version: 0.14.4
データの読み込み#
ここでは、scikit-learnに付属のCalifornia housing datasetを用いて線形回帰を実行してみます。こちらのデータは、1990年の米国国勢調査をもとにしたものであり、各行が国勢調査における1地区に対応しています。こちらのデータには下記の9カラムが存在しますが、今回は住宅価格の中央値（MedHouseVal）をそれ以外の8カラムから予測する回帰タスクを行います。
カラム |
説明 |
|---|---|
MedInc |
世帯所得の中央値（単位: 1万ドル） |
HouseAge |
住宅の築年数（単位: 年） |
AveRooms |
住宅の部屋数の平均 |
AveBedrms |
住宅の寝室数の平均 |
Population |
居住人数の合計 |
AveOccup |
世帯人数の平均 |
Latitude |
緯度 |
Longitude |
経度 |
MedHouseVal |
住宅価格の中央値（単位: 10万ドル） |
# データの読み込み
from sklearn.datasets import fetch_california_housing
dataset = fetch_california_housing(as_frame=True)
exp_cols = dataset.feature_names
tar_col = dataset.target_names[0]
df = dataset["frame"]
display(df)
| MedInc | HouseAge | AveRooms | AveBedrms | Population | AveOccup | Latitude | Longitude | MedHouseVal | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 8.3252 | 41.0 | 6.984127 | 1.023810 | 322.0 | 2.555556 | 37.88 | -122.23 | 4.526 |
| 1 | 8.3014 | 21.0 | 6.238137 | 0.971880 | 2401.0 | 2.109842 | 37.86 | -122.22 | 3.585 |
| 2 | 7.2574 | 52.0 | 8.288136 | 1.073446 | 496.0 | 2.802260 | 37.85 | -122.24 | 3.521 |
| 3 | 5.6431 | 52.0 | 5.817352 | 1.073059 | 558.0 | 2.547945 | 37.85 | -122.25 | 3.413 |
| 4 | 3.8462 | 52.0 | 6.281853 | 1.081081 | 565.0 | 2.181467 | 37.85 | -122.25 | 3.422 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 20635 | 1.5603 | 25.0 | 5.045455 | 1.133333 | 845.0 | 2.560606 | 39.48 | -121.09 | 0.781 |
| 20636 | 2.5568 | 18.0 | 6.114035 | 1.315789 | 356.0 | 3.122807 | 39.49 | -121.21 | 0.771 |
| 20637 | 1.7000 | 17.0 | 5.205543 | 1.120092 | 1007.0 | 2.325635 | 39.43 | -121.22 | 0.923 |
| 20638 | 1.8672 | 18.0 | 5.329513 | 1.171920 | 741.0 | 2.123209 | 39.43 | -121.32 | 0.847 |
| 20639 | 2.3886 | 16.0 | 5.254717 | 1.162264 | 1387.0 | 2.616981 | 39.37 | -121.24 | 0.894 |
20640 rows × 9 columns
ここでは全データのうち80%を訓練用（X_train, y_train）、残りの20%をテスト用（ X_test, y_test ）とします。
X_train, X_test, y_train, y_test = train_test_split(
df[exp_cols], df[tar_col], test_size=0.2, shuffle=True, random_state=0
)
scikit-learnによる係数推定#
scikit-learnでは、LinearRegressionクラスを用いて線形回帰を行うことができます。デフォルトでは定数項ありのモデルとなります。定数項なしのモデルを利用したい場合は、インスタンス生成時にfit_intercept=Falseを指定する必要があります。（デフォルトではfit_intercept=True）
lr = LinearRegression()
lr.fit(X_train, y_train) # 係数の推定
LinearRegression()In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
LinearRegression()
fitメソッドを実行すると係数の推定が行われ、係数および定数項の推定値がそれぞれcoef_, intercept_属性に格納されます。
# 推定された係数・定数項の確認
print(f"係数: {lr.coef_}")
print(f"定数項: {lr.intercept_}")
係数: [ 4.33333407e-01 9.29324337e-03 -9.86433739e-02 5.93215487e-01
-7.56192502e-06 -4.74516383e-03 -4.21449336e-01 -4.34166041e-01]
定数項: -36.85856910680128
推定した係数で予測を行いたい場合はpredictメソッドを使用します。下では、テストデータについて予測値を計算し、平均二乗誤差を計算しています：
# 予測の実行・MSEの算出
pred_sklearn = lr.predict(X_test)
mse_sklearn = mean_squared_error(y_test, pred_sklearn)
print(f"MSE (sklearn): {mse_sklearn}")
MSE (sklearn): 0.528984167036721
statsmodelsによる推定#
statsmodelsでは、statsmodels.api（実行例のsmに対応）内のOLSクラスを用いて線形回帰を行うことができます。この方法の場合、デフォルトでは定数項なしのモデルとなるので、定数項ありのモデルを利用する場合はsm.add_constant()を用いて定数項に対応するカラムを説明変数に追加する必要がある点にご注意ください。
# 定数項に対応するカラムを追加
X_train_aug = sm.add_constant(X_train, prepend=False)
model = sm.OLS(y_train, X_train_aug)
# 係数の推定
result_sm = model.fit()
上のコードではfitメソッドで推定を行い、推定結果をresult_sm（RegressionResults型）に格納しています。推定された係数はresult_smのparams属性にpandas.Series型で格納されています。
print(type(result_sm.params))
print(result_sm.params)
<class 'pandas.core.series.Series'>
MedInc 0.433333
HouseAge 0.009293
AveRooms -0.098643
AveBedrms 0.593215
Population -0.000008
AveOccup -0.004745
Latitude -0.421449
Longitude -0.434166
const -36.858569
dtype: float64
学習したモデルで予測を行う際は、predictメソッドを使用します。定数項ありのモデルを利用している場合は、学習時と同様にsm.add_constant()で定数項に対応するカラムを説明変数に追加する必要がある点にご注意ください。下では、scikit-learnの時と同様にテストデータについて予測値を計算し、平均二乗誤差を計算しています：
# 定数項に対応するカラムを追加
X_test_aug = sm.add_constant(X_test, prepend=False)
# 予測の実行・MSEの算出
pred_sm = result_sm.predict(X_test_aug)
mse_sm = mean_squared_error(y_test, pred_sm)
print(f"MSE (statsmodels): {mse_sm}")
MSE (statsmodels): 0.52898416703672
ここで、上のscikit-learnでの予測値・MSEを比較してみましょう。
print(f"sklearn: {pred_sklearn[:5]}")
print(f"statsmodels: {np.array(pred_sm[:5])}")
print(f"MSE (sklearn): {mse_sklearn}")
print(f"MSE (statsmodels): {mse_sm}")
sklearn: [2.28110738 2.79009128 1.90332794 1.01760331 2.94852425]
statsmodels: [2.28110738 2.79009128 1.90332794 1.01760331 2.94852425]
MSE (sklearn): 0.528984167036721
MSE (statsmodels): 0.52898416703672
数値計算によると思われる誤差はあるものの、両者の結果が一致していることが確認できます。
最後に、summaryメソッドを利用して推定結果の要約を出力してみます：
# 推定結果の要約の出力
print(result_sm.summary())
OLS Regression Results
==============================================================================
Dep. Variable: MedHouseVal R-squared: 0.609
Model: OLS Adj. R-squared: 0.609
Method: Least Squares F-statistic: 3212.
Date: Fri, 20 Dec 2024 Prob (F-statistic): 0.00
Time: 16:21:22 Log-Likelihood: -18085.
No. Observations: 16512 AIC: 3.619e+04
Df Residuals: 16503 BIC: 3.626e+04
Df Model: 8
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
MedInc 0.4333 0.005 92.651 0.000 0.424 0.443
HouseAge 0.0093 0.000 18.651 0.000 0.008 0.010
AveRooms -0.0986 0.007 -15.013 0.000 -0.112 -0.086
AveBedrms 0.5932 0.031 19.068 0.000 0.532 0.654
Population -7.562e-06 5.26e-06 -1.437 0.151 -1.79e-05 2.75e-06
AveOccup -0.0047 0.001 -5.404 0.000 -0.006 -0.003
Latitude -0.4214 0.008 -52.340 0.000 -0.437 -0.406
Longitude -0.4342 0.008 -51.498 0.000 -0.451 -0.418
const -36.8586 0.737 -49.994 0.000 -38.304 -35.413
==============================================================================
Omnibus: 3453.723 Durbin-Watson: 1.996
Prob(Omnibus): 0.000 Jarque-Bera (JB): 10941.189
Skew: 1.067 Prob(JB): 0.00
Kurtosis: 6.369 Cond. No. 2.40e+05
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.4e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
推定結果の要約にはサンプル数（No. Observations）や決定係数（R-squared）の他、t検定・F検定の検定統計量およびp値なども含まれています。次節以降では、これらの検定統計量の見方について解説します。
t検定#
線形回帰におけるt検定は、 \(p\) 番目の説明変数が目的変数の予測に寄与するかを検定します。すなわち帰無仮説 \(H_0\) 、対立仮説 \(H_1\) をそれぞれ
とする仮説検定です。
ノイズ \(e\) がサンプルごとに独立に平均 \(0\) 、分散 \(\sigma^2\) の正規分布 \(\mathcal{N}(0, \sigma^2)\) に従うと仮定し、 \((\mathbf{X}^T \mathbf{X})^{-1}\) の第 \((p, p)\) 成分を \(C_p\) とおくと検定統計量
は帰無仮説の下で自由度 \(N-P-1\) のt分布 \(t(N-P-1)\) に従うことが知られています。
上で出力したstatsmodelsの推定結果の要約では真ん中のブロックのtが各成分の検定統計量、P>|t|がそのp値となっています。例えば有意水準を5%としてMedInc（世帯所得）に注目すると、p値が0.000と5%を下回っていることから帰無仮説が棄却され、世帯所得が住宅価格の予測に寄与していると判定されます。一方で、Population（居住人数）に注目すると、p値が0.151と5%を上回っていることから、帰無仮説は棄却されず居住人数は住宅価格の予測に寄与しているとはいえない、ということになります。
検定の多重性
t検定に限らず、一般に統計的仮説検定では仮に帰無仮説が正しい場合であっても帰無仮説が棄却される第一種の過誤が有意水準 \(\alpha\) と同じ確率で発生します。統計的仮説検定を何度も行うと、少なくとも1つの検定で帰無仮説が棄却される確率は \(\alpha\) よりも大きくなり、特に検定を行う回数が増えれば増えるほどこの影響はより顕著になります。この問題は検定の多重性と呼ばれます。
例えば線形回帰におけるt検定では説明変数が多い場合、仮にどの説明変数も目的変数と全く関係なかったとしても、いずれかの変数についてt検定の帰無仮説が棄却されてしまい「実際には目的変数とは全く関係のない説明変数が、重要であると誤解される」ことにつながります。線形回帰で仮説検定を行う際はこういった検定の多重性の問題があることにご注意ください。また、実際の分析の際には統計的仮説検定の結果を過信せず、ドメイン知識などの他の要素とも合わせて複合的に判断することも重要です。
F検定#
線形回帰において単にF検定というと、一般には「少なくとも1つの説明変数が目的変数の予測に寄与するか」、すなわち帰無仮説 \(H_0\) 、対立仮説 \(H_1\) を
とする仮説検定のことを指します。なお、上の帰無仮説 \(H_0\) においては定数項 \(b_{P+1}\) が \(0\) であることまでは要請していない点にご注意ください。
この場合、各サンプルに対する予測値 \(\hat{y}_i = \sum_{p=1}^P b_p x_{i, p} + b_{P+1}\) および目的変数の平均値 \(\bar{y} = \frac{1}{N}\sum_{i=1}^N y_i\) を用いて検定統計量 \(F\) を
と定義すると、 \(F\) は帰無仮説 \(H_0\) の下で自由度 \((P, N-P-1)\) のF分布に従うことが知られています。
上で出力したstatsmodelsの推定結果の要約では上のブロックのF-statisticが検定統計量、Prob (F-statistic)がそのp値となっています。例えば有意水準を5%とすると、今回の結果ではp値が0.00と5%を下回っていることから帰無仮説が棄却され、いずれかの説明変数が予測に寄与していると判定されます。
参考
線形回帰におけるF検定では帰無仮説・対立仮説を(3)とすることが多いですが、より一般に \(\mathbf{C} \in \mathbb{R}^{K \times (P+1)}, ~ \mathbf{r} \in \mathbb{R}^K\) を用いて定義される帰無仮説
に対する仮説検定を行うことも可能です。この形のF検定の詳細については、付録をご参照ください。
小括#
本稿では、基本的な回帰モデルの一つである線形回帰について、係数推定の方法や各種仮説検定の方法をscikit-learnやstatsmodelsによる実行例も交えつつ解説しました。本文では解説を割愛した要素、例えば係数の推定値\(\hat{\mathbf{b}}\)の性質やその他の検定については下の付録にて補足していますので、そちらも適宜ご参照ください。
参考文献#
付録#
\(\hat{\mathbf{b}}\) の基本的な性質#
注意
本節で紹介する \(\hat{\mathbf{b}}\) の性質については、ノイズ \(e\) が正規分布に従わない場合でも成り立ちます。なお、次節で述べるF検定については \(e\) の正規性を仮定しているので、それぞれの仮定の違いにご留意ください。
データがモデル(1)に従い、かつノイズ \(\mathbf{e} \in \mathbb{R}^N\) の平均・分散について
が成り立つ場合、 \(\hat{\mathbf{b}}\) の期待値および分散共分散行列は下のようになります：
\(\mathbb{E}[\hat{\mathbf{b}}] = \mathbf{b}\)
\(\text{Var} [\hat{\mathbf{b}}] = \sigma^2 (\mathbf{X}^T \mathbf{X})^{-1}\)
実際、期待値については
と \(\mathbb{E}[\mathbf{e}] = \mathbf{0}\) から従います。
また、分散については
と \(\mathbb{E}[ \mathbf{e} \mathbf{e}^T] = \sigma^2 \mathbf{I}_N\) から従います。
期待値についての結果に注目すると、推定量 \(\hat{\mathbf{b}}\) の期待値が真値 \(\mathbf{b}\) と一致していることから、 \(\hat{\mathbf{b}}\) が \(\mathbf{b}\) の不偏推定量であることが確認できます。
さらに、下記のGauss-Markovの定理によれば\(\hat{\mathbf{b}}\)の分散は不偏な線形推定量の中では最小である、すなわち \(\hat{\mathbf{b}}\) が最良線形不偏推定量（BLUE） であることが知られています。
Gauss-Markovの定理
\(\mathbb{E}[\mathbf{e}] = \mathbf{0} \in \mathbb{R}^N\) かつ \(\text{Cov}[\mathbf{e}]=\sigma^2 \mathbf{I}_N\) であるならば、 \(\hat{\mathbf{b}}\) は最良線形不偏推定量（BLUE）である。
注釈
Gauss-Markovの定理の仮定は、ノイズの平均が0であり、等分散でかつ無相関であることとなっています。すなわち、ノイズの各成分が独立でない場合や、ノイズが正規分布に従わない場合でもGauss-Markovの定理の主張は成り立ちます。
なお、線形推定量・最良線形不偏推定量（BLUE）の定義は以下の通りです：
線形推定量：ある行列 \(C\) を用いて \(\tilde{\mathbf{b}}=\mathbf{C}\mathbf{y}\) の形で表される、すなわち \(\mathbf{y}\) について線形である推定量 \(\tilde{\mathbf{b}}\) （\(\mathbf{x}\) については非線形であってもよい）
最良線形不偏推定量（BLUE）：不偏な線形推定量のうち、分散が最小となるもの。すなわち、任意の線形不偏推定量 \(\tilde{\mathbf{b}}'\) について \(\text{Var}[\tilde{\mathbf{b}}'] - \text{Var}[\tilde{\mathbf{b}}]\) が半正定値対称行列であるとき、 \(\tilde{\mathbf{b}}\) は最良線形不偏推定量（BLUE）であるという
また、\(\hat{\mathbf{b}}\) は比較的緩い仮定の下で一致性を持つことが知られています
\(\lim_{n \to \infty} \hat{\mathbf{b}} = \mathbf{b}\)（一致性）
一致性は、 \(\mathbf{x}, \mathbf{e}\) の無相関性および大数の法則より \(\hat{\mathbf{b}} - \mathbf{b}=(\frac{1}{N} \mathbf{X}^T \mathbf{X})^{-1} \frac{1}{N} \mathbf{X}^T \mathbf{e} \to 0\) となることから従います。
一般化されたF検定#
実行例の節では係数の各成分についてのt検定と一般的なF検定について説明しましたが、ここではより一般化されたF検定について、その導出・証明も含めて解説します。
線形回帰におけるF検定では帰無仮説・対立仮説を(3)式のように設定することが多く、statsmodelsの推定結果の要約で表示されていたのも(3)についての検定統計量・p値でしたが、 \(\mathbf{C} \in \mathbb{R}^{K \times (P+1)}, ~ \mathbf{r} \in \mathbb{R}^{K}\) を用いた下のような帰無仮説・対立仮説に一般化することが可能です：
この場合、検定統計量 \(F\) を
により定義すると、 \(F\) は帰無仮説 \(H_0 : \mathbf{C} \mathbf{b} = \mathbf{r}\) の下で自由度 \((K, N-P-1)\) のF分布 \(F(K, N-P-1)\) に従うことが知られています。ここで、 \(s^2\) は \(\sigma^2\) のOLS推定量
です。
以降は上記の \(F\) が実際に自由度 \((K, N-P-1)\) のF分布 \(F(K, N-P-1)\) に従うことを示します。
\(\mathbf{P} := \mathbf{X}(\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T\) とおくと \(\mathbf{P}\) は \(\mathbf{X}\) の列ベクトルの張る部分空間への直交射影行列であり、
と表すことができます。ゆえに
が成り立ちます。
帰無仮説 \(\mathbf{C}\mathbf{b}=\mathbf{r}\) の下で
ですので、 \(\mathbf{W} = \mathbf{X} (\mathbf{X}^T \mathbf{X})^{-1}\mathbf{C}^T \in \mathbb{R}^{N \times K}\) とおくと
となります。今、
と定義するとこの \(\mathbf{Q}\) は \(\mathbf{W} \in \mathbb{R}^{N \times K}\) の列ベクトルの張る \(K\) 次元部分空間 \(V_\mathbf{Q} \subset \mathbb{R}^N\) への直交射影行列となっています。
\(\mathbf{W} = \mathbf{X} (\mathbf{X}^T \mathbf{X})^{-1}\mathbf{C}^T \in \mathbb{R}^{N \times K}\) の列ベクトルは \(\mathbf{X}\) の列ベクトルの線形結合となっているので、 \(V_\mathbf{Q}\) は \(V_\mathbf{P}\) の部分空間です。従って、 \(V_\mathbf{R} := V_\mathbf{P} \cap V_\mathbf{Q}^\perp\) への直交射影行列を \(\mathbf{R}\) とおくと
が成り立ちます。ゆえに下記のCochranの定理から \(\mathbf{e}^T\mathbf{Q}\mathbf{e}/\sigma^2, \mathbf{e}^T(\mathbf{I} - \mathbf{P})\mathbf{e} / \sigma^2\) はそれぞれ独立に自由度 \(K, N-P-1\) のカイ二乗分布に従うことが分かります。
Cochranの定理
\(\mathbf{x} \in \mathbf{R}^N\) が多変量正規分布 \(\mathcal{N}(\mathbf{0}, \mathbf{I}_N)\) に従っているとする。階数がそれぞれ \(r_1, r_2, \dots, r_m\) の行列 \(\mathbf{A}_1, \mathbf{A}_2, \dots, \mathbf{A}_m \in \mathbb{R}^{N \times N}\) が存在して \(\sum_{i=1}^m \mathbf{A}_i =\mathbf{I}_N\) かつ \(r_1+r_2+\cdots +r_m=N\) が成り立つならば、 \(\mathbf{x}^T \mathbf{A}_1 \mathbf{x}, \mathbf{x}^T \mathbf{A}_2 \mathbf{x}, \dots, \mathbf{x}^T \mathbf{A}_m \mathbf{x}\) はそれぞれ独立に自由度 \(r_i\) のカイ二乗分布 \(\chi^2(r_i)\) に従う。
従って
は自由度 \((K, N-P-1)\) のF分布 \(F(K, N-P-1)\) に従います。
次に、(5)で定義される \(F\) が上で述べた(3)の帰無仮説・対立仮説についての検定統計量(4)と一致することを確認します。まず、(3)の帰無仮説・対立仮説は
の場合に対応しています。ここで、 \(\mathbf{X} = \begin{pmatrix}\mathbf{X}_1 & \mathbf{X}_2\end{pmatrix}, \mathbf{X}_2 = \mathbf{1} \in \mathbb{R}^{N \times 1}, \mathbf{m}=\mathbf{X}_1^T \mathbf{1} \in \mathbb{R}^{P}\) とおくとブロック逆行列の公式から
となるので、このときのF統計量は \(\hat{\mathbf{b}}_{1:P} = (\hat{b}_1, \hat{b}_2, \dots, \hat{b}_P)^T \in \mathbb{R}^P\) を用いて
と表され、さらにこの \(F\) は(4)の \(F\) と一致します。このことは、下の式変形から確認することができます：
statsmodelsでの結果の詳細#
ここではstatsmodelsの推定結果の要約に記載されている情報のうち、本文では説明を省略した部分についていくつか補足します。
残差の正規性検定#
上のt検定やF検定では、ノイズが正規分布に従うことを仮定しています。そのため、ノイズが正規分布に従わない場合はt検定やF検定がうまく機能しなくなる恐れがあります。そこで正規性検定、すなわち「残差が正規分布に従う」という帰無仮説に対する仮説検定を行うことで、t検定やF検定の結果をどの程度信じてよいか確認することができます。statsmodelsの推定結果の要約では、残差の歪度（skewness）および尖度（kurtosis）に着目し、それらが正規分布の歪度・尖度と乖離しているか否かで検定を行う下記2つの手法についての結果が表示されています：
scipy.normaltestを用いるオムニバス検定Jarque-Bera検定
それぞれの概略を説明すると、まずオムニバス検定統計量は \(s^2+k^2\) であり、「確率変数（ここでは残差）が正規分布に従う」という帰無仮説の下でこれが自由度2のカイ二乗分布に従います。 \(s, k\) はそれぞれ
です。ここでいうZ値とは「確率変数が正規分布に従う」という帰無仮説の下で標準正規分布 \(\mathcal{N}(0, 1)\) に従うように調整した確率変数のことを指します。
もう一つの検定であるJarque-Bera検定は標本歪度 \(s=\frac{\frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^3}{\left[ \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2 \right]^{3/2}}\) および標本尖度 \(k=\frac{\frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^4}{\left[ \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2 \right]^{2}}\) を用いた
を検定統計量とする検定であり、標本が正規分布に従うという帰無仮説の下で上の \(JB\) が漸近的に自由度2のカイ二乗分布に従うことが知られています。
statsmodelsではオムニバス検定の検定統計量およびp値がそれぞれOmnibus, Prob(Omnibus)、Jarque-Bera検定の検定統計量およびp値がそれぞれJarque-Bera (JB), Prob(JB)に記載されています。今回の結果を見ると両者のp値がほぼ0となっていることから、残差が正規分布に従っているとは言い難く、t検定やF検定の結果の取り扱いには注意が必要そうだ、と判断できます。
残差の自己相関#
残差に自己相関があると線形回帰におけるノイズの独立性の仮定が満たされず、t検定やF検定がうまく機能しなくなる恐れがあります。なお、自己相関については自己相関関数と相互相関関数の記事もご参照ください。
そこで、残差の自己相関の有無を判定するために利用できる指標の一つがDurbin-Watson比です。残差が \(r_1, r_2, \dots, r_N\) であるとき、Durbin-Watson比 \(DW\) は
で定義されます。 \(DW\) は0以上4以下の値を取り、自己相関がない場合は2、正の自己相関をもつ場合は0、負の自己相関をもつ場合は4となります。
statsmodelsの推定結果の要約ではDurbin-Watson比がDurbin-Watsonに記載されています。今回の結果ではDurbin-Watson比が1.996と2に近い値なので、自己相関はなさそうであると判断できます。
説明変数の条件数#
説明変数の条件数は、多重共線性の有無を確認するために用いることができる指標の一つです。多重共線性が存在すると係数の推定値 \(\hat{\mathbf{b}}\) の分散が大きくなるなど、推定に悪影響が出てしまいます。（多重共線性の詳細や対処法については次元削減の記事もご参照ください。）
statsmodelsの推定結果の要約では計画行列 \(\mathbf{X}\) の条件数、すなわち \(\mathbf{X}\) の最大特異値 \(\sigma_\text{max}(\mathbf{X})\) と最小特異値 \(\sigma_\text{min}(\mathbf{X})\) の比
がCond. No.に記載されます。今回の結果では条件数が2.40e+05とかなり大きい値となっており、要約の末尾のNotesにも
[2] The condition number is large, 2.4e+05. This might indicate that there are strong multicollinearity or other numerical problems.
とあるように強い多重共線性が存在していたことが分かります。このような場合には、次元削減の記事にあるような変数選択をして多重共線性を回避することが推奨されます。
