wilcox_R: What the package does (short line)

ancovaWMW <-
function(x1,y1,x2,y2,fr1=1,fr2=1,alpha=05,
plotit=TRUE,pts=NA,xout=FALSE,outfun=out,LP=TRUE,sm=FALSE,est=hd,...){
#
# Compare two independent  groups using the ancova method in conjunction 
# with Cliff's improvement on the Wilcoxon-Mann-Whitney test.
# No parametric assumption is made about the form of
# the regression lines--a running interval smoother is used.
#
#  Assume data are in x1 y1 x2 and y2
#
#  sm=TRUE will use bootstrap bagging when plotting the regression lines
#  The plot is based on measure of location indicated by the argument
#  est. Default is the Harrell-Davis estimate of the median.
#
#   LP=TRUE: use running interval smoother followed by LOESS
#
if(ncol(as.matrix(x1))>1)stop('One covariate only is allowed with this function')
if(length(x1)!=length(y1))stop('x1 and y1 have different lengths')
if(length(x2)!=length(y2))stop('x2 and y2 have different lengths')
xy=elimna(cbind(x1,y1))
x1=xy[,1]
y1=xy[,2]
xy=elimna(cbind(x2,y2))
x2=xy[,1]
y2=xy[,2]
if(xout){
flag<-outfun(x1,...)$keep
x1<-x1[flag]
y1<-y1[flag]
flag<-outfun(x2,...)$keep
x2<-x2[flag]
y2<-y2[flag]
}
dv.sum=NULL
if(is.na(pts[1])){
npt<-5
CC=5
isub<-c(1:5)  # Initialize isub
test<-c(1:5)
xorder<-order(x1)
y1<-y1[xorder]
x1<-x1[xorder]
xorder<-order(x2)
y2<-y2[xorder]
x2<-x2[xorder]
n1<-1
n2<-1
vecn<-1
for(i in 1:length(x1))n1[i]<-length(y1[near(x1,x1[i],fr1)])
for(i in 1:length(x1))n2[i]<-length(y2[near(x2,x1[i],fr2)])
for(i in 1:length(x1))vecn[i]<-min(n1[i],n2[i])
sub<-c(1:length(x1))
isub[1]<-min(sub[vecn>=12])
isub[5]<-max(sub[vecn>=12])
isub[3]<-floor((isub[1]+isub[5])/2)
isub[2]<-floor((isub[1]+isub[3])/2)
isub[4]<-floor((isub[3]+isub[5])/2)
mat<-matrix(NA,5,8)
dimnames(mat)<-list(NULL,c('X','n1','n2','p.hat','ci.low','ci.hi','p.value','p.crit'))
for (i in 1:5){
g1<-y1[near(x1,x1[isub[i]],fr1)]
g2<-y2[near(x2,x1[isub[i]],fr2)]
g1<-g1[!is.na(g1)]
g2<-g2[!is.na(g2)]
test<-cidv2(g1,g2)
dv.sum=rbind(dv.sum,test$summary.dvals)
mat[i,1]<-x1[isub[i]]
mat[i,2]<-length(g1)
mat[i,3]<-length(g2)
mat[i,4]<-test$p.hat
mat[i,5]<-test$p.ci[1]
mat[i,6]<-test$p.ci[2]
mat[i,7]<-test$p.value
}}
if(!is.na(pts[1])){
CC=length(pts)
n1<-1
n2<-1
vecn<-1
for(i in 1:length(pts)){
n1[i]<-length(y1[near(x1,pts[i],fr1)])
n2[i]<-length(y2[near(x2,pts[i],fr2)])
}
mat<-matrix(NA,length(pts),8)
dimnames(mat)<-list(NULL,c('X','n1','n2','p.hat','ci.low','ci.hi','p.value','p.crit'))
for (i in 1:length(pts)){
g1<-y1[near(x1,pts[i],fr1)]
g2<-y2[near(x2,pts[i],fr2)]
g1<-g1[!is.na(g1)]
g2<-g2[!is.na(g2)]
test=cidv2(g1,g2)
dv.sum=rbind(dv.sum,test$summary.dvals)
mat[i,1]<-pts[i]
mat[i,2]<-length(g1)
mat[i,3]<-length(g2)
if(length(g1)<=5)print(paste('Warning, there are',length(g1),' points corresponding to the design point X=',pts[i]))
if(length(g2)<=5)print(paste('Warning, there are',length(g2),' points corresponding to the design point X=',pts[i]))
mat[i,4]<-test$p.hat
mat[i,5]<-test$p.ci[1]
mat[i,6]<-test$p.ci[2]
mat[i,7]<-test$p.value
}}
dvec<-alpha/c(1:CC)
temp2<-order(0-mat[,6])
mat[temp2,8]=dvec
if(plotit){
runmean2g(x1,y1,x2,y2,fr=fr1,est=est,sm=sm,xout=FALSE,LP=LP,...)
}
list(output=mat,summary=dv.sum)
}