wilcox_R: What the package does (short line)

bwrmcp <-
function(J,K,x,grp=NA,alpha=.05,bhop=TRUE){
#
# Do all pairwise comparisons of
# main effects for Factor A and B and all interactions
# using a rank-based method that tests for equal distributions.
#
#  A between by within subjects design is assumed.
#  Levels of Factor A are assumed to be independent and
#  levels of Factor B are dependent.
#
#  The data are assumed to be stored in x in list mode or in a matrix.
#  If grp is unspecified, it is assumed x[[1]] contains the data
#  for the first level of both factors: level 1,1.
#  x[[2]] is assumed to contain the data for level 1 of the
#  first factor and level 2 of the second factor: level 1,2
#  x[[j+1]] is the data for level 2,1, etc.
#  If the data are in wrong order, grp can be used to rearrange the
#  groups. For example, for a two by two design, grp<-c(2,4,3,1)
#  indicates that the second group corresponds to level 1,1;
#  group 4 corresponds to level 1,2; group 3 is level 2,1;
#  and group 1 is level 2,2.
#
#   Missing values are automatically removed.
#
 if(is.list(x))xrem=matl(x)
        JK <- J * K
        if(is.matrix(x)){
                xrem=x
                x <- listm(x)
}

        if(!is.na(grp[1])) {
                yy <- x
                x<-list()
                for(j in 1:length(grp))
                        x[[j]] <- yy[[grp[j]]]
        }
        if(!is.list(x))
                stop("Data must be stored in list mode or a matrix.")
#        for(j in 1:JK) {
#                xx <- x[[j]]
#                x[[j]] <- xx[!is.na(xx)] # Remove missing values
#        }
        #
if(JK != length(x))warning("The number of groups does not match the number of contrast coefficients.")
for(j in 1:JK){
temp<-x[[j]]
temp<-temp[!is.na(temp)] # Remove missing values.
x[[j]]<-temp
}
#
CC<-(J^2-J)/2
# Determine critical values
ncon<-CC*(K^2-K)/2
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
Fac.A<-matrix(0,CC,5)
dimnames(Fac.A)<-list(NULL,c("Level","Level","test.stat","p-value","sig.crit"))
mat<-matrix(c(1:JK),ncol=K,byrow=T)
ic<-0
for(j in 1:J){
for(jj in 1:J){
if(j < jj){
ic<-ic+1
Fac.A[ic,1]<-j
Fac.A[ic,2]<-jj
datsub=xrem[,c(mat[j,],mat[jj,])]
datsub=elimna(datsub)
#temp<-bwrank(2,K,elimna(x[,c(mat[j,],mat[jj,])]))
temp<-bwrank(2,K,datsub)
Fac.A[ic,3]<-temp$test.A
Fac.A[ic,4]<-temp$p.value.A
}}}
temp2<-order(0-Fac.A[,4])
Fac.A[temp2,5]<-dvec[1:length(temp2)]
CCB<-(K^2-K)/2
ic<-0
Fac.B<-matrix(0,CCB,5)
dimnames(Fac.B)<-list(NULL,c("Level","Level","test.stat","p-value","sig.crit"))
for(k in 1:K){
for(kk in 1:K){
if(k<kk){
ic<-ic+1
Fac.B[ic,1]<-k
Fac.B[ic,2]<-kk
mat1<-cbind(mat[,k],mat[,kk])
rv=c(mat1[,1],mat1[,2])
datsub=elimna(xrem[,sort(c(mat1[,1],mat1[,2]))])
temp<-bwrank(J,2,datsub)
Fac.B[ic,3]<-temp$test.B
Fac.B[ic,4]<-temp$p.value.B
}}}
temp2<-order(0-Fac.B[,4])
Fac.B[temp2,5]<-dvec[1:length(temp2)]
CCI<-CC*CCB
Fac.AB<-matrix(0,CCI,7)
dimnames(Fac.AB)<-list(NULL,c("Lev.A","Lev.A","Lev.B","Lev.B","test.stat","p-value","sig.crit"))
ic<-0
for(j in 1:J){
for(jj in 1:J){
if(j < jj){
for(k in 1:K){
for(kk in 1:K){
if(k<kk){
ic<-ic+1
Fac.AB[ic,1]<-j
Fac.AB[ic,2]<-jj
Fac.AB[ic,3]<-k
Fac.AB[ic,4]<-kk
val<-c(mat[j,k],mat[j,kk],mat[jj,k],mat[jj,kk])
val<-sort(val)
datsub=elimna(xrem[,val])
temp<-bwrank(2,2,datsub)
Fac.AB[ic,5]<-temp$test.AB
#Fac.AB[ic,6]<-temp$sig.AB
Fac.AB[ic,6]<-temp$p.value.AB
}}}}}}
temp2<-order(0-Fac.AB[,6])
Fac.AB[temp2,7]<-dvec[1:length(temp2)]
list(Factor.A=Fac.A,Factor.B=Fac.B,Factor.AB=Fac.AB)
}