線形回帰

1つの従属変数と1つの独立変数がある線形回帰の例。

線形回帰（せんけいかいき、英: linear regression）とは、説明変数（独立変数ともいう）に対して目的変数（従属変数、あるいは反応変数ともいう）が線形またはそれから近い値で表される状態。線形回帰は統計学における回帰分析の一種であり、非線形回帰と対比される。

線形回帰のうち、説明変数が1つの場合を線形単回帰（simple linear regression）や単純線形回帰や単変量線形回帰（univariate linear regression）、2つ以上の場合を線形重回帰（multiple linear regression）や多重線形回帰や多変量線形回帰（multivariate linear regression）と呼ぶ。単回帰と呼んだ場合、単変量の回帰のことであるが、多くの場合は非線形を含めずに線形単回帰の事を指す。

概要

線形回帰では，データから推定される線形予測関数を用いて関係性がモデル化される。このようなモデルは線形モデルと呼ばれる。 説明変数（または予測変数）に対して目的変数の条件付き期待値は、アフィン写像で与えられる。（通常は条件付き期待値だが、条件付メジアンまたは他の分位数を用いることもある。）

線形回帰が非線形回帰に比べて用いられる頻度が高いのは、未知のパラメータに線形に依存するモデルの方が、パラメータに非線形に依存するモデルよりもフィッティングが容易で、推定値の統計的性質を決定しやすいためである。

線形回帰が取り扱う範囲は、予測変数の値を与えられた応答の条件付き確率分布に限る。 全ての変数の同時確率分布は多変量解析の領域として、ここでは扱わない。

線形回帰の用途

線形回帰は多くの実用的な用途があり、大まかには以下の二種類の用途に分類される。

○予測、予想、またはエラーの削減を目的とする。 →線形回帰は、応答変数と説明変数の値の観測されたデータセットに予測モデルを適合させるために使用できる。 説明変数の追加値が収集された場合、このモデルから応答変数を予測できる。

○説明変数の変動に起因する応答変数の変動を説明することを目的とする。 →線形回帰分析を適用して、応答と説明変数の関係の強さを定量化できる。 これにより各説明変数が応答と全く線形関係を持たないかどうかを判断したり、説明変数のどのサブセットに応答に関する冗長な情報が含まれているかを特定できる。

線形モデルのフィッティング方法

線形回帰モデルは多くの場合、最小二乗法を用いてフィッティングされる。 それ以外のフィッティング方法としては、最小絶対値法や、リッジ回帰（L2ノルムペナルティ）やラッソ回帰（L1ノルムペナルティ）のように、最小二乗コスト関数のペナルティ付きバージョンを最小化する方法などがある。 逆に最小二乗法は、線形モデルではないモデルのフィットにも使用できる。 このように、「最小二乗法」と「線形モデル」という言葉は密接に関連しているが、同義ではない。

基本モデル

線形回帰モデルは、目的変数 Y と説明変数 X_i, i = 1, ..., p および擾乱項 ε の関係を以下のようにモデル化したものである。

$Y=\beta _{0}+\beta _{1}X_{1}+\beta _{2}X_{2}+\cdots +\beta _{p}X_{p}+\varepsilon \$

ここで β₀ は切片（「定数」項）、β_i は各々の説明変数の係数であり、p は説明変数の個数である。線形回帰においては、説明変数の係数および切片の組 {β_i}_i∈[0,p) をパラメタとするモデルを与える。また、擾乱項 ε は説明変数 X とは独立である。

ベクトル・行列記法を用いれば、線形回帰モデルは以下のように表せる。

$Y=X\beta +\varepsilon \$

線形とは

線形回帰が「線形」であるのは、目的変数 Y が説明変数 X の係数 β に対して線形であるためである。たとえば

$Y=\beta _{0}+\beta _{1}x+\beta _{2}x^{2}+\varepsilon$

という回帰は x に対して明らかに線形ではないが、係数 β に対して線形であるから、線形回帰の問題に分類される。

線形単回帰

線形単回帰や単純線形回帰や単変量線形回帰の場合、説明変数は1つだけであり回帰パラメタは2つである。上式は以下のようになる。

$y=a+bx+\varepsilon \$

最小二乗法を使用した場合、${\bar {x}}$ と ${\bar {y}}$ を $x_{i}$ と $y_{i}$ の平均としたとき、パラメータ $a$ と $b$ の推定量の ${\hat {a}}$ と ${\hat {b}}$ は以下のように求まる。

$${\begin{aligned}{\hat {a}}&={\bar {y}}-{\hat {b}}\,{\bar {x}}\\{\hat {b}}&={\frac {\sum _{i=1}^{n}(x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}}}\end{aligned}}$$

同等な定式化に、線形単回帰を条件付き期待値のモデルとして陽に表すものがある。

${\mbox{E}}(y|x)=\alpha +\beta x\$

ここで、所与の x に対する y の条件付き確率分布は擾乱項の確率分布に一致する。

線形回帰の種類

最小二乗モデル

最小二乗法はカール・フリードリッヒ・ガウスが1820年代に発展させた。本方法は、擾乱項 ε_i の振る舞いに次のような仮定をする（ガウス＝マルコフ仮定）。

• 擾乱 ε_i の期待値は 0 である

  $E[\varepsilon ]=0$

• 擾乱 ε_i は相互に無相関である（統計的な独立の仮定よりは弱い）

  $\operatorname {cov} (\varepsilon _{i},\varepsilon _{j})=0,\qquad i\neq j.$

• 擾乱 ε_i は等分散、すなわちみな等しい分散をもつ（ガウス＝マルコフの定理も参照）

  $V[\varepsilon _{i}]=\sigma ^{2},\qquad \forall i\in [n].$

以上の仮定は、最小二乗法がある意味で最適なパラメタの推定量を与えることを保証する。

説明変数の個数が p 個のモデルを考えると、線形回帰によって決定すべきパラメタは係数 β₁,..., β_p と切片 β₀ の p + 1 個である。目的変数と説明変数の測定結果の組 (y_k; x_k1,...,x_kp) を1つのデータとし、n 個のデータを用いた線形回帰は以下のように表すことができる。

${\begin{bmatrix}y_{1}\\y_{2}\\\vdots \\y_{n}\end{bmatrix}}={\begin{bmatrix}1&x_{11}&x_{12}&\dots &x_{1p}\\1&x_{21}&x_{22}&\dots &x_{2p}\\\vdots &\vdots &\vdots &&\vdots \\1&x_{n1}&x_{n2}&\dots &x_{np}\end{bmatrix}}{\begin{bmatrix}\beta _{0}\\\beta _{1}\\\vdots \\\beta _{p}\end{bmatrix}}+{\begin{bmatrix}\varepsilon _{1}\\\varepsilon _{2}\\\vdots \\\varepsilon _{n}\end{bmatrix}}$

上記の連立方程式は、目的変数の観測値を n 成分の列ベクトル Y、説明変数の観測値および切片 β₀ の係数 (=1) を n × (p + 1) 行列 X、回帰パラメタを (p + 1) 成分の列ベクトルβ、観測ごとの擾乱を n 成分の列ベクトル ε とすれば、行列の記法を用いて以下のように表せる。

$Y=\mathbf {X} \beta +\varepsilon$

n = p の場合、回帰パラメタの標準誤差は算出できない。n が p より小さい場合、パラメタは算出できない。

回帰パラメタの推定量は、

${\widehat {\beta }}=(\mathbf {X} ^{\top }\mathbf {X} )^{-1}\mathbf {X} ^{\top }{\vec {y}}$

で与えられ、ガウス＝マルコフの定理より推定量${\widehat {\beta }}$は最良線形不偏推定量になる。つまり、任意の線形不偏推定量$\beta$に対して

$V[\beta ]\geq V[{\widehat {\beta }}]$

が成立する。

回帰の二乗和 SSR は下式で与えられる。

${{\mathit {SSR}}=\sum {\left({{\hat {y_{i}}}-{\bar {y}}}\right)^{2}}={\hat {\beta }}^{\top }\mathbf {X} ^{\top }{\vec {y}}-{\frac {1}{n}}\left({{\vec {y}}^{\top }{\vec {u}}{\vec {u}}^{\top }{\vec {y}}}\right)}$

ここで ${\bar {y}}={\frac {1}{n}}\sum y_{i}$ であり ${\vec {u}}$ は n × 1 の1ベクトル（各要素が1）である。 項 ${\frac {1}{n}}y^{\top }uu^{\top }y$ は ${\frac {1}{n}}(\sum y_{i})^{2}$　と等価である。

誤差の二乗和 ESS は下式で与えられる。

${{\mathit {ESS}}=\sum {\left({y_{i}-{\hat {y_{i}}}}\right)^{2}}={\vec {y}}^{\top }{\vec {y}}-{\hat {\beta }}^{\top }\mathbf {X} ^{\top }{\vec {y}}}$

二乗和の全和 TSS' は下式で与えられる。

${{\mathit {TSS}}=\sum {\left({y_{i}-{\bar {y}}}\right)^{2}}={\vec {y}}^{\top }{\vec {y}}-{\frac {1}{n}}\left({{\vec {y}}^{\top }{\vec {u}}{\vec {u}}^{\top }{\vec {y}}}\right)={\mathit {SSR}}+{\mathit {ESS}}}$

決定係数, R² は下式で与えられる。

${R^{2}={\frac {\mathit {SSR}}{\mathit {TSS}}}=1-{\frac {\mathit {ESS}}{\mathit {TSS}}}}$

擾乱項が正規分布に従うモデル

以下では擾乱項ε_iが互いに独立な平均$0$, 分散$\sigma ^{2}$の正規分布に従うと仮定する。

残差は、観測値とモデルによる予測値の差を表し、以下のように決定される。

${\hat {\vec {\varepsilon }}}={\vec {y}}-\mathbf {X} {\hat {\beta }}\$

この時、統計量$S^{2}={\frac {{\hat {\vec {\varepsilon }}}{\;}^{\top }{\hat {\vec {\varepsilon }}}}{n-p-1}}$は分散$\sigma ^{2}$の不偏推定量（$E[S^{2}]=\sigma ^{2}$）になる。また、最小二乗推定量${\widehat {\beta }}$と統計量$S^{2}$について以下が成立することが知られている。証明は久保川(2017)や解説記事が詳しい。

1. ${\widehat {\beta }}$は多次元正規分布${\mathcal {N}}\left(\beta ,\sigma ^{2}({\boldsymbol {X}}^{\top }{\boldsymbol {X}})^{-1}\right)$に従う
2. ${\frac {(N-P-1)S^{2}}{\sigma ^{2}}}$は自由度$n-p-1$の$\chi _{n-p-1}^{2}$分布に従う
3. ${\widehat {\beta }}$と$S^{2}$は独立

上記の事実をもとに回帰係数の有意性検定、信頼区間や予測区間を構成できる。

回帰係数の有意性検定

回帰係数の推定量${\widehat {\beta }}_{i}$は正規分布${\mathcal {N}}\left(\beta _{i},\sigma ^{2}({\boldsymbol {X}}^{\top }{\boldsymbol {X}})_{ii}^{-1}\right)$に従うことから

$T={\dfrac {{\hat {\beta }}_{i}-\beta _{i}}{\sqrt {({\boldsymbol {X}}^{\top }{\boldsymbol {X}})_{ii}^{-1}S^{2}}}}$

は自由度$n-p-1$の$t$分布に従う。ここで$({\boldsymbol {X}}^{\top }{\boldsymbol {X}})_{ii}^{-1}$は行列${\boldsymbol {X}}^{\top }{\boldsymbol {X}}$の第$(i+1,i+1)$成分である。（添え字$i$は0から始まることに注意。）

これより適当な有意水準$\alpha$で

• 帰無仮説: $\beta _{i}=0$
• 対立仮説: $\beta _{i}\neq 0$

を検定することできる。

信頼区間と予測区間

値 ${\vec {x}}={\vec {x}}_{0}$における$100(1-\alpha )\%$ の信頼区間は下式で表される。

${{\vec {x_{0}}}{\widehat {\beta }}\pm t_{{\frac {\alpha }{2}},n-p-1}{\sqrt {{\vec {x_{0}}}(\mathbf {X} ^{\top }\mathbf {X} )_{}^{-1}{\vec {x_{0}}}^{\top }S^{2}}}}$

同様に値 ${\vec {x}}={\vec {x}}_{0}$における$100(1-\alpha )\%$ の予測区間は下式で表される。

${{\vec {x_{0}}}{\widehat {\beta }}\pm t_{{\frac {\alpha }{2}},n-p-1}{\sqrt {(1+{\vec {x_{0}}}(\mathbf {X} ^{\top }\mathbf {X} )_{}^{-1}{\vec {x_{0}}}^{\top })S^{2}}}}$