wilcox_R: What the package does (short line)

DregGOLS <-
function(x1,y1,x2,y2,xout=FALSE,outfun=outpro,SEED=TRUE,nboot=200,
STAND=TRUE,...){
#
#  Global test that two dependent (time 1 and time 2)
#  OLS regression lines are identical
#  
#  Use a variation of Hotelling's test coupled with a bootstrap 
#  estimate of the relevant covariance matrix associated with the differences
#  in the estimates of the parameters.
#
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){
opf=identical(outfun,outpro)
if(!opf){
flag1=outfun(x1)$out.id
flag2=outfun(x2)$out.id
}
if(opf){
flag1=outpro(x1,STAND=STAND)$out.id                               
flag2=outfun(x2,STAND=STAND)$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=lsfit,...)
# 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=lsfit,...)
dif=t(bvec1-bvec2)
S=cov(dif)
est1=lsfit(x1,y1)$coef
est2=lsfit(x2,y2)$coef
est=est1-est2
k <- (nk-p1)/((nk - 1)*p1)
        stat <- k * crossprod(est, solve(S, est))[1, ]
        pvalue <- 1 - pf(stat, p1, nk - p1)
list(test.statistic = stat, degrees_of_freedom = c(p1, nk - p1), p.value =
pvalue,est.1=est1,est.2=est2,estimate.dif = est)
}