第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
# 標準誤差
Sxx <- sum((data$x - mean(data$x))^2)
sqrt(sigma2_hat / Sxx)
## [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上でどのように検定結果が報告されているか確認しよう.

lm(y ~ x, data = data) |> summary()
## 
## 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()の出力が一致していることを確認せよ.組み込み関数の出力する数値は,すべて自分の手で再現できる.

8.4 本日の課題

仮想データを自ら作成して単回帰分析を実行し,summary()の出力について,Std. Errorと\(t\)値を自分で計算して再現せよ.そのうえで,検定結果を文章で報告せよ.Rスクリプトを実行した結果をWordにコンパイルしたファイルをLUNA提出