R/r1mcp.R

r1mcp <-
function(x,alpha=.05,bhop=FALSE){
#
# Do all pairwise comparisons using a modification of
# the Brunner, Dette and Munk (1997) rank-based method.
# FWE is controlled using Rom's technique.
#
#  Setting bhop=T, FWE is controlled using the
#  Benjamini-Hochberg Method.
#
#  The data are assumed to be stored in x in list mode or in a matrix.
#
#   Missing values are automatically removed.
#
        if(is.matrix(x))x <- listm(x)
        if(!is.list(x))
                stop("Data must be stored in list mode or a matrix.")
J<-length(x)
        for(j in 1:J) {
                xx <- x[[j]]
                x[[j]] <- xx[!is.na(xx)] # Remove missing values
        }
#
CC<-(J^2-J)/2
# Determine critical values
ncon<-CC
if(!bhop){
if(alpha==.05){
dvec<-c(.05,.025,.0169,.0127,.0102,.00851,.0073,.00639,.00568,.00511)
if(ncon > 10){
avec<-.05/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha==.01){
dvec<-c(.01,.005,.00334,.00251,.00201,.00167,.00143,.00126,.00112,.00101)
if(ncon > 10){
avec<-.01/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha != .05 && alpha != .01)dvec<-alpha/c(1:ncon)
}
if(bhop)dvec<-(ncon-c(1:ncon)+1)*alpha/ncon
output<-matrix(0,CC,5)
dimnames(output)<-list(NULL,c("Level","Level","test.stat","p.value","p.crit"))
ic<-0
for(j in 1:J){
for(jj in 1:J){
if(j < jj){
ic<-ic+1
output[ic,1]<-j
output[ic,2]<-jj
temp<-bdm(x[c(j,jj)])
output[ic,3]<-temp$output$F
output[ic,4]<-temp$output$sig
}}}
temp2<-order(0-output[,4])
output[temp2,5]<-dvec[1:length(temp2)]
list(output=output)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.