Почему оценки коэффициента регрессии rlm () отличаются от lm () в R?

Я использую rlm в пакете R MASS для регрессии многомерной линейной модели. Это хорошо работает для ряда образцов, но я получаю квазинулевые коэффициенты для конкретной модели:

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)

Для сравнения это коэффициенты, рассчитанные с помощью 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

График lm не показывает какой-либо особенно высокий выброс, измеренный расстоянием Кука:

лм диагностика

РЕДАКТИРОВАТЬ

Для справки и после подтверждения результатов, основанных на ответе, предоставленном макросом, команда R для установки параметра настройки kв оценке Хьюбера имеет значение ( k=100в данном случае):

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

r multiple-regression robust Роберт Кубрик
источник

Остаточные стандартные ошибки в сочетании с другой информацией создают впечатление, что rlmвесовая функция отбрасывает почти все наблюдения. Вы уверены, что это один и тот же Y в двух регрессиях? (Просто проверяю ...) Попробуйте method="MM"в вашем rlmвызове, затем попробуйте (если это не удастся) psi=psi.huber(k=2.5)(2.5 - произвольно, только больше, чем по умолчанию 1.345), что расширяет область, lmподобную весовой функции.

jbowman

@jbowman Y правильно. Добавлен метод ММ. Моя интуиция такая же, как вы упомянули. Остатки этой модели относительно компактны по сравнению с другими, которые я пробовал. Похоже, что методология отбрасывает большинство наблюдений.

Роберт Кубрик

@ РобертКубрик, ты понимаешь, что означает установка k на 100 , верно?

user603

Исходя из этого: кратный R-квадрат: 0,0182, скорректированный R-квадрат: 0,01812, вам следует еще раз изучить вашу модель. Выбросы, трансформация ответа или предикторы. Или вы должны рассмотреть нелинейную модель. Предиктор Х3 не имеет значения. То, что вы сделали, не является хорошей линейной моделью.

Мария Милоевич

Разница в том, что rlm()подходит для моделей, использующих выбранный вами ряд различных оценщиков, в то время как используются обычные наименьшие квадраты. $M$ lm()

В целом оценщик для коэффициента регрессии минимизирует $M$

Σ_{я знак равно 1}^{N} ρ (\frac{Y_{я} - {Икс}_{я} β}{σ})

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

как функция , где - это -й ответ, а - предикторы для индивидуума . Наименьшие квадраты - это частный случай этого случая, где Тем не менее, по умолчанию для параметра , который вы, похоже, используете, используется метод оценки Huber , который использует ${\boldsymbol \beta}$ $Y_i$ $i$ ${\bf X}_{i}$ $i$

ρ (Икс) знак равно {Икс}^{2}

$\rho(x) = x^2$ rlm()

M

$M$

ρ (Икс) знак равно {\begin{cases} \frac{1}{2} {Икс}^{2} & если | Икс | \leq К \\ К | Икс | - \frac{1}{2} К^{2} & если | Икс | > К, \end{cases}

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

где является константой. По умолчанию в это . Эти две оценки минимизируют разные критерии, поэтому неудивительно, что оценки разные. $k$ rlm() $k = 1.345$

Изменить: Из графика QQ, показанного выше, похоже, что у вас очень длинное хвостовое распределение ошибок. Это та ситуация, для которой предназначен М-оценщик Хубера, и в этой ситуации она может давать совершенно разные оценки:

Когда ошибки распределены нормально, то оценки будут очень похожи , так как, при нормальном распределении, большинство Huber функции будут подпадать под ситуации, которая эквивалентна наименьших квадратов. В ситуации длинного хвоста многие попадают в ситуацию , которая является отклонением от OLS, что объясняет расхождение. $ρ$ $|x|<k$ $|x|>k$

макрос
источник

Я пробовал несколько других моделей (такое же количество наблюдений, те же IV), и коэффициенты довольно похожи между RLM и LM. В этом конкретном наборе данных должно быть что-то, что создает большую разницу в коэффициентах.

Роберт Кубрик

Нет, не существует стандартизированных методов выбора - они являются параметрами настройки и обычно выбираются специальным образом. В оригинальной статье (Huber, 1964) он отмечает, что где-то между 1,0 и 2,0 дает приемлемые результаты и что выбор не имеет большого значения. В этой статье ( education.wayne.edu/jmasm/sawilowsky_lre.pdf ) авторы используют концепцию «Относительная эффективность местоположения» для выбора индексации. В любом случае, я не рекомендую рассматривать оценки наименьших квадратов как оценки максимального правдоподобия в ваших данных - ошибки очень длинные.

k

$k$

Макрос

Одна вещь, которую вы могли бы сделать, чтобы проверить (до некоторой степени) это, попробовать в функции и посмотреть, как изменяются остаточная стандартная ошибка и оценки параметров. Поскольку становится больше, должен быть некоторый подход к оценкам. Кроме того, возможно, что начальная оценка спреда (MAD) с этим набором данных очень и очень мала, что вы можете проверить, рассчитав MAD по остаткам от ; в этом случае выбрасывается все любой величины, потому что оценка разброса слишком мала, и варьирование k в некоторых не будет иметь значения.

k = 1.5, 2, 2.5, 3, 3.5, 4

$k=1.5, 2, 2.5, 3, 3.5, 4$ psi.huber

k

$k$ lmrlm

jbowman

Это для дополнительной информации, @jbowman - это полезные комментарии. Что касается вашего последнего комментария, то эти крупные наблюдения точно не отбрасываются - их влияние просто ослабляется (как, кажется, и должно быть), верно?

Макрос

@RobertKubrick, Huber (1964) показал, что это оценочное уравнение дает статистический вывод, который является правильным, несмотря на ошибки, которые представляют собой смесь между нормальными и длиннохвостыми ошибками, поэтому он устойчив в том смысле, что он может обрабатывать этот тип ненормальности , Re: ваш последний комментарий - это не так. Обратите внимание, что мы масштабируем по - плохо подходящая модель может иметь нормальные ошибки. Как только мы масштабируем по эти ошибки больше не будут "большими". Это, в некотором смысле, наблюдения с понижением веса с остатками, несовместимыми с нормальностью, хотя, как я уже сказал, метод не был получен.

σ

$\sigma$

σ

$\sigma$

Макрос

Почему оценки коэффициента регрессии rlm () отличаются от lm () в R?

Ответы: