wilcox_R: What the package does (short line)

DancovaV2 <-
function(x1=NULL,y1=NULL,x2=NULL,y2=NULL,xy=NULL,fr1=1,fr2=1,p.crit=NULL,
est=tmean,alpha=.05,plotit=TRUE,xlab='X',ylab='Y',qvals=c(.25,.5,.75),sm=FALSE,
xout=FALSE,eout=FALSE,outfun=out,DIF=FALSE,LP=TRUE,
nboot=500,SEED=TRUE,nreps=2000,MC=TRUE,cpp=FALSE,
SCAT=TRUE,pch1='*',pch2='+',
nmin=12,q=.5,...){
#
# Compare two dependent  groups using the ancova method.
# No parametric assumption is made about the form of
# the regression lines--a running interval smoother is used.
#
#  Like Dancova, only bootstrap samples are obtained by resampling
#  from c(x1,y1,x2,y2) rather than conditioning on the x value as done by Dancova.
#   This function tends to have more power than Dancova.
#
# One covariate only is allowed.
#
# To get critical p-value, need the following commands to get access to the software.
# library(`devtools')
# install_github( `WRScpp', `mrxiaohe')

#  Assume data are in xy having four columns: x1, y1, x2 and y2. 
#
#  Or can have the
#  data stored in four separate variables:
#   x1 y1 x2 and y2
#
#   x1 y1 are measures at time 1
#   x2 y2 are measures at time 2
#
#  LP=T, when plotting, running interval smoother is smoothed again using lplot.
#  sm=T will create smooths using bootstrap bagging.
#  pts can be used to specify the design points where the regression lines
#  are to be compared.
#
#  q=.5 means when est=hd (Harrell-Davis estimator), median is estimated.
#
#  eout=TRUE will eliminate all outliers when plotting. 
#
if(SEED)set.seed(2)
iter=nreps
if(!is.null(x1[1])){
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(x1)!=length(x2))stop('x1 and y2 have different lengths')
if(length(x2)!=length(y2))stop('x2 and y2 have different lengths')
if(length(y1)!=length(y2))stop('y1 and y2 have different lengths')
xy=cbind(x1,y1,x2,y2)
}
n=nrow(elimna(xy))
if(plotit){
ef=identical(est,hd)
if(!ef)runmean2g(xy[,1],xy[,2],xy[,3],xy[,4],fr=fr1,est=est,sm=sm,xout=xout,LP=LP,eout=eout,
xlab=xlab,ylab=ylab,SCAT=SCAT,pch1=pch1,pch2=pch2,...)
if(ef)runmean2g(xy[,1],xy[,2],xy[,3],xy[,4],fr=fr1,est=hd,sm=sm,xout=xout,LP=LP,q=q,eout=eout,
xlab=xlab,ylab=ylab,SCAT=SCAT,pch1=pch1,pch2=pch2,...)
}

if(is.null(p.crit)){
if(cpp)library(WRScpp)
p.crit=DancGLOB_pv(n,fr1=fr1,fr2=fr2,nboot=nboot,est=est,SEED=SEED,iter=iter,
nmin=nmin,MC=MC,alpha=alpha,qvals=qvals,cpp=cpp)$p.crit
#
}
pts=NULL
#if(is.null(pts)){
for(i in 1:length(qvals))pts=c(pts,qest(xy[,1],qvals[i]))
#}
if(SEED)set.seed(2)
ef=identical(est,hd)
n=nrow(xy)
est1=NA
est2=NA
J=length(pts)
est1=matrix(NA,nrow=nboot,ncol=J)
est2=matrix(NA,nrow=nboot,ncol=J)
#
data=matrix(sample(n,size=n*nboot,replace=TRUE),ncol=nboot,nrow=n)
if(!MC){
if(!ef){
est1=apply(data,2,DancGLOB_sub,xy=xy[,1:2],pts=pts,est=est,fr=fr1,nmin=nmin,...)
est2=apply(data,2,DancGLOB_sub,xy=xy[,3:4],pts=pts,est=est,fr=fr2,nmin=nmin,...)
}
if(ef){
est1=apply(data,2,DancGLOB_sub,xy=xy[,1:2],pts=pts,est=hd,fr=fr1,nmin=nmin,q=q,...)
est2=apply(data,2,DancGLOB_sub,xy=xy[,3:4],pts=pts,est=hd,fr=fr2,nmin=nmin,q=q,...)
}
est1=t(as.matrix(est1))
est2=t(as.matrix(est2))
}
if(MC){
library(parallel)
data=listm(data)
if(!ef){
est1=mclapply(data,DancGLOB_sub,xy=xy[,1:2],pts=pts,est=est,fr=fr1,nmin=nmin,...)
est2=mclapply(data,DancGLOB_sub,xy=xy[,3:4],pts=pts,est=est,fr=fr2,nmin=nmin,...)
}
if(ef){
est1=mclapply(data,DancGLOB_sub,xy=xy[,1:2],pts=pts,est=hd,fr=fr1,nmin=nmin,q=q,...)
est2=mclapply(data,DancGLOB_sub,xy=xy[,3:4],pts=pts,est=hd,fr=fr2,nmin=nmin,q=q,...)
}
est1=t(matl(est1))
est2=t(matl(est2))
}
pv=NA
for(j in 1:J){
pv[j]=mean(est1[,j]<est2[,j],na.rm=TRUE)+.5*mean(est1[,j]==est2[,j],na.rm=TRUE)
pv[j]=2*min(c(pv[j],1-pv[j]))
}
pvm=cbind(pts,pv)
dimnames(pvm)=list(NULL,c('X','p.values'))
list(output=pvm,n=n,p.crit=p.crit)
}