1\name{lmrob..M..fit}
2\alias{lmrob..M..fit}
3\title{Compute M-estimators of regression}
4\description{
5  This function performs RWLS iterations to find an
6  M-estimator of regression.  When started from an S-estimated
7  \code{beta.initial}, this results in an MM-estimator.
8}
9\usage{
10lmrob..M..fit(x, y, beta.initial, scale, control, obj,
11              mf = obj$model, method = obj$control$method)
12}
13\arguments{
14  \item{x}{design matrix (\eqn{n \times p}{n x p}) typically including a
15    column of \code{1}s for the intercept.}
16  \item{y}{numeric response vector (of length \eqn{n}).}
17  \item{beta.initial}{numeric vector (of length \eqn{p}) of initial
18    estimate.  Usually the result of an S-regression estimator.}
19  \item{scale}{robust residual scale estimate. Usually
20    an S-scale estimator.}
21  \item{control}{list of control parameters, as returned
22    by \code{\link{lmrob.control}}.  Currently, the components
23    \code{c("max.it", "rel.tol","trace.lev", "psi", "tuning.psi", "mts", "subsampling")}
24    are accessed.}
25  \item{obj}{an optional \code{lmrob}-object.  If specified, this is
26    typically used to set values for the other arguments.}
27  \item{mf}{unused and deprecated.}
28  \item{method}{optional; the \code{method} used for \emph{obj} computation.}
29}
30\details{
31  This function is used by \code{\link{lmrob.fit}} (and
32  \code{anova(<lmrob>, type = "Deviance")}) and typically not to be used
33  on its own.
34}
35\value{A list with the following elements:
36  \item{coef}{the M-estimator (or MM-estim.) of regression}
37  \item{control}{the \code{control} list input used}
38  \item{scale}{ The residual scale estimate}
39  \item{seed }{ The random number generator seed}
40  \item{converged}{ \code{TRUE} if the RWLS iterations converged,
41    \code{FALSE} otherwise}
42}
43\references{
44  Yohai, 1987
45}
46\seealso{
47  \code{\link{lmrob.fit}}, \code{\link{lmrob}};
48  \code{\link[MASS]{rlm}} from package \CRANpkg{MASS}.
49}
50\author{Matias Salibian-Barrera and Martin Maechler}
51\examples{
52data(stackloss)
53X <- model.matrix(stack.loss ~ . , data = stackloss)
54y <- stack.loss
55## Compute manual MM-estimate:
56## 1) initial LTS:
57m0 <- ltsReg(X[,-1], y)
58## 2) M-estimate started from LTS:
59m1 <- lmrob..M..fit(X, y, beta.initial = coef(m0), scale = m0$scale, method = "SM",
60                    control = lmrob.control(tuning.psi = 1.6, psi = 'bisquare'))
61## no 'method' (nor 'obj'):
62m1. <- lmrob..M..fit(X, y, beta.initial = coef(m0), scale = m0$scale,
63                     control = m1$control)
64stopifnot(all.equal(m1, m1., tol = 1e-15)) # identical {call *not* stored!}
65
66cbind(m0$coef, m1$coef)
67## the scale is kept fixed:
68stopifnot(identical(unname(m0$scale), m1$scale))
69
70##  robustness weights: are
71r.s <- with(m1, residuals/scale) # scaled residuals
72m1.wts <- Mpsi(r.s, cc = 1.6, psi="tukey") / r.s
73summarizeRobWeights(m1.wts)
74##--> outliers 1,3,4,13,21
75which(m0$lts.wt == 0) # 1,3,4,21 but not 13
76\dontshow{stopifnot(which(m0$lts.wt == 0) == c(1,3,4,21))
77}
78## Manually add M-step to SMD-estimate (=> equivalent to "SMDM"):
79m2 <- lmrob(stack.loss ~ ., data = stackloss, method = 'SMD')
80m3 <- lmrob..M..fit(obj = m2)
81
82## Simple function that allows custom initial estimates
83## (Deprecated; use init argument to lmrob() instead.) %% MM: why deprecated?
84lmrob.custom <- function(x, y, beta.initial, scale, terms) {
85  ## initialize object
86  obj <- list(control = lmrob.control("KS2011"),
87              terms = terms) ## terms is needed for summary()
88  ## M-step
89  obj <- lmrob..M..fit(x, y, beta.initial, scale, obj = obj)
90  ## D-step
91  obj <- lmrob..D..fit(obj, x)
92  ## Add some missing elements
93  obj$cov <- TRUE ## enables calculation of cov matrix
94  obj$p <- obj$qr$rank
95  obj$degree.freedom <- length(y) - obj$p
96  ## M-step
97  obj <- lmrob..M..fit(x, y, obj=obj)
98  obj$control$method <- ".MDM"
99  obj
100}
101
102m4 <- lmrob.custom(X, y, m2$init$init.S$coef,
103                   m2$init$scale, m2$terms)
104stopifnot(all.equal(m4$coef, m3$coef))
105
106## Start from ltsReg:
107m5 <- ltsReg(stack.loss ~ ., data = stackloss)
108m6 <- lmrob.custom(m5$X, m5$Y, coef(m5), m5$scale, m5$terms)
109}
110\keyword{robust}
111\keyword{regression}
112