別のブログでアイナちゃんの課題をミサワ君とサエコ先輩が助ける話を書いたのですが、検証用にその物語中のソースコードを公開しておきます。
1. データセットの作成
データセットは課題で出されたモノになりますが、ファイルを用意するのも面倒なので、今回は乱数から機械的に作成します。
obs <- 300
set.seed(1225)
dframe <- data.frame(
x1 = runif(obs, min=0, max=3),
x2 = runif(obs, min=-1, max=4)
)
attach(dframe)
dframe$y <- 1 + 2*x1 + 3*x2 + rnorm(obs, mean=0, sd=2)
detach(dframe)
2. ミサワ君のレクチャー
実は大学院レベルのテキストで無いとやらない気もする行列によるOLSの解法をミサワ君は説明しています。
# 1. データフレームの中身を確認
names(dframe)
# 2. データフレームをattachする
attach(dframe)
# 3. 方程式を立てる
frml <- y ~ x1 + x2
# 4. 説明変数の行列を取得
X <- model.matrix(frml)
# 5. 自由度を計算
df <- 300 - 3
# 6. 係数βを計算
b <- solve(t(X)%*%X)%*%(t(X)%*%y)
# 7. 残差平方和を計算
ssr <- sum((y - X %*% b)^2)
# 8. 分散を計算
s2 <- ssr/df
# 9. 被説明変数yの分散を計算
var_y <- sum( (y - mean(y))^2 )
# 10. 分散共分散行列を計算
vcov <- s2*solve(t(X)%*%X)
# 11. 標準誤差を計算
sqrt(diag(vcov))
# 12. マスク値の行列を作成（切片項ならc(1, 0, 0)、x2ならc(0, 0, 1)）
tm <- matrix(c(0, 1, 0), 1, 3)
# 13. x1の係数のt値を計算
tv <- (tm %*% b) / sqrt(tm %*% vcov %*% t(tm))
# 14. 両側検定でx1のP値を計算
if(0>tv) tp <- 2*pt(tv, df) else tp <- 2*(1 - pt(tv, df))
# 15. x1のt値とP値を表示
sprintf("%f [%0.5f%%]", tv, tp/100)
# 16. 決定係数を計算
1 - ssr/var_y
# 17. 自由度調整済み決定係数を計算
1 - (ssr/df)/(var_y/(obs-1))
# 18. データフレームをdetachする
detach(dframe)
