tukey <- function(fn)
{
d <- read.table(fn, header=TRUE)
attach(d)
fac <- factor(f)
a <- length(unique(fac))
phi.e <- sum(n.i <- tapply(x,fac,length))-a
m.i <- tapply(x,fac,mean)
v.i <- tapply(x,fac,var)
v.e <- sum((n.i-1)*v.i)/phi.e
t <- abs(outer(m.i, m.i, "-"))/outer(n.i, n.i, function(x, y) { sqrt(v.e*(1/x+1/y)) })
p <- ptukey(t*sqrt(2), a, phi.e, lower=F)
list(groups=a, n=n.i, mean=m.i, variance=v.i, phi=phi.e, v=v.e, t=t, p=p)
}
ありがたいことに，R にはスチューデント化した範囲の関数があります。
でもって，R だと非常に楽です。
p 値はこれでいいのでしょうね。
使用例の tykey.data は以下の通り。
f x
1 14
1 15
1 14
1 16
1 15
1 17
1 17
2 17
2 16
2 17
2 16
2 15
2 18
2 19
2 15
3 18
3 19
3 20
3 19
3 17
3 17
4 20
4 21
4 19
4 20
4 19
4 22
4 20
5 19
5 20
5 19
5 17
5 17
5 17
5 18
tukey("tukey.data")
結果のうち一番重要な p は以下のようになります。
$p
1 2 3 4 5
1 1.000000e+00 0.3630247954 0.001975521 6.970838e-07 0.002687670
2 3.630248e-01 1.0000000000 0.109056039 6.039780e-05 0.156759130
3 1.975521e-03 0.1090560392 1.000000000 9.414444e-02 0.998679526
4 6.970838e-07 0.0000603978 0.094144436 1.000000e+00 0.039780546
5 2.687670e-03 0.1567591303 0.998679526 3.978055e-02 1.000000000