wilcox_R: What the package does (short line)

ols1wayISO <-
function(x,y,xout=FALSE,outfun=outpro,STAND=TRUE,alpha=.05,pr=TRUE,BLO=FALSE,...){
#
#  Test hypothesis that for two or more independent groups, all slope parameters
#  are equal using OLS  estimator.
#
#  Use Johansen MANOVA 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')
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)
gmean=rep(0,p)
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]])$cov
ecovinv[[j]]=solve(ecov[[j]])  #W_j
gmean=gmean+ecovinv[[j]][2:K,2:K]%*%est[j,2:K]
W=W+ecovinv[[j]]
}
estall=solve(W[2:K,2:K])%*%gmean
estall=c(0,estall)
F=0
for(k in 2:K){
for(m in 2:K){
for(j in 1:J){
F=F+ecovinv[[j]][k,m]*(est[j,k]-estall[k])*(est[j,m]-estall[m])
}}}
pvalad=NULL
df=p*(J-1)
# Adjust critical value:
iden=diag(p)
Aw=0
for(j in 1:J){
temp=iden-solve(W[2:K,2:K])%*%ecovinv[[j]][2:K,2:K]
tempsq=temp%*%temp
Aw=Aw+(sum(diag(tempsq))+(sum(diag(temp)))^2)/(nv[j]-1)
}
Aw=Aw/2
crit=qchisq(alpha,df)
crit=crit+(crit/(2*df))*(Aw+3*Aw*crit/(df+2))
alval<-c(1:999)/1000
for(i in 1:999){
irem<-i
crit=qchisq(alval[i],df)
critad=crit+(crit/(2*df))*(Aw+3*Aw*crit/(df+2))
if(F<critad)break
pvalad=1-irem/1000
}
#
pval=1-pchisq(F,df)
crit=qchisq(1-alpha,df)  #UNADJUSTED CRITICAL VALUE
critad=NULL
critad=crit+(crit/(2*df))*(Aw+3*Aw*crit/(df+2))
est=data.frame(est)
list(n=nv,n.keep=nv.keep,test.stat=F,crit.value=crit,adjusted.crit=critad,p.value=pval,adjusted.p.value=pvalad,est=est)
}