Skip to content

[RsR] How does "rlm" in R decide its "w" weights for each IRLS iteration?

9 messages · Michael, Valentin Todorov, S Ellison +3 more

#
Hi all,

I am also confused about the manual:

           a.  The input arguments:

wt.method are the weights case weights (giving the relative importance of
case, so a weight of 2 means there are two of these) or the inverse of the
variances, so a weight of two means this error is half as variable?

w (optional) initial down-weighting for each case.

init (optional) initial values for the coefficients OR a method to find
initial values OR the result of a fit with a coef component. Known methods
are "ls" (the default) for an initial least-squares fit using weights
w*weights, and "lts" for an unweighted least-trimmed squares fit with 200
samples.



b. The returned values:

w the weights used in the IWLS process

wresid a working residual, weighted for "inv.var" weights only.



How to use these input arguments?

Anybody please shed some light?

Also, is my understanding below correct?

The input argument "w" is used for the initial values of the rlm IRLS
weighting and the output value "w" is the converged "w".

The "weights" input argument is actually what I want to apply - if I have
my own set of weights.
And the real/actual weights are the product of "weights"(I supplied) and
the converged output "w" (an output).


If my understanding above is correct, how does "rlm" decide its "w" for
each IRLS iteration then?

Any pointers/tutorials/notes to the calculation of these "w"'s in each IRLS
iteration?

Thanks to all.
1 day later
#
When you give rlm weights (called 'weights', not 'w' on input, though you can abbreviate to 'w'), you need to tell it which of these two possibilities you used. 
If you gave it case numbers, say wt.method="case"; if you gave it inverse variance weights, say wt.method="inv.var".
The default is "inv.var".
There is no input argument 'w' for rlm (see above). 
The output w are a  calculated using the psi function, so between 0 and 1.
The effective weights for the final estimate would then be something like w*weights, using the full name of the input argument (and if I haven't forgotten a square root somewhere). At least, that would work for a simple location estimate (eg rlm(x~1)).
It uses the given psi functions to calculate the iterative weights based on the scaled residuals.
Read the cited references for a detailed guide. Or, of course, MASS - the package is, after all, intended to support the book, not replace it.



S Ellison

*******************************************************************
This email and any attachments are confidential. Any use...{{dropped:8}}
#
Hi Michael, S Ellison,

I do not actually understand what you want to achieve with the M
estimates of rlm in MASS, but why you do not give a try of lmrob in
'robustbase'. Please have a llok in the references (?lmrob) about the
advantages of MM estimators over the M estimators.

Best regards,
Valentin
On Fri, Jul 20, 2012 at 5:11 PM, S Ellison <S.Ellison at lgcgroup.com> wrote:
2 days later
#
Hi Valentin,

If the contamination is mainly in the response direction, M-estimator
provides good estimates for parameters and rlm can be used.

Rohan

-----Original Message-----
From: r-sig-robust-bounces at r-project.org
[mailto:r-sig-robust-bounces at r-project.org] On Behalf Of Valentin
Todorov
Sent: Saturday, 21 July 2012 6:57 a.m.
To: S Ellison
Cc: r-sig-robust at r-project.org; r-help
Subject: Re: [RsR] How does "rlm" in R decide its "w" weights for each
IRLSiteration?

Hi Michael, S Ellison,

I do not actually understand what you want to achieve with the M
estimates of rlm in MASS, but why you do not give a try of lmrob in
'robustbase'. Please have a llok in the references (?lmrob) about the
advantages of MM estimators over the M estimators.

Best regards,
Valentin




On Fri, Jul 20, 2012 at 5:11 PM, S Ellison <S.Ellison at lgcgroup.com>
wrote:
you can abbreviate to 'w'), you need to tell it which of these two
possibilities you used.
inverse variance weights, say wt.method="inv.var".
and 1.
like w*weights, using the full name of the input argument (and if I
haven't forgotten a square root somewhere). At least, that would work
for a simple location estimate (eg rlm(x~1)).
based on the scaled residuals.
the package is, after all, intended to support the book, not replace it.
#
rlm includes an MM estimator.

S Ellison
*******************************************************************
This email and any attachments are confidential. Any use...{{dropped:8}}
#
To what it's worth:

Valentin very well knows about the MM implementation of rlm;

it is simply that lmrob is a more recent implementation which has
taken up _very_ recent current research results on this field and
also includes the combined case (i.e. factors  and numerics as
regressors) (see ?lmrob) or try
Error: lqs failed: all the samples were singular
Call:
lmrob(formula = form, method = "MM")

Coefficients:
(Intercept)           xn          xf1          xf2          xf3         
xf4 
     0.2260       0.8957       0.3633       1.9949       2.6660      
3.2966 
        xf5          xf6          xf7          xf8 
     4.9415       5.9208       6.4388       9.8285 

Hth, best, Peter

Am 23.07.2012 15:19, schrieb S Ellison:
#
Your code is running for M-estimator, but not for MM (I don't know the
reason).
See the results here.
Error: lqs failed: all the samples were singular
Call:
rlm(formula = form, method = "M")
Converged in 7 iterations

Coefficients:
(Intercept)          xn         xf1         xf2         xf3         xf4 
  0.2017169   0.8708498   0.4001550   2.0120847   2.7285940   3.2854743 
        xf5         xf6         xf7         xf8 
  4.9677143   6.0077015   6.4267159   9.8351416 

Degrees of freedom: 100 total; 90 residual
Scale estimate: 0.713

Regards
Rohan




-----Original Message-----
From: Dr. Peter Ruckdeschel
[mailto:peter.ruckdeschel at itwm.fraunhofer.de] 
Sent: Tuesday, 24 July 2012 2:07 a.m.
To: r-sig-robust at r-project.org
Cc: Valentin Todorov; Maheswaran Rohan; S.Ellison at lgcgroup.com
Subject: Re: [RsR] How does "rlm" in R decide its "w" weights for each
IRLSiteration?

To what it's worth:

Valentin very well knows about the MM implementation of rlm;

it is simply that lmrob is a more recent implementation which has
taken up _very_ recent current research results on this field and
also includes the combined case (i.e. factors  and numerics as
regressors) (see ?lmrob) or try

set.seed(123)
 xr=rpois(100,lambda=3)
 xn=rnorm(100)
 xf =as.factor(xr)
 yn=rnorm(100)+xn+xr
 form=yn~xn+xf
 rlm(form,method="MM")
Error: lqs failed: all the samples were singular
Call:
lmrob(formula = form, method = "MM")

Coefficients:
(Intercept)           xn          xf1          xf2          xf3         
xf4 
     0.2260       0.8957       0.3633       1.9949       2.6660      
3.2966 
        xf5          xf6          xf7          xf8 
     4.9415       5.9208       6.4388       9.8285 

Hth, best, Peter

Am 23.07.2012 15:19, schrieb S Ellison:
#
I am aware of the lmrob implementation; I've been following it with considerable interest. 
As a matter of fact I see robustbase all round as a better implementation of robust statistics than MASS, and lmrob as a good example of a specialised package with sound research overtaking the more general teaching package that MASS represents. It was just the brevity of Valentin's post - which pointed specifically to MM-estimates' performance - that led me to note for the OP's benefit that an MM estimate was also included in MASS. 

But thanks for the thought.

Steve Ellison,
Lab of the Government Chemist
Teddington
UK
#
Dear all,

We've put quite some work into lmrob lately. There are some
enhancements for small samples, where we advocate the SMDM estimate
which gives a better scale estimate, important when performing tests
(use setting="KS2011").

Also lmrob does not fail for binary/categorical designs like the
example Rohan posted below. Robustbase 0.9-3 features the "nonsingular
subsampling" (also called "contstrained s." during development)
algorithm which makes the distinction between categorical and
continuous predictors obsolete. The idea is, that it just works,
whatever your data might be. Note that the version 0.9-2 of robustbase
available on CRAN does most of it as well, we mainly improved
numerical properties in 0.9-3.

Best regards,

Manuel
On Mon, Jul 23, 2012 at 10:51 PM, Maheswaran Rohan <mrohan at doc.govt.nz> wrote: