rlm(x, …)# S3 method for 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)
# S3 method for default
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)
y ~ x1 + x2 + ….
formula are
preferentially to be taken.
x.
formula method only) find the model frame. MM-estimation
is M-estimation with Tukey's biweight initialized by a specific
S-estimator. See the ‘Details’ section.
lm.
coef component. Known
methods are "ls" (the default) for an initial least-squares fit
using weights w*weights, and "lts" for an unweighted
least-trimmed squares fit with 200 samples.
g(x, …, deriv) that for
deriv=0 returns psi(x)/x and for deriv=1 returns
psi'(x). Tuning constants will be passed in via ….
"Huber"
or "proposal 2").
rlm.default or to the psi
function.
lqs.
0 or 1: compute values of the psi function or of its
first derivative.
"rlm" inheriting from "lm".
Note that the df.residual component is deliberately set to
NA to avoid inappropriate estimation of the residual scale from
the residual mean square by "lm" methods. The additional components not in an lm object are "inv.var" weights only.
psi.huber, psi.hampel and
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 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 k0 = 1.548; this gives (for \(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 c > k0;
this is true for the default value of c that corresponds to
95% relative efficiency at the normal. Case weights are not
supported for method = "MM".lm, lqs.summary(rlm(stack.loss ~ ., stackloss))
rlm(stack.loss ~ ., stackloss, psi = psi.hampel, init = "lts")
rlm(stack.loss ~ ., stackloss, psi = psi.bisquare)
Run the code above in your browser using DataLab