Skip to Main Content

★手計算とRで学ぶ統計学: 重回帰分析 (重回帰モデル)

「なぜ統計学が必要か」という問いをひもとき，実践を通じて読者の方と統計学の心理的距離を縮めるガイドです。

基本・重回帰分析 (モデル式の推定)

(データ概要)
単回帰モデルで使用したデータに駐車場の有無を追加したデータを用いて重回帰モデルを作成し，駅から徒歩20分の物件は家賃がどれくらいになるかを予測します。

\begin{array}{cccc}
  \hline
  ID & time & parking & rent \\
  \hline
 1 &  5   & 1 &  55000   \\
  \hline
 2 &  6 & 1 & 57000 \\
  \hline
3 &  7  & 0 & 45000 \\
  \hline
4&   8  & 0 & 44000 \\
  \hline
5 &  9  & 1 & 40000 \\
  \hline
 6 & 10  & 1 & 41000 \\
  \hline
  7 & 11  & 0 & 38000 \\
  \hline
  8 &12  & 0 & 32000 \\
  \hline
 9 & 13 & 1 &35000\\
  \hline
 10 & 15 & 0 & 29000 \\
  \hline
\end{array}

①モデルと実際のデータの差 (残差) の平方和を求める
駅から物件までの所要時間と駐車場の有無，家賃のデータから推定された重回帰モデルは以下の式で示されます。

\begin{align}
y = \alpha + \beta_{1}x_{1} + \beta_{2} x_{2} \\
y: 目的変数 (家賃の推定値) \\
x_{1}: 説明変数1 (所要時間) \\
x_{2}: 説明変数2 (駐車場の有無) 
\end{align}

単回帰モデルと同様，重回帰モデルでもまずは実際の値 (観測値) とモデルによって予測できる値 (推定値) の間の差 (残差) の二乗和，つまり残差平方和を求めます。
重回帰モデルにおける残差平方和は以下の式のように求められます。
\begin{align}
Q (\hat{\alpha}, \hat{\beta_{1}},  \hat{\beta_{2}}) &= \sum_{i=1}^n{\left(y_{i}-\hat{y_{i}}\right)^2} \\
&=\sum_{i=1}^n{\left\{( y_{i} - \overline{y}) - ( \hat{y_{i}} - \overline{y} )\right\}^2} \\
\hat{y_{i}}: 推定値\\
\overline{y}: yの平均値
\end{align}

ここで，\(\overline{y}\) と\(\hat{y_{i}}\) が回帰式を通る性質を利用すると
\begin{align}
\overline{y} &= \hat{\alpha} + \hat{\beta_{1}} \overline{x_{1}} + \hat{\beta_{2}} \overline{x_{2}} \\
\hat{y_{i}} &=\hat{\alpha} + \hat{\beta_{1}} x_{1i} +\hat{\beta_{2}} x_{2i} \\
&\overline{x_{1}}: 説明変数1 (所要時間) の平均値\\
&\overline{x_{2}}: 説明変数2 (駐車場の有無) の平均値 \\
&x_{1i}: 説明変数1 (所要時間) のi番目のデータ\\
&x_{2i}: 説明変数2 (駐車場の有無) のi番目のデータ 
\end{align}
とおくことができます。

上記の\(\overline{y}\), \(\hat{y_{i}}\) を残差平方和の式に代入すると，以下のようになります。
\begin{align}
\sum_{i=1}^n{\left(y_{i}-\hat{y_{i}}\right)^2} &= \sum_{i=1}^{n}[(y_{i}-\overline{y}) - \{\hat{y_{i}}- (\hat{\alpha} + \hat{\beta_{1}} \overline{x_1}+ \hat{\beta_{2}} \overline{x_2})\}]^2\\
&= \sum_{i=1}^{n}[(y_{i}-\overline{y}) - \{(\hat{\alpha} + \hat{\beta_{1}} x_{1i} + \hat{\beta_{2}} x_{2i} ) - (\hat{\alpha} + \hat{\beta_{1}} \overline{x_1}+ \hat{\beta_{2}} \overline{x_2})\}]^2\\
&= \sum_{i=1}^{n}[(y_{i}-\overline{y}) - (\hat{\alpha} +\hat{\beta_{1}} x_{1i} + \hat{\beta_{2}} x_{2i}  - \hat{\alpha} - \hat{\beta_{1}}\overline{x_1}-\hat{\beta_{2}} \overline{x_2})]^2 \\
&= \sum_{i=1}^{n}[(y_{i}-\overline{y}) -  \{ \hat{\beta_{1}} (x_{1i}-\overline{x_1}) - +\hat{\beta_{2}}( x_{2i}  -  \overline{x_2})\}]^2 \\
&= \sum_{i=1}^{n}\{(y_{i}-\overline{y}) -   \hat{\beta_{1}} (x_{1i} - \overline{x_1}) -  \hat{\beta_{2}}( x_{2i}  -  \overline{x_2})\}^2 \\
\end{align}

ここで気づいた方もいるかもしれませんが，式を整理する過程で\(\hat{\alpha} \) が消え去ります。
今後の計算で\(\hat{\alpha} \)がないと少し不便なので，以下の式変形を利用して式の中に\(\hat{\alpha} \) を復活させます。
\begin{align}
\overline{y} = \hat{\alpha} + \hat{\beta_{1}}\overline{x_{1}} + \hat{\beta_{2}} \overline{x_{2}} \\
\overline{y} - (\hat{\alpha} + \hat{\beta_{1}} \overline{x_{1}}+ \hat{\beta_{2}} \overline{x_{2}}) = 0
\end{align}

0を代入しても計算結果は変わらないので，最終的に残差平方和は以下の式で表されます。
\begin{align}
Q (\hat{\alpha}, \hat{\beta_{1}},  \hat{\beta_{2}}) &=  \sum_{i=1}^{n}[(y_{i}-\overline{y}) - \{\hat{\beta_{1}} (x_{1i} - \overline{x_1}) - \hat{\beta_{2}}(x_{2i} - \overline{x_2}) + \{\overline{y} - (\hat{\alpha} + \hat{\beta_{1}} \overline{x_{1}} + \hat{\beta_{2}} \overline{x_{2}})\} \}]^2 \\

&= \sum_{i=1}^{n} \{ (y_{i} - \overline{y})^2 -2\hat{\beta_{1}}(y_{i} - \overline{y})(x_{1i} - \overline{x_{1}}) \\
&+ \hat{\beta_{1}}(x_{1i} - \overline{x_{1}})^2 +2 (y_{i} - \overline{y}) (\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}}) \\
&-2\hat{\beta_{2}} (y_{i} - \overline{y})(x_{2i} - \overline{x_{2}}) + 2 (x_{1i} - \overline{x_{1}}) (\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}}) \\
&+2\hat{\beta_{1}}\hat{\beta_{2}}(x_{1i} - \overline{x_{1}})(x_{2i} - \overline{x_{2}}) + (\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}})^2 \\
&-2\hat{\beta_{2}}(x_{2i} - \overline{x_{2}})(\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}}) +\hat{\beta_{2}}^2(x_{2i} - \overline{x_{2}})^2 \\
\}
\end{align}

② ①で求めた式を整理する
ここまででモデルの残差平方和は求められましたが，式があまりにも煩雑すぎます。そこで，展開式の中に含まれる式を以下のように整理します。

\begin{align}
\sum_{i=1}^{n} (y_{i} - \overline{y})^2 &= S_{yy}\\
\sum_{i=1}^{n} -2\hat{\beta_{1}}(y_{i} - \overline{y})(x_{1i} - \overline{x_{1}}) &= -2\hat{\beta_{1}}S_{1y}\\
\sum_{i=1}^{n} \hat{\beta_{1}}^2(x_{1i} - \overline{x_{1}})^2 &= \hat{\beta_{1}}^2S_{11}\\
\sum_{i=1}^{n} 2 (y_{i} - \overline{y}) (\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}}) &= 0\\
\sum_{i=1}^{n} -2\hat{\beta_{2}} (y_{i} - \overline{y})(x_{2i} - \overline{x_{2}}) &= -2\hat{\beta_{2}}S_{2y}\\
\sum_{i=1}^{n} 2 (x_{1i} - \overline{x_{1}}) (\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}}) &= 0\\
\sum_{i=1}^{n} 2\hat{\beta_{1}}\hat{\beta_{2}}(x_{1i} - \overline{x_{1}})(x_{2i} - \overline{x_{2}}) &= 2\hat{\beta_{1}}\hat{\beta_{2}}S_{12} \\
\sum_{i=1}^{n} (\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}})^2 &= n(\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}})^2 \\
\sum_{i=1}^{n} -2\hat{\beta_{2}}(x_{2i} - \overline{x_{2}})(\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}}) &= 0\\
\sum_{i=1}^{n} \hat{\beta_{2}}^2(x_{2i} - \overline{x_{2}})^2 &= \hat{\beta_{2}}^2S_{22}\\
\end{align}

上の式において，\(S_{12}\) は偏差の積和 (\(x_{1}\) の偏差と\(x_{2}\) の偏差をかけ，全て足し合わせたもの) を表しています。\(S_{1y}\) や\(S_{yy}\) についても同様です。
よって，重回帰モデルの残差平方和は以下の式で表されます。
\begin{align}
Q (\hat{\alpha}, \hat{\beta_{1}},  \hat{\beta_{2}}) = S_{yy} -2\hat{\beta_{1}}S_{1y} + \hat{\beta_{1}}^2S_{11} -2\hat{\beta_{2}}S_{2y} + 2\hat{\beta_{1}}\hat{\beta_{2}}S_{12} + n(\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}})^2 + \hat{\beta_{2}}^2S_{22} 
\end{align}

③ ②で求めた式を\(\hat{\beta_{1}}, \hat{\beta_{2}} \) でそれぞれ偏微分する
単回帰分析と同様に，\(Q (\hat{\alpha}, \hat{\beta_{1}},  \hat{\beta_{2}}) \)をそれぞれ\(\hat{\beta_{1}}\), \(\hat{\beta_{2}}\)で偏微分します。
なお，\(S_{yy}\), \(n(\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}})^2\) は定数項であるため (\((\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}}) = 0\) であることを利用して式変形を行いました)，微分すると0になります。そのため，偏微分の際は一旦ないものとして扱います (\(n(\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}})^2\) は後程登場します)。

まず，式 (1) を\(\hat{\beta_{1}} \) のみを変数として以下の通り微分します。
\begin{align}
\frac{\partial Q (\hat{\alpha}, \hat{\beta_{1}}, \hat{\beta_{2}}) }{\partial\hat{\beta_{1}}} = 2\hat{\beta_1}S_{11} -2S_{1y} +2\hat{\beta_2}S_{12}\;...(2)
\end{align}

次に，式 (1) を\(\hat{\beta_{2}}\) のみを変数として以下の通り微分します。
\begin{align}
\frac{\partial Q (\hat{\alpha}, \hat{\beta_{1}}, \hat{\beta_{2}}) }{\partial\hat{\beta_{2}}} = 2\hat{\beta_2}S_{22} -2S_{2y} +2\hat{\beta_1}S_{12}\;...(3)
\end{align}

④ ③で求めた式(2), (3) について正規方程式を求める
式(2)，式(3) について，それぞれ解が0となる方程式を作り (※2)，連立することで\(\hat{\alpha}\), \(\hat{\beta}\) の値を求めます。この連立方程式は正規方程式と呼ばれます。
式(2)，(3) について，以下のように方程式を整理します。
\begin{align}
2\hat{\beta_1}S_{11} -2S_{1y} +2\hat{\beta_2}S_{12} &= 0 \\
2\hat{\beta_1}S_{11} +2\hat{\beta_2}S_{12} &= 2S_{1y} \\
\hat{\beta_1}S_{11} +\hat{\beta_2}S_{12} &= S_{1y}\;...(2)'\\
\end{align}

\begin{align}
2\hat{\beta_2}S_{22} -2S_{2y} +2\hat{\beta_1}S_{12} &= 0 \\
2\hat{\beta_2}S_{22} +2\hat{\beta_1}S_{12} &= 2S_{2y} \\
\hat{\beta_2}S_{22} +\hat{\beta_1}S_{12} &= S_{2y}\;...(3)' \\
\end{align}

したがって，正規方程式は以下のように示されます。
\begin{align}
\begin{cases}\hat{\beta_1}S_{11} +\hat{\beta_2}S_{12} &= S_{1y}\\ \hat{\beta_1}S_{12} +\hat{\beta_2}S_{22} &= S_{2y}\end{cases}
\end{align}

⑤正規方程式にデータからわかる値を代入し，\(\hat{\beta_{1}}, \hat{\beta_{2}} \)の値を求める
データからわかる値 \(\left(S_{11},\; S_{22},\; S_{12},\; S_{1y},\; S_{2y} \right)\) を代入すると以下の連立方程式が得られます。
\begin{align}
\begin{cases}92.4\hat{\beta_1} -5\hat{\beta_2} &= -247600\\ -5\hat{\beta_1}+2.5\hat{\beta_2}  &= 20000\end{cases}
\end{align}

上の連立方程式を解くと，\(\hat{\beta_{1}}\fallingdotseq -2519.42,\; \hat{\beta_{2}}\fallingdotseq 2961.17\) という値が得られます。

⑥\(\hat{\alpha}\) の値を求め，重回帰モデルを完成させる
⑤で求めた\(\hat{\beta_{1}}, \hat{\beta_{2}} \)，データから求められる\(\overline{y}, \overline{x_1}, \overline{x_2} \) の値を方程式\(n(\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}})^2 = 0\) に代入し，\(\hat{\alpha}\)の値を求めます。
\begin{align}
n(\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}})^2 &= 0 \\
\overline{y} - \hat{\alpha} - \hat{\beta_{1}} \overline{x_{1}} - \hat{\beta_{2}} \overline{x_{2}} &= 0 \\
41600 - \hat{\alpha} - (-2519.4)\times9.6 -2961.17\times0.5 &= 0 \\
\hat{\alpha}&=41600 - (-2519.4)\times9.6 -2961.17\times0.5  \\
\hat{\alpha}&=41600 + 2519.4\times9.6 -2961.17\times0.5  \\
\hat{\alpha}&\fallingdotseq64305.83  \\
\end{align}

したがって，駅からの所要時間と駐車場の有無から家賃を推定する重回帰モデルは以下のようになります。
\begin{align}
\hat{y} = -2519.42x_{1} + 2961.17x_{2} + 64305.83
\end{align}

決定係数によるモデルの検討

重回帰分析においても，そのモデルがどの程度データに当てはまっているのか，そしてモデルの偏回帰係数に意味があるのかを検証する必要があります。加えて，モデル全体に統計的な意味があるのかについても検討する必要があります。
先ほど求めたモデル式を例に，決定係数を求め，モデル全体に意味があるかどうかを検討するための分散分析とそれぞれの偏回帰係数に意味があるかを検討するための\(t\) 検定を行います。

○決定係数の算出，検討
単回帰モデルの決定係数\(R^2\)は，以下の式によって求められます。式に出てくる\(\hat{y_i}\)は，「導出したモデル式の\(x\)に\(i\)番目のデータ\(x_i\)を代入したときの推定値」です。
\begin{align}
R^2 &= \frac{回帰変動 (回帰モデルで説明できた変動)}{全変動 (データの変動)} \\
&= 1 - \frac{\sum(y_{i} -\hat{y_{i}})^2 }{\sum(y_{i} -\overline{y_{i}})^2 }\\
&= 1 -\frac{\sum(y_{i} -\hat{\alpha} -\hat{\beta_{1}}x_{1i} -\hat{\beta_{2}}x_{2i})^2} {\sum(y_{i} -\overline{y})^2 } \\
&= 1 -\frac{(y_{1} -\hat{\alpha} -\hat{\beta_{1}}x_{11} -\hat{\beta_{2}}x_{21})^2 + (y_{2} -\hat{\alpha} -\hat{\beta_{1}}x_{12} -\hat{\beta_{2}}x_{22})^2 + ...  (y_{10} -\hat{\alpha} -\hat{\beta_{1}}x_{110} -\hat{\beta_{2}}x_{210})^2} 
{(y_{1} -\overline{y})^2 + (y_{2} -\overline{y})^2 + ... + (y_{10} -\overline{y})^2} \\
&= 0.9176 
\end{align}
この式は単回帰分析における決定係数の式を拡張 (説明変数の数\(k\)のぶんだけ項\(\hat{\beta_k}x_{ki}\) の数を増やす) した式です。

しかし，この式には説明変数の数が増えるほど\(\hat{y}\) が大きくなり，モデルとデータの残差が小さく，決定係数が大きくなりやすい という欠点があります。
そのため，説明変数の数\(k\)に応じて分母 (データ全体の変動) の自由度\(n-1\) と分子 (残差の二乗の総和) の自由度 \(n-k-1\) の比で修正した自由度修正済み決定係数\(\left(R_{adj}^2 \right) \)を用いる必要があります。
修正済み決定係数は，以下の式で表されます。
\begin{align}
R_{adj}^2 &=  1 - \frac{\sum(y_{i} -\hat{y})^2 / (n-k-1) }{\sum(y_{i} -\overline{y})^2 / (n-1) }\\
&= 1 - \frac{n-1}{n-k-1}(1-R^2) \\
& = 0.894 
\end{align}

Rでの分析結果を見ると，今回のデータにおける自由度修正済み決定係数\(R_{adj}^2\) は0.894だとわかります。補正をかけていない決定係数\(R^2 = 0.9176\) と比較してやや低い値ですが，それでもデータに対して十分な説明力を持つモデルだといえます。

○分散分析
重回帰分析では，モデル全体に統計的に意味があるかどうかを確認するために分散分析を行います。単回帰分析ではモデルの傾きが0でないか検討するために\(t\)検定を行いましたが，重回帰分析では回帰係数が0でないかどうかをまとめて検討します。
まず以下のように仮説を立てます。なお，有意水準は\(\alpha = 0.05\)とします。
帰無仮説\(H_{0}\): 全ての偏回帰係数が0である
対立仮説\(H_{1}\): 少なくとも1つの偏回帰係数は0でない

次に，以下の式を用いて検定統計量\(F\)を求めます。
\begin{align}
F &=   \frac{回帰変動/自由度}{残差平方和/自由度}\\
&=   \frac{\sum(\hat{y_{i}} -\overline{y})^2/ (k) }{\sum(y_{i} -\hat{y})^2 / (n-k-1) } \\
&=   \frac{\sum(\hat{y_{i}} -\overline{y})^2/ 2 }{\sum(y_{i} -\hat{y})^2 / (n-k-1) } \\
& =  \frac{(\hat{y_{1}} -\overline{y})^2 + (\hat{y_{2}} -\overline{y})^2 + ... + (\hat{y_{10}} -\overline{y})^2/ 2 }{\sum(y_{1} -\hat{y})^2 + (y_{2} -\hat{y})^2 + ... + (y_{10} -\hat{y})^2/ (10-2-1) } \\
&= 38.95
\end{align}

このモデルから計算された\(F\)値は\(F(2, 7) = 38.95\)でした。\(F\)分布表で\(F(2, 7)_{0.05}\)の数値を見ると19.37であるため，このモデルから計算された\(F\)値は棄却域の内側にあることがわかりました。よって，帰無仮説が棄却され，偏回帰係数のうち少なくとも1つは0ではない，つまりモデル全体は統計的に意味があることが示されました。

サイトポリシー／プライバシーポリシー