理学療法士がまとめたノート

統計学備忘録 since2016

Rを使って統計学を勉強するブログです

回帰係数の区間推定

赤文字をRで実行

 

サンプル)統計学入門(基礎統計学),東京大学出版会,1991,p258

次のデータをコピーしてください
東京 福岡
1019.4 1018.4
1005.7 1007.6
1002 1006.2
1006.7 1009.9
1005.1 1010.8
1010.1 1013.2
1016.7 1016.2
1011 1009.1
999.5 1003.1
1006.9 1012.5
1001.9 1006.4
1007.5 1006.3
1014.4 1012.2
1014.3 1015
1014.6 1017.4
1009 1016.5
1006.7 1012.1
1009.4 1008.7
1011.8 1009.2
1009.4 1009.2

x<-read.table("clipboard",header=T)

 

回帰分析
lm(x$東京~x$福岡)

 

Call:
lm(formula = x$東京 ~ x$福岡)

Coefficients:
(Intercept)       x$福岡 
    16.0876       0.9822 

標本回帰方程式
Y=16.09+0.98X

 

 

定義
yi   実測値   x$東京
xi   実測値   x$福岡
B 回帰係数(最小二乗法より) この信頼区間を求める
β 回帰係数の不偏推定量 
ei 誤差項
Ei   (回帰)残差
Txx    xの偏差平方和      ∑(xi-xi平均値)^2
Txy xyの偏差平方和    ∑(yi-yi平均値)(xi-xi平均値)
Tyy    yの偏差平方和      ∑(yi-yi平均値)^2
σ二乗 誤差の分散
s二乗    誤差の分散の不偏推定量

 


回帰残差
母回帰方程式 yi=α+βxi+Ei  Ei~(0,σ二乗)  (i=1,2,...,n)
回帰残差Ei=yi-α-βxi
Ei<-(x$東京-16.09-0.9822*x$福岡)

 

誤差の分散の不偏推定量
回帰残差の平方和を(n-2)で割ることで回帰方程式の当てはまりの良さがわかります
s二乗<-(sum(Ei^2))/(n-2)
s二乗<-(sum(Ei^2))/18

 

推定値の標準誤差
s<-sqrt(s二乗)

 

ほとんどのeiは0±2sの中に入入ります
残差の散布図で確認(言訳;グラフは初心者です)
plot(ei,type="p",ylim=c(-6,6))
abline(h=s, col='black' ,lty="twodash")
abline(h=-s, col='black' ,lty="twodash")
abline(h=2*s, col='black' ,lty="dashed")
abline(h=-2*s, col='black',lty="dashed")


graphics.off()

偏回帰係数の標本分布
σは未知なのでsに置き換えます
Z=(β-B)÷σ/√Txx ~N(0,σ二乗)
t=(β-B)÷s/√Txx ~t(n-2)

 

Bの信頼区間
β-t×s/√Txx ≦ B ≦ β+t×s/√Txx
Txx<-sum( (x$福岡-mean(x$福岡) )^2)
0.98-2.101*s/sqrt(Txx)
0.6189584
0.98+2.101*s/sqrt(Txx)
1.341042
     0.62 ≦ B ≦ 1.34

 

 

Rの関数で確認します
confint(lm(x$東京~x$福岡),level=0.95)

              
(Intercept) -348.9125481 381.087689
x$福岡         0.6211872   1.343239

 

単回帰分析

赤文字のみRで実行

 

線形回帰分析
説明変数と目的変数を直線関係で傾向を示す。
説明変数と目的変数との関係を直線でモデル化する回帰分析。

非線形回帰分析
非線形関係でモデル化する回帰分析。
説明変数と目的変数を非直線的関係で傾向を示す。

 

ここでは単回帰分析について解説します
以下のように定義します
観測値 ( xi , yi )
観測値の平均値 ( ya , xa )
予測式 y'= b0 + b1*xi

最小二乗法を用いて回帰係数を求めます
残差の二乗和(残差平方和)を最小にする回帰係数をもとめます

残差ei ( yi - y' )
残差平方和 Q( b0 , b1 )

Q(b0,b1) = ∑(yi-y')^2 = ∑{yi-(b0+b1*xi)}^2    最小二乗基準

Q(b0,b1) を最小にするb0とb1を求めることで回帰直線を得る

Q(b0,b1) はパラメータb0b1の二次関数で最小値が存在します
そこでb0b1偏微分し、0に等しいとします

 

f:id:yoshida931:20170317200530p:plain
Y=b0+b1*x1
この式をyのx上への回帰方程式、あるいは回帰直線といいます

b0 : y切片
b1 : 偏回帰係数

 

例)
統計学入門 ,東京大学出版会,p42,1991より
yi<-c(114,124,143,158,166)
xi<-c(35,45,55,65,75)

 

平均を求めます
(ya<-mean(yi))
141
(xa<-mean(xi))
55

 

偏回帰係数
B<-sum( (yi-ya)*(xi-xa) )/sum( (xi-xa)^2)
1.38

 

Y切片
A<-ya-1.38*55
65.1

y=65.1+1.38x

 

平方和の分解
次に平方和を分解してみます
観測値の平方和(または総平方和)をSyとします
Sy = ∑(yi-ya)^2
   = ∑{(yi-y')+(y'-ya)}^2
   = ∑(yi-y')^2+∑(y'-ya)^2+∑(yi-y')*(y'-ya)
となります
ここで∑(yi-y')*(y'-ya)を考えてみます
予測値と残差の相関は0になります(証明は省略)
ry'ei=Sy'ei/Sy'*Sei=0
Sy'ei=∑(yi-y')*(y'-ya)=0となります

したがって
     Sy     = ∑(yi-y')^2 + ∑(y'-ya)^2
     Sy     =    Sei      +    SR
観測値の平方和 =  残差平方和   +  回帰による平方和

 

決定係数
回帰直線で説明できる部分の割合
R^2 = SR / Sei
Se=0,R^2=1 完全に説明できる
SR=0,R^2=0 全く説明できない

さきほどの例の決定係数を求めてみます

決定係数<-sum( (65.1+1.38*xi-ya)^2)/sum( (yi-ya)^2)
0.9836777

 

Rの関数を使って確認してみます
summary(lm(formula=yi~xi) )

 

Call:
lm(formula = yi ~ xi)

Residuals:
   1    2    3    4    5
 0.6 -3.2  2.0  3.2 -2.6

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  65.1000     5.8284   11.17 0.001538 **
xi            1.3800     0.1026   13.45 0.000889 ***
---
Signif. codes: 
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.246 on 3 degrees of freedom
Multiple R-squared:  0.9837,    Adjusted R-squared:  0.9782
F-statistic: 180.8 on 1 and 3 DF,  p-value: 0.0008894

 

散布図と回帰直線
plot(xi,yi)
abline(lm(yi~xi))

f:id:yoshida931:20170613173429p:plain

参考) 豊田秀樹 (著, 編集);回帰分析入門 (Rで学ぶ最新データ解析) ,東京図書 ,2012

t検定のまとめ

修正日2017.6.15
修正日2017.6.19

Rでは次の式を使用します
t.test ( x, alternative=c("two.sided","less","greater") )
mu=0, paired=FALSE, var.equal=FALSE, conf.level=0.95)

引数
alternative:両側検定("two.sided")か片側検定("less","greater")か
mu=0:母平均が0であるかどうかを検定する
paired=TRUE:対応のあるt検定
var.equa=TRUE:等分散を仮定する
conf.level:信頼区間

統計量
t=(標本平均-母平均)/√(不偏分散/n) 
=(標本平均-母平均)/標準誤差  
自由度n-1のt分布t(n-1)に従う 

 

1標本問題(1つの母集団の母数に関する検定)

確率変数Xからのn個の標本x, μは母平均、S^2は不偏分散

     f:id:yoshida931:20170619141944p:plain
例)次のデータを図で概観し,母平均=160.0を有意水準5%で検定せよ
x = c(153.0,148.0,161.0,166.0,160.0,173.0,173.0,165.0,153.0,168.0,171.0)

箱ひげ図
boxplot(x)
ドット図
stripchart  (x, pch=16, at=0, method="stack")
ヒストグラム
hist  (x, right=FALSE, col="gray")

t検定を実行

t.test (x)

One Sample t-test
data: x
t = 62.79, df = 10, p-value = 2.554e-14
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
157.0405 168.5959
sample estimates:
mean of x
162.8182

 

標本xの平均は162.82
95%信頼区間[157.04, 168.60]で「160.0」が含まれている
したがって帰無仮説は5%有意水準で採択される(p値<0.01)

もし、この検定の有意水準が1%であった場合は

t.test(x,conf.level=0.01)

One Sample t-test
data: x
t = 62.79, df = 10, p-value = 2.554e-14
alternative hypothesis: true mean is not equal to 0
1 percent confidence interval:
162.7849 162.8515
sample estimates:
mean of x
162.8182

 

99%信頼区間[162.78, 162.85]で「160.0」が含まれていません
したがって帰無仮説は1%有意水準では棄却されることになります
研究前のαの設定は重要です

 

対応のある2標本の場合
1標本問題と同じ原理です
治療前後の比較を考えてみましょう
治療前のデータPRE,治療後のデータPOSTとします
治療の前後差を有意水準1%で検定せよ
PRE<-c(25,27,32,37,38,39,44,45,47,48,52,54,59,60,65,65,67,68,68)
POST<-c(41,28,32,60,50,48,44,55,95,48,61,67,60,83,66,77,70,81,87)
対応あるt検定の帰無仮説は「前後差=0」と設定します
def<-POST-PRE
つまり1標本となります

t.test(def)

One Sample t-test
data: def
t = 4.1537, df = 18, p-value = 0.0005965
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
5.540338 16.880715
sample estimates:
mean of x
11.21053

 

差の平均は11.21となり有意水準1%で差を認める
(当然この11.21という数値は検査内容により実質的な意味合いが異なります)

 

2標本問題(2つの母集団の母数に関する検定)
ここでは2標本の平均の差の検定のみを扱います
確率変数Xからのm個の標本x, 確率変数Yからのn個の標本y
標本x,y   標本数m,n   標本平均μx,μy   標本分散sx, sy  

 

母分散が既知の場合
標本x (m個) ,標本y (n個) は母分散が等しい正規分布N(μ1,σ^2), N(μ2,σ^2)に従う
標本xの平均μx~N(μ1,σ^2/m)、標本yの平均μy~N(μ2,σ^2/n)
μxμy も正規分布に従う
E[μxμy]=μ1-μ2、V[μxμy] = V(μx)+V(μy) = (σ^2)/m + (σ^2)/n

ここをすぐに忘れてしまうので、ちょいと書き加えておきます.
AとBを互いに独立な確率変数として、それぞれの期待値をμa、μbとします.
V(A-B) = E { (A-B) - (μa-μb) }^2 = E { (A-μa) - (A-μb) } ^2
            = V(A) + V(B) - 2E{ (A-μa)*(B-μb) }
共分散は確率変数X,Yが互いに独立な場合には0になるので、
E{ (Aa)*(Bb) }=0
したがって
V(A-B) = V(A) + V(B)

 母分散が既知の場合には、等しくなくても xme-yme を利用して検定可能

            Z = { (μxμy) - (μ1-μ2) } ÷ √(σ^2/m+σ^2/n)

 

母分散が未知で等しい場合

f:id:yoshida931:20170619134940p:plain   f:id:yoshida931:20170619140210p:plain

f:id:yoshida931:20170619135001p:plain


定量はプールした分散s^2(合併した分散)を使用

プールした分散
∑( x - μx )^2   +  ∑( y -μy )^2 ÷ ( m+n-2 )
標本分散sx, sy は
sx=( 1/m-1 )*∑( x - μx )^2   ,   sy=( 1/n - 1 )*∑( y - μy )^2

したがってプールした分散s^2
s^2 = { (m-1 )*sx + ( n-1 )*sy }  ÷  ( m+n-2 )  


t=  (μx - μy) - ( μ1 - μ2 )  / √( s^2/m + s^2/n )

t= (μx - μy) - ( μ1 - μ2 ) / {s*( √1/m + 1/n ) }
自由度 m+n-2のt分布t(m+n-2)にしたがう

 

例)x若年群の片足立ち(sec)、y高齢群の片足立ち(sec)
x<-c(39,47,48,59,65,68,74,79,80,82,87,88,89,94,99,
100,101,106,107,108,111,120,122,129,130,136,138,139,
147,157,158,179,189,197,213)
y<-c(25,27,32,37,38,44,45,52,54,60,65,67,68,69,71,73,
82,88,89,92,99,100,101,116,126,131,151,156,174)

t.test(x, y, var.equal=T)

Two Sample t-test
data: x and y
t = 2.9234, df = 62, p-value = 0.004829
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
9.671654 51.500760
sample estimates:
mean of x mean of y
111.00000 80.41379

有意水準1%で差があるという結果になるのだが、
95%信頼区間は[9.67,51.50]と範囲が大きい
約9秒の差、50秒の差をどのようにとらえるかは研究者が十分に考察しなければならない

 

検定力分析(事後の分析)
検定力1-βとは、帰無仮説が「偽」であるときに、有意差が検出される確率.
効果量(effect size)とは、母集団における効果の大きさ

標本効果量(cohensのd効果量)の求め方
cohensのd効果量=
  xの平均値とyの平均値の差の絶対値  ÷  2群の共通の標準偏差

cohens_d <- function(x, y) {    
lx <- length(x)- 1         #n‐1を算出
ly <- length(y)- 1         #n‐1を算出
md <- abs(mean(x) - mean(y))   # x群とy群の平均値の差の絶対値
var2 <- lx * var(x) + ly * var(y)   # ( n‐1) × 不偏分散
pov <- var2  / ( lx + ly )      
                #プールした分散 s^2 = { (m-1)*sx + (n-1)*sy } ÷ (m+n-2) 
cd <- md/ sqrt(pov)   
     # sqrt(pov) = 2群に共通の標準偏差
}
res <- cohens_d(x, y)
res
[1] 0.7340753

 http://tjo.hatenablog.com/entry/2014/02/24/192655

 
次にt値から標本効果量を算出
A=length(x)+length(y) 
B=length(x)*length(y) 

標本効果量(cohensのd効果量) = t値 × A / B            
2.9234*sqrt(64/(35*29))
[1] 0.7400945

 

検定力
次に5%水準で検定力を求めてみます
library(pwr)
pwr.t2n.test(n1=35,n2=29, d =0.7752781 , sig.level =0.05 , power = NULL,alternative = "greater")

t test power calculation

n1 = 35
n2 = 29
d = 0.7752781
sig.level = 0.05
power = 0.9205285
alternative = greater

検定力は0.92となり
検定力の高い検定であったといえます


母分散が未知で等しくない場合(ウェルチの検定)
まず等分散性を検定する
Leveneの検定
帰無仮説:xの分散=yの分散
library(car)

score<-c(100,100,101,101,106,107,108,111,116,
120,89,100,145,85,106,107,90,149,116,78)
group=factor(rep(c("x", "y"), c(10,10)))
leveneTest(score ~ group)

Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 6.3953 0.02101 *
18

帰無仮説は棄却され、母分散が等しくないことが判明したため
ウェルチの検定を行う

t=(標本xの平均-標本yの平均)/√( s1^2/m + s2^2/n )
このtは帰無仮説の下で近似的に次の自由度fをもつt分布に従う
g1=s1^2/m ,g2=s2^2/n
f=(g1+g2)^2/{(g1^2/(m-1))-(g2^2/(n-1))}
このtを用いる方法はウェルチの検定(Welch's test)とよばれる

ベクトルscoreよりデータを抽出する
x<-score[1:10]
y<-score[11:20]

ウェルチの補正による2標本のt検定
t.test(x,y,var=F)

Welch Two Sample t-test
data: x and y
t = 0.062759, df = 10.485, p-value = 0.9511
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-17.14089 18.14089
sample estimates:
mean of x mean of y
107.0 106.5

Rで簡単(適合度検定)

例1)次のデータはメンデルの法則 1:2:1 の確率に適合しているか?
AA Aa aa n
52 102 46 200
帰無仮説:標本はメンデルの法則に適合している

観測度数<-c(52,102,46)
確率<-c(1/4,1/2,1/4)      
chisq.test(観測度数, p=確率)

Chi-squared test for given probabilities
data: 観測度数
X-squared = 0.44, df = 2, p-value = 0.8025

結果:帰無仮説は棄却されず、標本はメンデルの法則によく適合していることが確認された。(完全に適合しているという意味ではありません)

 

例2)男性は女性よりも多いと言えるか
男性 女性
229  179
帰無仮説:男性と女性の数は同じである(1/2, 1/2に適合している)

観測度数<-c(229,189)
確率<-c(1/2,1/2)      
chisq.test(観測度数, p=確率)

Chi-squared test for given probabilities

data: 観測度数
X-squared = 3.8278, df = 1, p-value = 0.05041

結果:有意水準5%では帰無仮説は棄却されず、男女数に差があるとは言えない。
(「男女比が1:1に適合している」という結果になります.あくまで有意水準0.05での結果です.)

Rで簡単(独立性の検定)

赤文字のみをRで実行

次のクロス表を検定してみます
問)男性と女性で、yesの解答者数とnoの解答者数の比率に違いはあるか?

   Yes N0
男性 n1 n2
女性 n3 n4

帰無仮説:男性と女性では解答者数の比率に違いがない

例)仮想データ

n1<-15
n2<-11
n3<-4
n4<-16

観測度数を行列で作成して解いてみます
(x<-matrix(c(n1,n2,n3,n4),nrow=2,ncol=2, byrow=T))

カイ二乗検定
chisq.test(x, correct=F)

Pearson's Chi-squared test

data: x
X-squared = 6.6244, df = 1, p-value = 0.01006

結果:帰無仮説は有意水準5%で棄却され、男女間での解答の比率に違いが認められた。

全セル数の20%以上で期待値が5未満の場合には次のいずれかの補正を実施
Yates(イェ-ツ)補正
chisq.test(x)
フィッシャーの正確確率検定
fisher.test(x)


カイ二乗検定の結果を代入する
RES <- chisq.test(x)

期待値
RES$expected

[,1] [,2]
[1,] 10.73913 15.26087
[2,] 8.26087 11.73913

標準化残差(残差÷√期待値)
RES$residuals

[,1] [,2]
[1,] 1.300211 -1.090708
[2,] -1.482468 1.243599

調整済み(標準化)残差
RES$stdres

[,1] [,2]
[1,] 2.573799 -2.573799
[2,] -2.573799 2.573799

ファイ係数
(φ<-(n1*n4-n2*n3)/sqrt(sum(x[1,])*sum(x[2,]) *sum(x[,1])*sum(x[,2])))

 0.379486

結果:男性はYES、女性はNOと解答するという正の連関を認めた。

連関係数<|0.2| 連関はほとんどない
|0.2|≦連関係数<|0.4| やや連関がある
|0.4|≦連関係数<|0.7| 強い連関がある
|0.7|≦連関係数 かなり強い連関がある
   対馬栄輝; よくわかる医療統計 -「なぜ?」にこたえる道しるべ, 東京図書, 2015

 

二項分布のグラフ

こんな説明でどうでしょうか?
  f:id:yoshida931:20170524161843p:plain

コイン投げ(いかさまなし)
表の出る確率確率1/2
3回コイン投げを実施してみた場合・・・
表が出る回数(x)は0回、1回、2回、3回
xの確率は𝐵𝑖(3 , 0.5) に従う.
その確率分布を二項分布といいます.

    f:id:yoshida931:20170524161933p:plain

f:id:yoshida931:20170524161906p:plain

 

 

村上 正康,安田 正実; 統計学演習, 培風館,1989,p51のグラフをRで作ってみます
B(10,0.2),B(10,0.5),B(10,0.8)を重ねたグラフです
赤文字のみRで実行してください

 (foo.bar@baz.com様よりコメントいただき修正しております。簡単で理解しやすい方法です。)

t0.2 <- 0:6
P0.2 <- dbinom(t0.2, 10, 0.2)
t0.5 <- 1:9
P0.5 <- dbinom(t0.5, 10, 0.5)
t0.8 <- 4:10
P0.8 <- dbinom(t0.8, 10, 0.8)
plot(t0.2, P0.2, type="b", lty=2, xlab="", ylab="", xlim=c(0,10), ylim=c(0,0.4))
lines(t0.5, P0.5, type="b", lty=2)
lines(t0.8, P0.8, type="b", lty=2)

f:id:yoshida931:20170121213325p:plain

 

 

 

途中経過を忘れないために残しておきます

 
n=10 、確率0.2 
t0.2 <- 0:6;p<-0.2;n<-10
y<-factorial(n)*(p^t0.2)*(1-p)^(n-t0.2)
r<-factorial(t0.2)*factorial(n-t0.2)
y/r

n=10、確率0.5
t0.5 <- 1:9;p<-0.5;n<-10 
c<-factorial(n)*(p^t0.5)*(1-p)^(n-t0.5)
d<-factorial(t0.5)*factorial(n-t0.5)

n=10、確率0.8
t0.8 <- 4:10;p<-0.8;n<-10 
e<-factorial(n)*(p^t0.8)*(1-p)^(n-t0.8)
f<-factorial(t0.8)*factorial(n-t0.8)
e/f


p=0.2,0.5,0.8の結果をそれぞれ算出
間違わないように名前付け
P0.2<-(y/r)
P0.5<-(c/d)
P0.8<-(e/f)

linesでグラフを重ねる

plot(t0.2,P0.2,type="b",lty=2,xlab="",ylab="",xlim=c(0,10),ylim=c(0,0.4))
lines(t0.5,P0.5,type="b",lty=2)
lines(t0.8,P0.8,type="b",lty=2)

 

グラフに色をつける

library(RColorBrewer)
#RColorBrewerパッケージのサンプル
display.brewer.all()
 
ヒストグラムを塗ってみる
どの色セットを使用するかを指定する
cols <- brewer.pal(8,"Pastel1")   # brewer.pal(何色、パレット名)

y<-c(1,2,3,4,5,6,7)
p<-c(2,3,4,5,4,3,2)
q<-c(2,3,4,5,4,3,2)
par(mfrow = c(3,3),mar = c(5, 4, 1, 4)) #余白 底辺、左、上、右の順 
pos.x <- barplot(q,ylim=c(0,6),col=cols[1])
pos.x <- barplot(q,ylim=c(0,6),col=cols[2])
pos.x <- barplot(q,ylim=c(0,6),col=cols[3])
pos.x <- barplot(q,ylim=c(0,6),col=cols[4])
pos.x <- barplot(q,ylim=c(0,6),col=cols[5])
pos.x <- barplot(q,ylim=c(0,6),col=cols[6])
pos.x <- barplot(q,ylim=c(0,6),col=cols[7])
pos.x <- barplot(q,ylim=c(0,6),col=cols[8])
pos.x <- barplot(q,ylim=c(0,6),col=cols[9])

 

f:id:yoshida931:20170114232836p:plain