統計学 学習ノート

はじめに

変数を増やすと、決定係数は上がるが予測は悪くなりうる

「説明変数は、全部入れれば入れるほどいいんじゃないの？」

売上を予測するモデルを作るとします。気温、曜日、広告費、近くの人口、競合の数。手元の変数を全部入れると、決定係数はどんどん上がります。これで完成でしょうか？

実は、変数を増やすほど、新しいデータでの予測はかえって悪くなることがあります。似た変数が混ざると、係数が暴れて読めなくなることもあります。2級ではソフトの出力を読むところまででしたが、準1級で問われるのは、その中身と、どの変数を残すかの判断です。

そこで登場するのが、今回の主役 重回帰モデル の行列表現です。行列で書くと、推定も検定も1本の式で見通せるようになります。そのうえで、変数選択、残差分析、一般化最小二乗推定、多重共線性、L1正則化法・L2正則化法 と順に進みます。

前回の記事では、確率過程の基礎を学びました。今回の記事では、回帰モデルの実践を、導出と判断の基準まで含めて解説します。

行列で書く重回帰モデル

予測値はyを平面に落とした影。残差は平面に垂直

はじめに、重回帰モデルの行列表現について解説します。

2級で学んだとおり、説明変数が $p$ 個の 重回帰モデル は、$y_i = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_p x_{pi} + \varepsilon_i$ でした。$n$ 個の式を、行列でまとめて書きます。

$$\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}, \qquad E(\boldsymbol{\varepsilon}) = \boldsymbol{0}, \qquad V(\boldsymbol{\varepsilon}) = \sigma^2 I$$

$\boldsymbol{y}$ は $n$ 個の観測値を並べたベクトルです。$X$ は $n \times (p+1)$ の 計画行列 で、1列目はすべて1、残りの列が説明変数です。$\boldsymbol{\beta}$ は切片を含む $p+1$ 個の係数です。

最小二乗法は、残差平方和 $(\boldsymbol{y} - X\boldsymbol{\beta})^{\top}(\boldsymbol{y} - X\boldsymbol{\beta})$ を最小にします。$\boldsymbol{\beta}$ で微分して0とおくと、次の 正規方程式 が得られます。

$$X^{\top}X\hat{\boldsymbol{\beta}} = X^{\top}\boldsymbol{y}$$

$X^{\top}X$ に逆行列があれば、次の式で推定値が求まります。

$$\hat{\boldsymbol{\beta}} = (X^{\top}X)^{-1}X^{\top}\boldsymbol{y}$$

単回帰の $b = S_{xy} / S_{xx}$ が、そのまま行列に広がった形です。「$x$ と $y$ の関係」を「$x$ 自身の広がり」で割っている、と読めます。

行列で書く利点は、性質が一度に分かることです。

【最小二乗推定量の性質】

• 不偏性： $E(\hat{\boldsymbol{\beta}}) = \boldsymbol{\beta}$
• 分散： $V(\hat{\boldsymbol{\beta}}) = \sigma^2 (X^{\top}X)^{-1}$
• 誤差分散の推定： $s^2 = \dfrac{\sum e_i^2}{n - p - 1}$

点推定の記事で学んだガウス・マルコフの定理により、$\hat{\boldsymbol{\beta}}$ は線形不偏推定量の中で分散が最小です。

予測値は、$\hat{\boldsymbol{y}} = X\hat{\boldsymbol{\beta}} = H\boldsymbol{y}$ と書けます。$H = X(X^{\top}X)^{-1}X^{\top}$ を ハット行列 と呼びます。$\boldsymbol{y}$ に帽子（ハット）をかぶせる行列です。

図形で考えると、イメージがはっきりします。$X$ の列が張る平面に、$\boldsymbol{y}$ を真上から光を当てて落とした影が $\hat{\boldsymbol{y}}$ です。影と本体の差である残差は、平面に垂直です。これが $X^{\top}\boldsymbol{e} = \boldsymbol{0}$ の意味で、正規方程式そのものです。$H$ のトレース（対角成分の和）は $p+1$ で、これがモデルの自由度になります。

例題：偏差の行列から係数を求める

2つの説明変数の偏差平方和・偏差積和が、$S_{11} = 10$、$S_{22} = 20$、$S_{12} = 10$、$S_{1y} = 30$、$S_{2y} = 40$ でした。$n = 10$、$S_{yy} = 120$ です。偏回帰係数と、$b_1$ の $t$ 値を求めましょう。

切片を除くと、正規方程式は偏差の行列で書けます。

$$\begin{pmatrix} 10 & 10 \\ 10 & 20 \end{pmatrix}\begin{pmatrix} b_1 \\ b_2 \end{pmatrix} = \begin{pmatrix} 30 \\ 40 \end{pmatrix}$$

【解き方】

1. 逆行列を求める： 行列式は $200 - 100 = 100$ なので、対角が 0.2, 0.1、非対角が −0.1 の行列
2. 係数を求める： $b_1 = 0.2 \times 30 - 0.1 \times 40 = 2$、$b_2 = -0.1 \times 30 + 0.1 \times 40 = 1$
3. 回帰変動を求める： $b_1 S_{1y} + b_2 S_{2y} = 60 + 40 = 100$
4. 残差変動を求める： $120 - 100 = 20$ で、$s^2 = 20 \div 7 = 2.86$
5. 標準誤差を求める： $\sqrt{2.86 \times 0.2} = 0.756$ で、$t = 2 \div 0.756 = 2.65$

答えは $b_1 = 2$、$b_2 = 1$、$t = 2.65$ です。$t_{0.025}(7) = 2.365$ を超えるので、$b_1$ は有意です。$b_2$ の $t$ 値は1.87で、有意ではありません。

試験のポイント：係数の標準誤差は、$s^2$ に $(X^{\top}X)^{-1}$ の対角成分を掛けて平方根をとります。

以上が、重回帰モデルの行列表現についての解説です。

決定係数と回帰の検定

全変動のうち回帰の割合がR²。F検定で全体を調べる

次に、モデル全体の当てはまりの測り方を説明します。

決定係数 は、全変動のうち回帰で説明できた割合 $R^2 = 1 - S_e / S_{yy}$ でした。重回帰では、$y$ と $\hat{y}$ の相関係数を 重相関係数 $R$ と呼び、その2乗が決定係数に一致します。

例題では、$R^2 = 100 \div 120 = 0.833$ です。2級で学んだ 自由度調整済み決定係数 は、$\bar{R}^2 = 1 - (20 \div 7) \div (120 \div 9) = 0.786$ です。モデル全体のF検定は、$F = (100 \div 2) \div (20 \div 7) = 17.5$ で、$F_{0.05}(2, 7) = 4.74$ を超えます。

準1級で大切なのは、一部の係数をまとめて検定する 方法です。大きいモデルから $q$ 個の変数を除いた小さいモデルを考えて、$H_0$: 除いた変数の係数はすべて0、を調べます。

$$F = \frac{(S_{e0} - S_{e1})/q}{S_{e1}/(n - p - 1)}$$

$S_{e0}$ は小さいモデル、$S_{e1}$ は大きいモデルの残差変動です。$H_0$ のもとで、$F$ は自由度 $(q, n-p-1)$ のF分布に従います。言葉で言うと、「変数を足して減った残差変動は、たまたまの範囲か」を測っています。$q = 1$ なら、$t$ 検定の2乗と同じです。

以上が、決定係数と回帰の検定についての説明です。

変数選択とMallows' Cp規準

Cpが小さくkに近いモデルを選ぶ。荷物は多すぎず少なすぎず

続いて、変数選択 について解説します。

変数を増やすと、バイアスは減りますが、推定する係数が増えて分散が増えます。点推定の性質の記事で学んだ、バイアス・バリアンス分解の綱引きです。荷物を減らしすぎると必要なものが足りず、詰めすぎると歩けなくなる、登山の荷造りのようなものです。

選び方には、次の手順があります。

【変数選択の手順】

• 総当たり法： すべての組合せを調べる。変数が $p$ 個なら $2^p$ 通り
• 前進選択法： 変数なしから始めて、最もよい変数を1つずつ足す
• 後退消去法： 全部入りから始めて、最も不要な変数を1つずつ除く
• ステップワイズ法： 足すのと除くのを交互に行う

比べる物差しには、次の Mallows' Cp規準 がよく使われます。切片を含めて $k$ 個の係数を持つモデルの残差変動を $S_e(k)$、全部入りのモデルの誤差分散の推定値を $\hat{\sigma}^2$ として、次の式で定めます。

$$C_p = \frac{S_e(k)}{\hat{\sigma}^2} - n + 2k$$

なぜこの形なのでしょうか？ $C_p$ は、予測値の平均二乗誤差の合計を $\sigma^2$ で割った量の推定値です。モデルにバイアスがなければ $E(S_e(k)) = (n-k)\sigma^2$ なので、$C_p$ は $k$ の近くになります。バイアスがあると、その分だけ $C_p$ は $k$ より大きくなります。つまり、第1項の「当てはまり」と第3項の「係数の数の罰」がつり合う点を探す規準です。

$C_p$ が小さく、かつ $k$ に近いモデルを選びます。全部入りのモデルでは、定義から必ず $C_p = k$ になります。

例題：Cpで変数を選ぶ

$n = 30$、説明変数は5個、全変動は400です。全部入りのモデルの残差変動は96で、$\hat{\sigma}^2 = 96 \div 24 = 4$ です。候補のモデルは次のとおりです。

変数	$k$	残差変動	$R^2$	調整済み	$C_p$
$x_1$	2	180	0.55	0.534	19
$x_1, x_2$	3	112	0.72	0.699	4
$x_1 \sim x_3$	4	108	0.73	0.699	5
$x_1 \sim x_4$	5	100	0.75	0.710	5
$x_1 \sim x_5$	6	96	0.76	0.710	6

【解き方】

1. $C_p$ を計算する： 例えば、$x_1, x_2$ のモデルは $112 \div 4 - 30 + 6 = 4$
2. 最小の値を探す： $x_1, x_2$ のモデルの4
3. 部分F検定で確かめる： 全部入りと比べて $F = (16 \div 3) \div 4 = 1.33$ で、$F_{0.05}(3, 24) = 3.01$ より小さい

答えは「$x_1, x_2$ のモデルを選ぶ」です。

$R^2$ は変数を足すほど必ず増えるので、規準になりません。調整済み決定係数も、足した変数の $F$ が1を超えるだけで増えるので、変数の追加に甘い物差しです。

ただし、選んだ後のモデルの $t$ 値やp値には注意が必要です。データを見て変数を選んだので、見かけより有意に見えるのです。AICなどの情報量規準は、モデル選択の記事で学びます。

試験のポイント：$C_p$ の $k$ は切片を含む係数の数です。全部入りのモデルでは $C_p = k$ です。

以上が、変数選択についての解説です。

残差分析と一般化最小二乗推定

残差の形から不調を読む。分散が違えば重みをつける

さらに、モデルの前提の確かめ方と、前提が崩れたときの対処を見ていきましょう。

残差分析 は、残差のふるまいからモデルの前提を点検する作業です。健康診断の検査結果のように、残差のパターンから体の不調を読み取ります。

残差は $\boldsymbol{e} = (I - H)\boldsymbol{y}$ で、$V(\boldsymbol{e}) = \sigma^2(I - H)$ です。誤差 $\varepsilon_i$ と違って、$e_i$ の分散は等しくありません。そこで、$H$ の対角成分 $h_{ii}$ を使って、$r_i = e_i / (s\sqrt{1 - h_{ii}})$ と揃えます。これを 標準化残差 と呼びます。

【残差のパターンと疑うこと】

• U字などの曲線： 2乗の項の不足など、モデルの形が違う
• 右へ広がる漏斗形： 誤差の分散が等しくない（不均一分散）
• 時間の順に波打つ： 誤差に相関がある（系列相関）

図を使った診断の詳しい方法は、次回の記事で学びます。

では、誤差の分散が等しくないときは、どうすればよいのでしょうか？ 最小二乗推定量は不偏のままですが、最良ではなくなり、標準誤差の式も崩れます。

そこで、$V(\boldsymbol{\varepsilon}) = \sigma^2 \Omega$ という一般の形を考えます。$\Omega$ は既知の正定値行列です。$\Omega = CC^{\top}$ と分解して、モデルの両辺に左から $C^{-1}$ を掛けます。すると、誤差の分散が $\sigma^2 I$ に戻るので、変換後のデータにガウス・マルコフの定理が使えます。結果は次の式です。

$$\hat{\boldsymbol{\beta}}_{G} = (X^{\top}\Omega^{-1}X)^{-1}X^{\top}\Omega^{-1}\boldsymbol{y}, \qquad V(\hat{\boldsymbol{\beta}}_{G}) = \sigma^2(X^{\top}\Omega^{-1}X)^{-1}$$

これが 一般化最小二乗推定（GLS）です。誤差の大きな観測の発言力を下げる、というイメージです。声の大きさではなく、話の確かさで重みをつける会議のようなものです。

$\Omega$ が対角行列なら、各観測に重み $1/\omega_i$ を付ける 重み付き最小二乗法 になります。

具体的に見てみましょう。原点を通るモデル $y_i = \beta x_i + \varepsilon_i$ で、$V(\varepsilon_i) = \sigma^2 x_i$ とします。店の規模が大きいほど売上のばらつきも大きい、という場面です。重みは $1/x_i$ なので、推定値は次のようになります。

$$\hat{\beta}_{G} = \frac{\sum (1/x_i) x_i y_i}{\sum (1/x_i) x_i^2} = \frac{\sum y_i}{\sum x_i}$$

データが $(x, y) = (1, 3), (2, 5), (4, 12)$ なら、$\hat{\beta}_{G} = 20 \div 7 = 2.86$ です。普通の最小二乗法では $\sum x_i y_i / \sum x_i^2 = 61 \div 21 = 2.90$ で、大きい $x$ の影響が強く出ます。

試験のポイント：不均一分散や系列相関があっても、最小二乗推定量は不偏です。失われるのは効率で、標準誤差の式も使えなくなります。

以上が、残差分析と一般化最小二乗推定についての説明です。

多重共線性と正則化

ひし形の角で触れると、係数がちょうど0になる

最後に、多重共線性と、その対処である正則化を見ていきましょう。

多重共線性 は、説明変数どうしの強い相関で係数が不安定になる現象でした。理由は、分散の式を見ると分かります。$b_j$ の分散は次のように書けます。

$$V(b_j) = \frac{\sigma^2}{S_{jj}(1 - R_j^2)}$$

$R_j^2$ は、$x_j$ を他の説明変数で回帰したときの決定係数です。$1 / (1 - R_j^2)$ が2級の補足で触れたVIFです。例題では $x_1$ と $x_2$ の相関が0.707で、$R_j^2 = 0.5$ なので、VIFは2です。相関が0.95なら、VIFは約10.3まで増えます。

$X^{\top}X$ の逆行列を計算するときに、行列式がほぼ0で割っている状態です。データが少し動くだけで、係数が大きく振れます。

そこで、係数の大きさに罰を与えて、振れを抑えます。これが 正則化 です。あえて少しバイアスを受け入れて、分散を大きく減らします。1つずつ見ていきましょう。

1. L2正則化法（リッジ回帰）

L2正則化法 は、係数の2乗の和に罰を与えます。

$$\min_{\boldsymbol{\beta}} \left\{\sum_{i}(y_i - \hat{y}_i)^2 + \lambda\sum_{j=1}^{p}\beta_j^2\right\}, \qquad \hat{\boldsymbol{\beta}}_{R} = (X^{\top}X + \lambda I)^{-1}X^{\top}\boldsymbol{y}$$

$\lambda \geq 0$ は罰の強さです。対角に $\lambda$ を足すので、逆行列が安定します。ぐらつくテーブルの脚に、下駄を履かせるイメージです。

先ほどの例題の偏差の行列に $\lambda = 10$ を足して解くと、$(b_1, b_2) = (2, 1)$ が $(1, 1)$ になります。相関のある変数の間で、効果がならされるのです。

2. L1正則化法（ラッソ）

L1正則化法 は、係数の絶対値の和に罰を与えます。ラッソ（lasso）とも呼びます。

$$\min_{\boldsymbol{\beta}} \left\{\sum_{i}(y_i - \hat{y}_i)^2 + \lambda\sum_{j=1}^{p}|\beta_j|\right\}$$

ラッソの特徴は、係数がちょうど0になる ことです。つまり、推定と変数選択を同時に行えます。

違いがよく見えるのは、説明変数が互いに直交し、$X^{\top}X = I$ のときです。最小二乗推定値を $b_j$ として、次のようになります。

【直交のときの推定値】

• リッジ： $\dfrac{b_j}{1 + \lambda}$（一律に縮める）
• ラッソ： $\mathrm{sign}(b_j)\max(|b_j| - \lambda/2,\ 0)$（一定量を引き、足りなければ0にする）

例えば、$b = (3, 1, -0.4)$、$\lambda = 1$ なら、リッジは $(1.5, 0.5, -0.2)$ です。ラッソは $(2.5, 0.5, 0)$ で、小さい係数が消えました。リッジは全員の給料を半分にし、ラッソは全員から一定額を引いて、足りない人を解雇するようなものです。

形で考えると、理由が分かります。罰の制約の領域は、リッジでは円、ラッソでは軸の上に角を持つひし形です。残差平方和の等高線は、ひし形の角で触れやすく、角では一部の係数が0になります。

【正則化の注意点】

• 説明変数を標準化してから使う。単位で罰の効き方が変わるため
• 切片には罰を与えない
• $\lambda$ はクロスバリデーションで決める
• 正則化した推定量は不偏ではない

補足（ここは読み飛ばしてOK）： ベイズ法の言葉では、リッジは係数に正規分布の事前分布を、ラッソはラプラス分布の事前分布を置いたときの事後モードに当たります。2つの罰を混ぜた方法をエラスティックネットと呼びます。

試験のポイント：係数が0になるのはL1、一律に縮むのはL2です。$\lambda = 0$ なら最小二乗法に戻り、$\lambda$ を大きくするほど係数は0に近づきます。

以上が、多重共線性と正則化についての説明です。

細かい導出が分からなかったとしても問題ありません。ここでは、「$R^2$ だけでモデルを選ばず、バイアスと分散のつり合いで判断する」という点だけ押さえておきましょう。

今回のまとめ

平面への影、Cpでの選択、罰で係数を抑える

それでは、今回の内容を振り返ります。

今回は重回帰分析を、行列の形から変数選択と正則化まで紹介しました。推定の式はきれいでも、変数を増やせばよいわけではありません。当てはまりと安定性の綱引きをどう決着させるかが、実践の中心です。

今回のポイントを振り返りましょう。

振り返り

• 最小二乗推定は $\hat{\boldsymbol{\beta}} = (X^{\top}X)^{-1}X^{\top}\boldsymbol{y}$。予測値は $\boldsymbol{y}$ の平面への影。
• 係数の分散は $\sigma^2(X^{\top}X)^{-1}$。一部の係数の検定には部分F検定を使う。
• Mallows' Cp規準は、当てはまりと係数の数のつり合いを見る。小さくて $k$ に近いモデルを選ぶ。
• 誤差の分散が等しくないときは、重みをつける一般化最小二乗推定を使う。
• 多重共線性には正則化が効く。L2は一律に縮め、L1は係数を0にして変数も選ぶ。

キーワードのおさらい

• 重回帰モデル: 複数の説明変数で目的変数を説明する線形モデル
• 決定係数: 全変動のうち回帰で説明できた割合
• 自由度調整済み決定係数: 変数の数の分だけ罰を与えた決定係数
• 変数選択: モデルに残す説明変数を決める作業
• 残差分析: 残差のパターンからモデルの前提を点検する作業
• 一般化最小二乗推定: 誤差の分散・共分散の構造を重みに使う推定法
• 多重共線性: 説明変数どうしの強い相関で、係数が不安定になる現象
• L1正則化法: 係数の絶対値の和に罰を与える方法。ラッソ
• L2正則化法: 係数の2乗の和に罰を与える方法。リッジ回帰
• Mallows' Cp規準: 残差変動と係数の数からモデルを評価する規準

さて、今回は残差のパターンを大まかに見ました。では、1つの点がモデル全体を振り回しているとき、どうやって見つければよいのでしょうか？

次回は、外れ値や影響の大きい点を見つける「回帰診断法」について紹介します。

それでは、次の記事へ進んでいきましょう。