lmrob..M..fit {robustbase}R Documentation

Compute M-estimators of regression

Description

This function performs RWLS iterations to find an M-estimator of regression with Tukey's bisquare psi-function. When started from an S-estimated beta.initial, this results in an MM-estimator.

Usage

lmrob..M..fit(x, y, beta.initial, scale, c.psi, control)

Arguments

x design matrix (n x p) typically including a column of 1s for the intercept.
y numeric response vector (of length n).
beta.initial numeric vector (of length p) of initial estimate. Usually the result of an S-regression estimator.
scale robust residual scale estimate. Usually an S-scale estimator.
c.psi the tuning constant of the psi / psi-function of the M-estimator used.
control list of control parameters, as returned by lmrob.control. Currently, only the components c("max.it", "rel.tol","trace.lev") are accessed.

Details

This function is used by lmrob.fit.MM and typically not to be used on its own.

Value

A list with the following elements:

coef the M-estimator (or MM-estim.) of regression
control the control list input used
scale The residual scale estimate
seed The random number generator seed
converged TRUE if the RWLS iterations converged, FALSE otherwise

Author(s)

Matias Salibian-Barrera and Martin Maechler

References

Yohai, 1987

See Also

lmrob.fit.MM, lmrob; for M-estimators with more flexibility, e.g., choice of rho/psi: rlm from package MASS.

Examples

data(stackloss)
X <- model.matrix(stack.loss ~ . , data = stackloss)
y <- stack.loss
## Compute manual MM-estimate:
## 1) initial LTS:
m0 <- ltsReg(X[,-1], y)
## 2) M-estimate started from LTS:
m1 <- lmrob..M..fit(X, y, beta.initial = coef(m0),
                    scale = m0$scale, c.psi = 1.6,
                    control = lmrob.control())
cbind(m0$coef, m1$coef)
## the scale is kept fixed:
stopifnot(identical(unname(m0$scale), m1$scale))

##  robustness weights: are
r.s <- with(m1, resid/scale) # scaled residuals
m1.wts <- robustbase:::tukeyPsi1(r.s, cc = 1.6) / r.s
summarizeRobWeights(m1.wts)
##--> outliers 1,3,4,13,21
which(m0$lts.wt == 0) # 1,3,4,21 but not 13


[Package robustbase version 0.4-3 Index]