8
$\begingroup$

Context. I'd like to fit a regression line to study to relation between some response variable $y$ and some continuous covariate $x$. Because of the presence of bad leverage points, I have opted for an MM-estimator instead of the usual LS-estimator.

Methodology. Basically, MM-estimation is M-estimation initialised by an S-estimator. Hence, two loss functions have to be picked. I have chosen the widely used Tukey Biweight's loss function

$\rho ( u ) = \left\{ \begin{array}{ll} 1 - \left[ 1 - \left( \tfrac{u}{k} \right)^{2} \right]^{3} & \textrm{if } | u | \leq k \\ 1 & \textrm{if } | u | > k, \end{array} \right.$

with $k = 1.548$ at the preliminary S-estimator (which gives a breakdown point equal to $50 \%$), and with $k = 2.697$ at the M-estimation step (to guarantee $70\%$ Gaussian efficiency).

I'd like to use R to fit my robust regression line.

Question.

library(MASS)
rlm(y~x, 
    method="MM",
    k0=1.548, c=2.697,
    maxit=50)
  • Is my code consistent with the previous paragraph?
  • Would you use other optional arguments?

EDIT. Following my discussion with @Jason Morgan, I realise that my previous code is wrong. (@Jason Morgan: Thank you very much for this!) However, I am still not convinced by his proposal. Instead, here is what I propose now:

library(robustbase)
lmrob(y~x, 
      tuning.chi=1.548, tuning.psi=2.697)

I think it sticks to the methodology now. Do you agree?

Thanks!

$\endgroup$
5
$\begingroup$

By default, the documentation indicates that rlm uses psi=psi.huber weights. Thus, if you want to use Tukey's bisquare, you need to specify psi=psi.bisquare. The default settings are psi.bisquare(u, c = 4.685, deriv = 0), which you can change as desired. For instance, possibly something like

rlm(x ~ y, method="MM", psi=psi.bisquare, maxit=50)

You may also want to investigate whether you should use least-trimmed squares (init="lts") to initialize your starting values. The default is to use least squares.

$\endgroup$
3
  • $\begingroup$ @Janson Morgan: are you sure of what you put forward? Do you have any experience with that function? My documentation (R 2.13.1) actually indicates "The initial set of coefficients and the final scale are selected by an S-estimator with k0 = 1.548; this gives (for n >> p) breakdown point 0.5. The final estimator is an M-estimator with Tukey's biweight and fixed scale that will inherit this breakdown point provided c > k0; this is true for the default value of c that corresponds to 95% relative efficiency at the normal." $\endgroup$ – ocram Nov 1 '11 at 16:03
  • 1
    $\begingroup$ I have estimated these models in the past. As the documentation states, the first step in the MM estimation is done with Huber weights, the second with bisquared weights. My notes (from a couple years ago) state that in the first S-step, you can employ bisquare weights instead of Huber weights if you specify psi accordingly. I would probably leave c at its default to begin with (I will modify my answer accordingly). $\endgroup$ – Jason Morgan Nov 1 '11 at 16:28
  • 1
    $\begingroup$ I also use rlm, and use the bisquare psi function because of its redescending property. Sometimes there are convergence issues with it, though, especially with smaller samples. $\endgroup$ – jbowman Dec 2 '11 at 15:03

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.