ペリの方法による多重比較 Last modified: Mar 15, 2007
目的
ペリの方法による多重比較を行う。
使用法
Peritz(data, group)
引数
data 観察値ベクトル
group 群変数ベクトル
statistcs 使用する検定統計量の種類 "F"（デフォルト） か "Q"
ソース
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/Peritz.R", encoding="euc-jp")
# ペリの方法
Peritz <- function( data, # データベクトル
group, # 群別ベクトル
alpha=0.05, # 有意水準
statistics=c("F", "Q")) # 検定統計量の種類
{
# 書式付き print 関数
printf <- function(fmt, ...) # 書式と項目
{
cat(sprintf(fmt, ...))
}
# 一つずつの検定結果
judge <- function( P, # 検定対象となる群番号を表す数字の集合
p, # 検定対象となる群の個数（P の要素数）
stat.P) # 検定統計量（Q.P または F.P）
{
judge2 <- function(P, ns.list, stat.P, c.p) # 下請け関数
{
if (!is.null(ns.list) && # 保留された仮説のリストがあって，
any(mapply(function(tbl) # P がそれらの仮説の部分集合であるときは，
all(is.element(P, tbl)), ns.list))) {
return("NSI") # この仮説を誘導する仮説が保留されているのでこの帰無仮説も保留する
}
else if (stat.P >= c.p) { # 検定統計量が棄却限界値より大きければ，
return("S") # 帰無仮説を棄却する
}
else { # 棄却限界値未満ならば，
return("NS") # 帰無仮説を保留する
}
}
res <- judge2(P, NS.list, stat.P, c.p[p]) # まずは，テューキー・ウェルシュの方法により判定
if (res == "NS") { # 帰無仮説を保留するときは，
NS.list <<- append(NS.list, list(P)) # 群番号の集合を保留された仮説のリストに付け加える
}
res.NK <- judge2(P, NS.list.NK, stat.P, c.p.NK[p]) # 次に，ニューマン・コイルスの方法により判定
if (res.NK == "NS") { # 帰無仮説を保留するときは，
NS.list.NK <<- append(NS.list.NK, list(P)) # 群番号の集合を保留された仮説のリストに付け加える
}
n.res <<- n.res+1 # 結果を蓄積する
result[n.res,] <<- c(paste(P, collapse=","), p, stat.P, alpha.p[p], c.p[p], res, alpha, c.p.NK[p], res.NK)
}
# Q 統計量を使う検定関数
funcQ <- function(P) # 群番号を表す数字の集合
{
p <- length(P) # 群番号を表す数字の集合 P の要素数
pair <- combn(p, 2, function(i) P[i]) # P から 2 つの群を取り出す組み合わせ
Q.P <- max(apply(pair, 2,
function(ij) QP[ij[2], ij[1]])) # 検定統計量（事前に Q.P として計算してある）
judge(P, p, Q.P) # 検定結果の評価
}
# F 統計量を使う検定関数
funcF <- function(P) # 群番号を表す数字の集合
{
p <- length(P) # 群番号を表す数字の集合 P の要素数
np <- n[P] # P で表される群のサンプルサイズ
mp <- m[P] # P で表される群の平均値
nt <- sum(np) # 全体のサンプルサイズ
mean.p <- sum(np*mp)/nt # 全体の平均値
vb.p <- sum(np*(mp-mean.p)^2) # 群間分散
F.P <- vb.p/(p-1)/vw # F 統計量
judge(P, p, F.P) # 検定結果の評価
}
# テューキー・ウェルシュとニューマン・コイルスの組み合わせ
proc11 <- function(i) # 食い違う結果があるデータフレームの行番号
{
P <- as.integer(unlist(strsplit(result[i,1], split=","))) # P
p <- length(P)
if (!is.null(NS.list.Peritz) && # 保留された仮説のリストがあって，
any(mapply(function(tbl) # P がそれらの仮説の部分集合であるときは，
all(is.element(P, tbl)), NS.list.Peritz))) {
return("NS") # この仮説を誘導する仮説が保留されているのでこの帰無仮説も保留する
}
else if (result[i, 3] >= result[i, 5]) { # ZP ≧ cp なら帰無仮説を棄却する
return("S")
}
else {
comprimentP <- (1:a)[-P] # P の補集合
m <- length(comprimentP) # P の補集合の要素数
sufix <- NULL # チェックすべき結果のある行番号
for (i in m:2) {
temp <- combn(m, i, # P の補集合の全ての部分集合について，
function(arg) paste(comprimentP[arg], collapse=","))
sufix <- c(sufix, sapply(temp, # 文字列に一致する行番号を探す
function(temp) which(result[,1]==temp)))
}
temp <- result[sufix, 3] >= result[sufix, 5] # ZP ≧ cp か？
if (all(temp)) { # 全てが真なら，
return("S") # 帰無仮説を棄却する
}
}
NS.list.Peritz <<- append(NS.list.Peritz, list(P)) # 群番号の集合を保留された仮説のリストに付け加える
return("NS")
}
# 結果を出力する
print.out <- function(result) # データフレームに納めた結果
{
for (i in 1:nrow(result)) {
strP <- paste("{", result[i,1], "}", sep="")
if (i==1 || result[i, 2] != result[i-1, 2]) { # p が同じである P の最初の検定結果
printf(format1, strP, result[i,2], result[i,3], result[i,4], result[i,5],
result[i,6], result[i,7], result[i,8], result[i,9], result[i,10])
}
else { # 二番目以降の検定結果
printf(format2, strP, result[i,3], result[i,6], result[i,9], result[i,10])
}
}
result2 <- result[result[,2] == 2,c(1,10)] # 二群比較の部分
result3 <- matrix("--- ", a, a) # 結果表
sapply(1:nrow(result2),
function(i) {
temp <- ifelse(result2[i,2] == "S", "* ", "n.s.")
ij <- as.integer(unlist(strsplit(result2[i,1], split=",")))
result3[ij[1], ij[2]] <<- result3[ij[2], ij[1]] <<- temp
})
colnames(result3) <- rownames(result3) <- paste("A", 1:a, sep="")
printf("\n判定結果\n")
print(result3, quote=FALSE)
}
# 関数本体
statistics <- match.arg(statistics) # 引数の補完
OK <- complete.cases(data, group) # 欠損値を持つケースを除く
data <- data[OK]
group <- group[OK]
n <- tapply(data, group, length)
m <- tapply(data, group, mean)
s <- tapply(data, group, sd)
a <- length(n) # 群の数
nt <- sum(n) # 全体の標本サイズ
dfw <- nt-a # 群内平方和の自由度
vw <- sum((n-1)*s^2)/dfw # 群内分散
NS.list <- NULL # 検定結果を保留した仮説のリスト
NS.list.NK <- NULL # 検定結果を保留した仮説のリスト
NS.list.Peritz <- NULL # 検定結果を保留した仮説のリスト
result <- matrix("", 2^a-a-1, 9)
n.res <- 0
format0 <- sprintf("%%-%is %%s %%6s %%7s %%7s %%-4s %%7s %%7s %%-4s %%-6s\n", 2*a+1) # 結果の出力書式0
format1 <- sprintf("%%-%is %%i %%6.3f %%7.4f %%7.3f %%-4s %%7.4f %%7.3f %%-4s %%-6s\n", 2*a+1) # 結果の出力書式1
format2 <- sprintf("%%-%is %%6.3f %%-4s %%-4s %%-6s\n", 2*a+1) # 結果の出力書式2
alpha.p <- c(NA, 1-(1-alpha)^(2:(a-2)/a), alpha, alpha) # αp を計算
if (a == 3) alpha.p <- c(NA, alpha, alpha)
if (statistics == "Q") {
QP <- matrix(NA, a, a) # 二群の検定統計量を前もって計算しておく領域を確保
QP[lower.tri(QP)] <- combn(a, 2, function(ij) {
i <- ij[1]
j <- ij[2]
abs(m[i]-m[j])/sqrt(vw*(1/n[i]+1/n[j])) # 検定統計量
})
qTukey <- qtukey(alpha.p[1:a], 1:a, dfw, lower.tail=FALSE)
c.p <- c(NA, qTukey[2]/sqrt(2)) # cp を計算
sapply(3:a, function(i) c.p[i] <<- max(qTukey[i]/sqrt(2), c.p[i-1]))
qTukey.NK <- c(NA, qtukey(alpha, 2:a, dfw, lower.tail=FALSE))
c.p.NK <- c(NA, qTukey.NK[2]/sqrt(2)) # cp.NK を計算
sapply(3:a, function(i) c.p.NK[i] <<- max(qTukey.NK[i]/sqrt(2), c.p.NK[i-1]))
printf(format0, "P", "p", "QP", "αp", "cp(TW)", "判定", "αp", "cp(NK)", "判定", "ペリの判定")
sapply(a:2, function(i) combn(a, i, funcQ)) # a 群から a〜2 群を取り出す組み合わせについて検定を実行
}
else { # if (statistics == "F") {
c.p <- qf(alpha.p, 0:(a-1), dfw, lower.tail=FALSE)
c.p.NK <- c(NA, qf(alpha, 1:(a-1), dfw, lower.tail=FALSE))
printf(format0, "P", "p", "FP", "αp", "cp(TW)", "判定", "αp", "cp(NK)", "判定", "ペリの判定")
sapply(a:2, function(i) combn(a, i, funcF)) # a 群から a〜2 群を取り出す組み合わせについて検定を実行
}
result <- data.frame(result)
old <- options(width=256)
con <- textConnection(capture.output(result))
result <- read.table(con)
close(con)
options(old)
result[,1] <- as.character(result[,1])
result[,6] <- as.character(result[,6])
result[,9] <- as.character(result[,9])
result[,10] <- ifelse((result[,6] == "NS" | result[,6] == "NSI") & result[,9] == "S", "不一致", ifelse(result[,6] == "NSI", "NS", result[,6]))
colnames(result) <- c("P", "p", "ZP", "αp", "cp.TW", "判定.TW", "αp.NK", "cp.NK", "判定.NK", "判定.Peritz")
sapply(1:nrow(result), function(i) {
if (result[i, 10] == "不一致") {
result[i, 10] <<- proc11(i)
}
})
print.out(result)
invisible(result)
}
使用例
data <- c(
14, 15, 14, 16, 15, 17, 17, # 第 1 群のデータ，7 例
17, 16, 17, 16, 15, 18, 19, 15, # 第 2 群のデータ，8 例
18, 19, 20, 19, 17, 17, # 第 3 群のデータ，6 例
20, 21, 19, 20, 19, 22, 20, # 第 4 群のデータ，7 例
19, 20, 19, 17, 17, 17, 18 # 第 5 群のデータ，7 例
)
group <- rep(1:5, c(7, 8, 6, 7, 7)) # 群を表す数値
Peritz(data, group)
Peritz(data, group, statistics="Q")
出力結果例
> data <- c(14, 15, 14, 16, 15, 17, 17,
+ 17, 16, 17, 16, 15, 18, 19, 15,
+ 18, 19, 20, 19, 17, 17,
+ 20, 21, 19, 20, 19, 22, 20,
+ 19, 20, 19, 17, 17, 17, 18)
> group <-rep(1:5, c(7,8,6,7,7))
> Peritz(data, group)
P p FP αp cp(TW) 判定 αp cp(NK) 判定 ペリの判定
{1,2,3,4,5} 5 14.619 0.0500 2.690 S 0.0500 2.690 S S
{1,2,3,4} 4 19.100 0.0500 2.922 S 0.0500 2.922 S S
{1,2,3,5} 8.165 S S S
{1,2,4,5} 18.841 S S S
{1,3,4,5} 16.990 S S S
{2,3,4,5} 9.934 S S S
{1,2,3} 3 8.801 0.0303 3.937 S 0.0500 3.316 S S
{1,2,4} 27.225 S S S
{1,2,5} 8.336 S S S
{1,3,4} 25.424 S S S
{1,3,5} 11.529 S S S
{1,4,5} 25.210 S S S
{2,3,4} 14.866 S S S
{2,3,5} 4.145 S S S
{2,4,5} 14.883 S S S
{3,4,5} 5.388 S S S
{1,2} 2 3.438 0.0203 6.006 NS 0.0500 4.171 NS NS
{1,3} 17.536 S S S
{1,4} 50.037 S S S
{1,5} 16.587 S S S
{2,3} 6.437 S S S
{2,4} 29.720 S S S
{2,5} 5.533 NS S S
{3,4} 6.805 S S S
{3,5} 0.075 NS NS NS
{4,5} 9.006 S S S
判定結果
A1 A2 A3 A4 A5
A1 --- n.s. * * *
A2 n.s. --- * * *
A3 * * --- * n.s.
A4 * * * --- *
A5 * * n.s. * ---
> Peritz(data, group, statistics="Q")
P p QP αp cp(TW) 判定 αp cp(NK) 判定 ペリの判定
{1,2,3,4,5} 5 7.074 0.0500 2.901 S 0.0500 2.901 S S
{1,2,3,4} 4 7.074 0.0500 2.719 S 0.0500 2.719 S S
{1,2,3,5} 4.188 S S S
{1,2,4,5} 7.074 S S S
{1,3,4,5} 7.074 S S S
{2,3,4,5} 5.452 S S S
{1,2,3} 3 4.188 0.0303 2.688 S 0.0500 2.465 S S
{1,2,4} 7.074 S S S
{1,2,5} 4.073 S S S
{1,3,4} 7.074 S S S
{1,3,5} 4.188 S S S
{1,4,5} 7.074 S S S
{2,3,4} 5.452 S S S
{2,3,5} 2.537 NS S S
{2,4,5} 5.452 S S S
{3,4,5} 3.001 S S S
{1,2} 2 1.854 0.0203 2.451 NS 0.0500 2.042 NS NS
{1,3} 4.188 S S S
{1,4} 7.074 S S S
{1,5} 4.073 S S S
{2,3} 2.537 NSI S S
{2,4} 5.452 S S S
{2,5} 2.352 NSI S S
{3,4} 2.609 S S S
{3,5} 0.275 NSI NS NS
{4,5} 3.001 S S S
判定結果
A1 A2 A3 A4 A5
A1 --- n.s. * * *
A2 n.s. --- * * *
A3 * * --- * n.s.
A4 * * * --- *
A5 * * n.s. * ---
直前のページへ戻る E-mail to Shigenobu AOKI
