重回帰分析は重要と聞くので、念入りに勉強した記録です。
問16.2
問題[1]の要約
モデル2: Ozone = 
+ 
Solar.R + 
Wind + 
Temp + 
Month + 
Day + 
に対して重回帰分析を実施し、変数Dayに対する回帰係数の
を求め、有意性を論じよ。
解答
実際に重回帰分析を実行すると、
# 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列目は
であるので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] )
を求めておく。逆行列は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どおりだとうまくいかず(何故？)、こちらを見て
を用いてみたところ標準誤差たちが得られた。
> 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
のものを引き当てるために必要となる購入回数の期待値は
である(幾何分布)
が幾何分布
に従うとする。
とする。
回連続で外れを引いて
回目で初成功するような確率である
である
となる
を足して
となる。
