重回帰分析は重要と聞くので、念入りに勉強した記録です。

問16.2

問題[1]の要約

モデル2: Ozone = $\beta_0$ + $\beta_1$ Solar.R + $\beta_2$ Wind + $\beta_3$ Temp + $\beta_4$ Month + $\beta_5$ Day + $\varepsilon$ に対して重回帰分析を実施し、変数Dayに対する回帰係数の

最小二乗推定量
t統計量

を求め、有意性を論じよ。

解答

実際に重回帰分析を実行すると、

# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
#Day           0.27388    0.22967   1.192  0.23576

を得る。

従って、

変数Dayに対する回帰係数はEstimate列から0.27388
t統計量はt value列から1.192
- 問題文ではEstimateとStd.Errorしか提示されていないが、t value = Estimate÷Std.Error より 0.27388÷0.22967=1.192494を得る
p-値はPr(>|t|)列から0.23576である
- 問題文ではp-値は提示されていないので付表2. t分布のパーセント点から判断してみると、
- 自由度120で両側10%点が1.658、自由度120で両側20%点が1.289なのでp-値は20%よりも大きいことが推測できる
- このことから、両側10%でも20%でも有意ではないと言える

問題文のデータセットを使ってRで重回帰分析してみる

airqualityのデータはRに元から入っているデータセットなので重回帰分析してみる。

データ取得

airqualityをsummaryで見るとNAを含む行がいるのでomitする。

summary(airquality)

# NAのある行は除外する
air <- na.omit(airquality)

構造確認

ワークブックに従い、Ozoneを目的変数・それ以外を説明変数とする。

str(air)
> str(air)
'data.frame':   111 obs. of  6 variables:
 $ Ozone  : int  41 36 12 18 23 19 8 16 11 14 ...
 $ Solar.R: int  190 118 149 313 299 99 19 256 290 274 ...
 $ Wind   : num  7.4 8 12.6 11.5 8.6 13.8 20.1 9.7 9.2 10.9 ...
 $ Temp   : int  67 72 74 62 65 59 61 69 66 68 ...
 $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
 $ Day    : int  1 2 3 4 7 8 9 12 13 14 ...
 - attr(*, "na.action")= 'omit' Named int [1:42] 5 6 10 11 25 26 27 32 33 34 ...
  ..- attr(*, "names")= chr [1:42] "5" "6" "10" "11" ...

重回帰分析の実行

lm関数にて重回帰分析を実行する。

model <- lm(Ozone ~., air)
summary(model)
> summary(model)

Call:
lm(formula = Ozone ~ ., data = air)

Residuals:
    Min      1Q  Median      3Q     Max 
-37.014 -12.284  -3.302   8.454  95.348 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -64.11632   23.48249  -2.730  0.00742 ** 
Solar.R       0.05027    0.02342   2.147  0.03411 *  
Wind         -3.31844    0.64451  -5.149 1.23e-06 ***
Temp          1.89579    0.27389   6.922 3.66e-10 ***
Month        -3.03996    1.51346  -2.009  0.04714 *  
Day           0.27388    0.22967   1.192  0.23576    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 20.86 on 105 degrees of freedom
Multiple R-squared:  0.6249,    Adjusted R-squared:  0.6071 
F-statistic: 34.99 on 5 and 105 DF,  p-value: < 2.2e-16

ワークブックと見比べると、p.136と同じ結果が出ていることがわかる。

重回帰分析をスクラッチ実装する

せっかくワークブックで線形代数をゴリゴリ用いて重回帰分析を勉強したので、lm関数を用いずにスクラッチ実装(と言えるのか・・・？とりあえず標準的な関数だけで実装)してみる。

変数の準備

データ行列などを準備する。注意点として、

と、Xの1列目は $\vec{1}$ であるのでcbind関数で列をくっつけることを忘れずに。

# X: データ行列, n: 行数, k: パラメータ数
X_raw <- air[, -1]
n <- nrow(X_raw)
k <- ncol(X_raw)
X <- cbind(1, as.matrix(X_raw) )

# y: 目的変数
y <- as.numeric( air[, 1] )

例の‪逆行列

$(X^{T}X)^{-1}$ を求めておく。逆行列はsolve関数で求める。

# (X^{T}X)^{-1} 
XtXinv <- solve(t(X) %*% X)

各係数を求める

残差の標準誤差

残差の標準誤差は、

を用いて計算する。なお、mean関数だとn個足してnで割ってしまっているので、nを掛けて自由度n-k-1で割っている。

## 残差の標準誤差(Residual standard error)
s <-  sqrt( mean( (y - X %*% XtXinv %*% t(X) %*% y )^2 )*(n/(n-k-1)) )

値を確認すると、「esidual standard error: 20.86」と同じ値が得られていることがわかる。

> print(s)
[1] 20.85846

> sigma(model)
[1] 20.85846

回帰係数(Estimate)たち

を素直に計算すれば良い。最後の転置(t())は見栄えのために使用している。

> print( XtXinv %*% t(X) %*% y |> t() )
                  Solar.R      Wind     Temp     Month       Day
[1,] -64.11632 0.05027432 -3.318444 1.895786 -3.039957 0.2738775

これはlm関数で求めた値と等しい。

> coef(model)
 (Intercept)      Solar.R         Wind         Temp        Month          Day 
-64.11632110   0.05027432  -3.31844386   1.89578642  -3.03995664   0.27387752

標準誤差(Std. Error)たち

ワークブックp.130どおりだとうまくいかず(何故？)、こちらを見て $X^{T}X$ を用いてみたところ標準誤差たちが得られた。

> print( sqrt( s^2 * diag(XtXinv) ) )
              Solar.R       Wind       Temp      Month        Day 
23.4824887  0.0234186  0.6445095  0.2738869  1.5134568  0.2296708

> coef(summary(model))[,2]
(Intercept)     Solar.R        Wind        Temp       Month         Day 
 23.4824887   0.0234186   0.6445095   0.2738869   1.5134568   0.2296708

axjack's blog

### axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz ###

統計学実践ワークブック第16章重回帰分析問16.2を通じて・重回帰分析のスクラッチ実装