Contexto
Esta pergunta usa R, mas trata de questões estatísticas gerais.
Estou analisando os efeitos dos fatores de mortalidade (% de mortalidade por doenças e parasitismo) na taxa de crescimento populacional da mariposa ao longo do tempo, onde populações de larvas foram amostradas de 12 locais uma vez por ano durante 8 anos. Os dados da taxa de crescimento populacional exibem uma tendência cíclica clara, porém irregular, ao longo do tempo.
Os resíduos de um modelo linear generalizado simples (taxa de crescimento ~% de doença +% parasitismo + ano) apresentaram uma tendência cíclica igualmente clara, porém irregular, ao longo do tempo. Portanto, modelos de mínimos quadrados generalizados da mesma forma também foram ajustados aos dados com estruturas de correlação apropriadas para lidar com a autocorrelação temporal, por exemplo, simetria composta, ordem de processo autoregressiva 1 e estruturas de correlação de média móvel autoregressiva.
Todos os modelos continham os mesmos efeitos fixos, foram comparados usando o AIC e foram ajustados pelo REML (para permitir a comparação de diferentes estruturas de correlação pelo AIC). Estou usando o pacote R nlme e a função gls.
Questão 1
Os resíduos dos modelos GLS ainda exibem padrões cíclicos quase idênticos quando plotados contra o tempo. Esses padrões sempre permanecerão, mesmo em modelos que explicam com precisão a estrutura de autocorrelação?
Simulei alguns dados simplificados, mas similares, em R abaixo da minha segunda pergunta, que mostra o problema com base no meu entendimento atual dos métodos necessários para avaliar padrões autocorrelacionados temporalmente em resíduos de modelo , que agora sei que estão errados (ver resposta).
Questão 2
Eu ajustei os modelos GLS com todas as estruturas de correlação plausíveis possíveis para os meus dados, mas nenhum é realmente muito melhor do que o GLM sem nenhuma estrutura de correlação: apenas um modelo GLS é marginalmente melhor (pontuação AIC = 1,8 mais baixa), enquanto todo o resto tem valores mais altos de AIC. No entanto, esse é apenas o caso quando todos os modelos são ajustados pelo REML, não pelo ML, onde os modelos GLS são claramente muito melhores, mas eu entendo nos livros de estatísticas que você só deve usar o REML para comparar modelos com diferentes estruturas de correlação e os mesmos efeitos fixos por razões Não vou detalhar aqui.
Dada a natureza claramente correlacionada temporalmente dos dados, se nenhum modelo é moderadamente melhor que o GLM simples, qual é a maneira mais apropriada de decidir qual modelo usar para inferência, supondo que eu esteja usando um método apropriado (eu finalmente quero usar AIC para comparar diferentes combinações de variáveis)?
Q1 'simulação' explorando padrões residuais em modelos com e sem estruturas de correlação apropriadas
Gere variável de resposta simulada com um efeito cíclico de 'tempo' e um efeito linear positivo de 'x':
time <- 1:50
x <- sample(rep(1:25,each=2),50)
y <- rnorm(50,5,5) + (5 + 15*sin(2*pi*time/25)) + (x/1)
y deve exibir uma tendência cíclica ao longo do 'tempo' com variação aleatória:
plot(time,y)
E uma relação linear positiva com 'x' com variação aleatória:
plot(x,y)
Crie um modelo aditivo linear simples de "y ~ time + x":
require(nlme)
m1 <- gls(y ~ time + x, method="REML")
O modelo exibe padrões cíclicos claros nos resíduos quando plotados em relação ao 'tempo', como seria de esperar:
plot(time, m1$residuals)
E o que deve ser uma falta clara e agradável de qualquer padrão ou tendência nos resíduos quando plotados contra 'x':
plot(x, m1$residuals)
Um modelo simples de "y ~ time + x" que inclua uma estrutura de correlação autoregressiva de ordem 1 deve ajustar os dados muito melhor do que o modelo anterior devido à estrutura de autocorrelação, quando avaliada usando AIC:
m2 <- gls(y ~ time + x, correlation = corAR1(form=~time), method="REML")
AIC(m1,m2)
No entanto, o modelo ainda deve exibir resíduos idênticos 'temporalmente' autocorrelacionados:
plot(time, m2$residuals)
Muito obrigado por qualquer conselho.
Respostas:
Q1
Você está fazendo duas coisas erradas aqui. A primeira é uma coisa geralmente ruim; em geral, não mergulhe nos objetos de modelo e retire componentes. Aprenda a usar as funções do extrator, neste caso
resid()
. Nesse caso, você está obtendo algo útil, mas se você tivesse um tipo diferente de objeto de modelo, como um GLM deglm()
,mod$residuals
conteria resíduos de trabalho da última iteração IRLS e é algo que geralmente não deseja!A segunda coisa que você está fazendo de errado é algo que também me chamou a atenção. Os resíduos que você extraiu (e também teria extraído se você tivesse usado
resid()
) são os resíduos brutos ou de resposta. Essencialmente, essa é a diferença entre os valores ajustados e os valores observados da resposta, levando em consideração apenas os termos de efeitos fixos . Esses valores conterão a mesma autocorrelação residual quem1
a dos efeitos fixos (ou se você preferir, o preditor linear) são os mesmos nos dois modelos (~ time + x
).Para obter resíduos que incluem o termo de correlação especificado, você precisa dos resíduos normalizados . Você obtém estes fazendo:
Este (e outros tipos de resíduos disponíveis) é descrito em
?residuals.gls
:Por meio de comparação, aqui estão os ACFs dos resíduos brutos (resposta) e os resíduos normalizados
Para ver por que isso está acontecendo e onde os resíduos brutos não incluem o termo de correlação, considere o modelo que você ajustou
Onde
Os resíduos brutos, o padrão retornado por,
resid(m2)
são apenas da parte do preditor linear, portanto, deste bitQ2
Parece que você está tentando ajustar uma tendência não linear com uma função linear de
time
e explica a falta de ajuste à "tendência" com um AR (1) (ou outras estruturas). Se seus dados são parecidos com os dados de exemplo que você fornece aqui, eu ajustaria um GAM para permitir uma função suave das covariáveis. Esse modelo seriaonde
select = TRUE
aplica alguma contração extra para permitir que o modelo remova qualquer um dos termos do modelo.Este modelo dá
e possui termos suaves que são assim:
Os resíduos desse modelo também são mais bem comportados (resíduos brutos)
Agora uma palavra de cautela; há um problema com a suavização de séries temporais, na medida em que os métodos que decidem o quão suave ou distorcida são as funções assumem que os dados são independentes. O que isso significa em termos práticos é que a função suave de time (
s(time)
) pode caber em informações que são realmente erros aleatórios correlacionados automaticamente e não apenas na tendência subjacente. Portanto, você deve ter muito cuidado ao ajustar os smoothers aos dados de séries temporais.Há várias maneiras de contornar isso, mas uma maneira é mudar para ajustar o modelo através do
gamm()
qual chamalme()
internamente e que permite que você use ocorrelation
argumento usado para ogls()
modelo. Aqui está um exemplos(time)
s(time)
s(time)
O modelo com o AR (1) não representa uma melhoria significativa em relação ao modelo sem o AR (1):
Se olharmos para a estimativa de $ \ hat {\ rho}}, veremos
Phi
fonte