wilcox_R: What the package does (short line)

difregMC <-
function(x1,y1,x2,y2,regfun=tsreg,xout=FALSE,outfun=outpro,nboot=599,
alpha=.05,SEED=TRUE,plotit=FALSE,xlab='X',ylab='Y',pr=TRUE,...){
#
# regression data from two different times i.e., two dependent groups
#
#  compute confidence interval for the difference in the slopes and intercepts
#
library(parallel)
if(SEED)set.seed(2)
X=elimna(cbind(x1,y1,x2,y2))
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,...)$out.id
flag2=outfun(x2,...)$out.id
flag=unique(c(flag1,flag2))
X=X[-flag,]
x1=X[,1:p]
y1=X[,p1]
x2=X[,p2:p3]
y2=X[,p4]
}
flagF=identical(regfun,tsreg)
if(flagF){
if(pr){
if(sum(duplicated(y1)>0))print('Duplicate values detected; regfun=tshdreg might have more power than tsreg')
pr=FALSE
}
if(pr){
if(sum(duplicated(y2)>0))print('Duplicate values detected; regfun=tshdreg might have more power than tsreg')
}}
nk=length(y1)
x1=as.matrix(x1)
x2=as.matrix(x2)
data<-matrix(sample(length(y1),size=length(y1)*nboot,replace=TRUE),nrow=nboot)
data=listm(t(data))
bvec1<-mclapply(data,regboot,x1,y1,regfun,mc.preschedule=TRUE,xout=FALSE,...)
# bvec is a p+1 by nboot matrix. The first row
#                     contains the bootstrap intercepts, the second row
#                     contains the bootstrap values for first predictor, etc.
bvec2<-mclapply(data,regboot,x2,y2,regfun,mc.preschedule=TRUE,xout=FALSE,...)
bvec1=matl(bvec1)
bvec2=matl(bvec2)
dif=t(bvec1)-t(bvec2)
dif.sort=apply(dif,2,sort)
pvec=NA
for(i in 1:p1){
pvec[i]<-(sum(dif[,i]<0)+.5*sum(dif[,i]==0))/nboot
if(pvec[i]>.5)pvec[i]<-1-pvec[i]
}
pvec<-2*pvec
ilow<-round((alpha/2) * nboot)
ihi<-nboot - ilow
ilow<-ilow+1
ci=matrix(NA,nrow=p1,ncol=3)
ci[,1]=c(0:p)
for(i in 1:p1){
ci[i,2]=dif.sort[ilow,i]
ci[i,3]=dif.sort[ihi,i]
}
dimnames(ci)=list(NULL,c('Param','ci.low','ci.hi'))
if(plotit){
reg2plot(x1,y1,x2,y2,xlab=xlab,ylab=ylab,regfun=regfun,...)
}
lvec='Intercept'
for(j in 2:p1)lvec=c(lvec,paste('slope',j-1))
pvec=array(pvec,dimnames=lvec)
est1=regfun(x1,y1,xout=FALSE,...)$coef
est2=regfun(x2,y2,xout=FALSE,...)$coef
list(n=n,n.keep=nk,param=lvec,p.values=pvec,est.grp1=est1,est.grp2=est2,conf.intervals=ci)
}