wilcox_R: What the package does (short line)

hc4test <-
function(x,y,pval=c(1:ncol(x)),xout=FALSE,outfun=outpro,pr=TRUE,plotit=FALSE,xlab="X",ylab="Y",...){
#
# Perform omnibus test using OLS and HC4 estimator
# That is, test the hypothesis that all of the slope parameters
# are equal to 0 in a manner that allows heteroscedasticity.
#
# recommended by Cribari-Neto (2004).
# Seems to work well with p=1 but can be unsatisfactory wit p>4 predictors,
# Unknown how large n must be when p>1
#
x<-as.matrix(x)
if(ncol(x)>1 && pr)print("WARNING: more than 1 predictor, olstest might be better")
if(nrow(x) != length(y))stop("Length of y does not match number of x values")
m<-cbind(x,y)
m<-elimna(m)
p=ncol(x)
p1=p+1
y<-m[,p1]
x=m[,1:p]
nrem=length(y)
n=length(y)
n.keep=n
x<-as.matrix(x)
if(xout){
flag<-outfun(x,...)$keep
x<-as.matrix(x)
x<-x[flag,]
y<-y[flag]
n.keep=length(y)
x<-as.matrix(x)
}
n=n.keep
pvalp1<-pval+1
temp<-lsfit(x,y) # unrestricted
if(plotit){
if(p==1){
plot(x[,1],y,xlab=xlab,ylab=ylab)
abline(temp$coef)
}}
x<-cbind(rep(1,nrow(x)),x)
hval<-x%*%solve(t(x)%*%x)%*%t(x)
hval<-diag(hval)
hbar<-mean(hval)
delt<-cbind(rep(4,n),hval/hbar)
delt<-apply(delt,1,min)
aval<-(1-hval)^(0-delt)
x2<-x[,pvalp1]
pval<-0-pvalp1
x1<-x[,pval]
df<-length(pval)
x1<-as.matrix(x1)
imat<-diag(1,n)
M1<-imat-x1%*%solve(t(x1)%*%x1)%*%t(x1)
M<-imat-x%*%solve(t(x)%*%x)%*%t(x)
uval<-as.vector(M%*%y)
R2<-M1%*%x2
rtr<-solve(t(R2)%*%R2)
temp2<-aval*uval^2
S<-diag(aval*uval^2)
V<-n*rtr%*%t(R2)%*%S%*%R2%*%rtr
nvec<-as.matrix(temp$coef[pvalp1])
test<-n*t(nvec)%*%solve(V)%*%nvec
test<-test[1,1]
p.value<-1-pchisq(test,df)
list(n=nrem,n.keep=n.keep,test=test,p.value=p.value,coef=temp$coef)
}