Hace poco estuve comparando los resultados de regresiones LOWESS realizadas en R (y utilizando el módulo statsmodels de Python) y Stata. Me di cuenta de que algunos de los valores obtenidos por Stata parecen estar fuera de lugar; específicamente, son las colas las que parecen estar estimadas incorrectamente.
Me zambullí en el código fuente de la R lowess() (que parece estar basada en el código Fortran original de Cleveland que se encuentra en aquí ) y el código heredado de Stata para el ksm que se encuentra en el archivo de actualización de ado de Stata 7 (que se encuentra aquí ). Tenga en cuenta que los Statas posteriores a la v7 implementaron el lowess mediante el comando _LOWESS C que no está expuesta al usuario. Verifiqué que ejecutando la rutina ksm utilizando el comando opcional lowess en Stata 14 genera el mismo resultado ("incorrecto") que ejecutar lowess directamente. Eso es, ksm Y X, lowess es lo mismo que lowess Y X .
Después de inspeccionar los respectivos códigos fuente, creo que el problema con la implementación de Stata radica en el hecho de que el tamaño del subconjunto utilizado en cada regresión local depende de la posición ordinal de los valores X. Es decir, para valores extremos y casi extremos de X (en otras palabras, para valores cercanos a las colas), Stata utiliza un subconjunto más pequeño que para X más centrales. Intuitivamente, el problema se puede ilustrar utilizando un ejemplo simple con 100 puntos de datos donde el parámetro de ancho de banda se elige para ser 0.4 de modo que cada subconjunto es del tamaño 0.4*100=40. En R, el tamaño de los subconjuntos utilizados para estimar $Y_1$ , $Y_{10}$ , $Y_{20}$ , $Y_{30}$ y $Y_{40}$ quedaría así:
Por otra parte, Stata parece reducir el tamaño de los subconjuntos para todas las X de (aproximadamente) 1 a 20 y de 80 a 100:
Cuando reescribí el ksm de forma que el tamaño de los subconjuntos se mantuviera fijo para todas las X, obtuve los mismos resultados que en R (o Python). A continuación se muestran los pseudocódigos de ambas implementaciones:
# bwidth = bandwidth
# count = N
# i = index of local regression considered
# Cleveland's Fortran implementation
size = max(min(floor(bwidth*count), count), 2) #size of neighborhood
lower_bound = 1
upper_bound = size
for each i
x = X[i]
diff_left = x - X[lower_bound]
diff_right = X[upper_bound + 1] - x
if diff_left > diff_right and diff_right < count
increment lower_bound and upper_bound by 1
h = max(x - X[lower_bound], X[upper_bound] - x)
h9 = 0.999*h and h1 = 0.001*h
r = abs(X - x)
W = 1 - (abs(X - x)/h)^3)^3 if r > h1 and r <= h9, W = 1 if r <= h1, W = 0 otherwise
run WLS of Y on X using weights W if within lower_bound/upper_bound, get Y_hat[i]
# Stata implementation
k = int((bwidth*count - 0.5)/2) #half-bandwidth
for each i
x = X[i]
lower_bound = max(1, i - k)
upper_bound = min(i + k, count)
h = 1.0001*max(X[upper_bound] - x, x - X[lower_bound])
W = 1 - (abs(X - x)/h)^3)^3 if h != 0, W = 1 otherwise
run WLS of Y on X using weights W if within lower_bound/upper_bound, get Y_hat[i]Entiendo que a veces las implementaciones difieren entre herramientas de software (o incluso paquetes hechos para el mismo software). Dicho esto, yo esperaría encontrar al menos algunas referencias en línea que discutan tales enfoques no estándar. He intentado encontrar tantas otras implementaciones de la regresión LOWESS como he podido y todas parecen estar basadas en el código Fortran mencionado anteriormente y/o al menos siguen el mismo enfoque. No pude encontrar ninguna referencia que discutiera el enfoque específico utilizado por Stata.
Me gustaría saber si tal implementación de LOWESS es correcta, si se trata de algún método "propietario", cuasi-correcto o simplemente de un bug. En el primer caso, también agradecería que alguien me indicara una referencia (preferiblemente académica). Gracias.
EDITAR I : Me pregunto si esto se debe a que históricamente Stata lowess se basaba en el comando ksm que (cuando se utiliza con su configuración por defecto) no estaba destinado a estimar la regresión LOWESS. Parece que desde Stata 8, se tomó la decisión de abandonar ksm y en su lugar aplicar lowess donde se ha elegido por defecto la regresión LOWESS (mientras que ksm se implementaron como lowess 's no predeterminados). El tratamiento de las colas, sin embargo, sigue siendo el mismo que en ksm .
EDITAR II : Nótese que todas estas comparaciones se basaron en la coincidencia de todos los demás parámetros de las regresiones LOWESS. Como Stata no permite iteraciones múltiples y tampoco implementa la interpolación para reducir el número de regresiones locales requeridas (como lo hacen R y el código Fortran original), me aseguré de que el comando de R usara tanto iter=0 y delta=0 parámetros. En concreto, estaba comparando lowess(X, Y, f=0.4, delta=0, iter=0) (R) con lowess Y X, bwidth(0.4) (Stata) y statsmodels.api.nonparametric.lowess(Y, X, frac=0.4, it=0, delta=0) (Python). Tanto X como Y se portaron bien.
UPDATE : OK, parece que no fui el primero en darme cuenta de esto. Acabo de encontrar un viejo usuario escrito adjksm que "es idéntico a ksm excepto que el ancho de banda del suavizador es constante a lo largo del eje x". Véase aquí . También confirmé que adjksm Y X, bwidth(0.4) lowess genera el mismo resultado que R lowess() utilizando iter=0 y delta=0 hasta el cuarto decimal.
