「統計学関連なんでもあり」の過去ログ---027

★ ダネット検定 ★

1807.　ダネット検定　長崎　2004/01/15 (木) 15:18
└1810.　Re: ダネット検定　青木繁伸　2004/01/15 (木) 17:06
　└1813.　Re^2: ダネット検定　長崎　2004/01/15 (木) 17:24
　　├1816.　Re^3: ダネット検定　青木繁伸　2004/01/15 (木) 19:43
　　└1814.　Re^3: ダネット検定　俊坊　2004/01/15 (木) 17:48
　　　└1827.　Re^4: ダネット検定　青木繁伸　2004/01/16 (金) 13:32
　　　　├1861.　Re^5: ダネット検定　中澤　2004/01/17 (土) 23:22
　　　　│└1862.　Re^6: ダネット検定　青木繁伸　2004/01/17 (土) 23:44
　　　　│　└1867.　Re^7: ダネット検定　中澤　2004/01/18 (日) 17:58
　　　　└1830.　Re^5: ダネット検定　俊坊　2004/01/16 (金) 14:02
　　　　　├1833.　Re^6: ダネット検定　俊坊　2004/01/16 (金) 14:06
　　　　　└1832.　Re^6: ダネット検定　青木繁伸　2004/01/16 (金) 14:05
　　　　　　└1841.　Re^7: ダネット検定　青木繁伸　2004/01/16 (金) 15:08

1807.　ダネット検定　長崎　　2004/01/15 (木) 15:18

多重比較による検定において，有意水準0.001を求めたいのですが，一般的に市販されているソフトではほとんど0.01までで，0.001には対応してません。Tukey法に関しては，ネットでソフトを見つけた（スチューデント化された範囲を求めるもの）のですが，ダネットのほうはいまだにわかりません。何か求めるよい方法（簡単に求められるもの）はありませんか？もともと求められるものなんですか？？？

　　　　　[このページのトップへ]

1810.　Re: ダネット検定　青木繁伸　　2004/01/15 (木) 17:06

> 何か求めるよい方法（簡単に求められるもの）はありませんか？もともと求められるものなんですか？？？

個人的に使うなら以下を参照
http://www.r-project.org/nocvs/mail/r-help/2001/5714.html

一般に使ってもらえるように R なんかで書き直してもよいかもしれないが，著作権関係が面倒かな？

　　　　　[このページのトップへ]

1813.　Re^2: ダネット検定　長崎　　2004/01/15 (木) 17:24

> 個人的に使うなら以下を参照
> http://www.r-project.org/nocvs/mail/r-help/2001/5714.html
>
> 一般に使ってもらえるように R なんかで書き直してもよいかもしれないが，著作権関係が面

早速のご返事ありがとうございました。
正直なところ，FORTRANやらRやらまったくわかりません。
理論とかわからなくても使えるような（数値入力だけとか）ソフトはないんでしょうか？？
やっぱり，こういったところも著作権に引っかかるのかな

　　　　　[このページのトップへ]

1816.　Re^3: ダネット検定　青木繁伸　　2004/01/15 (木) 19:43

> 理論とかわからなくても使えるような（数値入力だけとか）ソフトはないんでしょうか？？

ちょっと調べたら，R にちゃんとありましたね。
mvtnorm をインストールして，pmvt 関数を使います。
そのままではちょっと使いにくいので，以下のように糖衣を着せて。

R はただですからね。

p <- function(a,df,x,r) { corr <- diag(a-1) corr[lower.tri(corr)] <- r delta <- prob <- pmvt(lower=-x, upper=x, delta=rep(0, a-1), df=df, corr=corr) print(prob) } p(5,17,2.604,0.7) > library(mvtnorm) > p <- function(a,df,x,r) + { + corr <- diag(a-1) + corr[lower.tri(corr)] <- r + delta <- + prob <- pmvt(lower=-x, upper=x, delta=rep(0, a-1), df=df, corr=corr) + print(prob) + } > p(5,17,2.604,0.7) [1] 0.9499043 attr(,"error") [1] 0.0007555838 attr(,"msg") [1] "Normal Completion"

　　　　　[このページのトップへ]

1814.　Re^3: ダネット検定　俊坊　　2004/01/15 (木) 17:48

p値を求める，ということなら，YUKMS社製のStat Light #4がやってくれます(#4は多群検定に特化)。最近バージョンアップしました。

　　　　　[このページのトップへ]

1827.　Re^4: ダネット検定　青木繁伸　　2004/01/16 (金) 13:32

もしYUKMS社製のStat Light #4をお持ちでしたら，申し訳ありませんが，お願いがあります。
以下のデータをダネットの両側検定をしたとき，結果（P値）が正しいか教えていただけませんでしょうか。R でプログラムを書いて，一応正しいと思うのですが，確認したく思っています。
A1 がコントロールです。データは縦方向が各群のデータです。
A1 は10ケース，他の3群は8ケースずつです。

A1 A2 A3 A4 7 11.59 11 13 9 9 12 12 8 10 12 12 6 9 10 11 9 9 11 14 8 9 13 12 11 10 9 11 10 12 10 10 8 8 結果 rho=0.44444 group:2 t=2.48520 p=0.04966 group:3 t=4.17208 p=0.00065 group:4 t=5.57615 p=0.00001

　　　　　[このページのトップへ]

1861.　Re^5: ダネット検定　中澤　　2004/01/17 (土) 23:22

Rのmultcompパッケージのsimtest()でtype="Dunnett"で計算させた結果と比べると，t値の絶対値は同じですが，pは違う値でした。
library(multcomp) y <- c(7,11.59,11,13,9,9,12,12,8,10,12,12,6,9,10,11,9,9,11,14,8,9,13,12,11,10,9,11,10,12,10,10,8,8) g <- factor(c(rep(c('A','B','C','D'),8),'A','A')) res <- simtest(y ~ g, type="Dunnett") summary(res) の結果， Simultaneous tests: Dunnett contrasts Call: simtest.formula(formula = y ~ g, type = "Dunnett") Dunnett contrasts for factor g Contrast matrix: gA gB gC gD gB-gA 0 -1 1 0 0 gC-gA 0 -1 0 1 0 gD-gA 0 -1 0 0 1 Absolute Error Tolerance: 0.001 Coefficients: Estimate t value Std.Err. p raw p Bonf p adj gD-gA 3.475 -5.576 0.623 0.000 0.000 0.000 gC-gA 2.600 -4.172 0.623 0.000 0.000 0.000 gB-gA 1.549 -2.485 0.623 0.019 0.019 0.019
となりました。

　　　　　[このページのトップへ]

1862.　Re^6: ダネット検定　青木繁伸　　2004/01/17 (土) 23:44

> Rのmultcompパッケージのsimtest()でtype="Dunnett"で計算させた結果と比べると，t値の絶対値は同じですが，pは違う値でした。

う。なんでだろ。
あのテストデータは，例の本の例解を変更して，ぎりぎり5%有意のデータにして，P値を確かめたのだけど(^_^;)

　　　　　[このページのトップへ]

1867.　Re^7: ダネット検定　中澤　　2004/01/18 (日) 17:58

> う。なんでだろ。
> あのテストデータは，例の本の例解を変更して，ぎりぎり5%有意のデータにして，P値を確かめたのだけど(^_^;)

よくみるとsimtest()の結果で，p rawとp Bonfとp adjとして同じ値が表示されているようなので，simtest()のバグかもしれません。

　　　　　[このページのトップへ]

1830.　Re^5: ダネット検定　俊坊　　2004/01/16 (金) 14:02

> rho=0.44444 > group:2 t=2.48520 p=0.04966 > group:3 t=4.17208 p=0.00065 > group:4 t=5.57615 p=0.00001 1 vs 2 p=0.0500 1 vs 3 p=0.000689117 1 vs 4 p=1.48972E-05
となりました。
rhoは表示されません。
以前Yukms社に問い合わせたところ，自由度が異なる場合の吉田らの方法で計算されているそうです。

　　　　　[このページのトップへ]

1833.　Re^6: ダネット検定　俊坊　　2004/01/16 (金) 14:06

ただし，Stat Lightはアドインソフトですので，計算しているのはビル・ゲイツの手下です。念のため。

　　　　　[このページのトップへ]

1832.　Re^6: ダネット検定　青木繁伸　　2004/01/16 (金) 14:05

> 以前Yukms社に問い合わせたところ，自由度が異なる場合の吉田らの方法で計算されているそうです。

早速，ありがとうございました。
この後の記事で，R で組んだダネット検定をアップします。

　　　　　[このページのトップへ]

1841.　Re^7: ダネット検定　青木繁伸　　2004/01/16 (金) 15:08

library(mvtnorm) Dunnett <- function(dat) { printf <- function(fmt, ...) cat(sprintf(fmt, ...)) get.rho <- function(ni) { k <- length(ni) rho <- outer(ni, ni, function(x, y) { sqrt(x/(x+ni[1])*y/(y+ni[1])) }) diag(rho) <- 0 sum(rho[-1,-1])/(k-2)/(k-1) } pdunnett <- function(a,df,x,r) { corr <- diag(a-1) corr[lower.tri(corr)] <- r prob <- 1-pmvt(lower=-x, upper=x, delta=rep(0, a-1), df=df, corr=corr, abseps=0.0001) } x <- dat[,2] ni <- table(g <- dat[,1]) gv <- as.numeric(names(ni)) a <- length(ni) n <- sum(ni) mi <- vi <- rep(0, a) for (i in 1:a) { mi[i] <- mean(xi <- x[g==gv[i]]) vi[i] <- var(xi) } Vw <- sum(vi*(ni-1))/(n-a) rho <- get.rho(ni) printf("rho=%5.3f\n", rho) for (i in 2:a) { ti <- abs(mi[i]-mi[1])/sqrt(Vw*(1/ni[i]+1/ni[1])) pi <- pdunnett(a, n-a, ti, rho) printf("group:%i t=%.3f p=%.3f\n", i, ti, pi[1]) } } Dunnett(dat) データは， dat <- matrix(c( 1,3.42, 1,3.84, 1,3.96, 1,3.76, 2,3.17, 2,3.63, 2,3.47, 2,3.44, 2,3.39, 3,3.64, 3,3.72, 3,3.91 ), ncol=2, byrow=TRUE)
のような感じで用意する。

　　　　　[このページのトップへ]