wilcox_R: What the package does (short line)

chregF <-
function(x,y,bend=1.345,SEED=FALSE,xout=FALSE,outfun=out,...){
#
# Compute Coakley Hettmansperger robust regression estimators
# JASA, 1993, 88, 872-880
#
# x is a n by p matrix containing the predictor values.
#
# No missing values are allowed
#
#  Comments in this function follow the notation used
#  by Coakley and Hettmansperger
#
library(MASS)
# with old version of R, need library(lqs) when using ltsreg
# as the initial estimate.
#
if(SEED)set.seed(12) # Set seed so that results are always duplicated.
x<-as.matrix(x)
p<-ncol(x)
m<-elimna(cbind(x,y))
x<-m[,1:p]
p1<-p+1
y<-m[,p1]
if(xout){
x<-as.matrix(x)
flag<-outfun(x,plotit=plotit,...)$keep
x<-x[flag,]
y<-y[flag]
x<-as.matrix(x)
}
x<-as.matrix(x)
cutoff<-bend
mve<-vector('list')
if(ncol(x)==1){
mve$center<-median(x)
mve$cov<-mad(x)^2
}
if(ncol(x)>=2)mve<-cov.mve(x)  # compute minimum volume ellipsoid measures of
                 # location and scale and store in mve.
reg0<-ltsgreg(x,y) # compute initial regression est using least trimmed
                 # squares.
# Next, compute the rob-md2(i) values and store in rob
rob<-1  # Initialize vector rob
mx<-mve$center
rob<-mahalanobis(x,mx,mve$cov)
k21<-qchisq(.95,p)
c62<-k21/rob
vecone<-c(rep(1,length(y))) # Initialize vector vecone to 1
c30<-pmin(vecone,c62)  # mallows weights put in c30
k81<-median(abs(reg0$residuals)) # median of absolute residuals
k72<-1.4826*(1+(5/(length(y)-p-1)))*k81 # lms scale
c60<-reg0$residuals/(k72*c30) # standardized residuals
#  compute psi and store in c27
cvec<-c(rep(cutoff,length(y))) # Initialize vector cvec to cutoff
c27<-pmin(cvec,c60)
c27<-pmax(-1*cutoff,c27)  #c27 contains psi values
#
# compute B matrix and put in c66.
#  Also, transform B so that i th diag elem = 0 if c27[i] is
#  between -cutoff and cutoff, 1 otherwise.
#
c66<-ifelse(abs(c27)<=bend,1,0) # Have derivative of psi in c66
m1<-cbind(1,x)  # X matrix with col of 1's added
m2<-t(m1)   #X transpose
m5<-diag(c30) # matrix W, diagonal contains weights
m4<-diag(c66) # B matrix
m6<-m4%*%m1   # BX
m7<-m2%*%m6   # X'BX (nD=X'BX)
m8<-solve(m7)  #m8 = (X'-B-X)inverse
m9<-m8%*%m2 #m9=X prime-B-X inverse X'
m9<-m9%*%m5 # m9=X prime-B-X inverse X'W
m10<-m9%*%c27
c20<-m10*k72
c21<-reg0$coef+c20 #update initial estimate of parameters.
res<-y-m1%*%c21
list(coef=t(c21),residuals=res)
}