wilcox_R: What the package does (short line)

ancdifplot <-
function(x1,y1,x2,y2,fr1=1,fr2=1,tr=.2,alpha=.05,pr=TRUE,xout=FALSE,outfun=out,LP=TRUE,
nmin=8,scat=TRUE,xlab='X',ylab='Y',report=FALSE,...){
#
# Compare two independent  groups using the ancova method
# 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
#
#  nmin indicates minimun number of values close to a point 
#
#  Similar to ancova, only compute a confidence band for the difference and plot it.
#
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){
flago<-outfun(x1,...)$keep
x1<-x1[flago]
y1<-y1[flago]
flag<-outfun(x2,...)$keep
x2<-x2[flago]
y2<-y2[flago]
}
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])
flag=vecn>=nmin
ptsum=sum(flag)
est=NA
low=NA
up=NA
ic=0
xp1=NA
xp2=NA
pv=NA
for (i in 1:length(x1)){
if(flag[i]){
g1<-y1[near(x1,x1[i],fr1)]
g2<-y2[near(x2,x2[i],fr2)]
test<-yuen(g1,g2,tr=tr)
ic=ic+1
xp1[ic]=x1[i]
xp2[ic]=x2[i]
est[ic]=test$dif
low[ic]=test$ci[1]
up[ic]=test$ci[2]
pv[ic]=test$p.value
}}
#print(length(pv))
#print(length(xp1))
if(LP){
xy=elimna(cbind(xp1,est,low,up,pv))
est=lplot(xy[,1],xy[,2],plotit=FALSE,pyhat=TRUE,pr=FALSE)$yhat
up=lplot(xy[,1],xy[,4],plotit=FALSE,pyhat=TRUE,pr=FALSE)$yhat
low=lplot(xy[,1],xy[,3],plotit=FALSE,pyhat=TRUE,pr=FALSE)$yhat
}
if(!report)output='DONE'
plot(c(x1,x2),c(y1,y2),xlab=xlab,ylab=ylab,type='n')
if(!LP){
lines(xp1,up,lty=2)
lines(xp1,low,lty=2)
lines(xp1,est)
if(scat)points(c(x1,x2),c(y1,y2))
if(report){
output=cbind(xp1,est,low,up,pv)
dimnames(output)=list(NULL,c(xlab,'est.dif','lower.ci','upper.ci','p-value'))
}}
if(LP){
lines(xy[,1],up,lty=2)                                                                   
lines(xy[,1],low,lty=2)                                        
lines(xy[,1],est)                                                       
 if(scat)points(c(x1,x2),c(y1,y2))
if(report){
output=cbind(xy[,1],est,low,up,xy[,5])
dimnames(output)=list(NULL,c(xlab,'est.dif','lower.ci','upper.ci','p-value'))
}
}
output
}