オッズ比とその信頼区間 - ぼやかないつもりのメモ（ブログ Ver）

Rでオッズ比とその信頼区間を求めるのには、少し手間がかかる。群馬大の青木先生のR -- オッズ比には青木先生の自作関数を含め3つの方法が記されていた。同じく群馬大の中澤先生の「Rによる保健医療データ解析演習」には、中澤先生の自作関数を含め5つの方法が記されている。以下の例を参考にいくつかの方法をメモしておく。

	疾患あり	疾患なし
曝露あり	120	40
曝露なし	80	60

青木先生の関数を使う

R -- オッズ比にあるodds.ratio()関数をコピペし、使用法に従う。この場合は以下のようにすればよい。

odds.ratio <- function( a, b, c, d,             # 4 つのセルの観察値
                        correct=FALSE)          # 補正をするとき TRUE にする
{
        cl <- function(x)
        {
                or*exp(c(1, -1)*qnorm(x)*sqrt(1/a+1/b+1/c+1/d))
        }
        if (correct || a*b*c*d == 0) {          # a,b,c,d のいずれかが 0 のときには必ず補正する
                a <- a+0.5
                b <- b+0.5
                c <- c+0.5
                d <- d+0.5
        }
        or <- a*d/(b*c)
        conf <- rbind(cl90=cl(0.05), cl95=cl(0.025), cl99=cl(0.005), cl999=cl(0.0005))
        conf <- data.frame(conf)
        colnames(conf) <- paste(c("下側","上側"), "信頼限界値", sep="")
        rownames(conf) <- paste(c(90, 95, 99, 99.9), "%信頼区間", sep="")
        list(or=or, conf=conf)
}

R -- オッズ比

odds.ratio(120, 40, 80, 60, correct=FALSE)

$or
[1] 2.25

$conf
              下側信頼限界値 上側信頼限界値
90%信頼区間        1.4914023       3.394456
95%信頼区間        1.3784218       3.672678
99%信頼区間        1.1817222       4.284002
99.9%信頼区間      0.9883667       5.122087

中澤先生のoddsratio2()関数を使う

中澤先生の場合は、.Rファイルにて公開してくれているので、source()関数を使うことで利用できる。ただし、青木先生とは入力の順番が違うので注意した方がよい*1。

source("http://phi.med.gunma-u.ac.jp/msb/c10-5.R")
oddsratio2(120, 80, 40, 60)

> oddsratio2(120, 80, 40, 60)
       疾病あり 疾病なし 合計
曝露群      120       40  160
対照群       80       60  140
合計        200      100  300
オッズ比の点推定量: 2.25  (p= 0.001083837 ) 95%信頼区間 = [ 1.378422 , 3.672678 ]

vcdライブラリのoddsratio()関数を使う

この方法からは、以下のように分割表を作成して使う。vcdライブラリのoddsratio()関数は、青木先生が書かれたころより少し変更が加わったようでsummary()しても信頼区間は表示されない。ヘルプを読むとconfint()してあげる必要があるようだ。

library(vcd)
x <- matrix(c(120, 40, 80, 60), byrow=TRUE, nc=2)
y <- oddsratio(x, log=F)
summary(y)
confint(y)

> summary(y)
     Odds Ratio
[1,]       2.25
> confint(y)
          lwr      upr
[1,] 1.381450 3.664627

fisher.test()関数を使う

これは以前から知っていたが、どうやらオッズ比は推定量で、信頼区間は少し広くなるようだ。

fisher.test(x)

        Fisher's Exact Test for Count Data

data:  x 
p-value = 0.001372
alternative hypothesis: true odds ratio is not equal to 1 
95 percent confidence interval:
 1.339755 3.788809 
sample estimates:
odds ratio 
  2.243767

*1:もちろん、同じ結果が返ってくるが。