統計ソフトRを使えば、回帰直線の分析、つまり単回帰モデルの分析は本当に簡単にできる。
基本に帰って、実際どんな計算をしているか改めて見てみるのもいいかなと思い立った。
単回帰モデルの計算を教科書に沿ってやってみた。
教科書
教科書は「医学への統計学」。6.2 回帰直線の検定と信頼区間
新版はこちら。
データ
例題6.5のデータを使う。
cd <- c(1.3,2.1,2.7,1.6,2.3,1.7,1.8,2.1,1.8,1.7,1.4,1.6,2.1,1.7,1.9,2.2,1.2) mg <- c(1.7,2.6,3.0,0.1,3.9,1.6,1.5,1.7,2.6,1.8,1.3,2.1,2.9,1.7,2.2,2.4,0.6)
散布図を描いてみる
どんなデータもグラフ表示は欠かせない。連続量 対 連続量の場合は散布図。
plot(cd,mg,las=1,xlim=c(0,4),ylim=c(0,4),xlab="Cd(microg/100g)", ylab="beta2-MG(mg/L)") title(" Linear Regression Model ")
散布図はこんな感じ。傾きはゼロではなさそうな感じがする。
まず統計ソフトRに用意されている関数で解析してみる
lm()が線形モデル（Linear Model）の関数。単回帰モデルは、線形モデルのもっとも単純なモデル。summary()が結果のサマリーを出力して、confint()が傾き（Slope）と切片（Intercept）の95%信頼区間を表示する。
lm1 <- lm(mg~cd) summary(lm1) confint(lm1)
結果はこちら。cdのEstimateが傾きの推定値で、InterceptのEstimateが切片の推定値。Residual standard errorが誤差分散の平方根。Multiple R-squaredが寄与率だ。
> summary(lm1) Call: lm(formula = mg ~ cd) Residuals: Min 1Q Median 3Q Max -1.47132 -0.27257 0.05399 0.45524 1.10586 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.2237 0.7664 -1.597 0.131205 cd 1.7469 0.4093 4.268 0.000673 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6286 on 15 degrees of freedom Multiple R-squared: 0.5484, Adjusted R-squared: 0.5183 F-statistic: 18.22 on 1 and 15 DF, p-value: 0.0006734 > confint(lm1) 2.5 % 97.5 % (Intercept) -2.8573187 0.4099371 cd 0.8745512 2.6192144
母回帰直線の推定（傾き、切片（それぞれ検定と信頼区間））、誤差分散、寄与率
function(){} を使って、回帰直線の傾き（b）と切片（a）の推定、ゼロでないかどうかの検定（pb, pa）、95%信頼区間（bl, bu, al, au）、誤差分散（VE）、Residual standard error（rse）、寄与率（R2）をいっぺんに計算させる。
linear.reg <- function(x,y){ xbar <- mean(x) ybar <- mean(y) n <- length(x) b <- sum((x-xbar)*(y-ybar))/sum((x-xbar)^2) a <- ybar-b*xbar yhat <- a+b*x R2 <- sum((yhat-ybar)^2)/sum((y-ybar)^2) VE <- 1/(n-2)*sum((y-(a+b*x))^2) rse <- sqrt(VE) SEb <- sqrt(VE/sum((x-xbar)^2)) Tb <- b/SEb pb <- 2*pt(abs(Tb), n-2, lower.tail=F) bl <- b-qt(1-0.05/2, n-2)*SEb bu <- b+qt(1-0.05/2, n-2)*SEb SEa <- sqrt(VE*(1/n+xbar^2/sum((x-xbar)^2))) Ta <- a/SEa pa <- 2*pt(abs(Ta), n-2, lower.tail=F) al <- a-qt(1-0.05/2, n-2)*SEa au <- a+qt(1-0.05/2, n-2)*SEa list(round( c(Slope=b,"t_Slope"=Tb, "p_Slope"=pb, "LL_Slope"=bl, "UL_Slope"=bu, Intercept=a,"t_Intcp"=Ta, "p_Intcp"=pa, "LL_Intcp"=al, "UL_Intcp"=au, "Error Var."=VE,"Resid. SE"=rse,"R^2"=R2),6)) } linear.reg(x=cd, y=mg)
結果はこちら。lm(), summary(), confint()で出力した結果と同じことが見て取れる。つまり統計ソフトRのlm()の中の計算は、単回帰モデルの場合、もっとも単純化すれば、上記 function(){} で書いた内容が実施されているわけである。傾きは統計学的に有意で（）、回帰直線は意味があるといえる。
> linear.reg(x=cd, y=mg) [[1]] Slope t_Slope p_Slope LL_Slope UL_Slope 1.746883 4.268322 0.000673 0.874551 2.619214 Intercept t_Intcp p_Intcp LL_Intcp UL_Intcp -1.223691 -1.596591 0.131205 -2.857319 0.409937 Error Var. Resid. SE R^2 0.395101 0.628571 0.548445
回帰直線を描き入れてみる
冒頭で描いた散布図に回帰直線を描き入れてみる。先ほどの出力 list() に対して、unlist()を実施するとリスト内の数値を指定することができる。
linreg1 <- unlist(linear.reg(x=cd, y=mg)) c(linreg1[6],linreg1[1]) abline(a=linreg1[6], b=linreg1[1])
6番目が切片（a）で、1番目が傾き（b）だ。
> c(linreg1[6],linreg1[1]) Intercept Slope -1.223691 1.746883
abline() は、a=で切片、b=で傾きを指定すると、直線を描いてくれる。ほかに、v=で垂直線、h=で水平線を描いてくれる。今回求めた回帰直線を描き入れた図は以下の通り。ピッタリ予測できるとは言えないが、ある程度予測に使えそうな回帰直線に見える。
まとめ
単回帰モデルの傾きと切片の推定、検定、95%信頼区間の算出、誤差分散、寄与率を改めて式のスクリプトを書いて計算してみた。改めてやってみると、傾きや切片の検定統計量が、どのように計算されているか、標準誤差はどのように計算されているかが鮮明にわかる。たまには初心に帰ってゼロから確認してみるのも楽しい。