多重比較の補正をRで

http://d.hatena.ne.jp/Tanakky/20081030
この前、話題にした多重比較時の補正をR言語で書いてみました。

#multiple experimental correction
#2008/11/04 1st version coded by Yoshiaki Tanaka

multi_correct = function(p_value, method, cutoff, ...){
	      k <- length(p_value)
	      if(method == "Bonferroni"){
			cutoff2 <- cutoff / k
			result <- which(p_value < cutoff2)
	      }else if(method == "Sidak"){
		        cutoff2 <- 1 - (1 - cutoff) ^ (1/k)
			result <- which(p_value < cutoff2)
	      }else if(method == "Holm"){
			p_value <- sort(p_value, decreasing=F)
			for(i in 1:k){
			      cutoff2 <- cutoff/(k+1-i)
			      if(p_value[i] > cutoff2){
					 if(i > 1){
					      result <- 1:(i-1)
					 }else{
					      result <- c()
					 }
					 break
			      }
			}		      
	      }else if(method == "BH" || method == "BY" || method == "AdaptiveBH"){
		       p_value <- sort(p_value, decreasing=F)

		       if(method == "BY"){
			     denomi <- 0
			     for(i in 1:k){
				   denomi <- denomi + 1/i 
			     } 
		       }

		       flag <- 0
		       for(i in k:1){
			     cutoff2 <- cutoff*i/k
			     if(method == "BY"){
				       cutoff2 <- cutoff2/denomi
			     }
			     if(p_value[i] < cutoff2){
					flag <- i
					break
			     }
			     #samp <- c(cutoff2, p_value[i])
			     #print(samp)
		       }
		       if(flag == 0){
		       	   result <- c()
		       }else{
			   result <- 1:flag
		       }
		       
		       if(method == "AdaptiveBH"){
			   if(is.null(result)==F){
				
				S <- numeric(k)
				for(i in 1:k){
				      S[i] = (1-p_value[i])/(k+1-i)
				      if(i >= 2){
					   if(S[i]-S[i-1] < 0){
						S_star <- S[i]
						break
					   }
				      }else if(i == k){
					    S_star <- S[i]
				      }
				}
				m <- min(c(floor(1/S_star + 1),k))

				result <- c()
				for(i in k:1){
				      cutoff2 <- cutoff*i/m
				      if(p_value[i] < cutoff2){
						    result <- 1:i
						    break
				      }
				}
			   }
		       }
	      }else{
		print("Bonferroni, Sidak, Holm, BH, AdaptiveBH and BY are available.")
		return(NULL)
	      }
	      if(length(result) == 0){
	        return(NULL)
	      }else{
	        #result <- p_value[result]
	        return(result)
	      }
}

一応この前とりあげた４つの方法に、AdaptiveBH法とBY法を加えました。
使い方は

> source("multi_correct.R")
> p_list <- c(0.33, 0.031,0.000001,0.00201,0.000061)
> multi_correct(p_list, "Bonferroni",0.05)
[1] 0.000001 0.002010 0.000061

もし間違いがあればご指摘くださいm(_ _)m