読者です 読者をやめる 読者になる 読者になる

野次馬エンジニア道

野次馬な気持ちでプログラミングをあれこれと綴ります

回帰分析 - 単回帰における回帰係数の推定

回帰分析

今日はモデリングの基本として単回帰から見ていく。Rを使った出力まで。

回帰モデル

母回帰方程式(population regression equation)と呼ばれる1次関数。{ \beta
_0, \beta_1} を母回帰係数、{ y_i} を従属変数 { x_i} を独立変数 { \varepsilon_i}を誤差項と呼ぶ。 他にも呼び方は色々ある。

{ \displaystyle
 y_i = \beta_0 + \beta_1x_i + \varepsilon_i \quad (i=1,\cdots,n)
}

それぞれの変数は下記の条件を満たす

  • 誤差の期待値は0。{ E[\varepsilon_i] = 0 }
  • 誤差の分散は一定。 { V[\varepsilon_i] = \sigma^2 }
  • 誤差は互いに独立。{ i\neq j } のとき、{ Cov[\varepsilon_i,\varepsilon_j] = E[\varepsilon_i\varepsilon_j] = 0 }

解釈としては、{y_i} が誤差を含む実現値で、その期待値は、 {\beta_0 + \beta_1x_i} である。

最小二乗法

当てはめた直線との誤差を二乗した

{ \displaystyle
S_e = \sum_{i=1}^ne_i^2 = \sum_{i=1}^n\left\{ y_i - \beta_0 + \beta_1x_i \right\}^2
}

指標を最小にする。その推定量は、 {\beta_0, \beta_1}偏微分して0と置いた式を解くと

{ \displaystyle
 \hat{\beta_1} = \frac{\sum_{i=1}^n(x_i - \overline{x})(y_i - \overline{y})}{ \sum_{i=1}^n (x_i - \overline{x} )^2 } = \frac{S_{xy}}{S_{xx}}
}

{ \displaystyle
 \hat{\beta_0} = \overline{y} - \hat{\beta_2}\overline{x}
}

と求まる。

予測値と残差

  • 予測値 { \hat{y} = \beta_1 + \beta_2x_i } の平均は観測値の平均と等しい
  • 残差 (residual)は { e_i = y_i - \hat{y_i} } の平均は0。{ \sum_{i=1}^n e_i = 0 }
  • 予測値の偏差と残差の積和は0。{ \sum_{i=1}^{n}(\hat{y_i}-\overline{y})(e_i-\overline{e}) = \sum_{i=1}^{n}(\hat{y_i}-\overline{y})e_i = 0 }

重相関係数 (multiple correlation coefficient)

観測値と予測値の相関係数は、

{ \displaystyle
r_{y\hat{y}} = \frac{S_{y\hat{y}}}{ \sqrt{S_{yy}S_{\hat{y}\hat{y}}}}
}

{ \displaystyle
S_{y\hat{y}} = \sum_{i=1}^{n}(y_i-\overline{y})(\hat{y_i}-\overline{y}) = \sum_{i=1}^{n}(\hat{y_i} + e_i-\overline{y})(\hat{y_i}-\overline{y})
}

{ \displaystyle
= \sum_{i=1}^{n}(\hat{y_i}-\overline{y})^2 + \sum_{i=1}^{n}(\hat{y_i}-\overline{y})e_i = S_{\hat{y}\hat{y}}
}

すると、

{ \displaystyle
r_{y\hat{y}} = \frac{ \sqrt{S_{\hat{y}\hat{y}}} }{ \sqrt{S_{yy}} } = \frac{ s_{\hat{y}} }{ s_{y} } = R
}

となる。{ y } {\hat{y}}の相関係数を重相関係数と呼ぶ。

(予測値の標準偏差)/(観測値の標準偏差)と等しくなる。

決定係数 (coefficient of determination)

当てはまりの良さ(どの程度 {x}{y} を説明できているか)を、説明できている変動とそうでない変動とに分けて評価する。{ y } のばらつきの総和(変動)を下記のようにすると

{ \displaystyle
 S_{yy} = \sum_{i=1}^n (y_i-\overline{y})^2 = \sum_{i=1}^n (\hat{y_i} + e_i -\overline{y})^2
}

{ \displaystyle
= \sum_{i=1}^n (\hat{y_i}-\overline{y})^2 + \sum_{i=1}^n e_i^2 + 2 \sum_{i=1}^{n}(\hat{y_i} - \overline{y})e_i
}

{ \displaystyle
= \sum_{i=1}^n (\hat{y_i}-\overline{y})^2 + \sum_{i=1}^n e_i^2
}

{ \displaystyle
= S_R + S_e
}

予測値と残差、基準化したいので、左辺で割ったものを決定数

{ \displaystyle
R^2 = 1 - \frac{\sum_{i=1}^n e_i^2}{\sum_{i=1}^n (y_i-\overline{y})^2 } = \frac{\sum_{i=1}^n (\hat{y_i}-\overline{y})^2}{\sum_{i=1}^n (y_i-\overline{y})^2 }
}

と定義する。また

{ \displaystyle
R^2 = 1 - \frac{S_e}{S_{yy}} = \frac{S_R}{S_{yy}}
}

説明ができれば1に近くなり、全く説明できなければ0となる。

{ \displaystyle
= \frac{\sum_{i=1}^n \frac{(\hat{y_i}-\overline{y})^2}{n-1}}{ \sum_{i=1}^n \frac{(y_i-\overline{y})^2}{n-1} } = \left(\frac{ s_{\hat{y}} }{ s_{y} } \right)^2
}

決定係数の平方根は重相関係数となる。

自由度調整済み決定係数

データの大きさ { n} が小さいときに、決定係数が大きくなる傾向がある。モデルの自由度を考慮して決定係数を調整する。独立変数の個数が { p } 個のとき、データの変動の自由度が {n-1} 残差の自由度は、 { n - 1 - p } となる。

{ \displaystyle
R^{*2} = 1 - \frac{ \frac{S_e}{n-1-p}} { \frac{S_{yy}}{n-1} }
}

推定値の標準誤差 (standard error of estimates)

モデルの誤差 {\varepsilon^2} を残差の平方和 {S_e} を用いて推定できる。

{ \displaystyle
s^2 =
\frac{ S_e }{n-2} =\sum_{i=1}^{n}\frac{e_i^2}{ n - 2} = \sum_{i=1}^{n}\frac{ ( y_i - \hat{y_i})^2}{ n - 2 }
}

{s} の値が小さければ当然回帰式の当てはまりはよい。

残差の平方和の計算

最小二乗の式に推定された係数を入れると

{
 S_e = S_{yy} - \hat{\beta_1} S_{xy}
}

となるので、{S_{yy} = S_{R} + S_{e} } の関係より

{
 S_R = \hat{\beta_1} S_{xy}
}

であることがわかる。{s^2} の算出に用いる残差の平方和は下記で計算可能。

{ \displaystyle
\sum_{i=1}^n (y_i-\hat{y_i})^2 = S_{e} = S_{yy} - S_{R} = S_{yy} - \hat{\beta_1} S_{xy}
}

回帰係数の推定

不偏推定量

最小二乗で求めた推定値は不偏推定量となる。{ \sum_{i=1}^n(x_i - \overline{x}) = 0 } となるので

{ \displaystyle
 \hat{\beta_1} = \frac{\sum_{i=1}^n(x_i - \overline{x})(y_i - \overline{y})}{ \sum_{i=1}^n (x_i - \overline{x} )^2 } = \frac{\sum_{i=1}^n(x_i - \overline{x}) y_i}{ \sum_{i=1}^n (x_i - \overline{x} )^2 }
}

{ \displaystyle
 = \frac{\sum_{i=1}^n(x_i - \overline{x})(\beta_0 + \beta_1x_i + \varepsilon_i)}{ \sum_{i=1}^n (x_i - \overline{x} )^2 }
}

{ \displaystyle
 = \beta_1\frac{\sum_{i=1}^n(x_i - \overline{x})x_i}{ \sum_{i=1}^n (x_i - \overline{x} )x_i } + \frac{\sum_{i=1}^n(x_i - \overline{x})\varepsilon_i}{ \sum_{i=1}^n (x_i - \overline{x} )^2 }
}

ここで、{ w_i = \frac{ x_i - \overline{x}}{\sum_{i=1}^n (x_i - \overline{x})^2}
} とすると、

{ \displaystyle
 \hat{\beta_1} = \beta_1 + \sum_{i=1}^n w_i\varepsilon_i
}

となる。期待値と分散は { E[\varepsilon_i] = 0, V[\varepsilon_i] = \sigma^2 } なので

{ \displaystyle
 E[\hat{\beta_1}] = \beta_1 + \sum_{i=1}^n w_i E[\varepsilon_i] = \beta_1
}

{ \displaystyle
 V[\hat{\beta_1}] = \sum_{i=1}^n w_i^2 V[\varepsilon_i] = \frac{\sigma^2}{\sum_{i=1}^n (x_i - \overline{x} )^2}
}

{ \beta_0 } の方は

{ \displaystyle
\hat{\beta_0} = \overline{x} - \hat{\beta_1}\overline{x} = (\beta_0 + \beta_1\overline{x} + \epsilon_i) -  \hat{\beta_1}\overline{x} = \beta_0 + \epsilon_i - \overline{x}( \hat{\beta_1} - \beta_1 )
}

{ \displaystyle
 = \beta_0 + \sum_{i=1}^n \left( \frac{1}{n} -  w_i \overline{x}\right)\epsilon_i
}

となる。

{ \displaystyle
 E[\hat{\beta_0}] = \beta_0 + \sum_{i=1}^n \left( \frac{1}{n} -  w_i \overline{x}\right) \epsilon_i
}

{ \displaystyle
 V[\hat{\beta_0}] = V[\beta_0] + \sum_{i=1}^n \left( \frac{1}{n} -  w_i \overline{x}\right)^2 V[\epsilon_i]
}

{ \displaystyle
= \sum_{i=1}^n \left( \frac{1}{n^2} -  2nw_i \overline{x} + w_i^2\overline{x}^2 \right) \sigma^2
}

{ \displaystyle
= \left( \frac{1}{n} + \frac{\overline{x}^2}{\sum_{i=1}^n(x_i - \overline{x})^2} \right) \sigma^2
}

{ \sum x_i^2 = \sum ( x_i - \overline{x} )^2 + n\overline{x}^2 } となるのを用いて

{ \displaystyle
 = \frac{\sigma^2\sum_{i=1}^n x_i^2}{n\sum_{i=1}^n (x_i - \overline{x} )^2}
}

{ \hat{\beta_0}, \hat{\beta_1}}ガウスマルコフの定理(Gauss-Markov's theorem)より最良線形不偏推定量(BLUE,best linear unbiased estimator)になっていることがわかる。

回帰係数の標本分布

回帰係数は、{ N(\beta_2, \frac{\sigma^2}{\sum(x_i-\overline{x})^2})} に従う。{ \sigma } を推定量 { s } を用いて

{ \displaystyle
 s.e. = \sqrt{s^2} = \sqrt{\sum_{i=1}^{n}\frac{e_i^2}{ n - 2}}
}

とすると、回帰係数の標準誤差は

{ \displaystyle
s.e.(\hat{\beta_1}) = \frac{s.e.}{\sqrt{\sum_{i=1}^{n}(x_i - \overline{x})^2}}
}

となる。標準化した

{ \displaystyle
  t = \frac{\hat{\beta_1}-\beta_1}{s.e.(\hat{\beta_1})}
}

は自由度 { n-2}{t} 分布に従う。

回帰分析の結果例

Rによる出力例をみてみる。

> summary(lm(y ~ x))

Call:
lm(formula = y ~ x)

Residuals:
   Min     1Q Median     3Q    Max 
-5.232 -3.132  1.132  2.450  5.026 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  11.9747     2.7856   4.299 0.000865 ***
x             0.3790     0.2774   1.366 0.194974    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.391 on 13 degrees of freedom
Multiple R-squared:  0.1256,    Adjusted R-squared:  0.05833 
F-statistic: 1.867 on 1 and 13 DF,  p-value: 0.195

{ y_i = \beta_0+ \beta_1x_i + \varepsilon_i
} とすると、推定値は、Estimate、Interceptは切片の推定値 {\beta_0} 、xの右が {\beta_1} { 11.9747 + 0.379 x_i + \varepsilon_i }となる。Residual standard error:が標準誤差の推定値。 Multiple R-squaredが決定係数。この例では、0.1256で変動を説明できていないことがわかる。

回帰結果の評価

長くなったが

  • 決定係数と自由度調整済み決定係数
  • 推定された係数および標準誤差

を用いて結果を評価する。長くなったので他の値は次回。