Hétéroscédasticité simultanée et queues lourdes dans un modèle de régression

8

J'essaie de créer un modèle de prédiction en utilisant la régression. Voici le tracé de diagnostic pour le modèle que j'obtiens en utilisant lm () dans R: parcelles de diagnostic de R

Ce que j'ai lu dans le graphique QQ, c'est que les résidus ont une distribution à queue lourde, et le graphique Residuals vs Fitted semble suggérer que la variance des résidus n'est pas constante. Je peux apprivoiser les queues lourdes des résidus en utilisant un modèle robuste:

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

Mais c'est là que les choses s'arrêtent. Le modèle robuste pèse plusieurs points 0. Après avoir supprimé ces points, voici à quoi ressemblent les résidus et les valeurs ajustées du modèle robuste:Résidus vs ajusté pour le modèle robuste

L'hétéroscédasticité semble être toujours là. En utilisant

logtrans(model, alpha) 

du paquet MASS, j'ai essayé de trouver un tel queα

rlm(formula, method = "MM") 

avec la formule étant a des résidus avec une variance constante. Une fois que j'ai trouvé le , le modèle robuste résultant obtenu pour la formule ci-dessus a le tracé résiduel vs ajusté suivant:Journal(Oui+α)X1++Xnα

Résidus vs ajustés pour une réponse transformée en journal

Il me semble que les résidus n'ont toujours pas de variance constante. J'ai essayé d'autres transformations de réponse (dont Box-Cox), mais elles ne semblent pas non plus être une amélioration. Je ne suis même pas sûr que la deuxième étape de ce que je fais (c'est-à-dire trouver une transformation de la réponse dans un modèle robuste) ne soit étayée par aucune théorie. J'apprécierais beaucoup tout commentaire, réflexion ou suggestion.

user765195
la source
2
Je pense que vous êtes un peu pointilleux sur la variance non constante. Ça me semble bien. Quel est le but de la régression? Test d'explication / hypothèse ou prédiction?
Probabilislogic
@probabilityislogic, merci pour votre commentaire. J'apprécie énormément. Mon objectif est la prédiction. Tu as raison. Je suis probablement trop pointilleux. Existe-t-il une mesure de l'hétéroscédasticité que je peux examiner? J'ai pensé à représenter la variance par rapport aux valeurs ajustées, mais il n'y a pas beaucoup de points pour chaque valeur prédite pour calculer la variance. Je suis également curieux de comprendre quelle est la solution à ce problème en général. Les transformations Box-Cox et log sont-elles également applicables aux modèles robustes?
user765195
Vous pouvez effectuer un test par paire pour l'égalité des variances en utilisant le test F pour un modèle avec des termes d'erreur gaussiens ou s'ils ont une distribution non gaussienne, il existe des tests robustes de dispersion tels que le test de Levene.
Michael R. Chernick,
Merci @MichaelChernick. J'apprécie beaucoup votre commentaire. J'ai finalement utilisé la généralisation par Koenker du test d'hétéroscédasticité de Breusch-Pagan tel qu'il est implémenté dans le package lmtest dans R ( hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/lmtest/html/… ).
user765195

Réponses:

3

L'hétéroscédasticité et la leptokurtose sont facilement confondues dans l'analyse des données. Prenez un modèle de données qui génère un terme d'erreur comme Cauchy. Cela répond aux critères d'homoscédastictie. La distribution de Cauchy a une variance infinie. Une erreur de Cauchy est la manière d'un simulateur d'inclure un processus d'échantillonnage de valeurs aberrantes.

Avec ces lourdes erreurs à queue, même lorsque vous ajustez le modèle moyen correct, la valeur aberrante conduit à un grand résidu. Un test d'hétéroscédasticité a considérablement gonflé l'erreur de type I dans ce modèle. Une distribution de Cauchy a également un paramètre d'échelle. La génération de termes d'erreur avec une augmentation d'échelle linéaire produit des données hétéroscédastiques, mais le pouvoir de détecter de tels effets est pratiquement nul, de sorte que l'erreur de type II est également gonflée.

Permettez-moi de suggérer alors que l'approche analytique des données appropriée n'est pas de s'embourber dans les tests. Les tests statistiques sont principalement trompeurs. Nulle part cela n'est plus évident que les tests destinés à vérifier les hypothèses de modélisation secondaires. Ils ne remplacent pas le bon sens. Pour vos données, vous pouvez clairement voir deux gros résidus. Leur effet sur la tendance est minime, car peu ou pas de résidus sont compensés par un écart linéaire par rapport à la ligne 0 dans le graphique des résidus par rapport aux ajustés. C'est tout ce que vous devez savoir.

Ce que l'on souhaite alors, c'est un moyen d'estimer un modèle de variance flexible qui vous permettra de créer des intervalles de prédiction sur une gamme de réponses ajustées. Fait intéressant, cette approche est capable de gérer la plupart des formes saines d'hétéroscédasticité et de kurtotis. Pourquoi ne pas alors utiliser une approche spline de lissage pour estimer l'erreur quadratique moyenne.

Prenons l'exemple suivant:

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)

Donne l'intervalle de prédiction suivant qui "s'élargit" pour s'adapter à la valeur aberrante. C'est toujours un estimateur cohérent de la variance et dit utilement aux gens: "Hé, il y a cette grande observation bancale autour de X = 4 et nous ne pouvons pas prédire les valeurs très utilement là-bas."

entrez la description de l'image ici

AdamO
la source
Est-ce que cela fonctionnerait pour d'autres types de lms, tels que les gls?
user2974951