Compute M-estimators of regression

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 = obj$x, y = obj$y,
              beta.initial = obj$coef, scale = obj$scale, control = obj$control,
              obj,
              mf,
              method = obj$control$method)

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, the components c("max.it", "rel.tol","trace.lev", "psi", "tuning.psi", "mts", "subsampling") are accessed.

obj

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

mf

defunct.

method

optional; the method used for obj computation.

Details

This function is used by lmrob.fit (and anova(<lmrob>, type = "Deviance")) 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

References

Yohai, 1987

See also

lmrob.fit, lmrob; rlm from package MASS.

Author

Matias Salibian-Barrera and Martin Maechler

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, method = "SM",
                    control = lmrob.control(tuning.psi = 1.6, psi = 'bisquare'))
## no 'method' (nor 'obj'):
m1. <- lmrob..M..fit(X, y, beta.initial = coef(m0), scale = m0$scale,
                     control = m1$control)
stopifnot(all.equal(m1, m1., tol = 1e-15)) # identical {call *not* stored!}

cbind(m0$coef, m1$coef)
#>                    [,1]         [,2]
#> Intercept  -37.65245890 -36.79032457
#> Air.Flow     0.79768556   0.84937712
#> Water.Temp   0.57734046   0.43210319
#> Acid.Conc.  -0.06706018  -0.07538038
## 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)
#> Robustness weights: 
#>  5 observations c(1,3,4,13,21)
#> 	 are outliers with |weight| <= 9.685e-06 ( < 0.004762); 
#>  2 weights are ~= 1. The remaining 14 ones are summarized as
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#> 0.4499621 0.7351054 0.8729619 0.8269551 0.9560550 0.9927093 
##--> outliers 1,3,4,13,21
which(m0$lts.wt == 0) # 1,3,4,21 but not 13
#> [1]  1  3  4 21
## Manually add M-step to SMD-estimate (=> equivalent to "SMDM"):
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)