robustbase (version 0.7-6)

lmrob..M..fit: Compute M-estimators of regression

Description

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

Usage

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

Arguments

x

design matrix ($n \times 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.

control

list of control parameters, as returned by lmrob.control. Currently, only the components c("max.it", "rel.tol","trace.lev", "psi", "tuning.psi", "method") are accessed.

obj

an optional lmrob-object. If given, this is used to set values for the other arguments.

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

Details

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

References

Yohai, 1987

See Also

lmrob.fit, lmrob; 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,
                    control = lmrob.control(tuning.psi = 1.6,
		                            psi = 'bisquare'))
cbind(m0$coef, m1$coef)
## the scale is kept fixed:
stopifnot(identical(unname(m0$scale), m1$scale))

##  robustness weights: are
r.s <- with(m1, residuals/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


## Add M-step to SMD-estimate
m2 <- lmrob(stack.loss ~ ., data = stackloss, method = 'SMD')
m3 <- lmrob..M..fit(obj = m2)