wilcox_R: What the package does (short line)

m2way <-
function(J,K,x,est=hd,alpha=.05,nboot=600,SEED=TRUE,grp=NA,pr=TRUE,...){
#
# Two-way ANOVA based on forming averages
#
#  By default
#  est=hd meaning that medians are used with the Harrell-Davis estimator.
#
#   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.data.frame(x))x=as.matrix(x)
xcen<-list()
        if(is.matrix(x))
                x <- listm(x)
        if(!is.list(x))
                stop("Data must be stored in list mode or a matrix.")
        if(!is.na(grp[1])) {
                yy <- x
                for(j in 1:length(grp))
                        x[[j]] <- yy[[grp[j]]]
        }
for(j in 1:JK){
temp<-x[[j]]
temp<-temp[!is.na(temp)] # Remove missing values.
x[[j]]<-temp
}
xx<-list()
mloc<-NA
for(i in 1:JK){
xx[[i]]<-x[[i]]
mloc[i]<-est(xx[[i]],...)
xcen[[i]]<-xx[[i]]-mloc[i]
}
x<-xx
mat<-matrix(mloc,nrow=J,ncol=K,byrow=T)
leva<-apply(mat,1,mean) # J averages over columns
levb<-apply(mat,2,mean)
gm<-mean(levb)
testa<-sum((leva-mean(leva))^2)
testb<-sum((levb-mean(levb))^2)
testab<-NA
tempab<-matrix(NA,nrow=J,ncol=K)
for(j in 1:J){
for(k in 1:K){
tempab[j,k]<-mat[j,k]-leva[j]-levb[k]+gm
}}
testab<-sum(tempab^2)
bvec<-matrix(NA,nrow=JK,ncol=nboot)
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
if(pr)print("Taking bootstrap samples. Please wait.")
for(j in 1:JK){
if(pr)print(paste("Working on group ",j))
data<-matrix(sample(xcen[[j]],size=length(xcen[[j]])*nboot,replace=TRUE),
nrow=nboot)
bvec[j,]<-apply(data,1,est,...) # JK by nboot matrix, jth row contains
#                          bootstrapped  estimates for jth group
}
boota<-NA
bootb<-NA
bootab<-NA
for(i in 1:nboot){
mat<-matrix(bvec[,i],nrow=J,ncol=K,byrow=T)
leva<-apply(mat,1,mean) # J averages over columns
levb<-apply(mat,2,mean)
gm<-mean(mat)
boota[i]<-sum((leva-mean(leva))^2)
bootb[i]<-sum((levb-mean(levb))^2)
for(j in 1:J){
for(k in 1:K){
tempab[j,k]<-mat[j,k]-leva[j]-levb[k]+gm
}}
bootab[i]<-sum(tempab^2)}
pvala<-1-sum(testa>=boota)/nboot
pvalb<-1-sum(testb>=bootb)/nboot
pvalab<-1-sum(testab>=bootab)/nboot
list(p.value.A=pvala,p.value.B=pvalb,p.value.AB=pvalab,
test.A=testa,test.B=testb,
test.AB=testab,est.loc=matrix(mloc,nrow=J,ncol=K,byrow=T))
}