第8章 単回帰分析(3)
前章では,回帰係数の推定量が \[ \begin{align} \hat{\beta}_1 \sim N(\beta_1,\sigma^2/S_{xx}) \end{align} \] という標本分布をもつことを見た.しかし,この分散には未知の母数\(\sigma^2\)が含まれているため,このままでは検定に使えない.本章では\(\sigma^2\)をデータから推定し,そのうえで回帰係数の検定を行う.
本章で用いるデータは前章と同じものである.
library(tidyverse)
set.seed(8931)
data <- tibble(
x = runif(100, 0, 50),
y = rnorm(100, x + 10, 10)
) 8.1 誤差分散の推定
誤差分散\(\sigma^2\)は,残差平方和\(S_e\)を利用して \[ \hat{\sigma}^2=\frac{S_e}{n-2} \] と推定する.自由度が\(n-2\)となるのは,\(\beta_0, \beta_1\)の2つを推定したぶんだけ自由に動ける残差の数が減るためである(「推測統計の復習」で標本分散の自由度が\(n-1\)であったのと同じ理屈である).
この推定量を用いると, \[ \begin{align} t_{\hat{\beta}_1}=\frac{\hat{\beta}_1-\beta_1}{\sqrt{\hat{\sigma}^2/S_{xx}}} \sim t(n-2) \end{align} \] が成り立つ.つまり,\(\hat{\beta}_1\)を標準化する際に\(\hat{\beta}_1\)の分散をデータから推定して得られた\(t\)統計量は,自由度\(n-2\)の\(t\)分布に従う.母分散が未知のときに正規分布ではなく\(t\)分布を用いる点は,「推測統計の復習」で扱った母平均の検定と同じ構造である.
Rで実際に計算してみよう.
lm(y ~ x, data) -> lm_res
# 残差平方和から誤差分散を推定
Se <- sum(residuals(lm_res)^2)
n <- nrow(data)
sigma2_hat <- Se / (n - 2)
sigma2_hat## [1] 109.822
## [1] 0.06951742
この標準誤差の値は,後で見るsummary()の出力のStd. Errorと一致する.
8.2 検定統計量
帰無仮説,対立仮説をそれぞれ \[ \begin{align} H_0: \beta_1=0,\qquad H_1:\beta_1\neq 0 \end{align} \] とする.このとき,帰無仮説が正しいという仮定の下で検定統計量は自由度\(n-2\)の\(t\)分布に従う.つまり, \[ \begin{align} t^*_{\hat{\beta}_1}=\frac{\hat{\beta}_1}{\sqrt{\hat{\sigma}^2/S_{xx}}} \sim t(n-2). \end{align} \] \(|t^*_{\hat{\beta}_1}|\geq t_{\alpha/2}(n-2)\)のとき,両側\(100\alpha\%\)の有意確率で帰無仮説を棄却し対立仮説を採用する.\(|t^*_{\hat{\beta}_1}|< t_{\alpha/2}(n-2)\)の場合帰無仮説は棄却できない.
8.3 検定の実際
実際にR上でどのように検定結果が報告されているか確認しよう.
##
## Call:
## lm(formula = y ~ x, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.6956 -7.0253 -0.8862 6.0086 21.5698
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.86720 2.02106 4.882 4.08e-06 ***
## x 0.96227 0.06952 13.842 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.48 on 98 degrees of freedom
## Multiple R-squared: 0.6616, Adjusted R-squared: 0.6582
## F-statistic: 191.6 on 1 and 98 DF, p-value: < 2.2e-16
データの分析結果から検定の仮定を再現する.結果出力におけるxのEstimateの値が\(\hat{\beta}_1\)に対応し,Std. Errorの値が\(\sqrt{\hat{\sigma}^2/S_{xx}}\)に対応する.ゆえに,\(t\)値は \[
\begin{align*}
t^*_{\hat{\beta}_1}=\frac{0.96227}{0.06952} = 13.842.
\end{align*}
\] となる.データサイズは\(n=100\)なので,自由度98の\(t\)分布より, \[\begin{align*}
P(|T|>13.842) = 2* 10^{-16}
\end{align*}\] が\(p\)値として得られる.
先ほど手で計算したStd. Errorの値と,summary()の出力が一致していることを確認せよ.組み込み関数の出力する数値は,すべて自分の手で再現できる.