Heteroscedasticidad simultánea y colas pesadas en un modelo de regresión.

Estoy tratando de crear un modelo de predicción usando la regresión. Este es el diagrama de diagnóstico para el modelo que obtengo al usar lm () en R:

Lo que leí en el gráfico QQ es que los residuos tienen una distribución de cola pesada, y el gráfico Residuals vs Fitted parece sugerir que la varianza de los residuos no es constante. Puedo domar las colas pesadas de los residuos utilizando un modelo robusto:

fitRobust = rlm(formula, method = "MM", data = myData)

Pero ahí es donde las cosas se detienen. El modelo robusto pesa varios puntos 0. Después de eliminar esos puntos, así es como se ven los residuos y los valores ajustados del modelo robusto:

La heterocedasticidad parece estar todavía allí. Utilizando

logtrans(model, alpha) 

del paquete MASS, intenté encontrar un $\alpha$ tal que

rlm(formula, method = "MM") 

con la fórmula siendo $\log(Y + \alpha) \sim X_1+\cdots+X_n$tiene residuos con varianza constante. Una vez que encuentre el$\alpha$, el modelo robusto resultante obtenido para la fórmula anterior tiene el siguiente gráfico Residuals vs Fitted:

Me parece que los residuos aún no tienen una variación constante. He intentado otras transformaciones de respuesta (incluida Box-Cox), pero tampoco parecen una mejora. Ni siquiera estoy seguro de que la segunda etapa de lo que estoy haciendo (es decir, encontrar una transformación de la respuesta en un modelo robusto) sea respaldada por alguna teoría. Agradecería mucho cualquier comentario, pensamiento o sugerencia.

La heterocedasticidad y la leptokurtosis se combinan fácilmente en el análisis de datos. Tome un modelo de datos que genere un término de error como Cauchy. Esto cumple con los criterios de homocedasticidad. La distribución de Cauchy tiene una varianza infinita. Un error de Cauchy es la forma en que un simulador incluye un proceso de muestreo atípico.

Con estos errores pesados, incluso cuando se ajusta al modelo medio correcto, el valor atípico conduce a un gran residuo. Una prueba de heteroscedasticidad ha inflado enormemente el error tipo I en este modelo. Una distribución de Cauchy también tiene un parámetro de escala. La generación de términos de error con un aumento lineal en la escala produce datos heteroscedasticos, pero el poder de detectar tales efectos es prácticamente nulo, por lo que el error tipo II también se infla.

Permítanme sugerir, entonces, que el enfoque analítico de datos adecuado no se enrede en las pruebas. Las pruebas estadísticas son principalmente engañosas. En ningún lugar es esto más obvio que las pruebas destinadas a verificar supuestos de modelado secundario. No sustituyen el sentido común. Para sus datos, puede ver claramente dos grandes residuos. Su efecto en la tendencia es mínimo, ya que pocos si los residuos se compensan en una desviación lineal de la línea 0 en la gráfica de residuos frente a ajustados. Eso es todo lo que necesitas saber.

Lo que se desea entonces es un medio para estimar un modelo de varianza flexible que le permitirá crear intervalos de predicción en un rango de respuestas ajustadas. Curiosamente, este enfoque es capaz de manejar la mayoría de las formas sensatas de heterocedasticidad y kurtotis. ¿Por qué no utilizar un enfoque de spline de suavizado para estimar el error cuadrático medio?

Tome el siguiente ejemplo:

set.seed(123)
x <- sort(rexp(100))
y <- rcauchy(100, 10*x)

f <- lm(y ~ x)
abline(f, col='red')
p <- predict(f)
r <- residuals(f)^2

s <- smooth.spline(x=p, y=r)

phi <- p + 1.96*sqrt(s$y)
plo <- p - 1.96*sqrt(s$y)

par(mfrow=c(2,1))
plot(p, r, xlab='Fitted', ylab='Squared-residuals')
lines(s, col='red')
legend('topleft', lty=1, col='red', "predicted variance")

plot(x,y, ylim=range(c(plo, phi), na.rm=T))
abline(f, col='red')
lines(x, plo, col='red', lty=2)
lines(x, phi, col='red', lty=2)

Da el siguiente intervalo de predicción que se "amplía" para acomodar los valores atípicos. Todavía es un estimador consistente de la varianza y útilmente le dice a la gente: "Oye, hay una observación grande y torpe alrededor de X = 4 y no podemos predecir valores muy útiles allí".