wilcox_R: What the package does (short line)

opregpbMC <-
function(x,y,nboot=1000,alpha=.05,om=TRUE,ADJ=TRUE,cop=3,
nullvec=rep(0,ncol(x)+1),plotit=TRUE,opdis=2,gval=sqrt(qchisq(.95,ncol(x)+1))){
#
#  Same as opregpb, only this function takes advantage of a multi-core
#  processor assuming one is availabe and that the R package
#  multicore has been installed.
#
# generate bootstrap estimates
# use projection-type outlier detection method followed by
# TS regression.
#
# om=T and ncol(x)>1, means an omnibus test is performed,
# otherwise only individual tests of parameters are performed.
#
# opdis=2, means that Mahalanobis distance is used
# opdis=1, means projection-type distance is used
#
# gval is critical value for projection-type outlier detection
# method
#
# ADJ=T, Adjust p-values as described in Section 11.1.5 of the text.
#
library(parallel)
x<-as.matrix(x)
m<-cbind(x,y)
p1<-ncol(x)+1
m<-elimna(m) # eliminate any rows with missing data
x<-m[,1:ncol(x)]
x<-as.matrix(x)
y<-m[,p1]
if(nrow(x)!=length(y))stop("Sample size of x differs from sample size of y")
if(!is.matrix(x))stop("Data should be stored in a matrix")
print("Taking bootstrap samples. Please wait.")
data<-matrix(sample(length(y),size=length(y)*nboot,replace=TRUE),nrow=nboot)
bvec<-apply(data,1,regboot,x,y,regfun=opregMC)
# 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.
# using Hochberg method
bvec<-t(bvec)
dvec<-alpha/(c(1:ncol(x)))
test<-NA
icl0<-round(alpha*nboot/2)
icl<-round(alpha*nboot/(2*ncol(x)))
icu0<-nboot-icl0
icu<-nboot-icl
output<-matrix(0,p1,6)
dimnames(output)<-list(NULL,c("Param.","sig.test","sig.crit",
"ci.lower","ci.upper","s.e."))
pval<-NA
for(i in 1:p1){
output[i,1]<-i-1
se.val<-var(bvec[,i])
temp<-sort(bvec[,i])
output[i,6]<-sqrt(se.val)
if(i==1){
output[i,4]<-temp[icl0+1]
output[i,5]<-temp[icu0]
}
if(i>1){
output[i,4]<-temp[icl+1]
output[i,5]<-temp[icu]
}
pval[i]<-sum((temp>nullvec[i]))/length(temp)
if(pval[i]>.5)pval[i]<-1-pval[i]
}
fac<-2
if(ADJ){
# Adjust p-value if n<60
nval<-length(y)
if(nval<20)nval<-20
if(nval>60)nval<-60
fac<-2-(60-nval)/40
}
pval[1]<-2*pval[1]
pval[2:p1]<-fac*pval[2:p1]
output[,2]<-pval
temp2<-order(0-pval[2:p1])
zvec<-dvec[1:ncol(x)]
sigvec<-(test[temp2]>=zvec)
output[temp2+1,3]<-zvec
output[1,3]<-NA
output[,2]<-pval
om.pval<-NA
temp<-opregMC(x,y)$coef
if(om && ncol(x)>1){
temp2<-rbind(bvec[,2:p1],nullvec[2:p1])
if(opdis==1)dis<-pdisMC(temp2,pr=FALSE,center=temp[2:p1])
if(opdis==2){
cmat<-var(bvec[,2:p1]-apply(bvec[,2:p1],2,mean)+temp[2:p1])
dis<-mahalanobis(temp2,temp[2:p1],cmat)
}
om.pval<-sum((dis[nboot+1]<=dis[1:nboot]))/nboot
}
# do adjusted p-value
nval<-length(y)
if(nval<20)nval<-20
if(nval>60)nval<-60
adj.pval<-om.pval/2+(om.pval-om.pval/2)*(nval-20)/40
if(ncol(x)==2 && plotit){
plot(bvec[,2],bvec[,3],xlab="Slope 1",ylab="Slope 2")
temp.dis<-order(dis[1:nboot])
ic<-round((1-alpha)*nboot)
xx<-bvec[temp.dis[1:ic],2:3]
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
lines(xx[temp,])
lines(xx[c(temp[1],temp[length(temp)]),])
}
list(output=output,om.pval=om.pval,adj.om.pval=adj.pval)
}