wilcox_R: What the package does (short line)

mlrreg <-
function(x,y,cov.fun=cov.mcd,ols.op=TRUE,mcd.op=TRUE,
quantile.used=floor(.75*n),RES=FALSE,...){
#
# Do Multivariate regression, using by default the method
#  in Rousseeuw, Van Aelst, Van Driessen Agullo
# Technometrics, 46, 293-305
#
# Note, to use the method recommended by Rousseeuw et al., the argument
# quantile.used=.75*n is used when calling cov.mcd.
#
#  RES=T, the residuals will be returned.
#
# y is assumed to be  multivariate with data stored in a matrix.
#
# an initial fit is found using the measures of scatter and location
# corresponding to cof.fun and mcd.op. If
# mcd.op=T, cov.mcd is used with quanitle.used=.75n
# mcd.op=F, cov.fun is used and defaults to cov.mcd with the
# default value usded by R for the argument quanitle.used
# But any function that returns location and scatter in $center and $cov
# can be used.
#
#  if ols.op=T, OLS is applied after points are removed based on iniital fit
#  if ols.op=F, Theil-Sen is used by calling the function mopreg
#
#  Early version of this function considered estimating
#  explanatory power in terms of the generalized variance
#  of the predicted y values and the observed y values
#  epow.cov determines which robust covariance matrix will be used.
#  This idea has not been explored enough
#  Some choices are:
# cov (the usual generalized variance)
# skipcov
# tbscov
# covout
# covogk
# mgvcov
# mvecov
# mcdcov
#
library(MASS)
if(!is.matrix(y))stop("y is not a matrix")
X<-cbind(x,y)
X<-elimna(X)
n<-nrow(X)
qy<-ncol(y)
qx<-ncol(x)
qxp1<-qx+1
tqyqx<-qy+qx
y<-X[,qxp1:tqyqx]
# compute initial estimate of slopes and intercept:
if(!mcd.op)locscat<-cov.fun(X,...)
if(mcd.op)locscat<-cov.mcd(X,quan=quantile.used)
sig<-locscat$cov
mu<-locscat$center
sigxx<-sig[1:qx,1:qx]
sigxy<-sig[1:qx,qxp1:tqyqx]
sigyy<-sig[qxp1:tqyqx,qxp1:tqyqx]
Bhat<-solve(sigxx)%*%sigxy
sige<-sigyy-t(Bhat)%*%sigxx%*%Bhat
sige.inv<-solve(sige)
Ahat<-t(mu[qxp1:tqyqx]-t(Bhat)%*%mu[1:qx])
resL<-matrix(nrow=nrow(X),ncol=qy)
for(i in 1:nrow(X))resL[i,]<-y[i,]-t(Bhat)%*%X[i,1:qx]
for(j in 1:qy)resL[,j]<-resL[,j]-Ahat[j]
drL<-NA
for(i in 1:nrow(X))drL[i]<-t(resL[i,])%*%sige.inv%*%resL[i,]
# In Rousseeuw notation, drL<- is d^2
w<-rep(0,nrow(X))
qdr<-qchisq(.99,qy)
iflag<-(drL<qdr)
w[iflag]<-1
term1<-0
vec<-c(1:nrow(X))
keep<-vec[iflag==1]
X<-X[keep,]
if(ols.op)output<-lsfit(X[,1:qx],X[,qxp1:tqyqx])
if(!ols.op)output<-mopreg(X[,1:qx],X[,qxp1:tqyqx],KEEP=T)
yhat=X[,qxp1:tqyqx]-output$residuals
res=NULL
if(RES)res=output$residuals
#epow=(gvarg(yhat,epow.cov)/gvarg(X[,qxp1:tqyqx],epow.cov))
#list(coef=output$coefficients,residuals=res,E.power=epow,Strength.Assoc=sqrt(epow))
list(coef=output$coefficients,residuals=res)
}