比率の差の多重比較（pairwise.prop.test の拡張） Last modified: Jan 05, 2010
目的
比率の差の多重比較を行うために，R に用意されている pairwise.prop.test を拡張する
1. pairwise.prop.test もとのままの機能を持つ（二群の比率の差の検定の多重比較）
2. pairwise.prop.test は prop.test を下請けに使っているが，fisher.test も使えるようにする
3. 分布の差の検定を行えるようにする（chisq.test と fisher.test）
使用法
paiwise.prop2.test(r, x, p.adjust.method=p.adjust.method, test.function=prop.test, ...)
引数
x 比率の差（2 カテゴリー）ならば，ベクトルまたは列数 2 の行列。または，列数 2 以上の行列
n 比率の差（2 カテゴリー）ならば，ベクトル。それ以外の場合は無視される
p.adjust.method P 値の調整法
test.function 下請けに使う関数名 prop.test の他に，fisher.test, chisq.test が使える
... pairwise.prop.test で使えるその他の引数
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/pairwise.prop2.test.R", encoding="euc-jp")
# 比率の差（分布の差）の多重比較（prop.test, fisher.test, chisq.test を用いる）
pairwise.prop2.test <- function ( x, # 比率の差（2 カテゴリー）ならば，ベクトルまたは列数 2 の行列。または，列数 2 以上の行列
n, # 比率の差（2 カテゴリー）ならば，ベクトル。それ以外の場合は無視される
p.adjust.method = p.adjust.methods, # P 値の調整法
test.function=prop.test, ...) # 下請けに使う関数名 prop.test の他に，fisher.test, chisq.test が使える
{
p.adjust.method <- match.arg(p.adjust.method)
METHOD <- deparse(substitute(test.function))
DNAME <- deparse(substitute(x))
if (is.matrix(x)) {
if (ncol(x) < 2)
stop("'x' must have at least 2 columns")
} else if (is.vector(x) && is.vector(n))
x <- cbind(x, n-x)
else
stop("'x' must be a matrix, or 'x', and 'n' must be a vector")
if (nrow(x) < 2)
stop("too few groups")
compare.levels <- function(i, j) {
test.function(x[c(i, j),], ...)$p.value # test.function で使用する検定を使い分ける
}
level.names <- names(x)
if (is.null(level.names))
level.names <- seq_along(1:nrow(x))
PVAL <- pairwise.table(compare.levels, level.names, p.adjust.method) # R の pairwise.*.test は有意水準ではなく P 値を調整するので解釈時に注意
ans <- list(method = METHOD, data.name = DNAME, p.value = PVAL,
p.adjust.method = p.adjust.method)
class(ans) <- "pairwise.htest"
ans
}
使用例
★ 比率の差の検定
# pairwise.prop.test の example にある例
#
smokers <- c( 83, 90, 129, 70 )
patients <- c( 86, 93, 136, 82 )
#
# このようなデータ
# group smokers non-smokers patients
# 1 83 3 86
# 2 90 3 93
# 3 129 7 136
# 4 70 8 82
#
# pairwise.prop.test を使う場合の指定法 -- 以下のような 2 つのベクトルを指定する。第 2 引数は，合計数
# c( 83, 90, 129, 70 ) と c( 86, 93, 136, 82 )
#
pairwise.prop.test(smokers, patients)
#
# pairwise.prop2.test を使う場合の指定法 -- pairwise.prop.test と同じ
#
pairwise.prop2.test(smokers, patients)
#
# pairwise.prop2.test を使う場合の別法 -- 以下のような行列，つまり，分割表の本体を列数 2 の行列として指定する
#
# smokers non-smokers
# 83 3
# 90 3
# 129 7
# 70 8
#
pairwise.prop2.test(cbind(smokers, patients-smokers))
#
# 検定関数として fisher.test を使う
#
pairwise.prop2.test(cbind(smokers, patients-smokers), test.function=fisher.test)
★ カテゴリーが 3 つ以上の場合（分布の差の検定）
# 以下のような集計表 4 群それぞれ，4 つのカテゴリーに該当するケース数が集計されている
# group Cat1 Cat2 Cat3 Cat4
# 1 22 41 52 31
# 2 12 31 25 14
# 3 12 25 14 17
# 4 23 14 16 28
#
# 行列の作り方は色々な方法がある（table 関数で集計したり，集計結果を matrix 関数で行列にするなど）
#
a <- rbind(c(22,41,52,31), c(12,31,25,14), c(12,25,14,17), c(23,14,16,28))
#
# 検定関数としては，chisq.test または fisher.test を使う
#
pairwise.prop2.test(a, test.function=chisq.test)
pairwise.prop2.test(a, test.function=fisher.test)
出力結果例
★ 比率の差の検定
> # pairwise.prop.test の example にある例
> #
> smokers <- c( 83, 90, 129, 70 )
> patients <- c( 86, 93, 136, 82 )
> #
> # このようなデータ
> # group smokers non-smokers patients
> # 1 83 3 86
> # 2 90 3 93
> # 3 129 7 136
> # 4 70 8 82
> #
> # pairwise.prop.test を使う場合の指定法 -- 以下のような 2 つのベクトルを指定する。第 2 引数は，合計数
> # c( 83, 90, 129, 70 ) と c( 86, 93, 136, 82 )
> #
> pairwise.prop.test(smokers, patients)
Pairwise comparisons using Pairwise comparison of proportions
data: smokers out of patients
1 2 3
2 1.000 - -
3 1.000 1.000 -
4 0.119 0.093 0.124
P value adjustment method: holm
警告メッセージ：
1: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
カイ自乗近似は不正確かもしれません
2: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
カイ自乗近似は不正確かもしれません
3: In prop.test(x[c(i, j)], n[c(i, j)], ...) :
カイ自乗近似は不正確かもしれません
> #
> # pairwise.prop2.test を使う場合の指定法 -- pairwise.prop.test と同じ
> #
> pairwise.prop2.test(smokers, patients)
Pairwise comparisons using prop.test
data: smokers
1 2 3
2 1.000 - -
3 1.000 1.000 -
4 0.119 0.093 0.124
P value adjustment method: holm
警告メッセージ：
1: In test.function(x[c(i, j), ], ...) :
Chi-squared approximation may be incorrect
2: In test.function(x[c(i, j), ], ...) :
Chi-squared approximation may be incorrect
3: In test.function(x[c(i, j), ], ...) :
Chi-squared approximation may be incorrect
> #
> # pairwise.prop2.test を使う場合の別法 -- 以下のような行列，つまり，分割表の本体を列数 2 の行列として指定する
> #
> # smokers non-smokers
> # 83 3
> # 90 3
> # 129 7
> # 70 8
> #
> pairwise.prop2.test(cbind(smokers, patients-smokers)) # test.function=chisq.test としても，同じ結果になる
Pairwise comparisons using prop.test
data: cbind(smokers, patients - smokers)
1 2 3
2 1.000 - -
3 1.000 1.000 -
4 0.119 0.093 0.124
P value adjustment method: holm
警告メッセージ：
1: In test.function(x[c(i, j), ], ...) :
Chi-squared approximation may be incorrect
2: In test.function(x[c(i, j), ], ...) :
Chi-squared approximation may be incorrect
3: In test.function(x[c(i, j), ], ...) :
Chi-squared approximation may be incorrect
> #
> # 検定関数として fisher.test を使う
> #
> pairwise.prop2.test(cbind(smokers, patients-smokers), test.function=fisher.test)
Pairwise comparisons using fisher.test
data: cbind(smokers, patients - smokers)
1 2 3
2 1.000 - -
3 1.000 1.000 -
4 0.075 0.075 0.096
P value adjustment method: holm
★ カテゴリーが 3 つ以上の場合（分布の差の検定）
> # 以下のような集計表 4 群それぞれ，4 つのカテゴリーに該当するケース数が集計されている
> # group Cat1 Cat2 Cat3 Cat4
> # 1 22 41 52 31
> # 2 12 31 25 14
> # 3 12 25 14 17
> # 4 23 14 16 28
> #
> # 行列の作り方は色々な方法がある（table 関数で集計したり，集計結果を matrix 関数で行列にするなど）
> #
> a <- rbind(c(22,41,52,31), c(12,31,25,14), c(12,25,14,17), c(23,14,16,28))
> #
> # 検定関数としては，chisq.test または fisher.test を使う
> #
> pairwise.prop2.test(a, test.function=chisq.test)
Pairwise comparisons using chisq.test
data: a
1 2 3
2 0.8626 - -
3 0.5112 0.8626 -
4 0.0086 0.0053 0.1600
P value adjustment method: holm
> pairwise.prop2.test(a, test.function=fisher.test)
Pairwise comparisons using fisher.test
data: a
1 2 3
2 0.8705 - -
3 0.4768 0.8705 -
4 0.0087 0.0051 0.1637
P value adjustment method: holm
直前のページへ戻る E-mail to Shigenobu AOKI
