\HeaderA{rlm}{Robust Fitting of Linear Models}{rlm}
\aliasA{predict.rlm}{rlm}{predict.rlm}
\aliasA{print.rlm}{rlm}{print.rlm}
\aliasA{psi.bisquare}{rlm}{psi.bisquare}
\aliasA{psi.hampel}{rlm}{psi.hampel}
\aliasA{psi.huber}{rlm}{psi.huber}
\methaliasA{rlm.default}{rlm}{rlm.default}
\methaliasA{rlm.formula}{rlm}{rlm.formula}
\keyword{models}{rlm}
\keyword{robust}{rlm}
\begin{Description}\relax
Fit a linear model by robust regression using an M estimator.
\end{Description}
\begin{Usage}
\begin{verbatim}
rlm(x, ...)

## S3 method for class 'formula':
rlm(formula, data, weights, ..., subset, na.action,
    method = c("M", "MM", "model.frame"),
    wt.method = c("inv.var", "case"),
    model = TRUE, x.ret = TRUE, y.ret = FALSE, contrasts = NULL)

## Default S3 method:
rlm(x, y, weights, ..., w = rep(1, nrow(x)),
    init = "ls", psi = psi.huber,
    scale.est = c("MAD", "Huber", "proposal 2"), k2 = 1.345,
    method = c("M", "MM"), wt.method = c("inv.var", "case"),
    maxit = 20, acc = 1e-4, test.vec = "resid", lqs.control = NULL)

psi.huber(u, k = 1.345, deriv = 0)
psi.hampel(u, a = 2, b = 4, c = 8, deriv = 0)
psi.bisquare(u, c = 4.685, deriv = 0)
\end{verbatim}
\end{Usage}
\begin{Arguments}
\begin{ldescription}
\item[\code{formula}] a formula of the form \code{y \textasciitilde{} x1 + x2 + ...}.

\item[\code{data}] data frame from which variables specified in \code{formula} are
preferentially to be taken.

\item[\code{weights}] a vector of prior weights for each case.

\item[\code{subset}] An index vector specifying the cases to be used in fitting.

\item[\code{na.action}] A function to specify the action to be taken if \code{NA}s are found.
The \sQuote{factory-fresh} default action in \R{} is
\code{\LinkA{na.omit}{na.omit}}, and can be changed by
\code{\LinkA{options}{options}(na.action=)}.

\item[\code{x}] a matrix or data frame containing the explanatory variables.

\item[\code{y}] the response: a vector of length the number of rows of \code{x}.

\item[\code{method}] currently either M-estimation or MM-estimation or (for the
\code{formula} method only) find the model frame.  MM-estimation
is M-estimation with Tukey's biweight initialized by a specific
S-estimator.  See the \sQuote{Details} section.

\item[\code{wt.method}] are the weights case weights (giving the relative importance of case,
so a weight of 2 means there are two of these) or the inverse of the
variances, so a weight of two means this error is half as variable?

\item[\code{model}] should the model frame be returned in the object?

\item[\code{x.ret}] should the model matrix be returned in the object?

\item[\code{y.ret}] should the response be returned in the object?

\item[\code{contrasts}] optional contrast specifications: se \code{\LinkA{lm}{lm}}.

\item[\code{w}] (optional) initial down-weighting for each case.

\item[\code{init}] (optional) initial values for the coefficients OR a method to find
initial values OR the result of a fit with a \code{coef} component.  Known
methods are \code{"ls"} (the default) for an initial least-squares fit
using weights \code{w*weights}, and \code{"lts"} for an unweighted
least-trimmed squares fit with 200 samples.

\item[\code{psi}] the psi function is specified by this argument.  It must give
(possibly by name) a function \code{g(x, ..., deriv)} that for
\code{deriv=0} returns psi(x)/x and for \code{deriv=1} returns
psi'(x).  Tuning constants will be passed in via \code{...}.

\item[\code{scale.est}] method of scale estimation: re-scaled MAD of the residuals (default)
or Huber's proposal 2 (which can be selected by either \code{"Huber"}
or \code{"proposal 2"}).

\item[\code{k2}] tuning constant used for Huber proposal 2 scale estimation.

\item[\code{maxit}] the limit on the number of IWLS iterations.

\item[\code{acc}] the accuracy for the stopping criterion.

\item[\code{test.vec}] the stopping criterion is based on changes in this vector.

\item[\code{...}] additional arguments to be passed to \code{rlm.default} or to the \code{psi}
function.

\item[\code{lqs.control}] An optional list of control values for \code{\LinkA{lqs}{lqs}}.

\item[\code{u}] numeric vector of evaluation points.

\item[\code{k, a, b, c}] tuning constants.

\item[\code{deriv}] \code{0} or \code{1}: compute values of the psi function or of its
first derivative.

\end{ldescription}
\end{Arguments}
\begin{Details}\relax
Fitting is done by iterated re-weighted least squares (IWLS).

Psi functions are supplied for the Huber, Hampel and Tukey bisquare
proposals as \code{psi.huber}, \code{psi.hampel} and
\code{psi.bisquare}. Huber's corresponds to a convex optimization
problem and gives a unique solution (up to collinearity). The other
two will have multiple local minima, and a good starting point is
desirable.

Selecting \code{method = "MM"} selects a specific set of options which
ensures that the estimator has a high breakdown point. The initial set
of coefficients and the final scale are selected by an S-estimator
with \code{k0 = 1.548}; this gives (for \eqn{n \gg p}{})
breakdown point 0.5.
The final estimator is an M-estimator with Tukey's biweight and fixed
scale that will inherit this breakdown point provided \code{c > k0};
this is true for the default value of \code{c} that corresponds to
95\% relative efficiency at the normal.  Case weights are not
supported for \code{method = "MM"}.
\end{Details}
\begin{Value}
An object of class \code{"rlm"} inheriting from \code{"lm"}.
Note that the \code{df.residual} component is deliberately set to
\code{NA} to avoid inappropriate estimation of the residual scale from
the residual mean square by \code{"lm"} methods.

The additional components not in an \code{lm} object are

\begin{ldescription}
\item[\code{s}] the robust scale estimate used

\item[\code{w}] the weights used in the IWLS process

\item[\code{psi}] the psi function with parameters substituted

\item[\code{conv}] the convergence criteria at each iteration

\item[\code{converged}] did the IWLS converge?

\item[\code{wresid}] a working residual, weighted for \code{"inv.var"} weights only.

\end{ldescription}
\end{Value}
\begin{References}\relax
P. J. Huber (1981)
\emph{Robust Statistics}.
Wiley.

F. R. Hampel, E. M. Ronchetti, P. J. Rousseeuw and W. A. Stahel (1986)
\emph{Robust Statistics: The Approach based on Influence Functions}.
Wiley.

A. Marazzi (1993)
\emph{Algorithms, Routines and S Functions for Robust Statistics}.
Wadsworth \& Brooks/Cole.

Venables, W. N. and Ripley, B. D. (2002)
\emph{Modern Applied Statistics with S.}  Fourth edition.  Springer.
\end{References}
\begin{SeeAlso}\relax
\code{\LinkA{lm}{lm}}, \code{\LinkA{lqs}{lqs}}.
\end{SeeAlso}
\begin{Examples}
\begin{ExampleCode}
summary(rlm(stack.loss ~ ., stackloss))
rlm(stack.loss ~ ., stackloss, psi = psi.hampel, init = "lts")
rlm(stack.loss ~ ., stackloss, psi = psi.bisquare)
\end{ExampleCode}
\end{Examples}

