��Ū ��Ψ�κ���¿����ӡ�Ryan ˡ��Tukey ˡ����Ԥ� ����ˡ p.multi.comp(n, r, alpha=0.05, method=c("ryan", "tukey")) ���� n �Ʒ������ r �Ʒ��������� alpha ͭ�տ�� method ��ˡ�ʾ�ά���� Ryan ˡ�� ������ ���ȡ���ϡ��ʲ��� 1 �Ԥԡ�����R ������˥ڡ����Ȥ��� source("http://aoki2.si.gunma-u.ac.jp/R/src/p_multi_comp.R", encoding="euc-jp") # �饤�������ˡ�ȥ��塼��������ˡ�ˤ����Ψ������� p.multi.comp <- function( n, # ɸ�ܥ����� r, # ���������� alpha=0.05, # ͭ�տ�� method=c("ryan", "tukey")) # ��ˡ { printf <- function(fmt, ...) { cat(sprintf(fmt, ...)) } check <- function(s, b) # ���ꤷ�褦�Ȥ��Ƥ�����������ޤǤ�ͭ�դǤʤ��Ȥ��줿�˶��ޤ�Ƥ��뤫 { if (ns.n > 1) { for (i in 1:ns.n) { if (ns.s[i] <= s && s <= ns.b[i] && ns.s[i] <= b && b <= ns.b[i]) { return(FALSE) # ���ꤹ��ޤǤ�ʤ�ͭ�դǤʤ��Ȥ��� } } } return(TRUE) # ���ꤷ�ʤ��ƤϤʤ�ʤ� } k <- length(n) # ���ο� stopifnot(k == length(r), k == length(n), n > 0, r >= 0, r <= n, floor(n) == n, floor(r) == r) method <- match.arg(method) # �������䴰 o <- order(r/n) # ɸ����Ψ���礭���ν�� sr <- r[o] sn <- n[o] sr.sn <- sr/sn # ��� num.significant <- ns.n <- 0 ns.s <- ns.b <- numeric(k*(k-1)/2) # ͭ�դǤʤ����ε�Ͽ�� for (m in k:2) { for (small in 1:(k-m+1)) { big <- small+m-1 if (check(small, big)) { prop <- sum(sr[small:big]) / sum(sn[small:big]) # ��������Ψ se <- sqrt(prop*(1-prop)*(1/sn[small]+1/sn[big])) # ɸ����� diff <- sr.sn[big]-sr.sn[small] # ��Ψ�κ� if (method == "ryan") { # Ryan ����ˡ nominal.alpha <- 2*alpha/(k*(m-1)) # ̾��Ūͭ�տ�� rd <- se*qnorm(nominal.alpha/2, lower.tail=FALSE) # ��������ȸ��ʤ��뺹���礭�� z <- if (se == 0) 0 else diff/se p <- pnorm(z, lower.tail=FALSE)*2 result <- p <= nominal.alpha } else { # Tukey ����ˡ if (m == k) { # ���緲�ȺǾ�������ӤΤȤ� qq <- q <- qtukey(alpha, k, Inf, lower.tail=FALSE) } else { # �嵭�ʳ��η�����ӤΤȤ� qq <- (q+qtukey(alpha, m, Inf, lower.tail=FALSE))/2 } WSD <- se*qq/sqrt(2) result <- diff != 0 && diff >= WSD } if (result) { # ͭ�դǤ���Ȥ� num.significant <- 1 printf("p[%2i]=%7.5f vs. p[%2i]=%7.5f : diff.= %7.5f, ", o[small], sr.sn[small], o[big], sr.sn[big], diff) if (method == "ryan") { printf("RD=%7.5f : P=%7.5f, alpha'=%7.5f\n", rd, p, nominal.alpha) } else { printf("WSD=%7.5f\n", WSD) } } else { # ͭ�դǤʤ��Ȥ� ns.n <- ns.n+1 ns.s[ns.n] <- small ns.b[ns.n] <- big } } } } if (num.significant == 0) { # ͭ�պ��Τ��뷲�ϰ�Ĥ�ʤ��ä� print("Not significant at all.") } } ������ p.multi.comp(c(30, 35, 47, 21, 45), c(2, 4, 14, 13, 39), method="ryan") p.multi.comp(c(30, 35, 47, 21, 45), c(2, 4, 14, 13, 39), method="tukey") ���Ϸ���� > p.multi.comp(c(30, 35, 47, 21, 45), c(2, 4, 14, 13, 39), method="ryan") p[ 1]=0.06667 vs. p[ 5]=0.86667 : diff.= 0.80000, RD=0.32472 : P=0.00000, alpha'=0.00500 p[ 1]=0.06667 vs. p[ 4]=0.61905 : diff.= 0.55238, RD=0.33341 : P=0.00001, alpha'=0.00667 p[ 2]=0.11429 vs. p[ 5]=0.86667 : diff.= 0.75238, RD=0.30528 : P=0.00000, alpha'=0.00667 p[ 1]=0.06667 vs. p[ 3]=0.29787 : diff.= 0.23121, RD=0.23054 : P=0.00979, alpha'=0.01000 p[ 2]=0.11429 vs. p[ 4]=0.61905 : diff.= 0.50476, RD=0.32612 : P=0.00007, alpha'=0.01000 p[ 3]=0.29787 vs. p[ 5]=0.86667 : diff.= 0.56879, RD=0.26479 : P=0.00000, alpha'=0.01000 p[ 3]=0.29787 vs. p[ 4]=0.61905 : diff.= 0.32118, RD=0.29877 : P=0.01239, alpha'=0.02000 > p.multi.comp(c(30, 35, 47, 21, 45), c(2, 4, 14, 13, 39), method="tukey") p[ 1]=0.06667 vs. p[ 5]=0.86667 : diff.= 0.80000, WSD=0.31555 p[ 1]=0.06667 vs. p[ 4]=0.61905 : diff.= 0.55238, WSD=0.32546 p[ 2]=0.11429 vs. p[ 5]=0.86667 : diff.= 0.75238, WSD=0.29800 p[ 1]=0.06667 vs. p[ 3]=0.29787 : diff.= 0.23121, WSD=0.22695 p[ 2]=0.11429 vs. p[ 4]=0.61905 : diff.= 0.50476, WSD=0.32104 p[ 3]=0.29787 vs. p[ 5]=0.86667 : diff.= 0.56879, WSD=0.26067 p[ 3]=0.29787 vs. p[ 4]=0.61905 : diff.= 0.32118, WSD=0.30102����ڡ���