wilcox_R: What the package does (short line)

difregOLS <-
function(x1,y1,x2,y2,regfun=lsfit,xout=FALSE,outfun=outpro,nboot=200,
alpha=.05,SEED=TRUE,plotit=FALSE,xlab='X',ylab='Y',...){
#
# OLS regression data from two different times i.e., two dependent groups
#
#  compute confidence interval for the difference between intercepts
#  and the slopes
#
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))
if(length(flag)>0)X=X[-flag,]
x1=X[,1:p]
y1=X[,p1]
x2=X[,p2:p3]
y2=X[,p4]
}
nk=length(y1)
x1=as.matrix(x1)
x2=as.matrix(x2)
data<-matrix(sample(length(y1),size=length(y1)*nboot,replace=TRUE),nrow=nboot)
bvec1<-apply(data,1,regboot,x1,y1,regfun,...)
# 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<-apply(data,1,regboot,x2,y2,regfun,...)
dif=t(bvec1)-t(bvec2)
est1=lsfit(x1,y1)$coef
est2=lsfit(x2,y2)$coef
estdif=est1-est2
se=apply(dif,2,sd)
pvec=NA
test=NA
test=estdif/se
df=nk-1
pvec=2*(1-pt(abs(test),df))
if(plotit){
reg2plot(x1,y1,x2,y2,xlab=xlab,ylab=ylab)
}
lvec='Intercept'
ci=matrix(NA,nrow=p1,ncol=3)
ci[,1]=c(0:p)
ci[,2]=estdif+qt(alpha/2,df)*se
ci[,3]=estdif-qt(alpha/2,df)*se
dimnames(ci)=list(NULL,c('Param','ci.low','ci.hi'))
for(j in 2:p1)lvec=c(lvec,paste('slope',j-1))
pvec=array(pvec,dimnames=lvec)
list(n=n,n.keep=nk,est.dif=estdif,est.1=est1,est.2=est2,
test.stat=test,standard.error=se,p.values=pvec,conf.intervals=ci)
}