15
$\begingroup$

I am using rlm in the R MASS package to regress a multivariate linear model. It works well for a number of samples but I am getting quasi-null coefficients for a particular model:

Call: rlm(formula = Y ~ X1 + X2 + X3 + X4, data = mymodel, maxit = 50, na.action = na.omit)
Residuals:
       Min         1Q     Median         3Q        Max 
-7.981e+01 -6.022e-03 -1.696e-04  8.458e-03  7.706e+01 

Coefficients:
             Value    Std. Error t value 
(Intercept)    0.0002   0.0001     1.8418
X1             0.0004   0.0000    13.4478
X2            -0.0004   0.0000   -23.1100
X3            -0.0001   0.0002    -0.5511
X4             0.0006   0.0001     8.1489

Residual standard error: 0.01086 on 49052 degrees of freedom
  (83 observations deleted due to missingness)

For comparison, these are the coefficients calculated by lm():

Call:
lm(formula = Y ~ X1 + X2 + X3 + X4, data = mymodel, na.action = na.omit)

Residuals:
    Min      1Q  Median      3Q     Max 
-76.784  -0.459   0.017   0.538  78.665 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.016633   0.011622  -1.431    0.152    
X1            0.046897   0.004172  11.240  < 2e-16 ***
X2           -0.054944   0.002184 -25.155  < 2e-16 ***
X3            0.022627   0.019496   1.161    0.246    
X4            0.051336   0.009952   5.159  2.5e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 2.574 on 49052 degrees of freedom
  (83 observations deleted due to missingness)
Multiple R-squared: 0.0182, Adjusted R-squared: 0.01812 
F-statistic: 227.3 on 4 and 49052 DF,  p-value: < 2.2e-16 

The lm plot doesn't show any particularly high outlier, as measured by Cook's distance:

lm diagnostic

EDIT

For reference and after confirming results based on the answer provided by Macro, the R command to set the tuning parameter, k, in the Huber estimator is (k=100 in this case):

rlm(y ~ x, psi = psi.huber, k = 100)
$\endgroup$
  • $\begingroup$ The residual standard errors, in combination with the other information, make it look as though the rlm weight function is throwing out almost all the observations. Are you sure it's the same Y in the two regressions? (Just checking...) Try method="MM" in your rlm call, then try (if that fails) psi=psi.huber(k=2.5) (2.5 is arbitrary, just bigger than the default 1.345) which spreads out the lm-like region of the weight function. $\endgroup$ – jbowman May 31 '12 at 13:44
  • $\begingroup$ @jbowman Y is correct. Added the MM method. My intuition is the same you mentioned. This model residuals are relatively compact compared to the others I have tried. It looks like the methodology is discarding most observations. $\endgroup$ – Robert Kubrick May 31 '12 at 13:56
  • 1
    $\begingroup$ @RobertKubrick you understand what setting k to 100 means, right? $\endgroup$ – user603 May 14 '13 at 13:41
  • $\begingroup$ Based on this: Multiple R-squared: 0.0182, Adjusted R-squared: 0.01812 you should examine your model one more time. Outliers, transformation of the response or predictors. Or you should consider nonlinear model. Predictor X3 is not significant. What you made is not good linear model. $\endgroup$ – Marija Milojevic Sep 20 '15 at 11:37
15
$\begingroup$

The difference is that rlm() fits models using your choice of a number of different $M$-estimators, while lm() uses ordinary least squares.

In general the $M$-estimator for a regression coefficient minimizes

$$ \sum_{i=1}^{n} \rho \left( \frac{Y_i - {\bf X}_{i} {\boldsymbol \beta}}{\sigma} \right) $$

as a function of ${\boldsymbol \beta}$, where $Y_i$ is the $i$'th response, and ${\bf X}_{i}$ is the predictors for individual $i$. Least squares is a special case of this where $$ \rho(x) = x^2 $$ However, the default setting for rlm(), which you appear to be using, is the Huber $M$-estimator, which uses

$$ \rho(x) = \begin{cases} \frac{1}{2} x^2 &\mbox{if } |x| \leq k\\ k |x| - \frac{1}{2} k^2 & \mbox{if } |x| > k. \end{cases} $$

where $k$ is a constant. The default in rlm() is $k = 1.345$. These two estimators are minimizing different criteria, so it is no surprise that the estimates are different.

Edit: From the QQ plot shown above, it looks like you have a very long tailed error distribution. This is the kind of situation the Huber M-estimator is designed for and, in that situation, can give quite different estimates:

When the errors are normally distributed, the estimates will be pretty similar since, under the normal distribution, most of the Huber $ρ$ function will fall under the $|x|<k$ situation, which is equivalent to least squares. In the long tailed situation you have, many fall into the $|x|>k$ situation, which is a departure from OLS, which would explain the discrepancy.

$\endgroup$
  • $\begingroup$ I have tried several other models (same number of observations, same IVs) and the coefficients are fairly similar between rlm and lm. There must be something in this particular data set that is producing the large difference in the coefficients. $\endgroup$ – Robert Kubrick May 31 '12 at 14:23
  • 1
    $\begingroup$ No, there are not standardized methods for choosing $k$ - they are tuning parameters and are typically chosen in an ad hoc way. In the seminal paper (Huber, 1964), he notes that anywhere between 1.0 and 2.0 gives acceptable results and that the choice doesn't matter a whole lot. In this paper (education.wayne.edu/jmasm/sawilowsky_lre.pdf) the authors using a concept called 'Location Relative Efficiency' to choose to index. In any case, I don't recommend treating the least squares estimates as maximum likelihood estimates in your data - the errors are very long tailed. $\endgroup$ – Macro May 31 '12 at 15:41
  • 1
    $\begingroup$ One thing you could do to validate (to some extent) this is to try $k=1.5, 2, 2.5, 3, 3.5, 4$ in the psi.huber function and see how the residual standard error and parameter estimates change. As $k$ gets larger, there should be some approach to the lm estimates. Also, it's possible that the starting estimate of spread (MAD) with this dataset is very, very small, which you can check by calculating MAD on the residuals from the rlm; in this case, everything of any magnitude is being thrown out because the estimate of spread is too small, and varying k some won't make a difference. $\endgroup$ – jbowman May 31 '12 at 16:27
  • 1
    $\begingroup$ That's for the added info, @jbowman - these are useful comments. Regarding in your last comment, those large observations aren't exactly being thrown out - their influence is just being dialed down (as it seems they should be), right? $\endgroup$ – Macro May 31 '12 at 16:43
  • 1
    $\begingroup$ @RobertKubrick, Huber (1964) showed that this estimating equation gives statistical inference that is correct in the face of errors that are a mixture between normal and long-tailed errors so it's robust in the sense that it can handle this type of non-normality. Re: your last comment - that's not true. Note that we scale by $\sigma$ - a poorly fitting model can have normal errors. Once we scale by $\sigma$ those errors will no longer be "large". This, in some sense, downweights observations with residuals inconsistent with normality although, as I said, this is not how the method was derived. $\endgroup$ – Macro May 31 '12 at 17:32

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.