Heterocedasticidade simultânea e caudas pesadas em um modelo de regressão

8

Estou tentando criar um modelo de previsão usando regressão. Este é o gráfico de diagnóstico para o modelo que recebo usando lm () no R: gráficos de diagnóstico de R

O que li do gráfico QQ é que os resíduos têm uma distribuição de cauda pesada, e o gráfico Residuais vs Ajustados parece sugerir que a variação dos resíduos não é constante. Eu posso domesticar as caudas pesadas dos resíduos usando um modelo robusto:

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

Mas é aí que as coisas param. O modelo robusto pesa vários pontos 0. Depois de remover esses pontos, é assim que os resíduos e os valores ajustados do modelo robusto se parecem:Resíduos vs Ajustados para o modelo robusto

A heterocedasticidade parece ainda estar lá. Usando

logtrans(model, alpha) 

α

rlm(formula, method = "MM") 

registro(Y+α)X1++Xnα

Residuais vs Ajustados para resposta transformada por log

Parece-me que os resíduos ainda não têm variação constante. Eu tentei outras transformações de resposta (incluindo Box-Cox), mas elas também não parecem melhorar. Não tenho certeza de que o segundo estágio do que estou fazendo (ou seja, encontrar uma transformação da resposta em um modelo robusto) seja suportado por qualquer teoria. Eu aprecio muito quaisquer comentários, pensamentos ou sugestões.

user765195
fonte
2
Eu acho que você está sendo um pouco exigente quanto à variação não constante. Parece-me bem. Qual é o objetivo da regressão? Explicação / teste de hipóteses ou previsão?
probabilityislogic
@probabilityislogic, obrigado pelo seu comentário. Eu aprecio muito isso. Meu objetivo é previsão. Você está certo. Eu provavelmente estou sendo muito exigente. Existe uma medida de heterocedasticidade que eu possa olhar? Pensei em plotar variação versus valores ajustados, mas não há muitos pontos para cada valor previsto para calcular a variação. Também estou curioso para entender qual é a solução para esse problema em geral. As transformações Box-Cox e log também são aplicáveis ​​a modelos robustos?
user765195
Você pode fazer um teste pareado para igualdade de variâncias usando o teste F para um modelo com termos de erro gaussiano ou se eles tiverem uma distribuição não gaussiana, existem testes robustos de dispersão, como o teste de Levene.
Michael R. Chernick 25/09/12
Obrigado @MichaelChernick. Eu aprecio muito o seu comentário. Finalmente, usei a generalização de Koenker do teste de Breusch-Pagan para heterocedasticidade, conforme implementado no pacote lmtest em R ( hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/lmtest/html/… ).
user765195

Respostas:

3

A heterocedasticidade e a leptokurtosis são facilmente confundidas na análise dos dados. Pegue um modelo de dados que gere um termo de erro como Cauchy. Isso atende aos critérios de homocedasticidade. A distribuição de Cauchy tem variação infinita. Um erro de Cauchy é a maneira de um simulador incluir um processo de amostragem externa.

Com esses erros de cauda pesados, mesmo quando você se encaixa no modelo médio correto, o outlier leva a um grande resíduo. Um teste de heterocedasticidade inflou bastante o erro do tipo I neste modelo. Uma distribuição Cauchy também possui um parâmetro de escala. A geração de termos de erro com um aumento linear na escala produz dados heterocedásticos, mas o poder de detectar esses efeitos é praticamente nulo, portanto o erro do tipo II também é inflado.

Deixe-me sugerir, então, que a abordagem analítica de dados adequada não seja atolada nos testes. Os testes estatísticos são principalmente enganosos. Onde isso é mais óbvio do que os testes destinados a verificar as suposições de modelagem secundárias. Eles não substituem o bom senso. Para seus dados, você pode ver claramente dois grandes resíduos. Seu efeito sobre a tendência é mínimo, se houver algum resíduo compensado em uma partida linear da linha 0 no gráfico de resíduos versus ajustado. Isso é tudo que você precisa saber.

O que se deseja, então, é um meio de estimar um modelo de variação flexível que permita criar intervalos de previsão em um intervalo de respostas ajustadas. Curiosamente, essa abordagem é capaz de lidar com a maioria das formas sãs de heterocedasticidade e kurtotis. Por que não usar uma abordagem de spline de suavização para estimar o erro quadrático médio.

Veja o seguinte exemplo:

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)

Fornece o seguinte intervalo de previsão que "aumenta" para acomodar os valores extremos. Ainda é um estimador consistente da variação e diz às pessoas: "Ei, essa observação grande e instável em torno de X = 4 e não podemos prever valores muito úteis lá".

insira a descrição da imagem aqui

AdamO
fonte
Isso funcionaria para outros tipos de lms, como gls?
user2974951 27/02