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

統計学備忘録 since2016

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

信頼区間のプロット

同じサイズのデータサンプルからt分布を利用した信頼区間の作図

まずは3×4の場合(サンプルサイズ3を4回実施する)

x <- matrix(NA,nrow=3,ncol=4)        #3×4の空セル

for (i in 1:4){                       #列数分乱数を代入
  x[,i] <- rnorm(3)              #標準正規分布の乱数を行数分繰り返す
}
mx <- rowMeans(x)    #行平均を一括計算
mean(x[1,])                #1行目のみ取り出して確認

vx <- apply(x,1,var)   #行分散を一括計算
var(x[1,])                  #1行目のみ取り出して確認

lx1 <- mx+qt(0.975,2)*sqrt(vx/3)    #95%信頼区間の上限(平均+t値*標準誤差), lower.tail=TRUE 
lx2 <- mx-qt(0.975,2)*sqrt(vx/3)     #95%信頼区間の下限(平均-t値*標準誤差),  (0.975,2, lower.tail=F )=qt(0.025,2,lower.tail=FALSE )

par(mfrow = c(1,2))       #1行2列に図を並べます
plot(lx1, ylim = c(-3,3), main = "上限のプロット", ylab="", xlab="")
plot(lx2, ylim = c(-3,3), main = "下限のプロット", ylab="", xlab="")
dev.off()          #作図を終了します

#t値はqt(0.025,2)~qt(0.975,2)になります, lower.tail=TRUE 
f:id:yoshida931:20180223160242p:plain:w600

行列数値を増やしていけば見えてきます 1000行10列でやってみます

x <- matrix(NA,nrow=1000,ncol=10)  
for (i in 1:10){   
  x[,i] <- rnorm(1000)   
}
mx <- rowMeans(x) 
vx <- apply(x,1,var) 
lx1 <- mx+qt(0.975,9)*sqrt(vx/10) 
lx2 <- mx-qt(0.975,9)*sqrt(vx/10)
par(mfrow = c(1,2))
plot(lx1, ylim = c(-3,3), main = "上限のプロット", ylab="", xlab="")
plot(lx2, ylim = c(-3,3), main = "下限のプロット", ylab="", xlab="")
dev.off()
f:id:yoshida931:20180223160331p:plain:w600

行数=r,列数=lとした場合

x <- matrix(NA,nrow=r,ncol=l)  
for (i in 1:l){   
  x[,i] <- rnorm(r)   
}
mx <- rowMeans(x) 
vx <- apply(x,1,var) 
lx1 <- mx+qt(0.975,l-1)*sqrt(vx/l) 
lx2 <- mx-qt(0.975,l-1)*sqrt(vx/l)
par(mfrow = c(1,2))
plot(lx1, ylim = c(-3,3), main = "上限のプロット", ylab="", xlab="")
plot(lx2, ylim = c(-3,3), main = "下限のプロット", ylab="", xlab="")
dev.off()

ベイズの定理でモンティ・ホール問題を考える

最終更新日2018-02-21
モンティ・ホール問題は不完全燃焼だったので、再々挑戦したいと思います.
今回は下記の文献をもとにベイズの定理を使って勉強していきます.

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

はじめての 統計データ分析 ―ベイズ的〈ポストp値時代〉の統計学―

参考web
yoshida931.hatenablog.com

 あるクイズ番組を想定して、あたり(○)とはずれ(×)で考えていきます.ABCの袋があり、それぞれの袋に〇か×のどちらかが入っていています.どれか一つに〇が入っており、司会者は答えを知っています(ポイント).ここで、あなたはAを選択したとします. 次に、あなたがAの袋を確かめる前に、司会者がBの袋を開けて見せてくれました.Bはハズレ(×)でした.「Aに当りのある確率は・・・? 」という問題です.残りA、Cのどちらかに入っているので、「確率は1/2」と考えるのが直観解とよばれている解答になります.

f:id:yoshida931:20180119115348p:plain:w250

ベイズの定理を用いてもうすこし深く考えてみたいと思います.まず、準備として以下のことを定義付けておきます.

P(A○)  # Aが当りの確率
P()  # Bがはすれの確率
P(C○)  # Cが当りの確率
P()  # Cがはずれの確率
求めたいのは、(B×)の場合の(A○)の確率です.

P(A○|B×) = \frac{P(B×|A○)*P(A○)}{P(B×)}



yoshida931.hatenablog.com

P(A○|)=P(|A○)*P(A○) / P()
# ここで分母のP(B×)を考えてみます
P() = P()*P(A○) + P()*P(C○) 
     = P(A○)*P(|A○) + P(C○)*P(|C○)
# 当然、当りは均等の確率なのでP(A○)=P(B○)=P(C○)=1/3です.
# P(B×|C○)は、Cが当たりの場合に司会者はCを開けません.必ずBを開けてB×(ハズレ)を示します.従って、P(B×|C○)の確率は1になります.
P()1/3 * P(|A○) + 1/3 * 1  
#したがって
P(A○|) =( P(|A○) * 1/3  )  /  ( 1/3 * P(|A○) + 1/3 * 1   )

つぎに P(B×|A○)について考えてみます
これは自分が選んだAがあたりの場合に、司会者がBを見せる確率になります.
Aがあたりの場合には、司会者はBまたはCを見せることになります.
つまり「Bを見せる確率」と「Cを見せる確率」の合計は1になります.

  P(|A○)    A○ →  B×
  P(|A○)    A○ →  C×

P(|A○) + P(|A○) = 1
P(|A○) = 1 - P(|A○)

P(B×|A○) は0~1の確立であることが理解できました.つまり次にようなことが考えられます.

# 例えば、P(C×|A○)=1/2の場合
P(|A○) = 1 - P(|A○)  = 1/2    求める解であるP(A○|) = 1/3 になり1/2(直観解)とは異なります

無情報的事前分布

直感的にはP(A○|B×)の確率は1/2(直観解)なのですが、「P(B×|A○)の確率が不明」と考えた方が妥当ということになります.
P(B×|A○)の確率, つまり「自分が選んだAが当たりの場合にBがハズレの確率」です.この確率を考えることが重要なポイントとなります.
このような場合にP(B×|A○)の確率として仮定するのが、無情報的事前分布となります.

Rを使って確率分布のシミュレーション
P(B×|A○)は確率なので0~1の一様分布を仮定します.求めたい確率P(A○|B×) はどのような分布になっているでしょうか.

P(A○|) = ( P(|A○) * 1/3  )  /  ( 1/3 * P(|A○) + 1/3 * 1   )
P(A○|) = y
P(|A○) = x

#Rを使用して一様分布[0,1]を1万個生成して
x <- runif(10000)
y <- x/(x+1)
hist(y)
P(A○|B×)の確率分布(事後分布
f:id:yoshida931:20180123100530p:plain:w300

確かに直観解の1/2が最も多い(MAP推定値)ので、そのまま1/2と答えても間違いではないのですが、大正解でもなさそうです.

summary(y)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000545 0.1994000 0.3325000 0.3073000 0.4298000 0.5000000 

P(A○|B×)の平均は0.307となりました(EAP推定値).
P(A○|B×)の中央値は1/3となりました(MED推定値).

問題の答え
Aに当りがある確率が1/2になる確率が一番多くなり、平均は1/3程度になります… といったところでしょうか?

まとめ

モンティ・ホール問題におけるP(A当たり|Bハズレ)の確率を求めるのですが、P(Bハズレ|A当たり)の確率には言及されていないません. P(Bハズレ|A当たり)の確率評価が重要になります.

f:id:yoshida931:20180221104529p:plain:w400

ロジスティック回帰分析(説明変数が単一かつ連続の場合)

今回は、この赤い線をどのようにして描いたか考えていきます

f:id:yoshida931:20180220163317p:plain:w400

ロジット関数とロジスティック関数

準備として関数の特徴を押さえておきます

ロジット関数 \ \  y=log\frac{P}{1-P} \ \  (0\verb|<|P\verb|<|1)\ \ \ →\ \ \ -∞\verb|<| f(P)\verb|<|∞

(標準)ロジスティック関数 \ \ ロジット関数の逆関数=P = \frac{1}{1+exp(-y)}

mathwords.net

サンプルirisより

品種"virginica=1"、"別の品種=0"という2値を目的変数、”Sepal.Length(がく片の長さ)(cm)”を説明変数とします. ロジスティック回帰で求めることは、”Sepal.Length(がく片の長さ)(cm)”から予測される、品種が"virginica=1"になる確率です.

class(iris$Species)  #属性を確認
is <- as.numeric(iris$Species)   #  iris$Speciesが要因になっているので数字に変換(as.integerでも同じ)
y <- gsub(1, 0, is)  # is の1を0に変換
y <- gsub(2, 0, y)  # is の2を0に変換
y <- gsub(3, 1, y)  # is の3を1に変換

x <- iris$Sepal.Length

#実際のデータは以下のようになっています
> y
  [1] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [24] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [47] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [70] "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0" "0"
 [93] "0" "0" "0" "0" "0" "0" "0" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
[116] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
[139] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
> x
  [1] 5.1 4.9 4.7 4.6 5.0 5.4 4.6 5.0 4.4 4.9 5.4 4.8 4.8 4.3 5.8 5.7 5.4 5.1 5.7 5.1 5.4 5.1 4.6
 [24] 5.1 4.8 5.0 5.0 5.2 5.2 4.7 4.8 5.4 5.2 5.5 4.9 5.0 5.5 4.9 4.4 5.1 5.0 4.5 4.4 5.0 5.1 4.8
 [47] 5.1 4.6 5.3 5.0 7.0 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6 6.7 5.6 5.8 6.2
 [70] 5.6 5.9 6.1 6.3 6.1 6.4 6.6 6.8 6.7 6.0 5.7 5.5 5.5 5.8 6.0 5.4 6.0 6.7 6.3 5.6 5.5 5.5 6.1
 [93] 5.8 5.0 5.6 5.7 5.7 6.2 5.1 5.7 6.3 5.8 7.1 6.3 6.5 7.6 4.9 7.3 6.7 7.2 6.5 6.4 6.8 5.7 5.8
[116] 6.4 6.5 7.7 7.7 6.0 6.9 5.6 7.7 6.3 6.7 7.2 6.2 6.1 6.4 7.2 7.4 7.9 6.4 6.3 6.1 7.7 6.3 6.4
[139] 6.0 6.9 6.7 6.9 5.8 6.8 6.7 6.7 6.3 6.5 6.2 5.9

ちなみに、そのまま単回帰分析を行ってみます.

summary(lm(formula=y~x) )
plot(x,y)
abline(lm(y~x)) #回帰直線
f:id:yoshida931:20180208184829p:plain:w300

目的変数は0、1なのですが、上図ではかなりはみ出してしまいます.0より少なくても、1より大きくても、良い予測とは言えません.したがって、この回帰直線は適当なモデルとは言えません.そこで指数を用いたロジスティックモデルを使用します.

ロジスティック変換 logistic transformation

従属変数が2値変数(品種=virginica=1、品種=別の品種=0)の回帰分析をおこなうためには、2値変数に与えられる上限と下限を取り除けば解決します.ロジスティック変換は、0~1の区間制限を取り除くことで、統計的推測に柔軟性を持たせた方法といえます.

品種が"virginica(y=1)"になる確率をPとした場合のロジスティック(ロジット)変換 を考えてみます.

オッズ=\frac{品種=virginica=1 である確率}{品種=他の品種=0 である確率}=\frac{p}{1-p}

ロジット(logit)とは、0から1の値をとる確率p に対し、そのオッズ\frac{p}{1-p}の対数から計算されます.

Pのロジスティック変換\ \  logit(p)=log\frac{P}{1-P}

f(P)=log\frac{P}{1-P} \ \  (0\verb|<|P\verb|<|1)\ \ \ →\ \ \ -∞\verb|<| f(P)\verb|<|∞

グラフで確認してみます

y <- function(p){ return(log(p/(1-p))) }
plot(y,xlab="P",ylab = "f(P)")
abline(0,0)

ロジットf(p)がどのような値になっても、0<p< 1となっています.

f:id:yoshida931:20180213133241p:plain:w300

ロジスティック回帰モデル

従属変数(virginica=1、別の品種=0)が1になる確率のロジットを従属変数にしたものが、ロジスティック回帰分析になります。

ロジスティック回帰モデル

\ \  f(P) = b0+b1*x

P = \frac{1}{1+exp(-(b0+b1*x))}

パラメータの推定(最尤推定法) P(yi=1|b0,b1)

尤度関数
L(θ)= \prod_{i=0}^{150} (\frac{1}{1+exp(-(b0+b1x))})^{yi})*(1-\frac{1}{1+exp(-(b0+b1x))})^{1-yi}

このままでは大変なので対数をとります
対数尤度関数 L(θ)= \sum_{i=0}^{150} log(\frac{1}{1+exp(-(b0+b1x))})^{yi}log(1-\frac{1}{1+exp(-(b0+b1x))})^{1-yi} 

yiは0か1です.例えば例題のy1と y150を見てみます.それぞれy1=1, y150=0となっています.
y1のときは
f\theta (y_1)=log(1-\frac{1}{1+exp(-(b0+b1x))})^{1-yi}
y150のときは
f\theta (y_{150})=log(\frac{1}{1+exp(-(b0+b1x))})^{yi}
となります.これを、1~150全て掛け合わせることになります.

ここからNewton法により最尤法を解いていくことになります.
私には解説が困難ので以下のサイトをご覧ください

qiita.com

Rを使ったロジスティック回帰分析

ここではRを使って楽させていただきます.一般化線形モデル(glm:Generalized Linear Models)は最尤法によりパラメータを推定する関数です! Rの関数はglmそのままです.

# glm関数で、応答変数の分布は二項分布なので引数をbinomialにすればOKのはず・・・
(ans <- glm(y~x,family=binomial))

Error in weights * y : non-numeric argument to binary operator
というエラーが出たので、例題で使用したxyの属性を調べます

typeof(x);typeof(y)
[1] "double"
[1] "character"

xは数で、yが文字になっていました.ということでyを数に置きなおして再トライ

y <- as.numeric(y)
(ans <- glm(y~x,family=binomial))

Call:  glm(formula = y ~ x, family = binomial)

Coefficients:
(Intercept)            x  
    -16.320        2.592  

Degrees of Freedom: 149 Total (i.e. Null);  148 Residual
Null Deviance:      191 
Residual Deviance: 117.3   AIC: 121.3

#やっぱりRは便利です

b0=-16.32, b1=2.596という解になりました

#対数尤度を求めてみます
x1 <- x[101:150]
x2 <- x[1:100]
e1 <- exp(-(-16.3198+2.592*x1))
e2 <- exp(-(-16.3198+2.592*x2))
(ly <- sum(log(1/(1+e1)))+sum(log(1-1/(1+e2))))

[1] -58.67273

視覚的にイメージしておきます.

plot(x,y,xlab = "",ylab = "")    #まずは散布図
abline(lm(y~x))    #回帰直線
par(new=T)   #図を重ねる命令
p <- function(x){ return(1/(1+exp(-(-16.320+2.592*x)))) } #pの関数を作ります
plot(p,ylim=c(0,1),xlim = c(4.3,7.9),col=2,
     ylab = "品種がvirginicaである確率",
     xlab="Sepal.Length(がく片の長さ)(cm)")  #赤線でp(確率)を挿入します
f:id:yoshida931:20180220114918p:plain

Rを使用して最尤法で推定したパラメータの有意性も確認することが可能です. 関数summaryを使用します.

summary(ans)

Call:
glm(formula = y ~ x, family = binomial)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9870  -0.5520  -0.2614   0.5832   2.7001  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -16.3198     2.6581  -6.140 8.27e-10 ***
x             2.5921     0.4316   6.006 1.90e-09 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 190.95  on 149  degrees of freedom
Residual deviance: 117.35  on 148  degrees of freedom
AIC: 121.35     # -2*対数尤度+2*パラメータ数

Number of Fisher Scoring iterations: 5

標準誤差の算出については、またいつか・・・

オッズについて

ロジスティック変換より

がく片の長さxcmのオッズ=\frac{P}{1-P} = exp(-16.3198+2.5921x)

がく片の長さxcmから1cm増えた場合のオッズ =\frac{P}{1-P} = exp(-16.3198+2.5921(x+1))

オッズ比= \frac{がく片の長さがx+1cmの場合に、品種がvirginicaである確率/品種がvirginica以外である確率}{ がく片の長さがxcmの場合に、品種がvirginicaである確率/品種がvirginica以外である確率}

\frac{\frac{P1}{1-P1}}{ \frac{P}{1-P}}= exp(-16.3198+2.5921(x+1))-exp(-16.3198+2.5921x)=exp(2.5921)=13.35646

xが1増えた場合の対数オッズ比は2.59、オッズ比は13.36になります.

 summary(x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.300   5.100   5.800   5.843   6.400   7.900 

xには大きな変動がないので、xが0.5増えた場合のオッズ比を求めてみます

\frac{\frac{P1}{1-P1}}{ \frac{P}{1-P}}=exp(2.5921×0.5)=3.654832

つまり、”Sepal.Length(がく片の長さ)”が5mm長くなると、品種が"virginica"である確率が約3.6倍高くなるということが言えます.オッズ比にはxが含まれていないので、どの幅(xcm)のに対しても同じことが言えます.


参考・引用

回帰分析入門 (Rで学ぶ最新データ解析)

回帰分析入門 (Rで学ぶ最新データ解析)

  • 作者: 豊田秀樹,尾崎幸謙,川端一光,岩間徳兼,鈴川由美,久保沙織,池原一哉,阿部昌利,大橋洸太郎,秋山 隆,大久保治信
  • 出版社/メーカー: 東京図書
  • 発売日: 2012/01/13
  • メディア: 単行本(ソフトカバー)
  • 購入: 2人 クリック: 36回
  • この商品を含むブログ (6件) を見る

LATEXでプレゼン

無料ソフトのみで統計からプレゼンまで!
spss, word, power pointを使用しないで以下のような2枚のスライドを作ってみました.

使ったもの

f:id:yoshida931:20180206114953p:plain:w500

f:id:yoshida931:20180205173355p:plain:w500
画像は オッズ比の信頼区間 - 統計学備忘録 since2016 からのコピペです

#LaTexの記載は以下の通りです.
\documentclass[slide,papersize]{jsarticle}
\usepackage[dvipdfmx]{color}
\usepackage{pxfonts}
\def \sheet #1{
\section*{\centering \large \bfseries #1}
}

\title{\LaTeX!でプレゼンに挑戦}

\usepackage[dvipdfmx]{graphicx}  %図
%http://www.wankuma.com/seminar/20090912nagoya09/5.pdf

\begin{document}

\maketitle

\begin{center}
   \includegraphics[width=60mm,scale=0.8]{WS000004.jpg} %絵も小さくしておきます
\end{center}
  
$p'A = \frac{a}{a+b}, p'B = \frac{c}{c+d}, \frac{\frac{p’A}{1-p’A}}{\frac{p’B}{1-p’B}}=\frac{p’A(1 - p'B)}{p’B(1 - p'A)}=\frac{\frac{a}{b}}{\frac{c}{d}}$

\end{document}  

全て、以下のページに書かれてます
http://www.wankuma.com/seminar/20090912nagoya09/5.pdf

オッズ比の信頼区間

オッズ比、見込み比(odds ratio)または交差積比(cross-product)

前提
比は分母が小さくなると、数値が大きくなりすぎて正規近似の精度が悪くなります.比の対数であれば高い精度で正規近似することが可能になります. したがって、比の対数を考えていくことになります.

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" 


 p'A = \frac{a}{a+b}    p'B = \frac{c}{c+d}

オッズ比=\frac{\frac{p’A}{1-p’A}}{\frac{p’B}{1-p’B}}=\frac{p’A(1 - p'B)}{p’B(1 - p'A)}=\frac{\frac{a}{b}}{\frac{c}{d}}

暴露なし(A群)が、暴露あり(B群)より疾病に罹患する可能性が大きいか、小さいかを表す尺度.
p'A, p'B が小さい場合には、リスク比と同じ値になります.


オッズ比の信頼区間

比の対数を考えていきます.
標本のオッズ比=標本オッズ比( \frac{\frac{a}{b}}{\frac{c}{d}} =\frac{ad}{cb}=or)

母集団のオッズ比を母オッズ比(OR)とします.

標本オッズ比の対数 log\ or、 母リスク比の対数 log OR

正規近似で信頼区間を求めていくためには、以下の式が必要になります.

Z = \frac{ log\ or - logOR}{ log\ or の標準誤差}

したがって、以下の式が95%信頼区間を求める式になります

log\ or - 1.96(log\ orの標準誤差) \leqq logOR \leqq  log\ or + 1.96(log\ orの標準誤差)

以下、log\ orの標準誤差 = SE とします

log\ or - 1.96SE \leqq logOR \leqq  log\ or + 1.96SE

or\times exp(-1.96\times SE) \leqq OR \leqq  or\times exp(1.96\times SE)


log\ orの分散

log\ or=log\ a + log\ d - log\ b - log\ c より

V(log\ or)=\frac{1}{a} + \frac{1}{d} + \frac{1}{b} + \frac{1}{c}  デルタ法(delta method)によって近似的に求めたもの(参考webより)

#Rで95%信頼区間を求めてみます
rr <- a*(c+d)/(a+b)/c      #  = (a/(a+b))/(c/(c+d)) 
exp(log(rr)+c(1, -1)*qnorm ( c(0.025,0.975) )*sqrt(b/a/(a+b)+d/c/(c+d)))

例)
クロス表イメージ
(x <- matrix(c(3, 5, 6, 8 ),2,byrow=T))
     [,1] [,2]
[1,]    3    5
[2,]    6    8

a<-3; b<-5; c<-6; d<-8      #参考webより
( or <-  a*d/(b*c) )        #オッズ比
[1] 0.8
or*exp(qnorm ( c(0.025, 0.975) )*sqrt( 1/a + 1/b + 1/c + 1/d) )
[1] 0.1348801  4.7449559

#Rの関数を使えば瞬間で色々出力してくれます(やっぱり凄い)
install.packages("Epi")  
library(Epi)            
twoby2(x)

2 by 2 table analysis: 
------------------------------------------------------ 
Outcome   : Col 1 
Comparing : Row 1 vs. Row 2 

      Col 1 Col 2    P(Col 1) 95% conf. interval
Row 1     3     5      0.3750    0.1254   0.7152
Row 2     6     8      0.4286    0.2065   0.6837

                                    95% conf. interval
             Relative Risk:  0.8750    0.2972   2.5763     #リスク比と95%信頼区間
         Sample Odds Ratio:  0.8000    0.1349   4.7450        #オッズ比と95%信頼区間
Conditional MLE Odds Ratio:  0.8081    0.0884   6.3780 #Fisherの正確検定によるオッズ比と95%信頼区間
    Probability difference: -0.0536   -0.3956   0.3312

             Exact P-value: 1       #Fisherの正確検定によるp値
        Asymptotic P-value: 0.8059   #正規近似によるp値
------------------------------------------------------

yoshida931.hatenablog.com


参考web1 統計学入門−第3章
参考web2 R による統計処理
参考文献)柳川 堯 ; 観察データの多変量解析―疫学データの因果分析,近代科学社 ,2016

リスク比の信頼区間

ポイント:比の対数をとり、正規近似する

リスク比は疫学における指標の1つです.一般的には相対危険度(相対リスク,relative risk,RR)として利用されています.

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}=p'A の母比率をPA

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

log\frac{a}{a+b}=log(p'A), log\frac{c}{c+d}=log(p'B)


リスク差(過剰絶対リスク) = p'A-p'B

リスク比 RR = \frac{p'A}{p'B} A群のB群における疾病の頻度の比

過剰リスク(過剰相対リスク) = \frac{p'A-p'B}{p'B} 

オッズ比=\frac{\frac{a}{b}}{\frac{c}{d}}


リスク比の信頼区間

比の対数を考えていきます.
標本のリスク比=標本リスク比( \frac{p'A}{p'B} =rr)、母集団のリスク比を母リスク比(\frac{PA}{PB} =RR)とします.

標本リスク比rrの対数 log\ rr、 母リスク比RRの対数 log RR

正規近似で信頼区間を求めていくためには、以下の式が必要になります.

Z = \frac{ log\ rr - logRR}{ log\ rr の標準誤差}

したがって、以下の式が95%信頼区間を求める式になります

log\ rr - 1.96(log\ rrの標準誤差) \leqq logRR \leqq  log\ rr + 1.96(log\ rrの標準誤差)

以下、log\ rrの標準誤差 = SE とします

log\ rr - 1.96SE \leqq logRR \leqq  log\ rr + 1.96SE

rr\times exp(-1.96\times SE) \leqq RR \leqq  rr\times exp(1.96\times SE)


log\ rrの分散

log\ rr=log\ a - log(a+b) - log\ c + log(c+d) より

V(log\ rr)=\frac{1}{a} - \frac{1}{a+b} - \frac{c}{c+d} + \frac{1}{c+d}  デルタ法(delta method)によって近似的に求めたもの(参考web1より)

#Rで95%信頼区間を求めてみます
rr <- a*(c+d)/(a+b)/c      #  = (a/(a+b))/(c/(c+d)) 
exp(log(rr)+c(1, -1)*qnorm ( c(0.025,0.975) )*sqrt(b/a/(a+b)+d/c/(c+d)))

例)
クロス表イメージ
matrix(c(76, 399, 129, 332),2,byrow=T)
     [,1] [,2]
[1,]   76  399
[2,]  129  332

a<-76; b<-399; c<-129; d<-332 #参考web2より
( rr <-  a*(c+d)/(a+b)/c )    #リスク比
[1] 0.5717829
rr*exp(qnorm ( c(0.025,0.975) )*sqrt(b/a/(a+b)+d/c/(c+d)))
[1] 0.4440632 0.7362370


参考web1 R -- 相対危険度(対応のない場合)
参考web2 R による統計処理
参考文献)柳川 堯 ; 観察データの多変量解析―疫学データの因果分析,近代科学社 ,2016

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

投稿日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