wilcox_R: What the package does (short line)

pbad2way <-
function(J,K,x,est=onestep,conall=TRUE,alpha=.05,nboot=2000,grp=NA,
op=FALSE,pro.dis=FALSE,MM=FALSE,...){
#
# This function is like the function pbadepth,
# only it is assumed that main effects and interactions for a
# two-way design are to be tested.
#
        #   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.
        #
        JK <- J * K
        if(is.matrix(x))
                x <- listm(x)
        if(!is.na(grp[1])) {
                yy <- x
                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)]
        }
        #
        # Create the three contrast matrices
        #
        if(!conall){
        ij <- matrix(c(rep(1, J)), 1, J)
        ik <- matrix(c(rep(1, K)), 1, K)
        jm1 <- J - 1
        cj <- diag(1, jm1, J)
        for(i in 1:jm1)
                cj[i, i + 1] <- 0 - 1
        km1 <- K - 1
        ck <- diag(1, km1, K)
        for(i in 1:km1)
                ck[i, i + 1] <- 0 - 1
        conA <- t(kron(cj, ik))
        conB <- t(kron(ij, ck))
        conAB <- t(kron(cj, ck))
        conAB <- t(kron(abs(cj), ck))
}
if(conall){
temp<-con2way(J,K)
conA<-temp$conA
conB<-temp$conB
conAB<-temp$conAB
}
        ncon <- max(nrow(conA), nrow(conB), nrow(conAB))
        if(JK != length(x))
                warning("The number of groups does not match the number of contrast coefficients.")
if(!is.na(grp[1])){  # Only analyze specified groups.
xx<-list()
for(i in 1:length(grp))xx[[i]]<-x[[grp[i]]]
x<-xx
}
mvec<-NA
for(j in 1:JK){
temp<-x[[j]]
temp<-temp[!is.na(temp)] # Remove missing values.
x[[j]]<-temp
mvec[j]<-est(temp,...)
}
bvec<-matrix(NA,nrow=JK,ncol=nboot)
set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
print("Taking bootstrap samples. Please wait.")
for(j in 1:JK){
data<-matrix(sample(x[[j]],size=length(x[[j]])*nboot,replace=TRUE),nrow=nboot)
bvec[j,]<-apply(data,1,est,...) # J by nboot matrix, jth row contains
#                          bootstrapped  estimates for jth group
}
bconA<-t(conA)%*%bvec #C by nboot matrix
tvecA<-t(conA)%*%mvec
tvecA<-tvecA[,1]
tempcenA<-apply(bconA,1,mean)
veczA<-rep(0,ncol(conA))
bconA<-t(bconA)
smatA<-var(bconA-tempcenA+tvecA)
bconA<-rbind(bconA,veczA)
if(!pro.dis){
if(!op)dv<-mahalanobis(bconA,tvecA,smatA)
if(op){
dv<-out(bconA)$dis
}}
if(pro.dis)dv=pdis(bconA,MM=MM)
bplus<-nboot+1
sig.levelA<-1-sum(dv[bplus]>=dv[1:nboot])/nboot
bconB<-t(conB)%*%bvec #C by nboot matrix
tvecB<-t(conB)%*%mvec
tvecB<-tvecB[,1]
tempcenB<-apply(bconB,1,mean)
veczB<-rep(0,ncol(conB))
bconB<-t(bconB)
smatB<-var(bconB-tempcenB+tvecB)
bconB<-rbind(bconB,veczB)
if(!pro.dis){
if(!op)dv<-mahalanobis(bconB,tvecB,smatB)
if(op){
dv<-out(bconA)$dis
}}
if(pro.dis)dv=pdis(bconB,MM=MM)
sig.levelB<-1-sum(dv[bplus]>=dv[1:nboot])/nboot
bconAB<-t(conAB)%*%bvec #C by nboot matrix
tvecAB<-t(conAB)%*%mvec
tvecAB<-tvecAB[,1]
tempcenAB<-apply(bconAB,1,mean)
veczAB<-rep(0,ncol(conAB))
bconAB<-t(bconAB)
smatAB<-var(bconAB-tempcenAB+tvecAB)
bconAB<-rbind(bconAB,veczAB)
if(!pro.dis){
if(!op)dv<-mahalanobis(bconAB,tvecAB,smatAB)
if(op){
dv<-out(bconAB)$dis
}}
if(pro.dis)dv=pdis(bconAB,MM=MM)
sig.levelAB<-1-sum(dv[bplus]>=dv[1:nboot])/nboot
list(sig.levelA=sig.levelA,sig.levelB=sig.levelB,sig.levelAB=sig.levelAB,conA=conA,conB=conB,conAB=conAB)

}