wilcox_R: What the package does (short line)

Dancols_sub <-
function(x1,y1,x2,y2,pts=NULL,fr1=1,fr2=1,alpha=.05,plotit=FALSE,xout=FALSE,outfun=out,nboot=100,SEED=TRUE,xlab='X',ylab='Y',...){
#
# Compare the OLS regression lines of two dependent (within) groups
# at specified design points
#
#  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 not specified, points are chosen for you.
#
x1=as.matrix(x1)
x2=as.matrix(x2)
if(ncol(x1)!=ncol(x2))stop('x1 and x2 have different number of columns')
if(SEED)set.seed(2)
FLAG=pts
X=elimna(cbind(x1,y1,x2,y2))
if(ncol(X)>4)stop('Only one covariate is allowed')
x1=as.matrix(x1)
x2=as.matrix(x2)
p=ncol(x1)
p1=p+1
p2=p+2
p3=p1+p
p4=p3+1
x1=X[,1:p]
y1=X[,p1]
x2=X[,p2:p3]
y2=X[,p4]
n=length(y1)
if(xout){
flag1=outfun(x1,SEED=SEED,...)$out.id
flag2=outfun(x2,SEED=SEED,...)$out.id
flag=unique(c(flag1,flag2))
if(length(flag)>0)X=X[-flag,]
x1=X[,1:p]
y1=X[,p1]
x2=X[,p2:p3]
y2=X[,p4]
}
n.keep=length(y1)
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
sqsd=difregYvar(x1,y1,x2,y2,regfun=lsfit,pts=pts,nboot=nboot,SEED=SEED)
est1=regYhat(x1,y1,xr=pts,regfun=lsfit) #Note: if xout=T, leverage points already removed
est2=regYhat(x2,y2,xr=pts,regfun=lsfit)
mat[,2]=est1
mat[,3]=est2
est=est1-est2
mat[,4]=est
sd=sqrt(sqsd)
mat[,6]=sd
tests=(est1-est2)/sd
mat[,5]=tests
df=length(y1)-1
pv=2*(1-pt(abs(tests),df))
mat[,9]=pv
crit<-smmcrit(df,5)
mat[,7]=est-crit*sd
mat[,8]=est+crit*sd
}
if(!is.null(FLAG)){
n1=1
n2=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),9)
dimnames(mat)<-list(NULL,c('X','Est1','Est2','DIF','TEST','se','ci.low','ci.hi','p.value'))
mat[,1]<-pts
sqsd=difregYvar(x1,y1,x2,y2,regfun=lsfit,pts=pts,nboot=nboot,SEED=SEED)
est1=regYhat(x1,y1,xr=pts,regfun=lsfit,,...)
est2=regYhat(x2,y2,xr=pts,regfun=lsfit,,...)
mat[,2]=est1
mat[,3]=est2
est=est1-est2
mat[,4]=est
sd=sqrt(sqsd)
mat[,6]=sd
tests=(est1-est2)/sd
mat[,5]=tests
df=length(y1)-1
pv=2*(1-pt(abs(tests),df))
mat[,9]=pv
crit<-smmcrit(df,length(pts))
mat[,7]=est-crit*sd
mat[,8]=est+crit*sd
}
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(lsfit(x1,y1)$coef)
abline(lsfit(x2,y2)$coef,lty=2)
}
list(n=n,n.keep=n.keep,output=mat)
}