第5章線形回帰（統計学） | Lemma

この章の目標

線形回帰モデルを行列で書き、最小二乗推定量を正規方程式から求められる

ハット行列が列空間への直交射影であることを証明し、残差・決定係数・てこ比を射影の言葉で説明できる

ガウス–マルコフの定理を、どの仮定を使うか（正規性は使わない）に注意して証明できる

正規線形モデルで係数の分布を導き、 $t$ 検定・ $F$ 検定・信頼区間・予測区間を正しく使える

多重共線性・外れ値・てこ比・ダミー変数を扱い、残差でモデルを診断できる

最小二乗問題を QR 分解で解く理由を、条件数を使って説明できる

前提：第1章、第2章、第4章、02 線形代数第7章。5.10 節では 02 線形代数第8章（特異値分解と作用素ノルム）を使う。

中古マンションの価格を面積と築年数で説明する、来週の需要を価格と広告費から予測する――ある量 $y$ をほかの量 $x_1, \dots, x_k$ の一次式で近似するのが線形回帰 (linear regression) である。統計学で最もよく使われる手法であり、次章の一般化線形モデルや機械学習の多くの手法の出発点でもある。

最小二乗法そのものは 02 第7章で直交射影の応用として学んだ。本章ではこれを確率モデルとして扱い、「推定値はどのくらいぶれるか」「この変数は本当に効いているか」「新しい物件の価格をどの程度の幅で予測できるか」に答える。鍵は幾何である。観測値のベクトルを説明変数の張る部分空間に直交射影したものが当てはめ値で、残差はその部分空間に直交する。標準誤差・ $t$ 検定・ $F$ 検定・決定係数は、すべてこの直交分解とピタゴラスの定理から出てくる。後半では、仮定が崩れたとき（強い相関、外れ値、一定でない誤差分散）に何が起こるかを確かめ、正規方程式を直接解いてはいけない数値計算上の理由を述べる。

記法. 第1章と同じく、転置を $A^{\top}$ （02 線形代数の ${}^tA$ と同じもの）、確率ベクトル $Z$ の共分散行列を $\operatorname{Cov}(Z) = E[(Z - E[Z])(Z - E[Z])^{\top}]$ と書く。第1章命題 1.10 の公式 $E[AZ + b] = AE[Z] + b$ 、 $\operatorname{Cov}(AZ + b) = A\operatorname{Cov}(Z)A^{\top}$ （ $A, b$ は定数）を「共分散行列の変換公式」と呼ぶ。 $\mathbf{1} = (1, \dots, 1)^{\top}$ 、 $u_i$ は第 $i$ 成分だけが $1$ の基本ベクトルである。対称行列 $A, B$ について、 $A - B$ が半正定値（02 第8章定義 8.12）であることを $A \succeq B$ と書く。

5.1 線形回帰モデル

定義 5.1（線形回帰モデル, linear regression model） $X$ を確率的でない $n \times p$ 行列、 $\beta \in \mathbb{R}^p$ を未知の定数ベクトル、 $\varepsilon$ を $n$ 次元確率ベクトルとして
$Y = X\beta + \varepsilon$
と表されるモデルを線形回帰モデルという。 $X$ を計画行列 (design matrix)、 $\beta$ を回帰係数 (regression coefficient)、 $\varepsilon$ を誤差 (error) という。誤差について次の仮定を考える。

(A1) $E[\varepsilon] = 0$ 。

(A2) $\operatorname{Cov}(\varepsilon) = \sigma^2 I_n$ （ $\sigma^2 > 0$ 。誤差は等分散で、互いに無相関）。

(A3) $\varepsilon \sim N_n(0, \sigma^2 I_n)$ 。

(A3) は (A1)(A2) を含む。(A3) を仮定したモデルを正規線形モデルという。

ふつう $X$ の第 1 列は $\mathbf{1}$ （定数項, intercept）で、残りの列が説明変数 $x_{i1}, \dots, x_{ik}$ である。このとき $p = k + 1$ 、 $\beta = (\beta_0, \dots, \beta_k)^{\top}$ で、 $Y_i = \beta_0 + \beta_1x_{i1} + \cdots + \beta_kx_{ik} + \varepsilon_i$ となる。「線形」とは $\beta$ について線形という意味で、 $x^2$ や $\log x$ の列、5.8 節のダミー変数の列を使ってよい。説明変数が確率変数であるときは、 $X$ を与えたときの条件付きの議論と読む。本章では断らない限り $\operatorname{rank} X = p$ （ $X$ の列が一次独立）を仮定する。

例 5.2（本章のデータ）ある地域の中古マンション 12 件の、専有面積 $x_1$ （m²）、築年数 $x_2$ （年）、価格 $y$ （百万円）である（説明用の架空のデータ）。

物件	1	2	3	4	5	6	7	8	9	10	11	12
面積 $x_1$	45	52	58	60	63	66	70	72	75	80	85	90
築年数 $x_2$	25	10	30	15	5	20	12	28	8	18	3	22
価格 $y$	16.2	27.7	21.1	25.6	36.2	30.2	34.0	30.2	40.4	38.5	47.5	42.7

モデル $Y_i = \beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \varepsilon_i$ の計画行列は、第 $i$ 行が $(1, x_{i1}, x_{i2})$ の $12 \times 3$ 行列である。

5.2 最小二乗推定量と正規方程式

定義 5.3（最小二乗推定量） $X$ の第 $i$ 行を $x_i^{\top}$ とする。 $\lVert y - X\beta \rVert^2 = \sum_i (y_i - x_i^{\top}\beta)^2$ を最小にする $\beta$ を最小二乗推定量 (ordinary least squares estimator, OLS) といい、 $\hat{\beta}$ と書く。 $\hat{y} = X\hat{\beta}$ を当てはめ値 (fitted value)、 $e = y - \hat{y}$ を残差 (residual)、 $\mathrm{RSS} = \lVert e \rVert^2$ を残差平方和という。

定理 5.4（正規方程式） $\beta$ が $\lVert y - X\beta \rVert^2$ を最小にするための必要十分条件は、正規方程式 $X^{\top}X\beta = X^{\top}y$ をみたすことである。 $\operatorname{rank} X = p$ ならば $X^{\top}X$ は正則で、最小点は $\hat{\beta} = (X^{\top}X)^{-1}X^{\top}y$ ただ一つである。

これは 02 第7章の定理 7.19 そのものだが、平方完成による短い証明を与えておく。

証明. 正規方程式の解 $\hat{\beta}$ を 1 つとる（存在は 02 の定理 7.19）。 $X^{\top}(y - X\hat{\beta}) = 0$ なので、任意の $\beta$ について交差項 $2(\hat{\beta} - \beta)^{\top}X^{\top}(y - X\hat{\beta})$ が消えて

\lVert y - X\beta \rVert^2 = \lVert (y - X\hat{\beta}) + X(\hat{\beta} - \beta) \rVert^2 = \lVert y - X\hat{\beta} \rVert^2 + \lVert X(\hat{\beta} - \beta) \rVert^2

となる。よって正規方程式の解は最小点である。逆に $\beta$ が最小点なら $X\beta = X\hat{\beta}$ なので $X^{\top}X\beta = X^{\top}X\hat{\beta} = X^{\top}y$ 。 $\operatorname{rank} X = p$ のとき、 $X^{\top}Xv = 0$ なら $\lVert Xv \rVert^2 = v^{\top}X^{\top}Xv = 0$ より $Xv = 0$ 、よって $v = 0$ なので $X^{\top}X$ は正則である。 $\square$

正規方程式 $X^{\top}e = 0$ は「残差がどの説明変数の列とも直交する」ことを表す。特に定数項があれば $\sum_i e_i = 0$ である。

例 5.5（単回帰） $X = (\mathbf{1}\ x)$ のとき、正規方程式は

n\hat{\beta}_0 + \Bigl(\sum_i x_i\Bigr)\hat{\beta}_1 = \sum_i y_i, \qquad \Bigl(\sum_i x_i\Bigr)\hat{\beta}_0 + \Bigl(\sum_i x_i^2\Bigr)\hat{\beta}_1 = \sum_i x_iy_i

である。第 1 式から $\hat{\beta}_0 = \bar{y} - \hat{\beta}_1\bar{x}$ （回帰直線は点 $(\bar{x}, \bar{y})$ を通る）。これを第 2 式に代入して整理すると、 $S_{xx} = \sum_i (x_i - \bar{x})^2$ 、 $S_{xy} = \sum_i (x_i - \bar{x})(y_i - \bar{y})$ として $\hat{\beta}_1 = S_{xy}/S_{xx}$ を得る。例 5.2 で面積だけを使うと、 $\bar{x} = 68$ 、 $\bar{y} = 32.525$ 、 $S_{xx} = 1964$ 、 $S_{xy} = 1207.5$ より $\hat{\beta}_1 = 0.6148$ 、 $\hat{\beta}_0 = -9.283$ 。

例 5.6（重回帰）例 5.2 の 2 変数モデルでは、正規方程式を解いて

\hat{y} = 2.054 + 0.5531x_1 - 0.4373x_2

を得る。「築年数が同じなら、面積が 1 m² 広いと価格は平均して約 0.553 百万円高い」「面積が同じなら、築年数が 1 年古いと約 0.437 百万円安い」と読む。面積の係数が単回帰の $0.6148$ より小さいのは、このデータでは広い物件ほど新しい傾向がある（ $x_1$ と $x_2$ の相関係数は $-0.21$ ）ため、単回帰の傾きに築年数の効果の一部が混ざっていたからである（問題 5.4）。「ほかの変数を一定にしたときの効果」の正確な意味は 5.7 節で与える。

5.3 ハット行列と射影の幾何

定義 5.7（ハット行列, hat matrix） $H = X(X^{\top}X)^{-1}X^{\top}$ をハット行列という。 $\hat{y} = Hy$ 、 $e = (I_n - H)y$ である。

$X$ の列空間を $\mathcal{C}(X) = \lbrace X\beta \mid \beta \in \mathbb{R}^p \rbrace$ と書く。

定理 5.8（ハット行列は直交射影） $\operatorname{rank} X = p$ とする。

$H^{\top} = H$ かつ $H^2 = H$ （対称かつ冪等）。

$H$ は $\mathcal{C}(X)$ への直交射影、 $I_n - H$ は $\mathcal{C}(X)^{\perp}$ への直交射影である。

$\operatorname{rank} H = \operatorname{tr} H = p$ 、 $\operatorname{tr}(I_n - H) = n - p$ 。

証明. (1) 対称行列 $X^{\top}X$ の逆行列は対称なので（ $((X^{\top}X)^{-1})^{\top} = ((X^{\top}X)^{\top})^{-1}$ ）、 $H^{\top} = H$ 。また $H^2 = X(X^{\top}X)^{-1}(X^{\top}X)(X^{\top}X)^{-1}X^{\top} = H$ 。

(2) 任意の $y$ について $Hy = X\bigl((X^{\top}X)^{-1}X^{\top}y\bigr) \in \mathcal{C}(X)$ であり、 $X^{\top}(y - Hy) = X^{\top}y - X^{\top}y = 0$ より $y - Hy \in \mathcal{C}(X)^{\perp}$ 。直交分解 $\mathbb{R}^n = \mathcal{C}(X) \oplus \mathcal{C}(X)^{\perp}$ （02 の定理 7.15）による分解はただ一通りなので、 $Hy$ は $y$ の $\mathcal{C}(X)$ 成分、すなわち直交射影による像であり、 $(I_n - H)y$ は $\mathcal{C}(X)^{\perp}$ 成分である（(1) と 02 の命題 7.24 からも従う）。

(3) $\operatorname{tr}(AB) = \operatorname{tr}(BA)$ （02 第1章問題 1.6）より $\operatorname{tr} H = \operatorname{tr}\bigl((X^{\top}X)^{-1}X^{\top}X\bigr) = \operatorname{tr} I_p = p$ 。 $H$ の像は $\mathcal{C}(X)$ で、その次元は $p$ 。 $\square$

定理 5.8 から次のことがわかる。

$y = \hat{y} + e$ は直交分解で、ピタゴラスの定理より $\lVert y \rVert^2 = \lVert \hat{y} \rVert^2 + \mathrm{RSS}$ 。
$\mathcal{C}(X)$ の正規直交基底 $q_1, \dots, q_p$ をとると $H = \sum_{j=1}^p q_jq_j^{\top}$ （02 の定理 7.15 の式 $P_W(v) = \sum_j \langle v, q_j \rangle q_j$ を行列で書いたもの）。
$H$ は $\mathcal{C}(X)$ だけで決まる。 $X$ を正則行列 $A$ で $XA$ に取り替えても（単位の変更、変数の中心化、ダミー変数の別のとり方など）、 $\hat{y}$ ・ $e$ ・RSS は変わらない。変わるのは係数の読み方だけである。

5.4 ガウス–マルコフの定理

命題 5.9 (A1)(A2) のもとで $E[\hat{\beta}] = \beta$ 、 $\operatorname{Cov}(\hat{\beta}) = \sigma^2(X^{\top}X)^{-1}$ 。

証明. $\hat{\beta} = (X^{\top}X)^{-1}X^{\top}(X\beta + \varepsilon) = \beta + (X^{\top}X)^{-1}X^{\top}\varepsilon$ に共分散行列の変換公式を使うと、 $E[\hat{\beta}] = \beta$ 、 $\operatorname{Cov}(\hat{\beta}) = (X^{\top}X)^{-1}X^{\top}(\sigma^2I_n)X(X^{\top}X)^{-1} = \sigma^2(X^{\top}X)^{-1}$ 。 $\square$

定理 5.10（ガウス–マルコフの定理, Gauss–Markov theorem） (A1)(A2) を仮定する。 $p \times n$ 定数行列 $C$ による推定量 $\tilde{\beta} = CY$ が、すべての $\beta \in \mathbb{R}^p$ について $E[\tilde{\beta}] = \beta$ をみたす（線形不偏推定量である）ならば、 $\operatorname{Cov}(\tilde{\beta}) \succeq \operatorname{Cov}(\hat{\beta})$ である。特に任意の $c \in \mathbb{R}^p$ について $\operatorname{Var}(c^{\top}\tilde{\beta}) \geq \operatorname{Var}(c^{\top}\hat{\beta})$ であり、すべての $c$ で等号が成り立つのは $\tilde{\beta} = \hat{\beta}$ のときに限る。すなわち $\hat{\beta}$ は最良線形不偏推定量 (best linear unbiased estimator, BLUE) である。

証明. $E[CY] = CX\beta$ がすべての $\beta$ で $\beta$ に等しいので $CX = I_p$ 。 $D = C - (X^{\top}X)^{-1}X^{\top}$ とおくと $DX = CX - I_p = O$ 。共分散行列の変換公式より

\begin{aligned} \operatorname{Cov}(CY) = \sigma^2CC^{\top} &= \sigma^2\bigl[(X^{\top}X)^{-1} + (X^{\top}X)^{-1}(DX)^{\top} + DX(X^{\top}X)^{-1} + DD^{\top}\bigr] \\ &= \sigma^2(X^{\top}X)^{-1} + \sigma^2DD^{\top} \end{aligned}

で、 $v^{\top}DD^{\top}v = \lVert D^{\top}v \rVert^2 \geq 0$ より $DD^{\top}$ は半正定値。 $\operatorname{Var}(c^{\top}\tilde{\beta}) = c^{\top}\operatorname{Cov}(\tilde{\beta})c$ なので後半の不等式が従う。すべての $c$ で等号なら、すべての $c$ で $\lVert D^{\top}c \rVert^2 = 0$ なので $D = O$ 、すなわち $C = (X^{\top}X)^{-1}X^{\top}$ 。 $\square$

注意 5.11（仮定の役割）

使ったのは (A1)(A2) だけで、誤差の正規性は使っていない。(A3) のもとではさらに強く、 $\hat{\beta}$ は線形に限らずすべての不偏推定量の中で分散が最小になる（第3章定理 3.25 のクラメール–ラオの不等式を多次元に拡張したものから従う。証明は省略する）。

最良なのは「線形かつ不偏」な推定量の中でである。偏りを許せば平均二乗誤差のより小さい推定量がありうる（リッジ推定量、第6章）。誤差の分布の裾が重いときは、最小絶対偏差法のような非線形の推定量の方が（漸近的に）分散が小さくなることがある。

(A2) が崩れて $\operatorname{Cov}(\varepsilon) = \Sigma$ となっても、(A1) だけで $\hat{\beta}$ は不偏である。しかし最良ではなくなり（ $\Sigma$ が既知の正定値行列なら、 $\Sigma^{-1/2}Y$ に定理 5.10 を適用して、一般化最小二乗推定量 $(X^{\top}\Sigma^{-1}X)^{-1}X^{\top}\Sigma^{-1}Y$ が BLUE になる）、何より共分散行列が $\operatorname{Cov}(\hat{\beta}) = (X^{\top}X)^{-1}X^{\top}\Sigma X(X^{\top}X)^{-1}$ に変わるので、 $\sigma^2(X^{\top}X)^{-1}$ にもとづく標準誤差は正しくなくなる。

ヒント

実務では 回帰係数そのものより、標準誤差の方が仮定に敏感である。売上の大きい店ほどばらつきも大きい、同じ顧客の観測が繰り返し入っている、時系列で誤差に自己相関がある――といった場合は (A2) が崩れ、通常の標準誤差は正しくなくなる（多くの場合は小さく出て、 $p$ 値が楽観的になる）。注意 5.11 の 3 の共分散行列の式で、中央の $X^{\top}\Sigma X$ を残差から推定したもの（不均一分散なら $\Sigma$ を $\operatorname{diag}(e_1^2, \dots, e_n^2)$ で置き換える。同じ顧客の観測が繰り返し入っているなら顧客ごとのまとまりで推定する）にもとづく頑健な標準誤差（サンドイッチ型）を併用し、通常の標準誤差と大きく違わないかを確かめるのが標準的である（ $\operatorname{diag}(e_1^2, \dots, e_n^2)$ 自体は $\Sigma$ のよい推定ではないが、 $n$ が大きければ、それを挟んだ $X^{\top}\operatorname{diag}(e_1^2, \dots, e_n^2)X$ は $X^{\top}\Sigma X$ のよい近似になる。正確な条件は省略する）。

5.5 誤差分散の推定と正規線形モデルでの推測

定理 5.12（ $\hat{\sigma}^2$ の不偏性） (A1)(A2) のもとで $E[\mathrm{RSS}] = (n - p)\sigma^2$ 。したがって $\hat{\sigma}^2 = \mathrm{RSS}/(n - p)$ は $\sigma^2$ の不偏推定量である。

証明. $(I_n - H)X = O$ より $e = (I_n - H)(X\beta + \varepsilon) = (I_n - H)\varepsilon$ 。 $I_n - H$ は対称かつ冪等なので $\mathrm{RSS} = \varepsilon^{\top}(I_n - H)\varepsilon$ 。スカラーはそのトレースに等しいので、 $\operatorname{tr}(AB) = \operatorname{tr}(BA)$ と期待値の線形性から

E[\mathrm{RSS}] = E\bigl[\operatorname{tr}\bigl((I_n - H)\varepsilon\varepsilon^{\top}\bigr)\bigr] = \operatorname{tr}\bigl((I_n - H)E[\varepsilon\varepsilon^{\top}]\bigr) = \sigma^2\operatorname{tr}(I_n - H) = (n - p)\sigma^2

（定理 5.8 の 3）。 $\square$

$n - p$ を残差の自由度という。 $\mathrm{RSS}/n$ （正規線形モデルでの $\sigma^2$ の最尤推定量）は $\sigma^2$ を過小に推定する。 $X = \mathbf{1}$ なら $\mathrm{RSS} = \sum_i (y_i - \bar{y})^2$ で、定理 5.12 は不偏分散の不偏性（第2章定理 2.3）そのものである。また同じ計算から $\operatorname{Cov}(e) = \sigma^2(I_n - H)$ で、誤差が無相関・等分散でも、残差は相関をもち、分散も等しくない（5.9 節）。

定理 5.13（正規線形モデルでの分布） (A3) を仮定し、 $p < n$ とする。

$\hat{\beta} \sim N_p\bigl(\beta, \sigma^2(X^{\top}X)^{-1}\bigr)$ 。

$\mathrm{RSS}/\sigma^2 \sim \chi^2(n - p)$ 。

$\hat{\beta}$ と $\mathrm{RSS}$ は独立である。

証明. (1) $\hat{\beta} = \beta + (X^{\top}X)^{-1}X^{\top}\varepsilon$ は $\varepsilon$ の一次式なので、多変量正規分布の線形変換の性質（第1章定理 1.22 の 2）より多変量正規分布に従い、平均と共分散行列は命題 5.9 で求めた。

(2)(3) $Y \sim N_n(X\beta, \sigma^2I_n)$ に、直交分解 $\mathbb{R}^n = \mathcal{C}(X) \oplus \mathcal{C}(X)^{\perp}$ について第2章定理 2.7（正規ベクトルの直交分解）を使う。2 つの部分空間への直交射影は $H$ と $I_n - H$ （定理 5.8）で、 $(I_n - H)X\beta = 0$ なので、 $\mathrm{RSS} = \lVert (I_n - H)Y \rVert^2$ について $\mathrm{RSS}/\sigma^2 \sim \chi^2(n - p)$ であり、 $HY$ と $(I_n - H)Y$ は独立である。 $H$ は対称で $HX = X$ なので $X^{\top}H = X^{\top}$ となり、 $\hat{\beta} = (X^{\top}X)^{-1}X^{\top}HY$ は $HY$ の関数、RSS は $(I_n - H)Y$ の関数である。よって $\hat{\beta}$ と RSS は独立である。 $\square$

$X = \mathbf{1}$ とすれば、これは第2章定理 2.8（正規標本の基本定理。標本平均と不偏分散の独立性）そのものである。

系 5.14（係数の $t$ 検定と信頼区間） (A3) を仮定し、 $p < n$ とする。 $(X^{\top}X)^{-1}$ の $\beta_j$ に対応する対角成分を $v_{jj}$ とし、 $\operatorname{SE}(\hat{\beta}_j) = \hat{\sigma}\sqrt{v_{jj}}$ （ $\hat{\beta}_j$ の標準誤差, standard error）とおくと
$T_j = \frac{\hat{\beta}_j - \beta_j}{\operatorname{SE}(\hat{\beta}_j)} \sim t(n - p)$
である。したがって $t(m)$ の上側 $\alpha/2$ 点を $t_{\alpha/2}(m)$ として、 $\hat{\beta}_j \pm t_{\alpha/2}(n - p)\operatorname{SE}(\hat{\beta}_j)$ は $\beta_j$ の信頼係数 $1 - \alpha$ の信頼区間であり、帰無仮説 $\beta_j = 0$ は $\lvert \hat{\beta}_j \rvert/\operatorname{SE}(\hat{\beta}_j) > t_{\alpha/2}(n - p)$ のとき有意水準 $\alpha$ で棄却される。

証明. 定理 5.13 の 1 より $U = (\hat{\beta}_j - \beta_j)/(\sigma\sqrt{v_{jj}}) \sim N(0, 1)$ 、2 より $V = \mathrm{RSS}/\sigma^2 \sim \chi^2(n - p)$ で、3 より $U$ と $V$ は独立。 $T_j = U/\sqrt{V/(n - p)}$ なので、 $t$ 分布の定義（第2章定義 2.10）より $T_j \sim t(n - p)$ 。区間と検定は第4章の構成による。 $\square$

同じ議論は予測にも使える。新しい対象の説明変数を $x_0$ 、 $Y_0 = x_0^{\top}\beta + \varepsilon_0$ （ $\varepsilon_0 \sim N(0, \sigma^2)$ は $\varepsilon$ と独立）とすると、 $Y_0 - x_0^{\top}\hat{\beta} = \varepsilon_0 - x_0^{\top}(\hat{\beta} - \beta)$ は平均 $0$ 、分散 $\sigma^2(1 + x_0^{\top}(X^{\top}X)^{-1}x_0)$ の正規分布に従い、RSS と独立である（ $\varepsilon_0$ は $Y$ と独立で、 $\hat{\beta}$ は定理 5.13 の証明の $HY$ の関数だから）。よって

x_0^{\top}\hat{\beta} \pm t_{\alpha/2}(n - p)\ \hat{\sigma}\sqrt{1 + x_0^{\top}(X^{\top}X)^{-1}x_0}

は $Y_0$ を確率 $1 - \alpha$ で含む予測区間 (prediction interval) である。根号の中の $1$ を除くと平均 $x_0^{\top}\beta$ の信頼区間になり、その幅は $n \to \infty$ で $0$ に近づくが、予測区間の幅は $0$ にならない。在庫を決めるときに必要なのは予測区間の方である。

例 5.15（例 5.2 の推測）残差の自由度は $12 - 3 = 9$ 、 $\mathrm{RSS} = 16.56$ 、 $\hat{\sigma} = 1.356$ で、

係数	推定値	標準誤差	$t$ 値	$p$ 値
定数項 $\beta_0$	2.054	2.439	0.84	0.42
面積 $\beta_1$	0.5531	0.0313	17.67	$2.7 \times 10^{-8}$
築年数 $\beta_2$	$-0.4373$	0.0467	$-9.36$	$6.2 \times 10^{-6}$

である。 $t_{0.025}(9) = 2.262$ より、 $\beta_1$ の 95% 信頼区間は $0.5531 \pm 2.262 \times 0.0313$ 、すなわち $[0.482, 0.624]$ 。面積 70 m²・築 10 年の物件では $\hat{y}_0 = 36.40$ で、平均価格の 95% 信頼区間は $[35.30, 37.50]$ 、価格そのものの 95% 予測区間は $[33.14, 39.66]$ とずっと広い。定数項の検定（ $p = 0.42$ ）は「面積 0 m²・築 0 年の物件の平均価格は 0 か」という意味のない問いで、これを理由に定数項を外してはいけない（定理 5.19 が使えなくなる）。

$F$ 検定

複数の係数をまとめて検定したいことがある（築年数と駅からの距離はどちらも不要か、カテゴリ変数のダミーはまとめて不要か、など）。

定理 5.16（ $F$ 検定, F-test） (A3) を仮定し、 $p < n$ とする。 $X_0$ を $\operatorname{rank} X_0 = r < p$ 、 $\mathcal{C}(X_0) \subset \mathcal{C}(X)$ をみたす $n \times r$ 行列とし、 $X_0$ による最小二乗法の残差平方和を $\mathrm{RSS}_0$ とする。帰無仮説「 $E[Y] \in \mathcal{C}(X_0)$ 」のもとで
$F = \frac{(\mathrm{RSS}_0 - \mathrm{RSS})/(p - r)}{\mathrm{RSS}/(n - p)} \sim F(p - r, n - p)$
である。

典型的には、 $X_0$ は $X$ から検定したい係数の列を除いたもので、帰無仮説は「それらの係数がすべて $0$ 」である。 $F$ が大きいときに棄却する。

証明. $\mathcal{C}(X_0)$ の正規直交基底 $q_1, \dots, q_r$ を $\mathcal{C}(X)$ の正規直交基底 $q_1, \dots, q_p$ に、さらに $\mathbb{R}^n$ の正規直交基底 $q_1, \dots, q_n$ に延長し（02 の系 7.10）、 $W_1, W_2, W_3$ をそれぞれ $q_1, \dots, q_r$ 、 $q_{r+1}, \dots, q_p$ 、 $q_{p+1}, \dots, q_n$ の張る部分空間とする。 $X_0$ のハット行列を $H_0$ とすると、5.3 節で見たように $H_0 = \sum_{i \leq r}q_iq_i^{\top}$ 、 $H = \sum_{i \leq p}q_iq_i^{\top}$ なので、 $W_1, W_2, W_3$ への直交射影は $H_0$ 、 $H - H_0$ 、 $I_n - H$ である。 $(I_n - H_0)Y = (H - H_0)Y + (I_n - H)Y$ は直交分解なので、ピタゴラスの定理より $\mathrm{RSS}_0 - \mathrm{RSS} = \lVert (H - H_0)Y \rVert^2$ 。帰無仮説のもとで $\mu = E[Y] \in \mathcal{C}(X_0)$ なので $(H - H_0)\mu = (I_n - H)\mu = 0$ であり、第2章定理 2.7（正規ベクトルの直交分解）より、 $(\mathrm{RSS}_0 - \mathrm{RSS})/\sigma^2 \sim \chi^2(p - r)$ と $\mathrm{RSS}/\sigma^2 \sim \chi^2(n - p)$ は独立である。 $F$ 分布の定義（第2章定義 2.14）より主張が従う。 $\square$

注意 5.17

$\mathrm{RSS}_0 - \mathrm{RSS} = \lVert (H - H_0)Y \rVert^2 \geq 0$ 。列を増やして残差平方和が増えることはない。

係数 1 つ（ $p - r = 1$ ）の $F$ 検定は、系 5.14 の両側 $t$ 検定と同じである（ $T_j = \hat{\beta}_j/\operatorname{SE}(\hat{\beta}_j)$ として $F = T_j^2$ 、問題 5.2）。 $X_0 = \mathbf{1}$ とした「説明変数はどれも効いていない」の検定を全体の $F$ 検定という。

(A3) のもとで対数尤度を $\beta, \sigma^2$ について最大化した値は $-\frac{n}{2}\log(2\pi\mathrm{RSS}/n) - \frac{n}{2}$ なので、尤度比検定（第4章定義 4.13）の統計量は $n\log(\mathrm{RSS}_0/\mathrm{RSS}) = n\log\bigl(1 + \frac{p - r}{n - p}F\bigr)$ で、 $F$ の単調増加関数である。正規線形モデルでは尤度比検定の帰無分布が、ウィルクスの定理（第4章定理 4.14）による近似でなく正確に求まっているわけである。

例 5.2 で「築年数は不要」（ $\beta_2 = 0$ ）を検定すると、面積だけのモデルの $\mathrm{RSS}_0 = 177.87$ より $F = (177.87 - 16.56)/(16.56/9) = 87.68$ で、 $T_2^2 = (-9.364)^2$ と一致する。全体の $F$ 検定は $F = 245.6$ （自由度 $(2, 9)$ 、 $p = 1.4 \times 10^{-8}$ ）である。

5.6 決定係数

定義 5.18（決定係数） $\mathrm{TSS} = \sum_i (y_i - \bar{y})^2$ （全平方和）、 $\mathrm{ESS} = \sum_i (\hat{y}_i - \bar{y})^2$ （回帰平方和）とし、 $R^2 = 1 - \mathrm{RSS}/\mathrm{TSS}$ を決定係数 (coefficient of determination) という。

定理 5.19（平方和の分解）モデルが定数項を含む（ $\mathbf{1} \in \mathcal{C}(X)$ ）とする。

$\sum_i e_i = 0$ 。したがって $\hat{y}_1, \dots, \hat{y}_n$ の平均は $\bar{y}$ に等しい。

$\mathrm{TSS} = \mathrm{ESS} + \mathrm{RSS}$ 。したがって $R^2 = \mathrm{ESS}/\mathrm{TSS}$ で、 $0 \leq R^2 \leq 1$ 。

$\mathrm{ESS} > 0$ ならば、 $R^2$ は $y$ と $\hat{y}$ の標本相関係数の 2 乗に等しい。

証明. (1) $e \in \mathcal{C}(X)^{\perp}$ かつ $\mathbf{1} \in \mathcal{C}(X)$ より $\sum_i e_i = \mathbf{1}^{\top}e = 0$ 。

(2) $y - \bar{y}\mathbf{1} = (\hat{y} - \bar{y}\mathbf{1}) + e$ で、 $\mathbf{1} \in \mathcal{C}(X)$ より $\hat{y} - \bar{y}\mathbf{1} \in \mathcal{C}(X)$ 、また $e \in \mathcal{C}(X)^{\perp}$ 。ピタゴラスの定理より $\lVert y - \bar{y}\mathbf{1} \rVert^2 = \lVert \hat{y} - \bar{y}\mathbf{1} \rVert^2 + \lVert e \rVert^2$ 。

(3) (1) より $\hat{y}$ の平均も $\bar{y}$ なので、標本相関係数は $r = \langle y - \bar{y}\mathbf{1}, \hat{y} - \bar{y}\mathbf{1} \rangle/\sqrt{\mathrm{TSS} \cdot \mathrm{ESS}}$ 。(2) の分解と $e \perp \hat{y} - \bar{y}\mathbf{1}$ より分子は ESS に等しいので、 $r^2 = \mathrm{ESS}^2/(\mathrm{TSS} \cdot \mathrm{ESS}) = R^2$ 。 $\square$

例 5.2 では $\mathrm{TSS} = 920.26$ 、 $\mathrm{ESS} = 903.70$ 、 $\mathrm{RSS} = 16.56$ で、 $R^2 = 0.982$ である。定数項がないと分解は一般に成り立たず、 $1 - \mathrm{RSS}/\mathrm{TSS}$ は負にもなりうる（定数項のないモデルで中心化しない $1 - \mathrm{RSS}/\lVert y \rVert^2$ を $R^2$ と表示するソフトもあり、定数項のあるモデルの値とは比べられない）。注意 5.17 の 1 より説明変数を増やすと $R^2$ は減らないので、自由度調整済み決定係数 $\bar{R}^2 = 1 - \frac{\mathrm{RSS}/(n - p)}{\mathrm{TSS}/(n - 1)}$ も使われる（例 5.2 では $0.978$ ）。モデルの選び方は第6章で扱う。定数項があるとき、全体の $F$ 統計量は $F = \frac{R^2/(p - 1)}{(1 - R^2)/(n - p)}$ と書ける（ $\mathrm{RSS}_0 = \mathrm{TSS}$ と定理 5.19 の 2 から）。

注意

$R^2$ が高いことは、モデルが正しいことも、係数が因果効果であることも意味しない。曲線的な関係に直線を当てはめても $R^2$ は高くなりうるし、時系列で 2 つの変数がともに時間とともに増えていれば、無関係でも $R^2$ は高くなる。逆に $R^2$ が低くても、 $n$ が大きければ係数は精度よく推定できる。目的変数を $\log y$ に変えたモデルとは $R^2$ を比べられない（TSS が違う）。

5.7 多重共線性と分散拡大要因

定理 5.20（フリッシュ–ウォー–ロヴェルの定理, Frisch–Waugh–Lovell theorem） $X$ の $\beta_j$ に対応する列を $x_j$ 、残りの列からなる行列を $X_{-j}$ 、そのハット行列を $H_{-j}$ とし、 $\tilde{x}_j = (I_n - H_{-j})x_j$ （ $x_j$ をほかの列に回帰した残差）とおく。 $\operatorname{rank} X = p$ ならば $\tilde{x}_j \neq 0$ で
$\hat{\beta}_j = \frac{\tilde{x}_j^{\top}y}{\tilde{x}_j^{\top}\tilde{x}_j}$

証明. $\tilde{x}_j = 0$ なら $x_j = H_{-j}x_j \in \mathcal{C}(X_{-j})$ となり、 $X$ の列が一次従属になる。 $\hat{y} = X_{-j}\hat{\beta}_{-j} + x_j\hat{\beta}_j$ （ $\hat{\beta}_{-j}$ は残りの係数）と $\tilde{x}_j$ の内積をとる。 $\tilde{x}_j \perp \mathcal{C}(X_{-j})$ より $\tilde{x}_j^{\top}\hat{y} = \tilde{x}_j^{\top}x_j\hat{\beta}_j$ 。 $x_j - \tilde{x}_j = H_{-j}x_j \in \mathcal{C}(X_{-j})$ も $\tilde{x}_j$ と直交するので $\tilde{x}_j^{\top}x_j = \tilde{x}_j^{\top}\tilde{x}_j$ 。一方 $\tilde{x}_j = x_j - H_{-j}x_j \in \mathcal{C}(X)$ で $y - \hat{y} \perp \mathcal{C}(X)$ なので $\tilde{x}_j^{\top}\hat{y} = \tilde{x}_j^{\top}y$ 。 $\square$

つまり $\hat{\beta}_j$ は、「 $x_j$ のうち、ほかの説明変数の一次式では説明できない部分 $\tilde{x}_j$ 」に $y$ を（定数項なしで）単回帰した傾きである。これが「ほかの変数を一定にしたときの効果」の正確な意味である。ただしこれは、手元のデータでほかの説明変数の一次式で説明できる部分を除いたうえでの関連の大きさであって、 $x_j$ を実際に動かしたときに $y$ がどれだけ変わるかという因果効果とは限らない。モデルに入っていない変数で $x_j$ とも $y$ とも関係するものがあれば、係数はその分だけずれる（問題 5.4、第8章）。

系 5.21（分散拡大要因） (A1)(A2) のもとで $\operatorname{Var}(\hat{\beta}_j) = \sigma^2/\lVert \tilde{x}_j \rVert^2$ 。さらにモデルが定数項を含み、 $x_j$ が定数項以外の列ならば、 $x_j$ をほかの列（定数項を含む）に回帰したときの決定係数を $R_j^2$ 、 $S_{jj} = \sum_i (x_{ij} - \bar{x}_j)^2$ として
$\operatorname{Var}(\hat{\beta}_j) = \frac{\sigma^2}{S_{jj}} \cdot \frac{1}{1 - R_j^2}$
である。 $\mathrm{VIF}_j = 1/(1 - R_j^2)$ を分散拡大要因 (variance inflation factor) という。

証明. 定理 5.20 より $\hat{\beta}_j = \tilde{x}_j^{\top}Y/\lVert \tilde{x}_j \rVert^2$ なので、共分散行列の変換公式より分散は $\sigma^2\lVert \tilde{x}_j \rVert^2/\lVert \tilde{x}_j \rVert^4$ 。 $\lVert \tilde{x}_j \rVert^2$ は $x_j$ を $X_{-j}$ に回帰した残差平方和で、 $X_{-j}$ が定数項を含むので定理 5.19 より $(1 - R_j^2)S_{jj}$ に等しい。 $\square$

$\sigma^2/S_{jj}$ は、 $x_j$ がほかの説明変数と（中心化したうえで）直交していた場合の分散で、 $\mathrm{VIF}_j$ はそこから何倍に増えたかを表す。 $R_j^2 = 0.9$ なら 10 倍、 $0.99$ なら 100 倍である。例 5.2 では説明変数が 2 つなので $R_1^2 = R_2^2$ は $x_1$ と $x_2$ の相関係数の 2 乗 $0.044$ で、 $\mathrm{VIF} = 1.05$ にすぎない。

説明変数どうしが強く相関している状態を多重共線性 (multicollinearity) という。系 5.21 のとおり、多重共線性は偏りを生じさせるのではなく（不偏性は (A1) だけで成り立つ）、分散を大きくする。そのため、個々の係数は有意でないのに全体の $F$ 検定は有意、データを少し変えると係数の符号まで変わる、といったことが起こる。一方、当てはめ値 $\hat{y} = Hy$ は $\mathcal{C}(X)$ だけで決まるので、データと同じ相関構造をもつ点での予測は悪くならない。問題は、相関した変数の効果を分離して解釈しようとするときに起こる。「 $\mathrm{VIF} > 10$ なら問題」といった目安は経験則であって定理ではない。対処には、変数を減らす・まとめる、相関の弱いデータを集める、リッジ回帰で分散を抑える（第6章）などがある。

5.8 ダミー変数と交互作用

カテゴリカルな説明変数（駅からの距離の区分、地域、曜日など）はダミー変数で表す。 $K$ 個の水準をもつ変数なら、ある水準を基準に選び、残りの $K - 1$ 個の水準 $k$ について「水準 $k$ なら $1$ 、そうでなければ $0$ 」という列を加える。その係数は、ほかの説明変数が同じときの、水準 $k$ と基準水準の平均の差である。定数項に加えて $K$ 個すべてのダミー変数を入れると、それらの和が $\mathbf{1}$ になって $\operatorname{rank} X < p$ となる（ダミー変数の罠、問題 5.5）。基準水準を変えても $\mathcal{C}(X)$ は同じなので、当てはめ値・RSS・ $R^2$ は変わらない。カテゴリ変数そのものが必要かどうかは、 $K - 1$ 個の係数の個別の $t$ 検定ではなく、まとめた $F$ 検定（定理 5.16）で調べる。

効果がほかの変数の値によって変わる場合は、積の列を加える。 $d$ をダミー変数として $Y = \beta_0 + \beta_1x + \beta_2d + \beta_3xd + \varepsilon$ とすると、 $d = 0$ の群では傾きが $\beta_1$ 、 $d = 1$ の群では $\beta_1 + \beta_3$ になる。 $\beta_3$ を交互作用 (interaction) の係数といい、 $\beta_3 = 0$ の検定は「傾きが群によって違うか」を調べる。このモデルの $\beta_2$ は「 $x = 0$ での群の差」で、 $x = 0$ がデータの範囲外なら単独では意味がない。 $x$ を中心化（ $x - \bar{x}$ に置き換え）しておくと、 $\beta_2$ は「 $x$ が平均値のときの群の差」になって解釈しやすい。

5.9 残差による診断、外れ値とてこ比

定義 5.22（てこ比, leverage） $H$ の対角成分 $h_{ii}$ を観測 $i$ のてこ比という。

命題 5.23 $\operatorname{rank} X = p$ とする。

$0 \leq h_{ii} \leq 1$ 、 $\sum_{i=1}^n h_{ii} = p$ 。

モデルが定数項を含むならば $h_{ii} \geq 1/n$ 。

$\hat{y}_i = \sum_k h_{ik}y_k$ であり、(A1)(A2) のもとで $\operatorname{Var}(e_i) = \sigma^2(1 - h_{ii})$ 。

証明. (1) $H = H^2 = HH^{\top}$ の $(i, i)$ 成分を比べると $h_{ii} = \sum_k h_{ik}^2 = h_{ii}^2 + \sum_{k \neq i}h_{ik}^2 \geq h_{ii}^2$ なので $0 \leq h_{ii} \leq 1$ 。和は $\operatorname{tr} H = p$ （定理 5.8）。

(2) $\mathbf{1}/\sqrt{n}$ を最初の元とする $\mathcal{C}(X)$ の正規直交基底 $q_1 = \mathbf{1}/\sqrt{n}, q_2, \dots, q_p$ をとる（02 の系 7.10）。 $q_j$ の第 $i$ 成分を $q_{j,i}$ と書くと、 $H = \sum_j q_jq_j^{\top}$ より $h_{ii} = \sum_j q_{j,i}^2 \geq q_{1,i}^2 = 1/n$ 。

(3) 前半は $\hat{y} = Hy$ 。後半は 5.5 節の $\operatorname{Cov}(e) = \sigma^2(I_n - H)$ の対角成分である。 $\square$

てこ比は、 $x_i$ が説明変数の「中心」からどれだけ離れているかを表す（単回帰では $h_{ii} = 1/n + (x_i - \bar{x})^2/S_{xx}$ 、問題 5.1）。 $h_{ii}$ は $\hat{y}_i$ が $y_i$ にどれだけ引っ張られるかでもあり、 $h_{ii}$ が $1$ に近い点では回帰式がその点の近くを通るので、 $y_i$ が異常でも残差は小さい。てこ比の平均は $p/n$ なので、 $h_{ii} > 2p/n$ を目安に注意する（経験則）。残差の分散がそろっていないので、外れ値の判定には標準化残差 $r_i = e_i/(\hat{\sigma}\sqrt{1 - h_{ii}})$ を使う。

定理 5.24（1 つ抜きの残差） $h_{ii} < 1$ ならば、 $X$ から第 $i$ 行を除いた行列の階数も $p$ である。観測 $i$ を除いて求めた最小二乗推定量を $\hat{\beta}_{(i)}$ とすると
$y_i - x_i^{\top}\hat{\beta}_{(i)} = \frac{e_i}{1 - h_{ii}}$

証明. 第 $i$ 行を除いた行列を $X_{(i)}$ とする。 $v \neq 0$ で $X_{(i)}v = 0$ なら、 $Xv$ は $0$ でなく第 $i$ 成分以外が $0$ なので $u_i \in \mathcal{C}(X)$ となり、 $Hu_i = u_i$ より $h_{ii} = u_i^{\top}Hu_i = 1$ となって仮定に反する。

次に $y^{\ast}$ を、 $y$ の第 $i$ 成分を $x_i^{\top}\hat{\beta}_{(i)}$ に置き換えたベクトルとする。任意の $\beta$ について

\lVert y^{\ast} - X\beta \rVert^2 \geq \sum_{k \neq i}(y_k - x_k^{\top}\beta)^2 \geq \sum_{k \neq i}(y_k - x_k^{\top}\hat{\beta}_{(i)})^2 = \lVert y^{\ast} - X\hat{\beta}_{(i)} \rVert^2

なので、 $\hat{\beta}_{(i)}$ はデータ $y^{\ast}$ に対する最小二乗推定量であり、 $X\hat{\beta}_{(i)} = Hy^{\ast}$ 。 $y - y^{\ast} = (y_i - x_i^{\top}\hat{\beta}_{(i)})u_i$ に注意して第 $i$ 成分を比べると

x_i^{\top}\hat{\beta}_{(i)} = (Hy)_i - h_{ii}(y_i - x_i^{\top}\hat{\beta}_{(i)}) = \hat{y}_i - h_{ii}(y_i - x_i^{\top}\hat{\beta}_{(i)})

両辺を $y_i$ から引くと $(1 - h_{ii})(y_i - x_i^{\top}\hat{\beta}_{(i)}) = y_i - \hat{y}_i = e_i$ 。 $\square$

1 つ抜きの残差は、観測 $i$ を使わずに $y_i$ を予測したときの誤差であり、 $n$ 回当てはめ直さなくても 1 回の当てはめの $e_i$ と $h_{ii}$ から求まる。第6章の交差検証で再び使う。

系 5.25（クックの距離, Cook's distance） $h_{ii} < 1$ のとき、 $D_i = \lVert X\hat{\beta} - X\hat{\beta}_{(i)} \rVert^2/(p\hat{\sigma}^2)$ を観測 $i$ のクックの距離という。
$D_i = \frac{r_i^2}{p} \cdot \frac{h_{ii}}{1 - h_{ii}}$

証明. 定理 5.24 の証明の記号で、 $X\hat{\beta} - X\hat{\beta}_{(i)} = H(y - y^{\ast}) = (y_i - x_i^{\top}\hat{\beta}_{(i)})Hu_i$ 。 $\lVert Hu_i \rVert^2 = u_i^{\top}H^{\top}Hu_i = h_{ii}$ と定理 5.24 より $\lVert X\hat{\beta} - X\hat{\beta}_{(i)} \rVert^2 = h_{ii}e_i^2/(1 - h_{ii})^2$ 。 $e_i^2 = r_i^2\hat{\sigma}^2(1 - h_{ii})$ を代入すればよい。 $\square$

クックの距離は、観測 $i$ を 1 つ除いたときに当てはめ値全体がどれだけ動くか（影響力, influence）を測る量で、標準化残差（外れているか）とてこ比（説明変数が端にあるか）の両方が大きいときに大きくなる。例 5.2 では、てこ比の最大は物件 12 の $0.421$ （目安 $2p/n = 0.5$ 未満）、標準化残差の絶対値とクックの距離の最大はどちらも物件 4（ $\lvert r_4 \rvert = 2.43$ 、 $D_4 = 0.28$ ）で、特に影響の大きい観測はない。

診断の基本は、残差と当てはめ値の散布図（曲がった傾向は非線形性、扇形の広がりは不均一分散を示す）、標準化残差の正規 Q-Q プロット（(A3) にもとづく検定・区間をどこまで信頼できるか）、残差と観測順の図（自己相関）を描くことである。外れ値を機械的に削除してはいけない。入力ミスなら直し、性質の違う対象が混ざっていたなら除外の理由を記録し、除いた場合と除かない場合の両方の結果を報告する。

5.10 数値計算：QR 分解と条件数

$\hat{\beta} = (X^{\top}X)^{-1}X^{\top}y$ は理論のための式であって、計算の手順ではない。計算機では逆行列を作らないのはもちろん、 $X^{\top}X$ を作ることも避け、QR 分解を使う。 $X = QR$ （02 第7章定理 7.13。 $Q$ は列が正規直交な $n \times p$ 行列、 $R$ は対角成分が正の $p \times p$ 上三角行列）とすると、 $X^{\top}X = R^{\top}R$ 、 $X^{\top}y = R^{\top}Q^{\top}y$ で $R^{\top}$ は正則なので、正規方程式は

R\hat{\beta} = Q^{\top}y

と同値であり、後退代入で解ける。 $Q$ の列は $\mathcal{C}(X)$ の正規直交基底なので $H = QQ^{\top}$ 、また $(X^{\top}X)^{-1} = R^{-1}(R^{-1})^{\top}$ で、てこ比と標準誤差も $Q, R$ から求まる。

定義 5.26（条件数, condition number）列の数と階数が等しい行列 $A$ の最大・最小の特異値（02 第8章定理 8.19）を $\sigma_{\max}, \sigma_{\min}$ とするとき、 $\kappa(A) = \sigma_{\max}/\sigma_{\min}$ を $A$ の条件数という。正則な正方行列については、作用素ノルム（02 の命題 8.21）を使って $\kappa(A) = \lVert A \rVert\lVert A^{-1} \rVert$ である。

命題 5.27

正則な $A$ と $b \neq 0$ について、 $Ax = b$ 、 $A(x + \delta x) = b + \delta b$ ならば $\lVert \delta x \rVert/\lVert x \rVert \leq \kappa(A)\lVert \delta b \rVert/\lVert b \rVert$ 。

$\operatorname{rank} X = p$ ならば $\kappa(X^{\top}X) = \kappa(X)^2$ 。また QR 分解 $X = QR$ について $\kappa(R) = \kappa(X)$ 。

証明. (1) $\delta x = A^{-1}\delta b$ より $\lVert \delta x \rVert \leq \lVert A^{-1} \rVert\lVert \delta b \rVert$ 、また $\lVert b \rVert = \lVert Ax \rVert \leq \lVert A \rVert\lVert x \rVert$ 。2 式を辺々かけて $\lVert x \rVert\lVert b \rVert$ （ $b \neq 0$ より $x \neq 0$ ）で割ればよい。

(2) 特異値分解 $X = U\Sigma V^{\top}$ より $X^{\top}X = V(\Sigma^{\top}\Sigma)V^{\top}$ は固有値 $\sigma_1^2, \dots, \sigma_p^2$ をもつ正定値対称行列で、対称行列の特異値は固有値の絶対値なので（02 の問題 8.6）、 $\kappa(X^{\top}X) = \sigma_1^2/\sigma_p^2 = \kappa(X)^2$ 。 $R^{\top}R = X^{\top}X$ より、 $R$ の特異値（ $R^{\top}R$ の固有値の平方根）は $X$ の特異値と一致する。 $\square$

倍精度の浮動小数点数では、データを丸めるだけで相対誤差 $u = 2^{-53} \approx 1.1 \times 10^{-16}$ 程度の誤差が入り、命題 5.27 の 1 によれば、方程式を解く過程でその誤差は最大で条件数倍に拡大されうる。正規方程式の係数行列の条件数は $\kappa(X)^2$ なので、 $\kappa(X) = 10^8$ なら $10^{16}$ となって有効数字がすべて失われうるが、 $R\hat{\beta} = Q^{\top}y$ の係数行列の条件数は $\kappa(X)$ のままである。ライブラリが使うハウスホルダー変換による QR 分解は後退安定（計算結果が、わずかに摂動したデータに対する厳密な結果になっている）であることが知られており、残差が小さい問題では誤差はおおむね $\kappa(X)u$ の程度に収まる。ただし、最小二乗問題そのものの感度には残差が大きいときに $\kappa(X)^2$ に比例する項が含まれ、これはどの解法でも避けられない（精密な誤差評価は数値線形代数の教科書に譲る）。なお、グラム–シュミットの直交化（02 の定理 7.9）を式のとおりに浮動小数点で計算すると、直交性が失われやすい。

次の例では、2001〜2024 年の年 $t$ について $1, t, t^2$ を列とする計画行列を作り、残差が $0$ になるようにデータを作って、真の係数との相対誤差を比べる。

import numpy as np

year = np.arange(2001, 2025, dtype=float)       # 2001〜2024 年
beta = np.array([1.0, 1.0, 1.0])
for t in (year, year - year.mean()):            # 西暦のまま／中心化
    X = np.column_stack([np.ones_like(t), t, t**2])
    y = X @ beta                                # 残差が 0 のデータ
    b_ne = np.linalg.solve(X.T @ X, X.T @ y)    # 正規方程式を直接解く
    Q, R = np.linalg.qr(X)
    b_qr = np.linalg.solve(R, Q.T @ y)          # QR 分解で解く
    err = lambda b: np.linalg.norm(b - beta) / np.linalg.norm(beta)
    print(f"cond(X) = {np.linalg.cond(X):.1e}  正規方程式 {err(b_ne):.1e}  QR {err(b_qr):.1e}")

出力：

cond(X) = 3.8e+11  正規方程式 2.5e-01  QR 1.9e-06
cond(X) = 9.6e+01  正規方程式 5.5e-15  QR 1.5e-15

西暦のままでは $\kappa(X) \approx 3.8 \times 10^{11}$ 、 $\kappa(X)^2 \approx 1.5 \times 10^{23}$ は $1/u$ を大きく超え、 $X^{\top}X$ は計算機の上では特異行列と区別できない。正規方程式の解は相対誤差 25% で使いものにならないが、QR 分解では $10^{-6}$ 程度で済む。年を中心化すると $\kappa(X) \approx 96$ となり、どちらの方法でも丸め誤差の程度に収まる（誤差の細かい値は計算環境の線形代数ライブラリによって変わりうる）。

ヒント

実務では 最小二乗法を inv(X.T @ X) @ X.T @ y と書かず、QR 分解や特異値分解にもとづくライブラリの関数（NumPy の numpy.linalg.lstsq は特異値分解を使う）で解く。西暦・円単位の金額・多項式の項のように桁の大きい変数や、互いにほぼ比例する変数は、中心化や単位の変更（標準化）をしてから当てはめる。5.3 節で見たとおり $X$ を $XA$ （ $A$ は正則）に取り替えても当てはめ値は変わらないので、係数を元の単位に換算し直せば同じモデルである。

まとめ

線形回帰 $Y = X\beta + \varepsilon$ の最小二乗推定量は正規方程式 $X^{\top}X\hat{\beta} = X^{\top}y$ の解で、 $\operatorname{rank} X = p$ なら $\hat{\beta} = (X^{\top}X)^{-1}X^{\top}y$ 。
ハット行列 $H = X(X^{\top}X)^{-1}X^{\top}$ は対称・冪等で、列空間 $\mathcal{C}(X)$ への直交射影であり、 $\operatorname{tr} H = p$ 。
(A1)(A2) だけで $\hat{\beta}$ は不偏、 $\operatorname{Cov}(\hat{\beta}) = \sigma^2(X^{\top}X)^{-1}$ で、線形不偏推定量の中で最良（ガウス–マルコフ）。正規性は不要だが、等分散・無相関が崩れると標準誤差の公式が誤りになる。
$\hat{\sigma}^2 = \mathrm{RSS}/(n - p)$ は不偏。正規線形モデルでは $\hat{\beta}$ は正規分布、 $\mathrm{RSS}/\sigma^2 \sim \chi^2(n - p)$ に従い、両者は独立。ここから $t$ 検定・信頼区間・予測区間・ $F$ 検定が正確に得られる。
定数項があれば $\mathrm{TSS} = \mathrm{ESS} + \mathrm{RSS}$ で、 $R^2$ は $y$ と $\hat{y}$ の相関係数の 2 乗。 $R^2$ は変数を増やすと減らず、モデルの正しさも因果も保証しない。
$\hat{\beta}_j$ は「 $x_j$ からほかの変数で説明できる部分を除いた残差」への回帰係数で（フリッシュ–ウォー–ロヴェル）、その分散は $\mathrm{VIF}_j$ 倍に拡大される。多重共線性は偏りではなく分散の問題である。
てこ比は $0 \leq h_{ii} \leq 1$ 、和は $p$ 。1 つ抜きの残差は $e_i/(1 - h_{ii})$ 、クックの距離は標準化残差とてこ比から決まる。
計算は QR 分解で行う。正規方程式は条件数を $\kappa(X)^2$ に悪化させる。変数の中心化・標準化も有効である。

演習問題

問題 5.1 ★ 単回帰（ $X = (\mathbf{1}\ x)$ ）で $h_{ii} = \frac{1}{n} + \frac{(x_i - \bar{x})^2}{S_{xx}}$ を示せ。また $x = (1, 2, 3, 4, 10)$ のときのてこ比を求め、和が $2$ であることを確かめよ。てこ比の最も大きい点について何がいえるか。

解答

$\mathcal{C}(X) = \operatorname{span}(\mathbf{1}, x - \bar{x}\mathbf{1})$ で、 $\langle \mathbf{1}, x - \bar{x}\mathbf{1} \rangle = 0$ 、 $\lVert x - \bar{x}\mathbf{1} \rVert^2 = S_{xx}$ 。正規直交基底 $\mathbf{1}/\sqrt{n}$ 、 $(x - \bar{x}\mathbf{1})/\sqrt{S_{xx}}$ を使うと（5.3 節） $H = \frac{1}{n}\mathbf{1}\mathbf{1}^{\top} + (x - \bar{x}\mathbf{1})(x - \bar{x}\mathbf{1})^{\top}/S_{xx}$ で、その対角成分が主張の式である。 $x = (1, 2, 3, 4, 10)$ では $\bar{x} = 4$ 、 $S_{xx} = 9 + 4 + 1 + 0 + 36 = 50$ より、てこ比は $0.2 + (9, 4, 1, 0, 36)/50 = (0.38, 0.28, 0.22, 0.20, 0.92)$ で、和は $2 = p$ 。 $x = 10$ の点のてこ比は $0.92$ で、回帰直線はほぼこの点を通る。残差の分散は $0.08\sigma^2$ しかないので、この点の $y$ が誤っていても残差はほとんど大きくならない。1 つ抜きの残差 $e_5/(1 - h_{55}) = 12.5e_5$ やクックの距離で調べる必要がある。

問題 5.2 ★★ 係数 1 つの $F$ 検定（ $X_0 = X_{-j}$ ）について、 $T_j = \hat{\beta}_j/\operatorname{SE}(\hat{\beta}_j)$ として $F = T_j^2$ を示せ（ヒント：定理 5.20 の $\tilde{x}_j$ を使って $H - H_{-j}$ を表せ）。

解答

$x_j = H_{-j}x_j + \tilde{x}_j$ より $\mathcal{C}(X) = \mathcal{C}(X_{-j}) + \operatorname{span}(\tilde{x}_j)$ で、 $\tilde{x}_j \perp \mathcal{C}(X_{-j})$ 。よって $\mathcal{C}(X_{-j})$ の正規直交基底に $\tilde{x}_j/\lVert \tilde{x}_j \rVert$ を加えると $\mathcal{C}(X)$ の正規直交基底になり、 $H = H_{-j} + \tilde{x}_j\tilde{x}_j^{\top}/\lVert \tilde{x}_j \rVert^2$ 。注意 5.17 の 1 と定理 5.20 より

\mathrm{RSS}_0 - \mathrm{RSS} = \lVert (H - H_{-j})y \rVert^2 = \frac{(\tilde{x}_j^{\top}y)^2}{\lVert \tilde{x}_j \rVert^2} = \hat{\beta}_j^2\lVert \tilde{x}_j \rVert^2

一方、命題 5.9 と系 5.21 を比べると $v_{jj} = 1/\lVert \tilde{x}_j \rVert^2$ なので、 $T_j^2 = \hat{\beta}_j^2/(\hat{\sigma}^2v_{jj}) = (\mathrm{RSS}_0 - \mathrm{RSS})/\hat{\sigma}^2 = F$ （ $p - r = 1$ ）。例 5.2 の $\beta_2$ では $F = 87.68 = (-9.364)^2$ だった。

問題 5.3 ★★（この結論は正しいか）ある小売チェーンの 40 店舗のデータで、月間売上をテレビ広告費とウェブ広告費に回帰したところ、2 つの係数の $p$ 値は $0.31$ と $0.44$ で有意でなかったが、全体の $F$ 検定の $p$ 値は $0.001$ 未満、 $R^2 = 0.85$ だった。2 つの広告費の相関係数は $0.97$ である。担当者は「どちらの広告も売上に効いていない」と結論した。この結論の問題点を述べよ。

解答

説明変数が 2 つなので $R_j^2 = 0.97^2 = 0.9409$ 、 $\mathrm{VIF} = 1/(1 - 0.9409) \approx 16.9$ 。各係数の標準誤差は、2 つの広告費が無相関だった場合の約 $\sqrt{16.9} \approx 4.1$ 倍に膨らんでいる（系 5.21）。個々の係数が有意でないのは、2 つの広告がほぼ同時に増減するため効果を分離できないことの表れであって、効果が $0$ であることの証拠ではない（「有意でない」は「効果がない」ではない）。
全体の $F$ 検定は「2 つの係数がともに $0$ 」を強く棄却している。少なくとも広告費全体と売上の関連は明らかで、結論はこれと矛盾する。2 つの係数の同時検定の結果を報告する、広告費の合計を 1 つの変数にする、2 つの広告を独立に動かす実験を行う、などが考えられる。
さらに、関連があっても因果とは限らない。売上の大きい店舗ほど広告費を多く割り当てているなら、広告費は売上の原因でなく結果を反映しているかもしれない（第8章）。

問題 5.4 ★★（欠落変数の偏り）真のモデルが $Y = X_1\beta_1 + X_2\beta_2 + \varepsilon$ （(A1) をみたす）であるのに、 $X_1$ だけで最小二乗法を行った推定量 $\tilde{\beta}_1 = (X_1^{\top}X_1)^{-1}X_1^{\top}Y$ について、 $E[\tilde{\beta}_1] = \beta_1 + (X_1^{\top}X_1)^{-1}X_1^{\top}X_2\beta_2$ を示せ。また例 5.2 で、面積の単回帰の傾き $0.6148$ と重回帰の係数の間に、標本の上で $0.6148 = \hat{\beta}_1 + \hat{\delta}\hat{\beta}_2$ （ $\hat{\delta}$ は $x_2$ を $x_1$ に単回帰した傾き）という恒等式が成り立つことを示し、数値で確かめよ。

解答

$E[\tilde{\beta}_1] = (X_1^{\top}X_1)^{-1}X_1^{\top}(X_1\beta_1 + X_2\beta_2) = \beta_1 + (X_1^{\top}X_1)^{-1}X_1^{\top}X_2\beta_2$ 。偏りは、省いた変数が入れた変数と相関し（ $X_1^{\top}X_2 \neq O$ ）、かつ $\beta_2 \neq 0$ のときに生じる。

恒等式： $X_1 = (\mathbf{1}\ x_1)$ とし、重回帰の結果を $y = X_1\hat{b} + x_2\hat{\beta}_2 + e$ （ $\hat{b} = (\hat{\beta}_0, \hat{\beta}_1)^{\top}$ ）と書く。両辺に $(X_1^{\top}X_1)^{-1}X_1^{\top}$ をかけると、 $X_1^{\top}e = 0$ より、単回帰の係数は $\hat{b} + (X_1^{\top}X_1)^{-1}X_1^{\top}x_2\ \hat{\beta}_2$ に等しい。 $(X_1^{\top}X_1)^{-1}X_1^{\top}x_2$ は $x_2$ を $x_1$ に単回帰した係数で、その傾きが $\hat{\delta}$ である。数値： $\bar{x}_2 = 196/12$ 、 $S_{12} = \sum_i (x_{i1} - 68)(x_{i2} - \bar{x}_2) = \sum_i x_{i1}x_{i2} - 12 \cdot 68 \cdot \bar{x}_2 = 13051 - 13328 = -277$ より $\hat{\delta} = -277/1964 = -0.1410$ で、 $0.5531 + (-0.1410)(-0.4373) = 0.5531 + 0.0617 = 0.6148$ 。広い物件ほど新しく（ $\hat{\delta} < 0$ ）、古いほど安い（ $\hat{\beta}_2 < 0$ ）ので、単回帰は「新しさ」の効果の一部を面積の効果として数えていた。

問題 5.5 ★★ $K$ 水準のカテゴリ変数だけを説明変数とするモデル $Y_i = \beta_0 + \sum_{k=2}^K \beta_kd_{ik} + \varepsilon_i$ （ $d_{ik}$ は観測 $i$ が水準 $k$ なら $1$ 、そうでなければ $0$ ）を考え、各水準に観測が少なくとも 1 つあるとする。(1) 当てはめ値は各観測が属する水準の標本平均であり、 $\hat{\beta}_0$ は水準 1 の平均、 $\hat{\beta}_k$ は水準 $k$ と水準 1 の平均の差であることを示せ。(2) 定数項に加えて水準 1 のダミー変数の列も入れると $\operatorname{rank} X < p$ となることを示せ。

解答

(1) 水準 $k$ の指示ベクトルを $d_k$ （ $k = 1, \dots, K$ ）、水準 $k$ の観測数を $n_k$ とする。 $\mathbf{1} = d_1 + \cdots + d_K$ より $\mathcal{C}(X) = \operatorname{span}(d_1, \dots, d_K)$ 。 $d_k$ どうしは $1$ の立つ位置が重ならないので互いに直交し、 $\lVert d_k \rVert^2 = n_k$ 。よって（5.3 節） $H = \sum_k d_kd_k^{\top}/n_k$ で、観測 $i$ が水準 $k$ に属するなら $\hat{y}_i = d_k^{\top}y/n_k = \bar{y}_k$ （水準 $k$ の平均）。 $\operatorname{rank} X = K = p$ なので係数は $X\hat{\beta} = \hat{y}$ から一意に決まり、水準 1 の観測から $\hat{\beta}_0 = \bar{y}_1$ 、水準 $k$ の観測から $\hat{\beta}_0 + \hat{\beta}_k = \bar{y}_k$ 、すなわち $\hat{\beta}_k = \bar{y}_k - \bar{y}_1$ 。

(2) 列 $\mathbf{1}, d_1, \dots, d_K$ は $d_1 + \cdots + d_K - \mathbf{1} = 0$ をみたすので一次従属であり、 $\operatorname{rank} X \leq K < K + 1 = p$ 。

問題 5.6 ★★（実装のどこが危ないか）2001〜2024 年の年次売上（円単位で $10^{10}$ 程度）を、年 $t$ 、 $t^2$ 、広告費（円単位）に回帰するため、次のように書いた（ t, ad, y は長さ 24 の配列）。どこが危ないか、どう直すべきかを述べよ。

X = np.column_stack([np.ones(24), t, t**2, ad])
beta = np.linalg.inv(X.T @ X) @ X.T @ y

解答

列の大きさの桁が大きく違い、しかも 2001〜2024 の範囲では $t$ と $t^2$ がほぼ比例するので、 $\kappa(X)$ が非常に大きい（ $1, t, t^2$ だけでも約 $3.8 \times 10^{11}$ 。5.10 節の例）。 $X^{\top}X$ の条件数は $\kappa(X)^2$ （命題 5.27）で $1/u \approx 9 \times 10^{15}$ を超えるので、 $X^{\top}X$ は数値的に特異であり、エラーも出ないまま意味のない係数が得られうる。直し方：(i) $t$ を中心化し（ $t - 2012.5$ など）、金額は百万円単位にするか標準化する。(ii) np.linalg.lstsq(X, y, rcond=None) や QR 分解で解く。(iii) np.linalg.cond(X) で条件数を確かめる。中心化や単位の変更は $X \mapsto XA$ の形なので、当てはめ値は変わらず、係数を換算し直せば同じモデルである。統計的にも、24 点の年次データの 2 次の傾向で先を予測するのは外挿の危険が大きく、誤差の自己相関で (A2) が崩れていないかも確かめる必要がある（5.4 節の TIP）。

問題 5.7 ★★★（不均一分散）(A1) をみたし $\operatorname{Cov}(\varepsilon) = \operatorname{diag}(\sigma_1^2, \dots, \sigma_n^2)$ であるとき、単回帰の傾きについて $\operatorname{Var}(\hat{\beta}_1) = \sum_i (x_i - \bar{x})^2\sigma_i^2/S_{xx}^2$ を示せ。さらに $\bar{\sigma}^2 = \frac{1}{n}\sum_i \sigma_i^2$ とし、 $(x_i - \bar{x})^2$ が大きい観測ほど $\sigma_i^2$ が大きい（ $(x_i - \bar{x})^2 > (x_k - \bar{x})^2$ なら $\sigma_i^2 \geq \sigma_k^2$ ）ならば $\operatorname{Var}(\hat{\beta}_1) \geq \bar{\sigma}^2/S_{xx}$ であることを示せ。

解答

$\sum_i (x_i - \bar{x})\bar{y} = 0$ より $\hat{\beta}_1 = \sum_i (x_i - \bar{x})Y_i/S_{xx}$ は $Y$ の一次式で、 $Y_i$ は無相関だから $\operatorname{Var}(\hat{\beta}_1) = \sum_i (x_i - \bar{x})^2\sigma_i^2/S_{xx}^2$ 。 $a_i = (x_i - \bar{x})^2$ 、 $b_i = \sigma_i^2$ とおくと、 $\sum_i a_i = S_{xx}$ より

\operatorname{Var}(\hat{\beta}_1) - \frac{\bar{\sigma}^2}{S_{xx}} = \frac{1}{S_{xx}^2}\Bigl(\sum_i a_ib_i - \frac{1}{n}\sum_i a_i\sum_k b_k\Bigr) = \frac{1}{2nS_{xx}^2}\sum_i\sum_k (a_i - a_k)(b_i - b_k)

（最後の等号は右辺を展開すれば確かめられる）。仮定より $a_i > a_k$ なら $b_i \geq b_k$ なので各項は $0$ 以上で、主張が従う。

つまり、等分散を仮定した公式 $\sigma^2/S_{xx}$ に平均的な分散を入れると、真の分散を過小評価する。実際、通常の $\hat{\sigma}^2$ の期待値は $\sum_i (1 - h_{ii})\sigma_i^2/(n - 2)$ （定理 5.12 の証明と同じ計算）で、和が $1$ の重みによる加重平均だが、てこ比の大きい観測ほど重みが小さいので、この状況では過小評価はむしろ強まる。頑健な標準誤差（5.4 節の TIP）はこの問題を避けるためのものである。