最小二乗法からp値まで

線形回帰の基本は最小二乗法を実行して決定係数やp値を確認しながらモデルをブラッシュアップしていくことです。

Rにはlm関数などの便利なものがあるので内部ロジックがどうなっているのか、この関数を使っている限りは気になりません。しかし、いざ人に説明しようとすると内部ロジックに精通する必要が出てきます。というわけで、今日は最小二乗法かp値算出までの線形回帰の基本を確認したいと思います。結構、このあたりをきちんと理解していると次のステップに進みやすいでしょう。

■データ
・今回は次のデータを使います。

x <- c(45,33,39,40,42,30,53,45,42,31) #説明変数
y <- c(41,29,36,36,40,26,46,33,35,28) #目的変数
X <- matrix(1,10,2) #行列作成
X[,2] <- x  #1列目 = 切片, 2列目=説明変数

■最小二乗法によるパラメータの推定
・線形回帰の場合には最小二乗法=最尤法で次のシンプルな式を解くだけでパラメータを求めることができます。

B <- t(X) %*% X; c <- t(X) %*% y; beta <- solve(B, c)
a <- beta[1,]; b <- beta[2,]
yhat <- a + b*x

・カンタンカンタン！

■決定係数
・修正済み決定係数: γ^2

#試行錯誤そのままなのでコードが冗長
n <- length(x)
sy2 <- sum((y-mean(y))^2)/n
sy <- sqrt(sy2)
syx2 <- sum((y-a-b*x)^2)/n
syx <- sqrt(syx2)

sig2 <- n*syx2/(n-2)
sig <- sqrt(sig2)
sigy2 <- n*sy2/(n-1)

rbar2 <- 1 - (sig2/sigy2)

■t値
・βの仮説検定を行うにはt値を算出しておく必要があります。

sx2 <- sum((x-mean(x))^2)/n
sx <- sqrt(sx2)

siga2 <- sig2 * sum(x^2) /(n^2*sx2)
siga <- sqrt(siga2)

sigb2 <- sig2/(n*sx2)
sigb <- sqrt(sigb2)

ta <- a / siga
tb <- b / sigb

■p値
・最後に上記で求めた値からt分布を使ってp値を求めます。

pa <- (1 - pt(ta, n-2))*2
pb <- (1 - pt(tb, n-2))*2

■参考文献

基本統計学

• 作者: 宮川 公男
• 出版社/メーカー: 有斐閣
• 発売日: 1999/03
• メディア: 単行本

・宮川先生のこの本では線形回帰がとても分かりやすく書かれています。上記のコーディングもこの本を参考にしました。おススメの一冊です。初学者には下記の本よりも分かりやすいことでしょう。（入門書としては下記が紹介されることが多いですが、こちらの方が記載が優しいので下記本に挫折した人でもこちらでリカバリーできる可能性があります。）

統計学入門 (基礎統計学)

• 作者:
• 出版社/メーカー: 東京大学出版会
• 発売日: 1991/07/09
• メディア: 単行本

・この本も分かりやすいですが、入門書を数冊読んだ後の方がよいですね。宮川先生の本を先に一読されることをおススメします。

■全コード

####DATA
x <- c(45,33,39,40,42,30,53,45,42,31) #説明変数
y <- c(41,29,36,36,40,26,46,33,35,28) #目的変数
X <- matrix(1,10,2) #行列作成
X[,2] <- x  #1列目 = 切片, 2列目=説明変数

####Estimation
B <- t(X) %*% X; c <- t(X) %*% y; beta <- solve(B, c)
a <- beta[1,]; b <- beta[2,]
yhat <- a + b*x

####R-squared
n <- length(x)
sy2 <- sum((y-mean(y))^2)/n
sy <- sqrt(sy2)
syx2 <- sum((y-a-b*x)^2)/n
syx <- sqrt(syx2)

sig2 <- n*syx2/(n-2)
sig <- sqrt(sig2)
sigy2 <- n*sy2/(n-1)

rbar2 <- 1 - (sig2/sigy2)
rbar <- a / (sig*sqrt(n)*1)

####t-value
sx2 <- sum((x-mean(x))^2)/n
sx <- sqrt(sx2)

siga2 <- sig2 * sum(x^2) /(n^2*sx2)
siga <- sqrt(siga2)

sigb2 <- sig2/(n*sx2)
sigb <- sqrt(sigb2)

ta <- a / siga
tb <- b / sigb

####p-value
pa <- (1 - pt(ta, n-2))*2
pb <- (1 - pt(tb, n-2))*2

####Print
print(paste("a=", a))
print(paste("b=", b))
print(paste("rbar2=", rbar2))
print(paste("ta=", ta))
print(paste("tb=", tb))
print(paste("pa=", pa))
print(paste("pb=", pb))

タグ：最小二乗法 p値