回帰分析のうちの1つである線形回帰を解説する。
線形回帰（重回帰分析）のモデルを与え、回帰係数の最小二乗推定量の導出や回帰係数に関する検定統計量の導出を行っていく。
R言語での線形回帰の実行方法は以下の記事を参照。
-
【R言語】関数lmによる線形回帰 単回帰分析・重回帰分析 回帰分析で予測する
R言語の回帰分析の実行方法を紹介していきます。回帰分析の実行方法だけでなく、回帰分析により作ったモデルでテストデータの予測も行っていきます。 実行方法では、回帰分析の結果の見方やその抽出方法を詳しく解 ...
続きを見る
線形回帰
線形回帰とは回帰係数と説明変数の線形結合で表されるモデルを用いた回帰分析のことをいう。目的変数と説明変数が比例関係にあるときに有効である分析方法である。線形回帰の用途として予測や目的変数に影響を及ぼす因子の特定などがある。
線形回帰
線形回帰は、回帰係数\(\beta_j,\ i= 0,1 \ldots, k\)、説明変数\(x_{ij}, \ i = 1, 2, \ldots, n, j = 1, 2, \ldots, k\)を用いて次のように表される。
\begin{align}y_i&= \beta_0 + \beta_1x_{i1} + \beta_k x_{ik} + \varepsilon_i,\quad i = 1, 2, \ldots, n\end{align}
ここに、\(\varepsilon_i\ i.i.d. \sim N(0, \sigma^2)\)である。また、目的変数のベクトルを\(\boldsymbol{y} = (y_1, y_2, \ldots, y_n)^T\)、回帰係数ベクトルを\(\boldsymbol{\beta} = (\beta_0, \beta_1, \ldots, \beta_k)^T\)、説明変数のベクトルを
\(\boldsymbol{X}_j = (1, x_{1j},x_{2j} \ldots . x_{nj})^T,\ j = 1,2 ,\ldots , n \)
とすると次のようにも表せる。
\begin{align}\boldsymbol{y} &= \begin{pmatrix}\boldsymbol{X}_1 \boldsymbol{\beta} + \varepsilon_1\\ \boldsymbol{X}_2 \boldsymbol{\beta} + \varepsilon_2 \\\vdots \\ \boldsymbol{X}_n \boldsymbol{\beta} + \varepsilon_n \end{pmatrix} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} , \end{align}
ここに、
\begin{gather}\boldsymbol{X} &= \begin{pmatrix} \boldsymbol{X}_1^T\\ \boldsymbol{X}_2^T \\ \vdots \\ \boldsymbol{X}_n^T\end{pmatrix}, \quad \boldsymbol{\varepsilon} &= \begin{pmatrix} \varepsilon_1\\ \varepsilon_2 \\ \vdots\\ \varepsilon_n \end{pmatrix}.\end{gather}
また線形予測子は
\begin{align}\hat{\boldsymbol{y}} = \boldsymbol{X}\hat{\boldsymbol{\beta}}\end{align}
で与えられる。ここに\(\hat{\boldsymbol{\beta}}\)は次で与えられる回帰係数ベクトルの最小二乗推定量である。
\begin{align}\label{eq1} \hat{\boldsymbol{\beta}} = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}\tag{1} \end{align}
また、次の「回帰係数\(\beta_j\)はモデルに影響を与えないという」仮説検定を考える。
\begin{align} &H_0:\ \beta_j = 0\\ & H_1:\ \beta_j \neq 0\end{align}
検定統計量として次を用いる。
\begin{align}\label{eq2}T = \cfrac{\hat{\beta}_j - \beta_j}{\sqrt{(\boldsymbol{X}^T\boldsymbol{X})^{i+1, i+1}S^2} } \sim t_{n-k-1} ,\tag{2}\end{align}
ここに\((\boldsymbol{X}^T\boldsymbol{X})^{i + 1, i + 1}\)は\((\boldsymbol{X}^T\boldsymbol{X})^{-1}\)の\(i+1\)行目、\( i+1\)列目の成分であり、\(S^2 = \boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon} / (n-k-1)\)。さらにこの検定の有意水準\(\alpha\)の棄却域は次で与えられる。
\begin{align}(-\infty , -t_{n-k-1}(\alpha / 2) )\cup (t_{n-k-1}(\alpha / 2) , \infty).\end{align}
線形予測子\(\boldsymbol{X}\hat{\boldsymbol{\beta}}\)の\(\hat{\boldsymbol{\beta}}\)を求めることがいわゆる最小二乗法である。目的変数と説明変数の観測値\(\boldsymbol{y}\)、\(\boldsymbol{X}\)が与えられ、\(\boldsymbol{y}\)を\(\boldsymbol{XB}\)で予測したいとする。\(\boldsymbol{y}\)とのずれが最小になるような\(\boldsymbol{B}\)が\(\hat{\boldsymbol{\beta}}\)となる。
回帰係数の最小二乗推定量と検定
回帰係数の推定量
回帰係数ベクトル\(\boldsymbol{\beta}\)の最小二乗推定量\eqref{eq1}の導出を行う。\(\boldsymbol{\beta}\)の最小二乗推定量は目的変数の真値\(\boldsymbol{y}\)と予測値\(\hat{\boldsymbol{y}} = \boldsymbol{XB}\)の差の二乗を最小にする\(\boldsymbol{\beta}\)で定義される。すなわち
\begin{align}\label{eq3} (\boldsymbol{y} - \hat{\boldsymbol{y}})^T(\boldsymbol{y} - \hat{\boldsymbol{y}}) &= (\boldsymbol{y} - \boldsymbol{XB})(\boldsymbol{y} - \boldsymbol{XB}) \tag{3}\end{align}
を最小にする\(\boldsymbol{B}\)が\(\boldsymbol{\beta}\)の最小二乗推定量である。\eqref{eq3}を\(\boldsymbol{B}\)に関して偏微分すると
\begin{align} \cfrac{\partial }{\partial \boldsymbol{B}} (\boldsymbol{y} - \boldsymbol{XB})(\boldsymbol{y} - \boldsymbol{XB}) &= \cfrac{\partial }{\partial \boldsymbol{B}}\boldsymbol{y}^T\boldsymbol{y} -2 \cfrac{\partial }{\partial \boldsymbol{B}} \boldsymbol{y}^T\boldsymbol{XB} + \cfrac{\partial }{\partial \boldsymbol{B}} (\boldsymbol{XB})^T\boldsymbol{XB} \\ &= \boldsymbol{0} -2 \boldsymbol{y}^T\boldsymbol{X} + 2\boldsymbol{X}^T\boldsymbol{XB} \\ &= -2 \boldsymbol{X}^T\boldsymbol{y} + 2\boldsymbol{X}^T\boldsymbol{XB}.\end{align}
上式を\(\boldsymbol{0}\)と置いて方程式を解くことで次の最小二乗解を得る。
\begin{align}& -2 \boldsymbol{X}^T\boldsymbol{y} + 2\boldsymbol{X}^T\boldsymbol{XB} = \boldsymbol{0} \\ &\Leftrightarrow \boldsymbol{X}^T\boldsymbol{XB} = \boldsymbol{X}^T \boldsymbol{y} \\ &\Leftrightarrow \boldsymbol{B} = (\boldsymbol{X}^T\boldsymbol{X})^{-1} \boldsymbol{X}^T\boldsymbol{y}.\end{align}
したがって、\eqref{eq1}の\(\boldsymbol{\beta}\)の最小二乗推定量\(\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^T\boldsymbol{X})^{-1} \boldsymbol{X}^T\boldsymbol{y}\)を導出することができた。□
回帰係数の検定
次に回帰係数\(\beta_j,\ j = 0,1, \ldots, k\)の検定統計量の導出を行う。回帰係数の推定量\(\hat{\boldsymbol{\beta}}\)は多変量正規分布の平均ベクトル・共分散行列の性質より、次の平均ベクトルと共分散行列をもつ。
\begin{align}\mathrm{E}[\hat{\boldsymbol{\beta}}] &= \mathrm{E}[(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}]\\ &= \mathrm{E}[(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T(\boldsymbol{X \beta} + \boldsymbol{\varepsilon})] \\&= \boldsymbol{\beta} + (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T \boldsymbol{0} \\\label{eq4} &= \boldsymbol{\beta}, \tag{4}\\ \mathrm{Var}[\hat{\boldsymbol{\beta}}] &= \mathrm{Var}[(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T(\boldsymbol{X \beta} + \boldsymbol{\varepsilon})]\\ &= \boldsymbol{0} + (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\mathrm{Var}[\boldsymbol{\epsilon}]\boldsymbol{X}(\boldsymbol{X}^T\boldsymbol{X})^{-1}\\ \label{eq5} &= \sigma^2(\boldsymbol{X}^T\boldsymbol{X})^{-1}.\tag{5}\end{align}
故に、\(\hat{\beta}_j\sim N(\beta_j, \sigma^2(\boldsymbol{X}^T\boldsymbol{X})^{jj})\)である。ここで、\(\boldsymbol{P} = \boldsymbol{X}(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\)とおく。\(\mathrm{rank}(\boldsymbol{X}) = k+1\)であることから、\(\mathrm{rank}(\boldsymbol{I}_n - \boldsymbol{P}) = n-k-1\)である。\(\boldsymbol{P}\)に関して次の性質が成り立つ。
\begin{align}\boldsymbol{P}^2 &= (\boldsymbol{X}(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T)^2\\&= \boldsymbol{X}(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T,\\ (\boldsymbol{I} - \boldsymbol{P})^2 &= \boldsymbol{I}^2 -2\boldsymbol{P} + \boldsymbol{P}^2\\&= \boldsymbol{I} - \boldsymbol{P}.\end{align}
\(\boldsymbol{P}\)と\(\boldsymbol{I} - \boldsymbol{P}\)はべき乗をしても同じになり冪等行列であることがいえる。残差\(\boldsymbol{e}\)とすると
\begin{align}\boldsymbol{e} &= \boldsymbol{y} - \boldsymbol{X}\hat{\boldsymbol{\beta}}\\ &= \bigl[\boldsymbol{I} - \boldsymbol{X}(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\bigr] \boldsymbol{y}\\&= (\boldsymbol{I} - \boldsymbol{P})\boldsymbol{y}\end{align}
で与えられる。今
\begin{align}(\boldsymbol{I} - \boldsymbol{P}) \boldsymbol{X} = \boldsymbol{X} - \boldsymbol{X}(\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{X} = \boldsymbol{0}\end{align}
より、
\begin{align} \boldsymbol{e} &= (\boldsymbol{I} - \boldsymbol{P})\boldsymbol{y} \\&= (\boldsymbol{I} - \boldsymbol{P}) (\boldsymbol{X\beta} + \boldsymbol{\varepsilon}) \\&= \boldsymbol{0} + (\boldsymbol{I} - \boldsymbol{P})\boldsymbol{\varepsilon}\\ &= (\boldsymbol{I} - \boldsymbol{P})\boldsymbol{\varepsilon}.\end{align}故に\begin{align}\mathrm{E}[\boldsymbol{e}] &= \boldsymbol{0}.\end{align}
ここで冪等行列\((\boldsymbol{I} - \boldsymbol{P})\)について、
\begin{align}\boldsymbol{U}(\boldsymbol{I} - \boldsymbol{P})\boldsymbol{U}^T = \begin{pmatrix}\boldsymbol{I}_{n-k-1} & \boldsymbol{0} \\\boldsymbol{0} &\boldsymbol{0} \end{pmatrix}\end{align}
となるような直交行列\(\boldsymbol{U}\)が存在する。したがって、残差平方和\(\boldsymbol{e}^T\boldsymbol{e}\)は次のように表現できる。
\begin{align} \boldsymbol{e}^T\boldsymbol{e} &= \boldsymbol{\varepsilon}^T(\boldsymbol{I} - \boldsymbol{P})^T(\boldsymbol{I} - \boldsymbol{P})\boldsymbol{\varepsilon} \\ &= \boldsymbol{\varepsilon}^T(\boldsymbol{I} - \boldsymbol{P})^2\boldsymbol{\varepsilon} \\ &= \varepsilon^T(\boldsymbol{I} - \boldsymbol{P}) \boldsymbol{\varepsilon} \\ &= \boldsymbol{\varepsilon}\boldsymbol{U}^T\boldsymbol{U}(\boldsymbol{I} - \boldsymbol{P})\boldsymbol{U}^T\boldsymbol{U}\boldsymbol{\varepsilon}\\ &= \boldsymbol{z}^T\begin{pmatrix}\boldsymbol{I}_{n-k-1} & \boldsymbol{0} \\ \boldsymbol{0} &\boldsymbol{0}\end{pmatrix}\boldsymbol{z} \\&= \sum_{i= 1}^{n-k-1} Z_i^2, \end{align}
ここに、\(Z_i ,\ i.i.d.\sim N(0, \sigma^2),\ i = 1, \ldots, n\)である。 故に\( \boldsymbol{e}^T\boldsymbol{e} / \sigma^2 \sim \chi_{n-k-1}^2\)であり、\(\boldsymbol{e}^T\boldsymbol{e} \)は\(\hat{\boldsymbol{\beta}}\)に依らず、独立な正規変数の二乗和で表されることから、これらの分布は独立である。よって\( S^2=\boldsymbol{e}^T\boldsymbol{e} / (n-k-1) \)とすると、次の帰無仮説\(H_0:\ \beta_i = 0\ i = 0, 1, \ldots, k\)の検定統計量\eqref{eq2}を得る。
\begin{align}T &= \cfrac{(\hat{\beta}_i - \beta_i) / \sqrt{\sigma^2(\boldsymbol{X}^T\boldsymbol{X})^{i+1. i+1}} }{\sqrt{S^2 / \sigma^2 }} \\ &= \cfrac{\hat{\beta}_i - \beta_i}{ \sqrt{(\boldsymbol{X}^T\boldsymbol{X})^{i+1, i+1} S^2 }} \sim T_{n-k-1}.\end{align}
