一元配置分散分析（exact test）　　　　　Last modified: Nov 15, 2011

目的
　正確な P 値を計算する一元配置分散分析（正確確率検定；並べ替え検定 permutation test）である。データによっては計算量が多くなり実用的な時間内で計算が終了できないこともあるので，そのような場合にはモンテカルロ法による近似計算もできる。
　なお，一元配置分散分析は，各群の分散が等しいことを仮定しないものである（ウェルチの方法）。 
使用法

exact.onway.test(x, exact=TRUE, hybrid=FALSE, loop=10000)

引数

x      各群のデータをベクトル，それらをリスト要素とするリスト（使用例参照）
exact  正確な P 値を求める場合，またはシミュレーションにより近似的な P 値を求めるときには TRUE（デフォルト）
       FALSE の場合には標本に対して一元配置分散分析のみを行う
hybrid TRUE を指定すれば，シミュレーションにより近似的な P 値を計算する
loop   シミュレーションの回数

ソース



以下のプログラムは短いが時間がかかりすぎる。（注を参照）
インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/permutation-oneway-test.R", encoding="euc-jp")

#usage: permutation.oneway.test(list(c(36.7, 52.4, 65.8), c(45.7, 61.9, 65.3), c(52.6, 76.6, 81.3))))
permutation.oneway.test <- function(x)
{
        n <- sapply(x, length)
        g <- factor(rep(1:length(x), n))
        z <- unlist(x)
        library(e1071)
        perm <- permutations(sum(n))
        p0 <- oneway.test(z~g)$p.value
        ps <- apply(perm, 1, function(i) oneway.test(z[i]~g)$p.value)
        invisible(list(p0=p0, permutation.p=mean(ps < p0+1e-10), ps=ps))
}


使用例

exact.oneway.test(list(c(36.7, 52.4, 65.8), c(45.7, 61.9, 65.3), c(52.6, 76.6, 81.3)))

出力結果例

観察値による一元配置分散分析の P 値 = 0.440037
並べ替え検定による P 値 = 0.4107142857
査察した分割表の個数は 1680 個
n個の中に，p, q, r, ...個同じものがあるとき（n=p+q+r+…），これらを並べてできる順列のリストは，e1071 の permutations が速いので，以下のようにしてやれば，一応はできる。ただし，permutations の戻り値を求める順列に変換することと，重複を除くことに時間がかかるので，ベストではない。
comb2 <- function(n.i)
{
	n <- sum(n.i)
	x <- rep(1:length(n.i), n.i)
	library(e1071)
	a <- permutations(n)
	b <- unique(data.frame(t(apply(a, 1, function(i) x[i]))))
	rownames(b) <- 1:nrow(b)
	return(b)
}
これを実行すると以下のような結果が得られる。
> ans <- comb2(c(3,1,2)) # 順列の個数は (3+1+2)! / 3! / 1! / 2! = 60
> ans
   X1 X2 X3 X4 X5 X6
1   1  1  1  2  3  3
2   1  1  2  1  3  3
3   1  2  1  1  3  3
    途中省略
57  3  3  1  1  1  2
58  3  3  1  1  2  1
59  3  3  1  2  1  1
60  3  3  2  1  1  1

・ 直前のページへ戻る　　・ E-mail to Shigenobu AOKI

Made with Macintosh