wilcox_R: What the package does (short line)

ancovaV2 <-
function(x1=NULL,y1=NULL,x2=NULL,y2=NULL,fr1=1,fr2=1,p.crit=NULL,padj=FALSE, pr=TRUE,
method='hochberg',FAST=TRUE,
est=tmean,alpha=.05,plotit=TRUE,xlab='X',ylab='Y',qpts=FALSE,qvals=c(.25,.5,.75),sm=FALSE,
xout=FALSE,eout=FALSE,outfun=out,LP=TRUE,nboot=500,SEED=TRUE,nreps=2000,MC=FALSE,
nmin=12,q=.5,SCAT=TRUE,pch1='*',pch2='+',...){
#
# 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.
#
# One covariate only is allowed.
#
#  Like ancova, only bootstrap samples are obtained by resampling
#  from c(x1,y1) and c(x2,y2) rather than conditioning on the x value as done by ancova.
#   This function tends to have more power than ancova.
#
#
#  padj=TRUE, p-values are adjusted using the method indicated by the argument
#  method. By default, Hochberg's method is used Setting method='hommel' will use Hommel's method.
# By default, adjusted p-values are not reported; by default the critical p-value, p_c, is used to get better power.
#
#   x1 y1 are measures for group 1
#   x2 y2 are measures for group 2
#
#  LP=T, when plotting, running interval smoother is smoothed again using lplot.
#  sm=T will create smooths using bootstrap bagging.
#
#  qvals: Determines which quantiles of the covariate x1 will be used when 
#   qpts=TRUE
#  qpts=FALSE means covariate chosen as done by the function ancova
#
#  q=.5 means when est=hd (Harrell-Davis estimator), median is estimated.
#
#  eout=TRUE will eliminate all outliers when plotting. 
#
#  nreps: indicates number of replications used to determine a critical value and p-value.
#
if(SEED)set.seed(2)
if(pr){
if(!FAST){
if(!MC & is.null(p.crit))print('For faster execution time, set MC=TRUE, assuming a multi-core processor is available')
}}
if(padj)p.crit=0
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')
xy1=elimna(cbind(x1,y1))
xy2=elimna=cbind(x2,y2)
n1=nrow(xy1)
n2=nrow(xy2)
if(plotit){
ef=identical(est,hd)
if(!ef)runmean2g(xy1[,1],xy1[,2],xy2[,1],xy2[,2],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(xy1[,1],xy1[,2],xy2[,1],xy2[,2],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(FAST){
if(alpha==.05){
nm=max(c(n1,n2))
if(nm<=800){
nv=c(50,60,80,100,200,300,500,800)
if(qpts){
pv=c(.02709,.0283,.0306,.02842,.02779,.02410,.02683,.01868,.02122)
p.crit=lplot.pred(1/nv,pv,1/n1)$yhat
}
if(!qpts){
pv=c(.020831,.017812,.015796,.014773,.012589,.015664,.011803,.012479)
p.crit=lplot.pred(1/nv,pv,1/n1)$yhat
}}
}}}
if(is.null(p.crit)){
p.crit=ancovaV2.pv(n1,n2,nreps=nreps,MC=MC,qpts=qpts,est=est,qvals=qvals,SEED=SEED,
alpha=alpha,nboot=nboot)$p.crit
}
pts=NULL
if(qpts)for(i in 1:length(qvals))pts=c(pts,qest(xy1[,1],qvals[i]))
if(!qpts)pts=ancova(x1,y1,x2,y2,pr=FALSE,plotit=FALSE)$output[,1]
if(SEED)set.seed(2)
ef=identical(est,hd)
est1=NA
est2=NA
J=length(pts)
est1=matrix(NA,nrow=nboot,ncol=J)
est2=matrix(NA,nrow=nboot,ncol=J)
#
data1=matrix(sample(n1,size=n1*nboot,replace=TRUE),ncol=nboot,nrow=n1)
data2=matrix(sample(n2,size=n2*nboot,replace=TRUE),ncol=nboot,nrow=n2)
if(!MC){
if(!ef){
est1=apply(data1,2,DancGLOB_sub,xy=xy1[,1:2],pts=pts,est=est,fr=fr1,nmin=nmin,...)
est2=apply(data2,2,DancGLOB_sub,xy=xy2[,1:2],pts=pts,est=est,fr=fr2,nmin=nmin,...)
}
if(ef){
est1=apply(data1,2,DancGLOB_sub,xy=xy1[,1:2],pts=pts,est=hd,fr=fr1,nmin=nmin,q=q,...)
est2=apply(data2,2,DancGLOB_sub,xy=xy2[,1:2],pts=pts,est=hd,fr=fr2,nmin=nmin,q=q,...)
}
est1=t(as.matrix(est1))
est2=t(as.matrix(est2))
}
if(MC){
library(parallel)
data1=listm(data1)
data2=listm(data2)
if(!ef){
est1=mclapply(data1,DancGLOB_sub,xy=xy1[,1:2],pts=pts,est=est,fr=fr1,nmin=nmin,...)
est2=mclapply(data2,DancGLOB_sub,xy=xy2[,1:2],pts=pts,est=est,fr=fr2,nmin=nmin,...)
}
if(ef){
est1=mclapply(data1,DancGLOB_sub,xy=xy1[,1:2],pts=pts,est=hd,fr=fr1,nmin=nmin,q=q,...)
est2=mclapply(data2,DancGLOB_sub,xy=xy2[,1:2],pts=pts,est=hd,fr=fr2,nmin=nmin,q=q,...)
}
est1=t(matl(est1))
est2=t(matl(est2))
}
pv=NULL
if(J==1){
est1=t(as.matrix(est1))
est2=t(as.matrix(est2))
}
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]))
}
pvadj=rep(NA,length(pts))
if(padj)pvadj=p.adjust(pv,method=method)
pvm=cbind(pts,pv,pvadj)
dimnames(pvm)=list(NULL,c('X','p.values','p.adjusted'))
list(output=pvm,n=c(n1,n2),p.crit=p.crit)
}