# 関数定義
# 　引数：　　file データファイル名
# 　　　　　　cochran Cochran の方法を採用するとき TRUE に設定
# 　　　　　　　　　　デフォールトでは Mantel-Haenszel の方法により計算
# 

Cochran.Mantel.Haenszel <- function(file, cochran = FALSE)
{
	data <- read.table(file, header=TRUE)
	attach(data)
	n1 <- n11+n12
	n2 <- n21+n22
	p1 <- n11/n1
	p2 <- n21/n2
	n <- n1+n2
	p.bar <- (n11+n21)/n
	if (cochran) {
		d <- (p1-p2)/(p.bar*(1-p.bar))
		w <- p.bar*(1-p.bar)*n1*n2/n
	}
	else {
		d <- (p1-p2)/(p.bar*(1-p.bar))*(n-1)/n
		w <- p.bar*(1-p.bar)*n1*n2/(n-1)
	}
	s.w <- sum(w)
	s.wd <- sum(w*d)
	d.bar <- s.wd/s.w
	SE.d.bar <- 1/sqrt(s.w)
	cat("標準差（standardized difference）の推定値 =", d.bar, "\n")
	cat("推定値の標準誤差 =", SE.d.bar, "\n\n")
	if (cochran == FALSE) {
		cat("見込み比の推定値 =", sum(n1*n2*p1*(1-p2)/n) / sum(n1*n2*p2*(1-p1)/n), "\n\n")
	}

	chi.sq.total <- sum(w*d*d)
	chi.sq.homog <- chi.sq.total-s.wd^2/s.w
	chi.sq.assoc <- (d.bar / SE.d.bar)^2
	name <- c("total", "homogenity", "association")
	value <- c(chi.sq.total, chi.sq.homog, chi.sq.assoc)
	df <- c(length(d), length(d)-1, 1)
	p <- pchisq(value, df, lower = FALSE)
	data.frame(source=name, chi.sq.=value, d.f.=df, p.value=p)
}



データファイル

各研究ごとに一行ずつ，以下の4つの数値を入力
最初の行は“n11	n12	n21	n22”であること

　　　　陽性数　陰性数
群１　　ｎ１１　ｎ１２
群２　　ｎ２１　ｎ２２

===== データファイル（cmh.data）開始 =====
n11	n12	n21	n22
38	67	33	72
56	136	61	113
43	102	33	112
===== データファイル（cmh.data）終了 =====



使用例

> Cochran.Mantel.Haenszel("cmh.data")
標準差（standardized difference）の推定値 = 0.04637281 
推定値の標準誤差 = 0.1477957 

見込み比の推定値 = 1.047385 

       source   chi.sq. d.f.   p.value
1       total 3.7588369    3 0.2887105
2  mohogenity 3.6603897    2 0.1603823
3 association 0.0984472    1 0.7537011


> Cochran.Mantel.Haenszel("cmh.data",cochran=TRUE)
標準差（standardized difference）の推定値 = 0.04653455 
推定値の標準誤差 = 0.1480532 

       source    chi.sq. d.f.   p.value
1       total 3.77149801    3 0.2872188
2  mohogenity 3.67270745    2 0.1593976
3 association 0.09879056    1 0.7532859