統計学実践ワークブック 第16章 重回帰分析 pp.135-136

tjwkbk_ch16_pp135-136.R

# ------------------------------
# 統計学実践ワークブック
# 第16章 重回帰分析 pp.135-136
# ------------------------------

rm(list=ls())
# NAのある行は除外する
# library(tidyverse)
# air <-  drop_na(airquality,everything())
air <- na.omit(airquality)

air |> dim()
air2 |> dim()

object.size(air)
object.size(air2)

identical(air, air2)
all.equal(air,air2)

str(air)
str(air2)

air  |> rownames()
air2 |> rownames()

attr(air2,"na.action")
attr(air2,"na.action") <- NULL

rownames(air2)
rownames(air2) <- NULL

object.size(air)
object.size(air2)

identical(air, air2)
all.equal(air,air2)


# データの準備
air <- na.omit(airquality)

# 構造確認
str(air)

# 重回帰分析の実行
model <- lm(Ozone ~., air)
summary(model)



# 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} 
XtXinv <- solve(t(X) %*% X)
# 残差の標準誤差(Residual standard error)
s <-  sqrt(mean( (y - X %*% XtXinv %*% t(X) %*% y )^2 ) *(n/(n-k-1)))


# 残差の標準誤差(Residual standard error)
s

# 回帰係数(Estimate)たち
XtXinv %*% t(X) %*% y |> t()

# 標準誤差(Std. Error)たち
sqrt( s^2 * diag(XtXinv) )

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

# Dayのt統計量=1.192, 自由度=105に対するp-値(Pr(>|t|) of Day)
pt(q = 1.192, df = 105, lower.tail = F) * 2

# 各種t値
qt(c(0.025, 0.05, 0.1, 0.2359484/2, 0.2), df = 105, lower.tail = F)