wilcox_R: What the package does (short line)

mulwmw <-
function(m1,m2,plotit=TRUE,cop=3,alpha=.05,nboot=1000,pop=4,fr=.8,pr=FALSE){
#
#
# Determine center correpsonding to two
# independent groups, project all  points onto line
# connecting the centers,
# then based on the projected distances,
# estimate p=probability that a randomly sampled
# point from group 1 is less than a point from group 2
# based on the projected distances.
#
# plotit=TRUE creates a plot of the projected data
# pop=1 plot two dotplots based on projected distances
# pop=2 boxplots
# pop=3 expected frequency curve.
# pop=4 adaptive kernel density
#
#  There are three options for computing the center of the
#  cloud of points when computing projections:
#  cop=1 uses Donoho-Gasko median
#  cop=2 uses MCD center
#  cop=3 uses median of the marginal distributions.
#
#  When using cop=2 or 3, default critical value for outliers
#  is square root of the .975 quantile of a
#  chi-squared distribution with p degrees
#  of freedom.
#
#  Donoho-Gasko (Tukey) median is marked with a cross, +.
#
if(is.null(dim(m1))||dim(m1)[2]<2){print("Data are assumed to be stored in")
print(" a matrix or data frame having two or more columns.")
stop(" For univariate data, use the function outbox or out")
}
m1<-elimna(m1) # Remove missing values
m2<-elimna(m2)
n1=nrow(m1)
n2=nrow(m2)
if(cop==1){
if(ncol(m1)>2){
center1<-dmean(m1,tr=.5)
center2<-dmean(m2,tr=.5)
}
if(ncol(m1)==2){
tempd<-NA
for(i in 1:nrow(m1))
tempd[i]<-depth(m1[i,1],m1[i,2],m1)
mdep<-max(tempd)
flag<-(tempd==mdep)
if(sum(flag)==1)center1<-m1[flag,]
if(sum(flag)>1)center1<-apply(m1[flag,],2,mean)
for(i in 1:nrow(m2))
tempd[i]<-depth(m2[i,1],m2[i,2],m2)
mdep<-max(tempd)
flag<-(tempd==mdep)
if(sum(flag)==1)center2<-m2[flag,]
if(sum(flag)>1)center2<-apply(m2[flag,],2,mean)
}}
if(cop==2){
center1<-cov.mcd(m1)$center
center2<-cov.mcd(m2)$center
}
if(cop==3){
center1<-apply(m1,2,median)
center2<-apply(m2,2,median)
}
if(cop==4){
center1<-smean(m1)
center2<-smean(m2)
}
center<-(center1+center2)/2
B<-center1-center2
if(sum(center1^2)<sum(center2^2))B<-(0-1)*B
BB<-B^2
bot<-sum(BB)
disx<-NA
disy<-NA
if(bot!=0){
for (j in 1:nrow(m1)){
AX<-m1[j,]-center
tempx<-sum(AX*B)*B/bot
disx[j]<-sign(sum(AX*B))*sqrt(sum(tempx^2))
}
for (j in 1:nrow(m2)){
AY<-m2[j,]-center
tempy<-sum(AY*B)*B/bot
disy[j]<-sign(sum(AY*B))*sqrt(sum(tempy^2))
}
}
if(plotit){
if(pop==1){
par(yaxt="n")
xv<-rep(2,length(disx))
yv<-rep(1,length(disy))
plot(c(disx,disy),c(xv,yv),type="n",xlab="",ylab="")
xv<-rep(1.6,length(disx))
yv<-rep(1.4,length(disy))
points(disx,xv)
points(disy,yv)
}
if(pop==2)boxplot(disx,disy)
if(pop==3)rd2plot(disx,disy,fr=fr)
if(pop==4)g2plot(disx,disy,fr=fr)
}
m<-outer(disx,disy,FUN="-")
m<-sign(m)
phat<-(1-mean(m))/2
if(bot==0)phat<-.5
print("Computing critical values")
m1<-t(t(m1)-center1)
m2<-t(t(m2)-center2)
v1<-mulwmwcrit(m1,m2,cop=cop,alpha=alpha,iter=nboot,pr=pr)
list(phat=phat,lower.crit=v1[1],upper.crit=v1[2],n1=n1,n2=n2)
}