R

前回に引き続き、統計の備忘録。

オッズ比とP値はエクセルで計算可。

X1 X2 Total
Y1 a b a+b
Y2 c d c+d
Total a+c b+d a+b+c+d

上記のようなデータがあった場合、期待値は、

期待値 X1 X2 Total
Y1 (a+c)*(a+b) /(a+b+c+d) (a+b)*(b+d) /(a+b+c+d) a+b
Y2 (a+c)*(c+d) /(a+b+c+d) (b+d)*(c+d) /(a+b+c+d) c+d
Total a+c b+d a+b+c+d

で求められる。

カイ二乗値は、(実測値ー期待値)^2/期待値 でそれぞれの欄を計算。

カイ二乗値 X1 X2 Total
Y1 実測値ー期待値)^2/期待値 実測値ー期待値)^2/期待値
Y2 実測値ー期待値)^2/期待値 実測値ー期待値)^2/期待値
Total

5%棄却限界は、

=CHISQ.INV.RT(0.05,自由度)

自由度は(列ー1x行ー1)なので、上記の場合は1。

★の数値が5%棄却限界より大きかどうかで判断する。

P値は、

=CHISQ.DIST.RT(★,1)

で出せる。

 

信頼区間(とオッズ比)は、青木先生のサイトで完結。

まず、下記プログラムをRコンソールにペースト

# オッズ比を計算する
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)
}

使用例:> odds.ratio(a, b, c, d)

> odds.ratio(9, 184, 4, 222)

出力例:

$or
[1] 2.71467
$conf 下側信頼限界値
90%信頼区間 0.99673
95%信頼区間 0.82265
99%信頼区間 0.56532
99.9%信頼区間 0.36578
上側信頼限界値
90%信頼区間 7.39361
95%信頼区間 8.95814
99%信頼区間 13.0359
99.9%信頼区間 20.1471

 

瞬殺。

青木先生、いつもありがとうございます。

 

にほんブログ村 科学ブログ 脳科学へ
にほんブログ村 子育てブログ 海外ワーキングマザー育児へ
PVアクセスランキング にほんブログ村