3
$\begingroup$

Could you please point me to tutorial/notes that can help me understand "rlm" better? Here is an example: summary(rlm(stack.loss ~ ., stackloss, weights=myweights))

My questions are:

a. How does the iterative procedure work?
b. What's the relation between the weights I supplied above "myweights" and the final used weights?

My understanding is that the final weights used are actually the Huber loss weights? So are myweights not used at all? Ultimately, I want to use myweights... I have a feeling that actually I should put "myweights" into the "w" argument not the "weights" argument...

c. Could you please give an example showing whether the final used weights are actually the Huber weights or my weights?

I am also confused about the manual:

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.

The returned values:

w the weights used in the IWLS process

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

Anybody please shed some light?

Thank you!

|cite|improve this question
$\endgroup$
  • $\begingroup$ What do you mean by the "Huber loss weights"? If you're user the Huber $M$-estimator (which I believe is the default in rlm()), then you're optimizing the function seen in this answer - stats.stackexchange.com/questions/29563/…. $\endgroup$ – Macro Jul 17 '12 at 14:44
  • $\begingroup$ Thanks a lot Macro for that very informative pointer! However my question is a bit different, more on the operational side and more on the "weights" and "w" arguments ... Thank you! $\endgroup$ – Luna Jul 17 '12 at 15:16
  • $\begingroup$ I also wanted to take one Y observation as an example, and see what weight it had applied and converged to, etc. $\endgroup$ – Luna Jul 17 '12 at 15:19
  • $\begingroup$ Is this the rlm function in the MASS package for R? (that is what it looks like, but there could easily be other functions with that name). If that is the case then the whole MASS package is a suport package to go along with a book (who's intials happen to be 'MASS'). So that book would seem to be the obvious place to start when looking for tutorials or notes to help with understanding the function and the science behind it. $\endgroup$ – Greg Snow Jul 17 '12 at 20:50
  • $\begingroup$ I actually read the book and still couldn't figure out the w's ... Thank you! $\endgroup$ – Luna Jul 17 '12 at 22:27
2
$\begingroup$
  • a. How does the iterative procedure work?

It's best explained here (first two equations). The key elements are the starting $\beta^0\in\mathbb{R}^{p+1}$'s. In this case (rlm), the algorithm uses a numbers ($M$) of starting $\beta^0$'s $\{\beta^{0,1},...,\beta^{0,M}\}$ and carries the iterative procedure on each of them until convergence.

Normally, each of these $\beta^{0,m}$ should be the slope of the regression line passing by $p+1$ observations drawn at random from the $n$ you have (denoting $n\times p$ the dimensions of your design matrix). But in the MASS implementation this is only true if init is set to lts. Using the default init=ls uses as starting point a single ($M=1$) weighted lm fit (with weights w or if none is supplied a single, un-weighted, lm fit). From a robustness perspective the choice of this default behavior is hard to comprehend but consistent with the generally low quality of the implementation of the robust functions in MASS.

--if you are looking for good implementation of robust procedures, i would advise having a look at the dedicated task view--

When $M>1$, the final $\beta^{F,*}$ shown is chosen among all the $\beta^{F,m}$ as the one which has smallest value of the sum of the $n/2$ smallest normalized residuals.

  • b. What's the relation between the weights I supplied above "myweights" and the final used weights?

At each iteration (and therefore also at the final one) the weights used in the "myweights"-weighted IRLS are equal to "myweights" times the normal IRLS weights.

  • c. Could you please give an example showing whether the final used weights are actually the Huber weights or my weights?

Neither: a product of both.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ So which robust lm would you recommend? $\endgroup$ – Luna Jul 17 '12 at 18:33
  • $\begingroup$ Also, could you please explain the relation between the input arguments "w", "weights" and the output values "w" (and wresid)? If I want to have my own weights as used in the lm case, shall I use "w" or "weights"? $\endgroup$ – Luna Jul 17 '12 at 18:35
  • 2
    $\begingroup$ weights is what you want; w is a vector of starting values for the IRLS weights. The ability to use the latter helps sometimes especially with the bisquare psi-function, which doesn't have a unique minimum and can wander off into fail-land. I personally like the bisquare, as it can completely remove extreme outliers. My thought process goes like this: Are my errors likely to have tails more like a t distribution or an exponential distribution? If the former, bisquare; the latter, Huber. (The MLE for a t dist'n location is redescending, like the bisquare.) $\endgroup$ – jbowman Jul 17 '12 at 19:25
  • $\begingroup$ Thanks Jbowman! How does "rlm" decide its w for each IRLS iteration then? Thanks again! $\endgroup$ – Luna Jul 17 '12 at 20:06
  • 1
    $\begingroup$ @Luna: i think this warrants another question... $\endgroup$ – user603 Jul 17 '12 at 22:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.