統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

回帰直線の検定と信頼区間を改めて計算してみる

ブログランキングに参加しています。
まずはぽちぽちっとお願いします。
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
にほんブログ村 科学ブログ 数学へ

統計ソフト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 ")

散布図はこんな感じ。傾きはゼロではなさそうな感じがする。

f:id:toukeier:20190105130353p:plain

まず統計ソフト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(){} で書いた内容が実施されているわけである。傾きは統計学的に有意で( p_{Slope}=0.000673)、回帰直線は意味があるといえる。

> 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=で水平線を描いてくれる。今回求めた回帰直線を描き入れた図は以下の通り。ピッタリ予測できるとは言えないが、ある程度予測に使えそうな回帰直線に見える。

f:id:toukeier:20190105131748p:plain

まとめ

単回帰モデルの傾きと切片の推定、検定、95%信頼区間の算出、誤差分散、寄与率を改めて式のスクリプトを書いて計算してみた。改めてやってみると、傾きや切片の検定統計量が、どのように計算されているか、標準誤差はどのように計算されているかが鮮明にわかる。たまには初心に帰ってゼロから確認してみるのも楽しい。