¿Por qué las estimaciones de los coeficientes de regresión de rlm() son diferentes a las de lm() en R?

Question

Estoy utilizando rlm en el paquete MASS de R para hacer una regresión de un modelo lineal multivariante. Funciona bien para una serie de muestras, pero estoy obteniendo coeficientes casi nulos para un modelo concreto:

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)

Para comparar, estos son los coeficientes calculados por 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 

El gráfico lm no muestra ningún valor atípico particularmente alto, medido por la distancia de Cook:

EDITAR

Como referencia y después de confirmar los resultados basados en la respuesta proporcionada por Macro, el comando R para establecer el parámetro de sintonía, k en el estimador de Huber es ( k=100 en este caso):

rlm(y ~ x, psi = psi.huber, k = 100)

Answer 1

La diferencia es que rlm() se ajusta a los modelos usando su elección de un número de diferentes $M$ -estimadores, mientras que lm() utiliza los mínimos cuadrados ordinarios.

En general, el $M$ -El estimador de un coeficiente de regresión minimiza

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

en función de ${\boldsymbol \beta}$ , donde $Y_i$ es el $i$ respuesta, y ${\bf X}_{i}$ son los predictores del individuo $i$ . Los mínimos cuadrados son un caso especial en el que $$ \rho(x) = x^2 $$ Sin embargo, la configuración por defecto para rlm() que parece estar utilizando, es el Huber $M$ -estimador, que utiliza

$$ \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} $$

donde $k$ es una constante. El valor por defecto en rlm() es $k = 1.345$ . Estos dos estimadores minimizan criterios diferentes, por lo que no es de extrañar que las estimaciones sean diferentes.

Editar: A partir del gráfico QQ mostrado arriba, parece que tienes una distribución de errores de cola muy larga. Este es el tipo de situación para la que está diseñado el estimador M de Huber y, en esa situación, puede dar estimaciones bastante diferentes:

Cuando los errores se distribuyen normalmente, las estimaciones serán bastante similares ya que, bajo la distribución normal, la mayoría de los 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$, que se aleja de OLS, lo que explicaría la discrepancia.