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

統計学備忘録 since2016

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

母比率の推定、母比率の差の検定

投稿日2017.6.19
更新日2018.1.30

比率の信頼区間

例)有権者から2,400人を無作為抽出した結果、1,250人は支持していたことがわかった.有権者の支持率の95%信頼区間を求めよ.
出典 日本統計学会 (編集); 日本統計学会公式認定 統計検定2級対応 統計学基礎, 東京図書, 2012, p149

二項分布の問題です
標本を試行回数n=2,400回、成功回数x=1,250回、成功確率p'=1,250÷2,400 として考えてみます
標本の成功確率p'から母比率(母集団の成功確率)Pの推定を目標とします.
nが十分に大きいと考えた場合、確率変数である成功回数 x は二項分布 B(n, P) に従うこととします.

確立変数X期待値 \ E(X)=nP,\ \ 分散\ V(X)=nP(1-P)

標本比率\ p' =\, \frac{x}{n}

標本比率p'の期待値 \ E(p')=P(母比率)

標本比率p'の分散V(p') = V(x/n) = \frac{V(x)}{n^{2}}=\frac{P(1-P)}{n}

二項分布の正規近似

二項分布B(n, p)は nが十分大きい場合,平均 np,分散 np(1-p)正規分布に近づきます(ラプラスの定理).したがって次のZは近似的に標準正規分布  N(0, 1) に従うことになります.

Z=\frac{p'-P}{\sqrt{\frac{P(1-P)}{n}}}

正規分布に従い、95%信頼区間を求めてみます

n=2400
x=1250
(x/n)+qnorm(c(0.025,0.975))*sqrt(((x/n)*(1-x/n))/n)

[1] 0.5008469 0.5408198


Rの関数 binom.test( )を使えば、簡単に算出できます

binom.test(1250, 2400)

    Exact binomial test

data:  1250 and 2400
number of successes = 1250, number of trials = 2400, p-value =
0.04328
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.5006234 0.5409923         # 二項分布からもとめた信頼区間
sample estimates:
probability of success 
             0.5208333 


比率の差の信頼区間

リスク差は疫学で用いられる指標の1つです.一般的には、寄与危険度(寄与リスク,絶対リスク)として利用されています. ここではリスク比の信頼区間、つまり母比率の差の信頼区間を推定してみます.

xm <- matrix(c("a","b","c","d"), nrow=2, byrow=T)
name <- list("暴露"=c("あり(A)","なし(B)"),"疾病"=c("あり","なし"))
dimnames(xm)<-name
xm

         疾病
暴露      あり なし
  あり(A) "a"  "b" 
  なし(B) "c"  "d" 

リスク差 = {\frac{a}{a+b}}-{\frac{c}{c+d}}

暴露あり群 = A群、暴露なし群 = B群と定義します.

A群の標本比率\ p'A=\frac{a}{a+b} 、母比率PA

B群の標本比率\ p'B=\frac{c}{c+d} 、母比率PB

標本比率piは 、近似的に正規分布N(pi, \frac{pi(1-pi)}{ni})にしたがいます.
それらの差も次の正規分布に従います
p'A-p'B ~ N(PA-PB, \frac{PA(1-PA)}{nA}+\frac{PB(1-PB)}{nB})

参考)分散の加法性
toukeigaku-jouhou.info

したがって

Z=\frac{(p'A-p'B)-(PA-PB)}{\sqrt{\frac{PA(1-PA)}{nA}+\frac{PB(1-PB)}{nB}}}

このZ値を使用して、95%区間を求めてみます.

xm2 <- matrix(c("116","76","192","244","44","288"), nrow=2, byrow=T)
name2 <- list("暴露"=c("あり(A)","なし(B)"),"疾病"=c("あり","なし","計"))
dimnames(xm2)<-name2
xm2
         疾病
暴露      あり  なし 計   
  あり(A) "116" "76" "192"
  なし(B) "244" "44" "288"

出典 日本統計学会 (編集); 日本統計学会公式認定 統計検定2級対応 統計学基礎, 東京図書, 2012, p156

# 正規分布からの95%信頼区間
va <- ((116/192)*(1-116/192)/192)+((244/288)*(1-244/288)/288)      #分散
se <- sqrt(va)      # 標準偏差
(116/192)-(244/288) + se*qnorm ( c(0.025,0.975) )
[1] -0.3237481 -0.1623630

# P値
帰無仮説:PA=PBより
z<- ( (116/192)-(244/288) ) / se       # Z値を求めます
pnorm(z)*2                                      #Rの関数でP値を求めます
  = 3.555521e-09   

Rの関数を使用してみます

#パッケージ "PropCIs"
install.packages("PropCIs")
library(PropCIs)
diffscoreci(116, 192, 244, 288, 0.95)

data:  

95 percent confidence interval:
 -0.3238230 -0.1627735

#パッケージ "epiR"
install.packages("epiR")
library(epiR)
grp1 <- matrix(cbind(116, 76), ncol = 2)
grp2 <- matrix(cbind(244, 44), ncol = 2)
dat <- as.matrix(cbind(grp1, grp2))
epi.conf(dat, ctype = "prop.unpaired")
         est         se      lower      upper
1 -0.2430556 0.04117041 -0.3227125 -0.1621579

# 関数BinomCIも試しましたが、計算が長かったので諦めました

A群とB群の支持率の差の95%信頼区間は、[ -0.32 ~ -0.16 ]  となりました.
有意水準5%では差があるのですが・・・AとBの支持率に差はあるのでしょうか?
Rの関数で差の検定を行ってみます.

prop.test(c(116, 244), c(192,288), alternative = c("two.sided"),correct=T)

    2-sample test for equality of proportions with continuity
    correction

data:  c(116, 244) out of c(192, 288)
X-squared = 35.012, df = 1, p-value = 3.278e-09
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.3280883 -0.1580228  #カイ二乗分布から求めた信頼区間
sample estimates:
   prop 1    prop 2 
0.6041667 0.8472222 

参考web
比率の差の信頼区間の手法色々 - r-statistics-fanの日記


P値はかなり小さな値になり、有意水準0.01でも差があるという結果になります.この信頼区間とP値から差を考察しなければいけません.統計学的な差はあるのですが、実質的な差の有無については研究者の解釈次第ということになります.

#検定結果は独立性の検定と同じ結果になります
chisq.test(matrix(c(116, 244, 76, 44), nrow=2))

    Pearson's Chi-squared test with Yates' continuity correction

data:  matrix(c(116, 244, 76, 44), nrow = 2)
X-squared = 35.012, df = 1, p-value = 3.278e-09

参考図書)柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010, p142