Ajuste una línea de regresión robusta usando un estimador MM en R

8

Contexto. Me gustaría ajustar una línea de regresión para estudiar la relación entre alguna variable de respuesta y alguna covariable continua . Debido a la presencia de puntos de apalancamiento incorrectos, he optado por un estimador MM en lugar del estimador LS habitual.yx

Metodología. Básicamente, la estimación MM es la estimación M inicializada por un estimador S. Por lo tanto, se deben elegir dos funciones de pérdida. Elegí la función de pérdida ampliamente utilizada de Tukey Biweight

ρ(tu)={1-[1-(tuk)2]3Si El |tuEl |k1Si El |tuEl |>k,

con en el estimador S preliminar (que da un punto de ruptura igual al ), y con en el paso de estimación M (para garantizar una eficiencia gaussiana del ).k=1.54850%k=2,69770%

Me gustaría usar R para ajustar mi robusta línea de regresión.

Pregunta.

library(MASS)
rlm(y~x, 
    method="MM",
    k0=1.548, c=2.697,
    maxit=50)
  • ¿Es mi código consistente con el párrafo anterior?
  • ¿Usarías otros argumentos opcionales?

EDITAR. Después de mi discusión con @Jason Morgan, me doy cuenta de que mi código anterior está equivocado. (@Jason Morgan: ¡Muchas gracias por esto!) Sin embargo, todavía no estoy convencido por su propuesta. En cambio, esto es lo que propongo ahora:

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

Creo que se apega a la metodología ahora. ¿Estás de acuerdo?

¡Gracias!

ocram
fuente

Respuestas:

5

Por defecto, la documentación indica que rlmusa psi=psi.huberpesos. Por lo tanto, si desea utilizar el bisquare de Tukey, debe especificarlo psi=psi.bisquare. La configuración predeterminada es psi.bisquare(u, c = 4.685, deriv = 0), que puede cambiar según lo desee. Por ejemplo, posiblemente algo como

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

También es posible que desee investigar si debe usar cuadrados menos recortados ( init="lts") para inicializar sus valores iniciales. El valor predeterminado es usar mínimos cuadrados.

Jason Morgan
fuente
@Janson Morgan: ¿estás seguro de lo que propones? ¿Tienes alguna experiencia con esa función? Mi documentación (R 2.13.1) en realidad indica "El conjunto inicial de coeficientes y la escala final son seleccionados por un estimador S con k0 = 1.548; esto da (para n >> p) punto de ruptura 0.5. El estimador final es un El estimador M con el biweight de Tukey y la escala fija que heredarán este punto de ruptura proporcionaron c> k0; esto es cierto para el valor predeterminado de c que corresponde al 95% de eficiencia relativa en la normalidad ".
ocram
1
He estimado estos modelos en el pasado. Como dice la documentación, el primer paso en la estimación de MM se realiza con los pesos de Huber, el segundo con los pesos bicares. Mis notas (de hace un par de años) afirman que en el primer paso S, puede emplear pesos bisquare en lugar de pesos Huber si así lo especifica psi. cPara empezar , probablemente dejaría su valor predeterminado (modificaré mi respuesta en consecuencia).
Jason Morgan el
1
También uso rlm, y uso la función psi bisquare debido a su propiedad de redescending. Sin embargo, a veces hay problemas de convergencia, especialmente con muestras más pequeñas.
jbowman