wilcox_R: What the package does (short line)

anctspb <-
function(x1,y1,x2,y2,pts=NULL,regfun=tshdreg,fr1=1,fr2=1,alpha=.05,plotit=TRUE,xout=FALSE,outfun=outpro,nboot=500,SEED=TRUE,xlab='X',ylab='Y',...){
#
# Compare the regression lines of two independent groups 
# at specified design points using a robust regression estimator.
#
#  Like ancts but uses 
#   a percentile bootstrap method is used. 
# This might help when there are tied values among the dependent variable. 
#
#  Assume data are in x1 y1 x2 and y2
#
#  pts can be used to specify the design points where the regression lines
#  are to be compared.
#
if(SEED)set.seed(2)
FLAG=pts
xy=elimna(cbind(x1,y1))
if(ncol(xy)>2)stop('Only one covariate is allowed')
x1=xy[,1]
y1=xy[,2]
xy=elimna(cbind(x2,y2))
if(ncol(xy)>2)stop('Only one covariate is allowed')
x2=xy[,1]
y2=xy[,2]
if(is.null(pts[1])){
npt<-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,length(pts),9)
dimnames(mat)<-list(NULL,c('X','Est1','Est2','DIF','TEST','se','ci.low','ci.hi','p.value'))
pts=x1[isub]
}
mat<-matrix(NA,length(pts),7)
dimnames(mat)<-list(NULL,c('X','Est1','Est2','DIF','ci.low','ci.hi','p.value'))
mat[,1]<-pts
bvec1=matrix(NA,nrow=nboot,ncol=length(pts))
bvec2=matrix(NA,nrow=nboot,ncol=length(pts))
x1=as.matrix(x1)
x2=as.matrix(x2)
p1=ncol(x1)+1
data<-matrix(sample(length(y1),size=length(y1)*nboot,replace=TRUE),nrow=nboot)
for(ib in 1:nboot){
bvec1[ib,]=regYsub(x1[data[ib,],],y1[data[ib,]],pts,p1=p1,regfun=regfun,...)
}
data<-matrix(sample(length(y2),size=length(y2)*nboot,replace=TRUE),nrow=nboot)
for(ib in 1:nboot){
bvec2[ib,]=regYsub(x2[data[ib,],],y2[data[ib,]],pts,p1=p1,regfun=regfun,...)
}
dif=bvec1<bvec2
L=apply(dif,2,mean)
E=bvec1==bvec2
T=apply(E,2,mean)
pvec=L+.5*T
pop=1-pvec
pb=cbind(pvec,pop)
pv=2*apply(pb,1,min)
ilow<-round((alpha/2) * nboot)
ihi<-nboot - ilow
ilow<-ilow+1
ciL=NA
ciU=NA
difb=bvec1-bvec2
for(i in 1:length(pts)){
bs=sort(difb[,i])
ciL[i]=bs[ilow]
ciU[i]=bs[ihi]
}
est1=regYhat(x1,y1,xr=pts,xout=xout,outfun=outfun,regfun=regfun,...)
est2=regYhat(x2,y2,xr=pts,xout=xout,outfun=outfun,regfun=regfun,...)
mat[,2]=est1
mat[,3]=est2
est=est1-est2
mat[,4]=est
mat[,5]=ciL
mat[,6]=ciU
mat[,7]=pv
if(!is.null(FLAG[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)])
}
est1=regYhat(x1,y1,regfun=regfun,xr=pts,xout=xout,outfun=outfun,...)
est2=regYhat(x2,y2,regfun=regfun,xr=pts,xout=xout,outfun=outfun,...)
mat[,2]=est1
mat[,3]=est2
est=est1-est2
mat[,4]=est
mat[,7]=pv
}
if(plotit){
if(xout){
flag<-outfun(x1,...)$keep
x1<-x1[flag]
y1<-y1[flag]
flag<-outfun(x2,...)$keep
x2<-x2[flag]
y2<-y2[flag]
}
plot(c(x1,x2),c(y1,y2),type='n',xlab=xlab,ylab=ylab)
points(x1,y1,pch='o')
points(x2,y2,pch='+')
abline(regfun(x1,y1)$coef)
abline(regfun(x2,y2)$coef,lty=2)
}
list(output=mat)
}