wilcox_R: What the package does (short line)

olsci <-
function(x,y,alpha=.05,CN=FALSE,xout=FALSE,outfun=outpro,HC3=FALSE,plotit=FALSE,xlab = "X", ylab = "Y", zlab = "Z",...){
#
# Compute confidence intervals via least squares
# regression using heteroscedastic method
# recommended by Cribari-Neto (2004).
# CN=F, degrees of freedom are n-p
# CN=T  degrees of freedom are infinite, as done by Cribari-Neto (2004)
# All indications are that CN=F is best for general use.
#
#  HC3=TRUE, will replace the HC4 estimator with the HC3 estimator.
#
x<-as.matrix(x)
pnum=ncol(x)
if(nrow(x) != length(y))stop("Length of y does not match number of x values")
m<-cbind(x,y)
m<-elimna(m)
y<-m[,ncol(x)+1]
x=m[,1:ncol(x)]
n=length(y)
nrem=n
n.keep=length(y)
x<-as.matrix(x)
if(xout){
flag<-outfun(x,plotit=FALSE,...)$keep
x<-as.matrix(x)
x<-x[flag,]
y<-y[flag]
n.keep=length(y)
x<-as.matrix(x)
}
temp<-lsfit(x,y)
x<-cbind(rep(1,nrow(x)),x)
xtx<-solve(t(x)%*%x)
h<-diag(x%*%xtx%*%t(x))
n<-length(h)
d<-(n*h)/sum(h)
for(i in 1:length(d)){
        d[i]<-min(4, d[i])
}
if(HC3)d=2
hc4<-xtx%*%t(x)%*%diag(temp$res^2/(1-h)^d)%*%x%*%xtx
df<-nrow(x)-ncol(x)
crit<-qt(1-alpha/2,df)
if(CN)crit=qnorm(1-alpha/2)
al<-ncol(x)
p=al-1
ci<-matrix(NA,nrow=al,ncol=6)
lab.out=rep("Slope",p)
dimnames(ci)<-list(c("(Intercept)",lab.out),c("Coef.","Estimates",
"ci.lower","ci.upper","p-value","Std.Error"))
for(j in 1:al){
ci[j,1]<-j-1
ci[j,2]<-temp$coef[j]
ci[j,3]<-temp$coef[j]-crit*sqrt(hc4[j,j])
ci[j,4]<-temp$coef[j]+crit*sqrt(hc4[j,j])
test<-temp$coef[j]/sqrt(hc4[j,j])
ci[j,5]<-2*(1-pt(abs(test),df))
if(CN)ci[j,5]<-2*(1-pnorm(abs(test),df))
}
ci[,6]=sqrt(diag(hc4))
if(plotit){
if(pnum==1){
plot(x[,-1],y,xlab=xlab,ylab=ylab)
abline(ci[,2])
}
if(pnum==2){
regp2plot(x[,-1],y,regfun=ols,xlab=xlab,ylab=ylab,zlab=zlab)
}}
list(n=nrem,n.keep=n.keep,ci=ci, cov=hc4)
}