目的

Breslow-Day 検定を行う
他のサイトに，Tarone の補正を行うプログラムがある。
metafor ライブラリの rma.mh でも行える。

使用法

BD.test(m)

引数

m	2×2×k の配列

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/BD_test.R", encoding="euc-jp")

# Breslow-Day 検定
BD.test <- function(m)               # 2×2×k の配列
{
        method <- "Breslow-Day 検定"
        data.name <- deparse(substitute(m))
        k <- dim(m)[3]
        Nk <- apply(m, 3, sum)
        psiMH <- sum(m[1, 1,]*m[2, 2,]/Nk)/sum(m[2, 1,]*m[1, 2,]/Nk)
        nk1 <- m[1, 1,]+m[1, 2,]
        nk2 <- m[2, 1,]+m[2, 2,]
        xkt <- m[1, 1,]+m[2, 1,]
        a1 <- psiMH-1
        b1 <- -psiMH*(nk1+xkt)-nk2+xkt
        c1 <- psiMH*nk1*xkt
        e <- (-b1-sqrt(b1^2-4*a1*c1))/(2*a1)
        v <- 1/(1/e+1/(nk1-e)+1/(xkt-e)+1/(nk2-xkt+e))
        chisqBD <- sum((m[1, 1,]-e)^2/v)
        df <- k-1
        p <- pchisq(chisqBD, df, lower.tail=FALSE)
        return(structure(list(statistic=c("chi sq."=chisqBD), parameter=c(df=df), p.value=p,
                method=method, data.name=data.name), class="htest"))
}


使用例

> m <- array(c(
+ 	1, 8, 0, 10,
+ 	2, 47, 0, 60,
+ 	3, 29, 1, 13,
+ 	3, 17, 0, 7,
+ 	13, 45, 1, 45,
+ 	13, 26, 0, 18,
+ 	8, 15, 0, 10,
+ 	11, 4, 0, 4
+ 	), dim=c(2, 2, 8))

> m            # 引数の構造
, , 1

     [,1] [,2]
[1,]    1    0
[2,]    8   10

, , 2

     [,1] [,2]
[1,]    2    0
[2,]   47   60

　　途中省略

, , 8

     [,1] [,2]
[1,]   11    0
[2,]    4    4

> BD.test(m)

	Breslow-Day 検定

data:  m 
chi sq. = 8.9432, df = 7, p-value = 0.2568

参考文献

1) Breslow NE and Day NE (1980) Statistical Methods in Cancer Research:
   Volume 1-The Ananlysis of Case-Control Studies, IARC Scientific Pub.

2) 宮原英夫, 丹後俊郎　編 (1995) 医学統計学ハンドブック, 朝倉書店.