¿Los patrones residuales autocorrelacionados permanecen incluso en modelos con estructuras de correlación apropiadas y cómo seleccionar los mejores modelos?

Contexto

Esta pregunta usa R, pero trata sobre cuestiones estadísticas generales.

Estoy analizando los efectos de los factores de mortalidad (% de mortalidad por enfermedad y parasitismo) en la tasa de crecimiento de la población de polillas a lo largo del tiempo, donde se tomaron muestras de las poblaciones de larvas de 12 sitios una vez al año durante 8 años. Los datos de la tasa de crecimiento de la población muestran una tendencia cíclica clara pero irregular a lo largo del tiempo.

Los residuos de un modelo lineal generalizado simple (tasa de crecimiento ~% enfermedad +% parasitismo + año) mostraron una tendencia cíclica igualmente clara pero irregular a lo largo del tiempo. Por lo tanto, los modelos de mínimos cuadrados generalizados de la misma forma también se ajustaron a los datos con estructuras de correlación apropiadas para tratar la autocorrelación temporal, por ejemplo, simetría compuesta, orden de proceso autorregresivo 1 y estructuras de correlación de media móvil autorregresiva.

Todos los modelos contenían los mismos efectos fijos, se compararon usando AIC y se ajustaron mediante REML (para permitir la comparación de diferentes estructuras de correlación por AIC). Estoy usando el paquete R nlme y la función gls.

Pregunta 1

Los residuos de los modelos GLS todavía muestran patrones cíclicos casi idénticos cuando se trazan contra el tiempo. ¿Permanecerán siempre tales patrones, incluso en modelos que representan con precisión la estructura de autocorrelación?

He simulado algunos datos simplificados pero similares en R debajo de mi segunda pregunta, que muestra el problema basado en mi comprensión actual de los métodos necesarios para evaluar patrones temporalmente autocorrelacionados en los residuos del modelo , que ahora sé que están equivocados (ver respuesta).

Pregunta 2

He ajustado modelos GLS con todas las posibles estructuras de correlación plausibles a mis datos, pero en realidad ninguno es sustancialmente mejor que el GLM sin ninguna estructura de correlación: solo un modelo GLS es marginalmente mejor (puntaje AIC = 1.8 menor), mientras que el resto tiene valores de AIC más altos. Sin embargo, este es solo el caso cuando REML ajusta todos los modelos, no ML donde los modelos GLS son claramente mucho mejores, pero entiendo de los libros de estadísticas que solo debe usar REML para comparar modelos con diferentes estructuras de correlación y los mismos efectos fijos por razones No detallaré aquí.

Dada la naturaleza claramente autocorrelacionada temporalmente de los datos, si ningún modelo es incluso moderadamente mejor que el GLM simple, ¿cuál es la forma más adecuada de decidir qué modelo usar para la inferencia, suponiendo que esté usando un método apropiado (eventualmente quiero usar AIC para comparar diferentes combinaciones de variables)?

Q1 'simulación' que explora patrones residuales en modelos con y sin estructuras de correlación apropiadas

Genere una variable de respuesta simulada con un efecto cíclico de 'tiempo' y un efecto lineal 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 debe mostrar una tendencia cíclica a lo largo del 'tiempo' con variación aleatoria:

plot(time,y)

Y una relación lineal positiva con 'x' con variación aleatoria:

plot(x,y)

Cree un modelo aditivo lineal simple de "y ~ time + x":

require(nlme)
m1 <- gls(y ~ time + x, method="REML")

El modelo muestra patrones cíclicos claros en los residuos cuando se grafica contra el 'tiempo', como se esperaría:

plot(time, m1$residuals)

Y lo que debería ser una falta clara y agradable de cualquier patrón o tendencia en los residuos cuando se traza contra 'x':

plot(x, m1$residuals)

Un modelo simple de "y ~ tiempo + x" que incluye una estructura de correlación autorregresiva de orden 1 debería ajustarse a los datos mucho mejor que el modelo anterior debido a la estructura de autocorrelación, cuando se evalúa usando AIC:

m2 <- gls(y ~ time + x, correlation = corAR1(form=~time), method="REML")
AIC(m1,m2)

Sin embargo, el modelo aún debería mostrar residuos autocorrelacionados casi idénticamente 'temporalmente':

plot(time, m2$residuals)

Muchas gracias por cualquier consejo.

Q1

Estás haciendo dos cosas mal aquí. El primero es algo generalmente malo; hacer no en ahondar en general en los objetos del modelo y arrancar componentes. Aprenda a usar las funciones del extractor, en este caso resid(). En este caso, está obteniendo algo útil, pero si tuviera un tipo diferente de objeto modelo, como un GLM glm(), ¡entonces mod$residualscontendría residuos de trabajo de la última iteración de IRLS y es algo que generalmente no desea!

La segunda cosa que estás haciendo mal es algo que también me ha pillado. Los residuos que extrajo (y también habría extraído si los hubiera usado resid()) son los residuos sin procesar o de respuesta. Esencialmente, esta es la diferencia entre los valores ajustados y los valores observados de la respuesta, teniendo en cuenta únicamente los términos de efectos fijos . Estos valores contendrán la misma autocorrelación residual que la de m1porque los efectos fijos (o si lo prefiere, el predictor lineal) son los mismos en los dos modelos ( ~ time + x).

Para obtener los residuos que incluyen el término de correlación que especificó, necesita los residuos normalizados . Obtienes estos haciendo:

resid(m1, type = "normalized")

Este (y otros tipos de residuos disponibles) se describen en ?residuals.gls:

type: an optional character string specifying the type of residuals
      to be used. If ‘"response"’, the "raw" residuals (observed -
      fitted) are used; else, if ‘"pearson"’, the standardized
      residuals (raw residuals divided by the corresponding
      standard errors) are used; else, if ‘"normalized"’, the
      normalized residuals (standardized residuals pre-multiplied
      by the inverse square-root factor of the estimated error
      correlation matrix) are used. Partial matching of arguments
      is used, so only the first character needs to be provided.
      Defaults to ‘"response"’.

A modo de comparación, aquí están los ACF de la materia prima (respuesta) y los residuos normalizados.

layout(matrix(1:2))
acf(resid(m2))
acf(resid(m2, type = "normalized"))
layout(1)

Para ver por qué sucede esto y dónde los residuos brutos no incluyen el término de correlación, considere el modelo que ajustó

dónde

$\Lambda$$\hat{\rho}$$\rho^{|d|}$$d$

Los residuos brutos, los valores predeterminados devueltos por, resid(m2)son solo de la parte del predictor lineal, por lo tanto, de este bit

$\Lambda$

Q2

Parece que está tratando de ajustar una tendencia no lineal con una función lineal timey explicar la falta de ajuste a la "tendencia" con un AR (1) (u otras estructuras). Si sus datos son similares a los datos de ejemplo que da aquí, ajustaría un GAM para permitir una función fluida de las covariables. Este modelo seria

$\Lambda = \mathbf{I}$

library("mgcv")
m3 <- gam(y ~ s(time) + s(x), select = TRUE, method = "REML")

donde se select = TRUEaplica una contracción adicional para permitir que el modelo elimine cualquiera de los términos del modelo.

Este modelo da

> summary(m3)

Family: gaussian 
Link function: identity 

Formula:
y ~ s(time) + s(x)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  23.1532     0.7104   32.59   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
          edf Ref.df      F  p-value    
s(time) 8.041      9 26.364  < 2e-16 ***
s(x)    1.922      9  9.749 1.09e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

y tiene términos suaves que se ven así:

Los residuos de este modelo también se comportan mejor (residuos brutos)

acf(resid(m3))

Ahora una palabra de precaución; Hay un problema con el suavizado de series de tiempo en que los métodos que deciden qué tan suaves o onduladas son las funciones suponen que los datos son independientes. Lo que esto significa en términos prácticos es que la función fluida del tiempo ( s(time)) podría ajustarse a información que es realmente un error autocorrelacionado aleatorio y no solo la tendencia subyacente. Por lo tanto, debe tener mucho cuidado al ajustar suavizadores a datos de series temporales.

Hay varias formas de evitar esto, pero una es cambiar para ajustar el modelo a través de gamm()las llamadas lme()internas y que le permite utilizar el correlationargumento que utilizó para el gls()modelo. Aquí hay un ejemplo

mm1 <- gamm(y ~ s(time, k = 6, fx = TRUE) + s(x), select = TRUE,
            method = "REML")
mm2 <- gamm(y ~ s(time, k = 6, fx = TRUE) + s(x), select = TRUE,
            method = "REML", correlation = corAR1(form = ~ time))

s(time)s(time)$\rho = 0$s(time)$\rho >> .5$

El modelo con el AR (1) no representa una mejora significativa sobre el modelo sin el AR (1):

> anova(mm1$lme, mm2$lme)
        Model df      AIC      BIC    logLik   Test   L.Ratio p-value
mm1$lme     1  9 301.5986 317.4494 -141.7993                         
mm2$lme     2 10 303.4168 321.0288 -141.7084 1 vs 2 0.1817652  0.6699

Si miramos la estimación de $ \ hat {\ rho}} vemos

> intervals(mm2$lme)
....

 Correlation structure:
         lower      est.     upper
Phi -0.2696671 0.0756494 0.4037265
attr(,"label")
[1] "Correlation structure:"

Phi$\rho$$\rho$