wilcox_R: What the package does (short line)

ancJNmp <-
function(x1,y1,x2,y2,regfun=qreg,p.crit=NULL,DEEP=FALSE,
plotit=TRUE,xlab='X1',ylab='X2',null.value=0,FRAC=.5,cov1=FALSE,SMM=TRUE,
alpha=.05,nreps=1000, MC=FALSE, pts=NULL,SEED=TRUE,nboot=100,xout=FALSE,outfun=out,...){
#
# Compare two independent  groups using a generalization of the ancts function that
#  allows more than one covariate.
#
# DEEP=FALSE: If pts=NULL, design points are chosen to be deepest  point in 
# x1  plus points on the .5 depth contour.
#
# DEEP=TRUE, choose deepest half of c(x1,x2) and use critical p-value indicated by
# p.crit, the critical p-value,which defaults to .015 when alpha=.05. 
# If alpha!=.05, p.crit must be computed, which can require high execution time.
# MC=TRUE will reduce execution time considerably.
#
#  cov1=TRUE: the covariates that are used are taken to be the points in x1. If
#
#  plotit=TRUE: if p=2 covariates, plot covariate points with non-significant points indicated by * and
#  significant points by +

# (This function replaces anctsmp, which does not have an option for using the deepest half of covariate points.)
#
if(SEED)set.seed(2) # now cov.mve always returns same result

x1=as.matrix(x1)
x2=as.matrix(x2)
if(ncol(x1)!=ncol(x2))stop('x1 and x2 have a different number of columns')
p=ncol(x1)
p1=p+1
m1=elimna(cbind(x1,y1))
x1=m1[,1:p]
y1=m1[,p1]
m2=elimna(cbind(x2,y2))
x2=m2[,1:p]
y2=m2[,p1]
#
if(xout){
m<-cbind(x1,y1)
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]
}
n1=length(y1)
n2=length(y2)
nv=c(n1,n2)
if(cov1){
pts=x1
if(alpha==.05)p.crit=0.00676
}
if(DEEP)pts=NULL
if(!is.null(pts[1])){
p.crit=NULL
DEEP=FALSE
}
if(is.null(pts[1])){
if(!DEEP){
x1<-as.matrix(x1)
pts<-ancdes(unique(rbind(x1,x2)))
p.crit=NULL
}
if(DEEP){
pts=ancov2COV(x1,y1,x2,y2,DETAILS=TRUE,cr=.27,pr=FALSE,FRAC=FRAC)$all.points.used[,1:2]
}}
pts<-as.matrix(pts)
ntests=nrow(pts)
mat<-matrix(NA,ntests,8)
dimnames(mat)<-list(NULL,c('Est 1', 'Est 2','DIF','TEST','se','ci.low','ci.hi','p.value'))
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,...)
#  xout=F because leverage points have already been removed.
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[,1]=est1
mat[,2]=est2
est=est1-est2
mat[,3]=est
sd=sqrt(sqsd1+sqsd2)
mat[,5]=sd
tests=(est1-est2)/sd
mat[,4]=tests
pv=2*(1-pnorm(abs(tests)))
mat[,8]=pv
crit=NULL
if(!cov1){
if(!DEEP){
if(ntests==1)crit=qnorm(1-alpha/2)
if(length(pts)>1){
if(ntests<=28){
if(alpha==.05)crit<-smmcrit(Inf,ntests)
if(alpha==.01)crit<-smmcrit01(Inf,ntests)
}
if(ntests>28 || is.null(crit))crit=smmvalv2(dfvec=rep(Inf,nrow(pts)),alpha=alpha)
}}}
if(cov1){
if(!DEEP){
if(alpha==.05)p.crit=0.00676
if(alpha!=.05)p.crit=ancJNmpcp(n1,n2,alpha=alpha,regfun=regfun,nreps=nreps,MC=MC,cov1=cov1,
DEEP=FALSE)$pc.est
crit=qnorm(1-p.crit/2)
}}
if(DEEP){
if(p==2){
p.crit=.012
if(alpha!=.05)p.crit=ancJNmpcp(n1,n2,alpha=alpha,regfun=regfun,nreps=nreps,MC=MC,cov1=cov1,
DEEP=TRUE)$pc.est
crit=qnorm(1-p.crit/2)
}
if(p>2){
if(length(pts)>1){
if(SMM){
if(ntests<=28){
if(alpha==.05)crit<-smmcrit(Inf,ntests)
if(alpha==.01)crit<-smmcrit01(Inf,ntests)
}
if(ntests>28 || is.null(crit))crit=smmvalv2(dfvec=rep(Inf,nrow(pts)),alpha=alpha)
}
if(!SMM){
p.crit=ancJNmpcp(n1,n2,alpha=alpha,regfun=regfun,nreps=nreps,MC=MC,cov1=FALSE)
crit=qnorm(1-p.crit/2)
}}}}
mat[,6]=est-crit*sd
mat[,7]=est+crit*sd
flag=rep(FALSE,nrow(mat))
flag.chk1=as.logical(mat[,6]>null.value)
flag.chk2=(mat[,7]<null.value)
flag.chk=(flag.chk1+flag.chk2>0)
num.sig=sum(flag.chk)
if(p==2){
if(plotit){
plot(pts[,1],pts[,2],xlab=xlab,ylab=ylab,type='n')
flag[flag.chk]=TRUE
points(pts[!flag,1],pts[!flag,2],pch='*')
points(pts[flag,1],pts[flag,2],pch='+')
}}
list(n=nv,num.sig=num.sig,p.crit=p.crit,points=pts,output=mat)
}