wilcox_R: What the package does (short line)

bwmarpb <-
function(J,K,x,est=hd,JK=J*K,grp=c(1:JK),nboot=599,
SEED=TRUE,na.rm=TRUE,...){
#
#  Multiple comparisons using
#  a percentile bootstrap for performing a split-plot design
#  using marginal measures of location
#  By default, 20% trimming is used with B=599 bootstrap samples.
#
#  The R variable x is assumed to contain the raw
#  data stored in list mode or in a matrix or a data frame
#  If in list mode, 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: level 1,2
#  x[[K]] is the data for level 1,K
#  x[[K+1]] is the data for level 2,1, x[[2K]] is level 2,K, etc.
#
#  If the data are in a matrix, column 1 is assumed to
#  correspond to x[[1]], column 2 to x[[2]], etc.
#
#  When in list mode x is assumed to have length JK, the total number of
#  groups being tested, but a subset of the data can be analyzed
#  using grp
#
if(is.matrix(x))
{
if(ncol(x)!=JK)print("WARNING: number of groups is not equal to JK")
}
if(is.list(x)){
if(length(x)!=JK)print("WARNING: number of groups is not equal to JK")
}
if(is.data.frame(x)){
if(ncol(x)!=JK)print("WARNING: number of groups is not equal to JK")
}
if(SEED)set.seed(2)
if(is.data.frame(x) || is.matrix(x)) {
y <- list()
ik=0
il=c(1:K)-K
for(j in 1:J){
il=il+K
zz=x[,il]
if(na.rm)zz=elimna(zz)
for(k in 1:K){
ik=ik+1
y[[ik]]=zz[,k]
}}
                x <- y
}

data<-list()
for(j in 1:length(x)){
data[[j]]<-x[[grp[j]]] # Now have the groups in proper order.
}
x<-data

#
con=con2way(J,K)
estA=psihat(x,est,con=con$conA,...)
estB=psihat(x,est,con=con$conB,...)
estAB=psihat(x,est,con=con$conAB,...)
set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
# Next determine the n_j values
nvec<-NA
jp<-1-K
for(j in 1:J){
jp<-jp+K
nvec[j]<-length(x[[j]])
}
blist<-list()
testmatA<-matrix(NA,ncol=ncol(con$conA),nrow=nboot)
testmatB<-matrix(NA,ncol=ncol(con$conB),nrow=nboot)
testmatAB<-matrix(NA,ncol=ncol(con$conAB),nrow=nboot)
for(iboot in 1:nboot){
iv<-0
for(j in 1:J){
temp<-sample(nvec[j],replace = T)
for(k in 1:K){
iv<-iv+1
tempx<-x[[iv]]
blist[[iv]]<-tempx[temp]
}}
#
# Now do all linear contrasts on bootstrap samples
testmatA[iboot,]<-psihat(blist,est,con=con$conA,...)
testmatB[iboot,]<-psihat(blist,est,con=con$conB,...)
testmatAB[iboot,]<-psihat(blist,est,con=con$conAB,...)
}
pbA=NA
pbB=NA
pbAB=NA
for(j in 1:ncol(con$conA))pbA[j]=mean(testmatA[,j]<0)+.5*mean(testmatA[,j]==0)
for(j in 1:ncol(con$conB))pbB[j]=mean(testmatB[,j]<0)+.5*mean(testmatB[,j]==0)
for(j in 1:ncol(con$conAB))pbAB[j]=mean(testmatAB[,j]<0)+.5*mean(testmatAB[,j]==0)
for(j in 1:ncol(con$conA))pbA[j]=2*min(c(pbA[j],1-pbA[j]))
for(j in 1:ncol(con$conB))pbB[j]=2*min(c(pbB[j],1-pbB[j]))
for(j in 1:ncol(con$conAB))pbAB[j]=2*min(c(pbAB[j],1-pbAB[j]))
p.valueA=pbA
p.valueB=pbB
p.valueAB=pbAB
pbA=cbind(estA,p.valueA)
pbB=cbind(estB,p.valueB)
pbAB=cbind(estAB,p.valueAB)
list(FacA=pbA,FacB=pbB,p.FacAB=pbAB,conA=con$conA,conB=con$conB,conAB=con$conAB)
}