wilcox_R: What the package does (short line)

ancts <-
function(x1,y1,x2,y2,pts=NULL,regfun=tsreg,fr1=1,fr2=1,SCAT=TRUE,pch1='*',pch2='+',
alpha=.05,plotit=TRUE,xout=FALSE,outfun=out,nboot=100,SEED=TRUE,xlab="X",ylab="Y",...){
#
# Compare the regression lines of two independent groups at specified design points
# using a robust regression estimator. 
# By default, use the Theil--Sen estimator
#
#  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(identical(outfun,boxplot))stop('Use outfun=outbox')
if(SEED)set.seed(2)
FLAG=pts
xy=elimna(cbind(x1,y1))
if(ncol(xy)>2)stop("Only one covariate is allowed. Use ancJNmp")
x1=xy[,1]
y1=xy[,2]
xy=elimna(cbind(x2,y2))
if(ncol(xy)>2)stop("Only one covariate is allowed. Use ancJNmp")
x2=xy[,1]
y2=xy[,2]
if(xout){
m<-cbind(x1,y1)
p1=ncol(m)
p=p1-1
flag<-outfun(x1,plotit=FALSE,...)$keep
m<-m[flag,]
x1<-m[,1:p]
y1<-m[,p1]
m<-cbind(x2,y2)
flag<-outfun(x2,plotit=FALSE,...)$keep
m<-m[flag,]
x2<-m[,1:p]
y2<-m[,p1]
}
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,5,9)
dimnames(mat)<-list(NULL,c("X","Est1","Est2","DIF","TEST","se","ci.low","ci.hi","p.value"))
pts=x1[isub]
mat[,1]=pts
sqsd1=regYvar(x1,y1,pts=pts,regfun=regfun,nboot=nboot,SEED=FALSE,xout=FALSE,outfun=outfun,...)
sqsd2=regYvar(x2,y2,pts=pts,regfun=regfun,nboot=nboot,SEED=FALSE,xout=FALSE,outfun=outfun,...)
est1=regYhat(x1,y1,xr=pts,regfun=regfun,xout=FALSE,outfun=outfun,...)
est2=regYhat(x2,y2,xr=pts,regfun=regfun,xout=FALSE,outfun=outfun,...)
mat[,2]=est1
mat[,3]=est2
est=est1-est2
mat[,4]=est
sd=sqrt(sqsd1+sqsd2)
mat[,6]=sd
tests=(est1-est2)/sd
mat[,5]=tests
pv=2*(1-pnorm(abs(tests)))
mat[,9]=pv
crit<-smmcrit(Inf,5)
mat[,7]=est-crit*sd
mat[,8]=est+crit*sd
}
if(!is.null(FLAG)){
n1=NA
n2=NA
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),9)
dimnames(mat)<-list(NULL,c("X","Est1","Est2","DIF","TEST","se","ci.low","ci.hi","p.value"))
mat[,1]<-pts
sqsd1=regYvar(x1,y1,regfun=regfun,pts=pts,nboot=nboot,SEED=FALSE,xout=FALSE,outfun=outfun,...)
sqsd2=regYvar(x2,y2,regfun=regfun,pts=pts,nboot=nboot,SEED=FALSE,xout=FALSE,outfun=outfun,...)
est1=regYhat(x1,y1,regfun=regfun,xr=pts,xout=FALSE,outfun=outfun,...)
est2=regYhat(x2,y2,regfun=regfun,xr=pts,xout=FALSE,outfun=outfun,...)
mat[,2]=est1
mat[,3]=est2
est=est1-est2
mat[,4]=est
sd=sqrt(sqsd1+sqsd2)
mat[,6]=sd
tests=(est1-est2)/sd
mat[,5]=tests
pv=2*(1-pnorm(abs(tests)))
mat[,9]=pv
crit<-smmcrit(Inf,length(pts))
mat[,7]=est-crit*sd
mat[,8]=est+crit*sd
}
reg1=regfun(x1,y1,...)$coef
reg2=regfun(x2,y2,...)$coef
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)
if(SCAT){
points(x1,y1,pch=pch1)
points(x2,y2,pch=pch2)
}
abline(reg1)
abline(reg2,lty=2)
}
list(intercept.slope.group1=reg1,intercept.slope.group2=reg2,output=mat)
}