Pesos para residuos en regresión robusta

Pesos para residuos en regresión robusta


Uno de los primeros métodos para la regresión robusta fue la regresión iterativa de mínimos cuadrados reponderados (Huber, 1964). Este es un proceso iterativo que asigna un peso a cada observación. Inicialmente, todos los pesos son 1. El método ajusta un modelo de mínimos cuadrados a los datos ponderados y utiliza la magnitud de los residuos para determinar un nuevo conjunto de pesos para cada observación. En la mayoría de las implementaciones, las observaciones con residuos grandes se ponderan hacia abajo (se les asigna pesos inferiores a 1), mientras que las observaciones con residuos pequeños reciben pesos cercanos a 1. Los residuos son mucho más grandes que los demás.

Esta forma de regresión robusta también se conoce como estimación M. La «M» significa que el método es similar METROEstimación de máxima verosimilitud. El método intenta minimizar una función de los residuos ponderados. Este método no es robusto para puntos de alto apalancamiento, ya que tienen un impacto desproporcionado en el ajuste inicial, lo que a menudo significa que tienen pequeños residuos y nunca se ponderan.

Los resultados de la regresión dependen de la función que utilice para ponderar los residuos. El procedimiento ROBUSTREG en SAS admite 10 funciones de «peso» diferentes. Algunos de estos asignan un peso cero a las observaciones con residuos grandes, esencialmente excluyendo esa observación del ajuste. Otros usan una función geométrica o exponencialmente decreciente para asignar pesos pequeños (pero distintos de cero) a valores residuales grandes. Este artículo visualiza las 10 funciones de ponderación disponibles en el software SAS para realizar la estimación M en el método ROBUSTREG. Puede usar la subopción WEIGHTFUNCTION= de la opción METHOD=M en la declaración PROC ROBUSTREG para seleccionar una función de peso.

Las fórmulas y los valores de parámetros predeterminados para las 10 funciones de ponderación se pueden encontrar en la documentación del procedimiento ROBUSTREG. Los residuos se estandarizan utilizando una estimación de escala robusta. En el resto de este artículo, «el residuo» en realidad se refiere a un residuo estandarizado, no al residuo sin procesar. Las gráficas de las funciones de ponderación se muestran en el intervalo[-6, 6] y mostrar cómo las funciones asignan pesos en función del tamaño de los residuos estandarizados.

Funciones de ponderación diferenciables

Si usa mínimos cuadrados reponderados iterativamente para calcular las estimaciones, no importa si las funciones de ponderación son diferenciables. Sin embargo, si desea utilizar métodos de optimización para resolver las estimaciones de los parámetros, muchos algoritmos de optimización comunes funcionan mejor cuando la función objetivo es diferenciable. La función objetivo a menudo se relaciona con una suma que contiene los residuos ponderados, así que primero veamos las funciones de peso que son diferenciables. Las siguientes cuatro funciones son transformaciones diferenciables que asignan pesos a las observaciones en función del tamaño del residuo estandarizado.

Todas estas funciones tienen un aspecto similar. Todas son funciones simétricas y se parecen a la función cuadrática 1 – a Correcto2 cerca Correcto = 0. Sin embargo, las funciones difieren en la velocidad a la que disminuyen Correcto = 0:

  • La función biquad (también llamada función Tukey) es una función cuártica en el dominio [-c, c] y fuera de este intervalo es igual a 0. Esta es la función de ponderación predeterminada para la estimación M en PROC ROBUSTREG. La función se elige de tal manera que la derivada en ±C 0, por lo que la función es derivable. Una observación se excluye del cálculo si el tamaño del residual estandarizado excede c. Por defecto, C=4.685.
  • Las funciones de Cauchy y logística son funciones suaves que decrecen exponencialmente. Cuando se utilizan estas funciones, el peso de una observación nunca es 0, sino que el peso disminuye rápidamente con el tamaño de la observación.
    Correcto aumenta
  • La función de Welsch es una función gaussiana reescalada. El peso de una observación nunca es 0, pero el peso cae rápidamente (como exp(-son2)) de la Orden de
    Correcto aumenta

Funciones de peso no diferenciables

Las siguientes seis funciones no son diferenciables en uno o más lugares.

  • La función de Andrews es la función C*Pecado(control remoto)/Correcto dentro del intervalo [-c π, c π] y 0 fuera. una función cuártica en el rango [-c, c] y fuera de este intervalo es igual a 0. Es similar a la función biquad (ver sección anterior) pero no es diferenciable en ±C π.
  • Las funciones justa y mediana son similares. Asignas un peso que es asintóticamente proporcional a 1/| esCorrecto|. Utilice la función justa cuando desee reducir la ponderación de las observaciones según la magnitud inversa de su valor residual. No utilice la función mediana, que asigna grandes pesos a las observaciones con residuos pequeños.
  • Las funciones de Hampel, Huber y Talworth son similares en el sentido de que todas asignan pesos unitarios a residuos de pequeña magnitud. La función de Talworth es la función de peso más draconiana: asigna un peso de 1 o 0 a cada observación, dependiendo de si el tamaño del residuo supera un umbral. Las funciones de Hampel y Huber asignan residuos pequeños con un peso de 1 y luego asignan pesos decrecientes. Si usa los valores de parámetro predeterminados para la función de salto, los pesos fuera del intervalo son exactamente 0 [-8, 8].

Cuando uso M-estimate, normalmente uso una de las siguientes tres funciones:

  • La función bicuadrada que excluye las observaciones cuyos residuos tienen una magnitud superior a 4,685.
  • La función de Welsch que reduce suavemente los residuos grandes.
  • La función de Huber, que asigna 1 al peso de los residuos pequeños y asigna pesos a los residuos más grandes que son proporcionales a la inversa de la magnitud.

resumen

M-estimate es una técnica de regresión robusta que asigna pesos a cada observación. La idea es asignar pesos pequeños (o incluso cero) a las observaciones con residuos grandes, y asignar pesos cercanos a uno a los residuos pequeños. (En la práctica, el algoritmo utiliza residuos estandarizados para asignar ponderaciones). Hay muchas formas de asignar ponderaciones y SAS ofrece 10 funciones de ponderación diferentes. Este artículo visualiza las funciones de ponderación. Los resultados de la regresión dependen de la función que utilice para ponderar los residuos.

Apéndice: Programa SAS para visualizar las funciones de ponderación en PROC ROBUSTREG

Usé SAS/IML para definir las funciones de ponderación porque el lenguaje admite argumentos estándar para funciones mientras que PROC FCMP no lo hace. Cada función evalúa un vector de valores de entrada (x) y devuelve un vector de valores para la función de ponderación. Los valores de los parámetros por defecto se especifican con ARG=valor Notación en la firma de la función. Si omite el parámetro opcional, se utiliza el valor predeterminado.

Las siguientes declaraciones producen un panel de gráficos de 2 x 2 que muestra las formas de las funciones bicuadrática, de Cauchy, logística y de Welsch:

/***********************************/
/* Differentiable weight functions */
/* Assume x is a column vector     */
/***********************************/
 
proc iml;
/* The derivarive of Bisquare is f'(x)=4x(x^2-c^2)/c^4, so f'( +/-c )= 0 
   However, second derivative is f''(x) = 4(3x^2-c^2)/c^4, so f''(c) ^= 0 
*/
start Bisquare(x, c=4.685);   /* Tukey's function */
   w = j(nrow(x), 1, 0);   /* initialize output */
   idx = loc( abs(x) <= c );
   if ncol(idx)>0 then do;
      y = x[idx]/c;
      w[idx] = (1 - y##2)##2;
   end;
   return ( w );
finish;
 
start Cauchy(x, c=2.385);
   w = 1 / (1 + (abs(x)/c)##2);
   return ( w );
finish;
 
start Logistic(x, c=1.205);
   w = j(nrow(x), 1, 1);   /* initialize output */
   idx = loc( abs(x) > constant('maceps') ); /* lim x->0 f(x) = 1 */
   if ncol(idx)>0 then do;
      y = x[idx]/c;
      w[idx] = tanh(y) / y;
   end;
   return ( w );
finish;
 
start Welsch(x, c=2.985);
   w = exp( -(x/c)##2 / 2 );
   return ( w );
finish;
 
x = T(do(-6, 6, 0.005));
Function = j(nrow(x), 1, BlankStr(12));
 
create WtFuncs var {'Function' 'x' 'y'};
y = Bisquare(x);  Function[,] = "Bisquare"; append;
y = Cauchy(x);    Function[,] = "Cauchy";   append;
y = Logistic(x);  Function[,] = "Logistic"; append;
y = Welsch(x);    Function[,] = "Welsch";   append;
close;
 
quit;
 
ods graphics/ width=480px height=480px;
title "Differentiable Weight Functions for M Estimation";
title2 "The ROBUSTREG Procedure in SAS";
proc sgpanel data=WtFuncs;
   label x = "Scaled Residual" y="Weight";
   panelby Function / columns=2;
   series x=x y=y;
run;

Las siguientes declaraciones producen un panel de gráficos de 3 x 2 que muestra las formas de las funciones de Andrews, Fair, Hampel, Huber, Median y Talworth:

/*******************************/
/* Non-smooth weight functions */
/* Assume x is a column vector */
/*******************************/
proc iml;
 
start Andrews(x, c=1.339);
   w = j(nrow(x), 1, 0);   /* initialize output */
   idx = loc( abs(x) <= constant('pi')*c );
   if ncol(idx)>0 then do;
      y = x[idx]/c;
      w[idx] = sin(y) / y;
   end;
   return ( w );
finish;
 
start Fair(x, c=1.4);
   w = 1 / (1 + abs(x)/c);
   return ( w );
finish;
 
start Hampel(x, a=2, b=4, c=8);
   w = j(nrow(x), 1, 0);   /* initialize output */
   idx = loc( abs(x) <= a );
   if ncol(idx)>0 then
     w[idx] = 1;
   idx = loc( a < abs(x) & abs(x) <= b );
   if ncol(idx)>0 then
     w[idx] = a / abs(x[idx]);
   idx = loc( b < abs(x) & abs(x) <= c );
   if ncol(idx)>0 then
     w[idx] = (a / abs(x[idx])) # (c - abs(x[idx])) / (c-b);
   return ( w );
finish;
 
start Huber(x, c=1.345);
   w = j(nrow(x), 1, 1);   /* initialize output */
   idx = loc( abs(x) >= c );
   if ncol(idx)>0 then do;
      w[idx] = c / abs(x[idx]);
   end;
   return ( w );
finish;
 
start MedianWt(x, c=0.01);
   w = j(nrow(x), 1, 1/c);   /* initialize output */
   idx = loc( x ^= c );
   if ncol(idx)>0 then do;
      w[idx] = 1 / abs(x[idx]);
   end;
   return ( w );
finish;
 
start Talworth(x, c=2.795);
   w = j(nrow(x), 1, 0);   /* initialize output */
   idx = loc( abs(x) < c );
   if ncol(idx)>0 then do;
      w[idx] = 1;
   end;
   return ( w );
finish;
 
x = T(do(-6, 6, 0.01));
Function = j(nrow(x), 1, BlankStr(12));
create WtFuncs var {'Function' 'x' 'y'};
 
y = Andrews(x);  Function[,] = "Andrews"; append;
y = Fair(x);     Function[,] = "Fair";    append;
y = Hampel(x);   Function[,] = "Hampel";  append;
y = Huber(x);    Function[,] = "Huber";   append;
y = MedianWt(x); Function[,] = "Median";  append;
y = Talworth(x); Function[,] = "Talworth";append;
close;
quit;
 
ods graphics/ width=480px height=640px;
title "Non-Differentiable Weight Functions for M Estimation";
title2 "The ROBUSTREG Procedure in SAS";
proc sgpanel data=WtFuncs;
   label x = "Scaled Residual" y="Weight";
   panelby Function / columns=2;
   series x=x y=y;
   rowaxis max=1;
run;

Related post

Directores de omnicanalidad: cada vez más importantes para las empresas

Directores de omnicanalidad: cada vez más importantes para las…

Hace algunas semanas llamó la atención un post en Linkedin en que una de las empresas más reconocidas en el mundo…
Acelere el análisis para todos – Blog de Cloudera

Acelere el análisis para todos – Blog de Cloudera

 ¿Qué pasaría si pudiera acceder rápidamente a todos sus datos y ejecutar todos sus análisis en un flujo de trabajo…
Materiales 2D para la computación de última generación

Materiales 2D para la computación de última generación

En un comentario compacto publicado en Nature Communications, los investigadores de Graphene Flagship y 2D-EPL describen los campos de aplicación más…

Leave a Reply

Tu dirección de correo electrónico no será publicada.