wilcox_R: What the package does (short line)

olsW2g <-
function(x,y,grp=c(1,2),iv=1,xout=FALSE,outfun=outpro,STAND=TRUE,alpha=.05,
pr=TRUE,BLO=FALSE,HC3=FALSE,...){
#
#  Test hypothesis that for two or more independent groups, all slope parameters
#  are equal using OLS  estimator.
#
#  Use Welch type test statistic in conjunction with HC4 estimate of covariances.
#
#  x and y are assumed to have list mode having length J equal to the number of groups
#  For example, x[[1]] and y[[1]] contain the data for group 1.
#
#   xout=T will eliminate leverage points using the function outfun, 
#   which defaults to the MVE method.
#
#   BLO=TRUE, only bad leverage points are removed.
#
#  OUTPUT: 
#   n is sample size after missing values are removed
#   nv.keep is sample size after leverage points are removed.
#
if(!is.list(x))stop('Argument x should have list mode')
iv1=iv+1
J=length(x) # number of groups
x=lapply(x,as.matrix)
pchk=lapply(x,ncol)
temp=matl(pchk)
if(var(as.vector(temp))!=0)stop('Something is wrong. Number of covariates differs among the groups being compared')
nv=NULL
p=ncol(x[[1]])
p1=p+1
for(j in 1:J){
xy=elimna(cbind(x[[j]],y[[j]]))
x[[j]]=xy[,1:p]
y[[j]]=xy[,p1]
x[[j]]=as.matrix(x[[j]])
nv=c(nv,nrow(x[[j]]))
}
nv.keep=nv
critrad=NULL
if(xout){
temp=lapply(x,outfun,plotit=FALSE,STAND=STAND,...)
for(j in 1:J){
x[[j]]=x[[j]][temp[[j]]$keep,]
y[[j]]=y[[j]][temp[[j]]$keep]
}
if(BLO){
for(j in 1:J){
temp=reglev(x[[j]],y[[j]],plotit=FALSE)
ad1=c(temp1[[j]]$out.id,temp$regout)
flag1=duplicated(ad1)
if(sum(flag1)>0){
flag1=ad1[flag1]
x[[j]]=as.matrix(x[[j]])
x[[j]]=x[[j]][-flag1,]
y[[j]]=y[[j]][-flag1]
}}}}
x=lapply(x,as.matrix)
K=p1
est=matrix(NA,nrow=J,ncol=p1)
grpnum=NULL
for(j in 1:J)grpnum[j]=paste("Group",j)
vlabs="Intercept"
for(j in 2:p1)vlabs[j]=paste("Slope",j-1)
dimnames(est)=list(grpnum,vlabs)
ecov=list()
ecovinv=list()
W=rep(0,p1)
for(j in 1:J){
est[j,]=ols(x[[j]],y[[j]],xout=FALSE,plotit=FALSE,...)$coef
nv.keep[j]=nrow(x[[j]])
ecov[[j]]=olshc4(x[[j]],y[[j]],HC3=HC3)$cov
}
q1=ecov[[grp[1]]][iv1,iv1]
q2=ecov[[grp[2]]][iv1,iv1]
top=est[grp[1]]-est[grp[2]]
F=(est[grp[1],iv1]-est[grp[2],iv1])/sqrt(q1+q2)
df=(q1+q2)^2/(q1^2/(nv[grp[1]]-1)+q2^2/(nv[grp[2]]-1))
pv=2*(1-pt(abs(F),df))
crit=qt(1-alpha/2,df)
ci=est[grp[1],iv1]-est[grp[2],iv1]-crit*sqrt(q1+q2)
ci[2]=est[grp[1],iv1]-est[grp[2],iv1]+crit*sqrt(q1+q2)
list(n=nv,n.keep=nv.keep,test.stat=F,conf.interval=ci,
est=c(est[grp[1],iv1],est[grp[2],iv1]),est.dif=est[grp[1],iv1]-est[grp[2],iv1],p.value=pv)
}