Skip to Main Content

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

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

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

(データ概要)
以下に示すような駅から物件まで徒歩でかかる時間と家賃のデータがあります。このデータを用いて単回帰モデルを作成し,駅から徒歩20分の物件は家賃がどれくらいになるかを予測します。

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

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

\begin{align}
y = \alpha + \beta x \\
y: 目的変数 (家賃の推定値) \\
x: 説明変数 (所要時間)
\end{align}

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

今回は一旦シグマの中身を放置して次のステップに進みます。ここで以下のようにシグマを展開してしまうと,その後非常に煩雑な計算を強いられます。数学アレルギーが発症しそうですね。
\begin{align}
\sum_{i=1}^n{\left\{y_{i}-\left(\hat{\alpha} +\hat{\beta}\;x_{i}\right)\right\}^2} = \left\{y_{1}-\left(\hat{\alpha} +\hat{\beta}\;x_{1}\right)\right\}^2 + \left\{y_{2}-\left(\hat{\alpha} +\hat{\beta}\;x_{2}\right)\right\}^2 + ... +\left\{y_{10}-\left(\hat{\alpha} +\hat{\beta}\;x_{10}\right)\right\}^2 \\
\end{align}

② ①で求めた式を\(\hat{\alpha}\), \(\hat{\beta}\) でそれぞれ微分する 
線形モデリングでは,観測値と理論値の差が最も小さくなるようなモデルを選択します。そのため,先程求めた残差平方和が最も小さくなるような\(\hat{\alpha}\), \(\hat{\beta}\) を求める必要があります。
この際,式に含まれる変数のうちどちらか一方のみを変数とみなして微分を行う,偏微分という方法を用います。

まず,\(\sum_{i=1}^n{\hat{u_{i}}^2}  = \sum_{i=1}^n{\left\{y_{i}-\left(\hat{\alpha} +\hat{\beta}\;x_{i}\right)\right\}^2}\ \) という式のうち,\(\hat{\alpha}\) のみを変数として以下の式の通り微分します (※1)
\begin{align}
\frac{\partial \sum{\hat{u_{i}}^2}}{\partial \hat{\alpha}} &= \sum\left[2(y_{i} - \hat{\alpha} - \hat{\beta}x_{i})'\; (y_{i} - \hat{\alpha} - \hat{\beta}x_{i})\right] \\
&=-2\sum(y_{i} - \hat{\alpha} - \hat{\beta}x_{i})...\;(1)
\end{align}

次に,\(\sum{\hat{u_{i}}^2}\) について,\(\hat{\beta}\) のみを変数として以下の式の通り微分します。
\begin{align}
\frac{\partial \sum{\hat{u_{i}}^2}}{\partial \hat{\beta}} &= \sum\left[2(y_{i} - \hat{\alpha} - \hat{\beta}x_{i})'\; (y_{i} - \hat{\alpha} - \hat{\beta}x_{i})\right] \\
&=-2x_{i}\sum(y_{i} - \hat{\alpha} - \hat{\beta}x_{i})...\;(2)
\end{align}

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

\begin{align}
-2x_{i}\sum(y_{i} - \hat{\alpha} - \hat{\beta}x_{i}) &= 0 \\
\sum(x_{i}y_{i} - x_{i}\hat{\alpha} - \hat{\beta}x_{i}^2) &= 0 \\
\sum x_{i}y_{i} - \hat{\alpha}\sum x_{i} - \hat{\beta}\sum x_{i}^2 &= 0 \\
\hat{\alpha}\sum x_{i} + \hat{\beta}\sum x_{i}^2 &= \sum x_{i}y_{i} ...(2)' 
\end{align}

したがって,正規方程式は以下のように示されます。
\(\begin{cases}n\hat{\alpha} + \hat{\beta}\sum x_{i} &=\sum y_{i}\\\hat{\alpha}\sum x_{i} + \hat{\beta}\sum x_{i}^2 &= \sum x_{i}y_{i}\end{cases}\)

④ \(\hat{\alpha}\), \(\hat{\beta}\) の値を求め,回帰モデル式を求める
式(1)',(2)'を連立させ,データからわかる値 \(\left(n,\;\sum x_{i},\; \sum y_{i},\; \sum x_{i}^2,\; \sum x_{i}y_{i} \right)\) を代入すると以下の連立方程式が得られます。
\(\begin{cases}10\hat{\alpha} + 96\hat{\beta} &=416000\\96\hat{\alpha} + 1014\hat{\beta} &= 3746000\end{cases}\)

上の連立方程式を解くと,\(\hat{\alpha}\fallingdotseq67324.7,\; \hat{\beta}\fallingdotseq-2679.7\) という値が得られます。よって,家賃と駅からの所要時間から推定された単回帰モデルは\(\hat{y} = 67324.7 - 2679.7x\) です。

(※1) 式が煩雑になるのを防ぐために\(\sum{\hat{u_{i}}^2}\) と表記していますが,意味としては\(\sum_{i=1}^n{\hat{u_{i}}^2} \)と全く同じだと考えてください。
(※2) 詳しく記述すると,まずは式(1),(2) を二次関数の接線の傾きを表す式 (導関数) として考えます。そして,接線の傾きが0になる値をそれぞれ求めることで誤差が最小となるような\(\hat{\alpha}\), \(\hat{\beta}\)を求めます。

基本・単回帰分析 (決定係数によるモデルの検討)

回帰分析 (線形モデリング) はモデル式だけを求めて終わり,ではなく,そのモデルがどの程度データに当てはまっているのか,そしてモデルの傾きに意味があるのかを検証する必要があります。
先ほど求めたモデル式を例に,モデルの当てはまり度合いを示す決定係数の計算とモデルの傾き\(\hat{\beta}\) についての\(t\)検定を行ってみましょう。

○決定係数の算出
単回帰モデルの決定係数\(R^2\)は,以下の式によって求められます。式に出てくる\(\hat{y_i}\)は,「導出したモデル式の\(x\)に\(i\)番目のデータ\(x_i\)を代入したときの値」です。
\begin{align}
R^2 &= \frac{回帰変動 (回帰モデルで説明できた部分)}{全変動 (データの変動)} \\
&= \frac{\sum(\hat{y_i} -\overline{y})^2 }{\sum(y_{i} -\overline{y})^2 }\\
&\hat{y_i} = 67324.7 -2679.7x_{i} \\
& (ただし,0 \leq R^2 \leq 1) 
\end{align}
\(R^2\)は観測値の変動に占める「回帰モデルで説明できる変動」の比を示しています。そのため,1に近いほど観測値の変動をモデルでよく説明できている (モデルとデータの当てはまりが良い) ことを,0に近いほど観測値の変動をモデルで説明できていない (モデルとデータの当てはまりが悪い) ことを示します。

今回のデータにおける決定係数は以下のように計算されます。決定係数から,データとモデルの適合度がかなり高いことがわかります。
\begin{align}
R^2 &= \frac{\sum(\hat{y_i} -\overline{y})^2 }{\sum(y_{i} -\overline{y})^2 }\\
&= \frac{(53926.2 - 41600)^2 +  (51246.5 - 41600)^2 + ... + (27129.2 - 41600)^2}{(13400 - 41600)^2 +  (15400 - 41600)^2 + ... + (-12600 - 41600)^2}\\
&= \frac{663505190.88}{744400000} \\
&\fallingdotseq 0.891
\end{align}

ひとつ留意しておくべきことは,決定係数\(R^2\)の大きさは必ずしもモデルの良し悪しを決定しないということです。分野や調査方法によっても当てはまりがよいとする目安は異なりますし,あまりにも1に近いとオーバーフィッティング (モデルがデータに対して過剰に当てはまっている状態。モデルを作ったデータだけにはよく当てはまるが,それ以外のデータはうまく説明できないモデルである可能性がある) を起こしている可能性もあります。

○モデルの傾きについての\(t\)検定
回帰モデルの傾きについての\(t\)検定は,まず以下のような仮説を立てます。真の傾きが0かそうでないかに着目しているので,両側検定を採用します。なお,有意水準は\(\alpha = 0.05\)とします。
帰無仮説\(H_{0}\): 真のモデルの傾きは\(\beta = 0\)である
対立仮説\(H_{1}\): 真のモデルの傾きは\(\beta \neq 0\)である

次に,以下の式を用いてデータから推定された回帰モデルの傾き\(\hat{\beta}\)の検定統計量\(t_{\hat{\beta}}\) の値を求めます。
\begin{align}
t_{\hat{\beta}} &= \frac{\hat{\beta} - \beta}{SE_{\hat{\beta}}} \\
&= \frac{\hat{\beta} }{\sqrt{\frac{\sum{\hat{u_{i}}^2} / n - k - 1} {\sum{(x_{i} - \overline{x})^2}}}} \\
&= \frac{-2679.7}{\sqrt{\frac{80917750.88/ 10 - 1 - 1} {92.4}}} \\
&\fallingdotseq -8.10 \\
&\sum{\hat{u_{i}}^2}: 残差平方和 (上のボックスの式にデータの値を代入することで求められる)\\
&SE_{\hat{\beta}}: 回帰モデルの傾き{\hat{\beta}}の標準誤差\\
&n: サンプルサイズ \\
&k: 説明変数の数 \\
& \overline{x}: xの平均値
\end{align}

回帰モデルの傾きの自由度は\(n - k - 1 = 8\) であるため,棄却域は \(2.306 < x <\infty \) または  \(- \infty < x < -2.306 \) となります。
今回得られた\(t_{\hat{\beta}} = -8.10\) と棄却域を比較すると,\(t_{\hat{\beta}} \) は棄却域に含まれています。よって,帰無仮説を棄却し,「真のモデルの傾きは\(\beta \neq 0\)である」と結論づけます。