| lmrob..M..fit {robustbase} | R Documentation |
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.
lmrob..M..fit(x, y, beta.initial, scale, control, obj, mf = obj$model)
x |
design matrix (n x p) typically including a
column of |
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 |
obj |
an optional |
mf |
(optional) a model frame as returned by |
This function is used by lmrob.fit
and typically not to be used on its own.
A list with the following elements:
coef |
the M-estimator (or MM-estim.) of regression |
control |
the |
scale |
The residual scale estimate |
seed |
The random number generator seed |
converged |
|
Matias Salibian-Barrera and Martin Maechler
Yohai, 1987
lmrob.fit, lmrob;
rlm from package MASS.
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 <- Mpsi(r.s, cc = 1.6, psi="tukey") / 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)
## Simple function that allows custom initial estimates
## (Deprecated; use init argument to lmrob() instead.) %% MM: why deprecated?
lmrob.custom <- function(x, y, beta.initial, scale, terms) {
## initialize object
obj <- list(control = lmrob.control("KS2011"),
terms = terms) ## terms is needed for summary()
## M-step
obj <- lmrob..M..fit(x, y, beta.initial, scale, obj = obj)
## D-step
obj <- lmrob..D..fit(obj, x)
## Add some missing elements
obj$cov <- TRUE ## enables calculation of cov matrix
obj$p <- obj$qr$rank
obj$degree.freedom <- length(y) - obj$p
## M-step
obj <- lmrob..M..fit(x, y, obj=obj)
obj$control$method <- ".MDM"
obj
}
m4 <- lmrob.custom(X, y, m2$init$init.S$coef,
m2$init$scale,m2$terms)
stopifnot(all.equal(m4$coef, m3$coef))
## Start from ltsReg:
m5 <- ltsReg(stack.loss ~ ., data = stackloss)
m6 <- lmrob.custom(m5$X, m5$Y, coef(m5), m5$scale, m5$terms)