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

Question

Contexto

Esta pregunta se usa R, pero se acerca general de las cuestiones estadísticas.

Estoy analizando los efectos de los factores de mortalidad (% de mortalidad debido a las enfermedades y parasitismo) sobre la polilla de la tasa de crecimiento de población a lo largo del tiempo, donde las poblaciones de larvas fueron muestreados de 12 sitios una vez al año durante 8 años. La tasa de crecimiento poblacional de los datos muestra una clara, pero de forma irregular tendencia cíclica a lo largo del tiempo.

Los residuos de un simple modelo lineal generalizado (tasa de crecimiento ~ %enfermedad + %parasitismo + año) muestran igualmente claro, pero de forma irregular tendencia cíclica a lo largo del tiempo. Por lo tanto, la generalización de los mínimos cuadrados de los modelos de la misma forma también fueron equipados con los datos con la adecuada correlación de estructuras para lidiar con la autocorrelación temporal, por ejemplo, el compuesto de la simetría, proceso autorregresivo de orden 1 y autorregresivo de media móvil de correlación de las estructuras.

Los modelos de todos los contenidos de la misma de efectos fijos, se compararon mediante la AIC, y fueron equipados por REML (para permitir la comparación de diferentes estructuras de correlación por AIC). Estoy utilizando el paquete R nlme y el gls función.

Pregunta 1

El GLS modelos de los residuos sigue siendo la pantalla casi idénticos patrones cíclicos, cuando conspiró contra el tiempo. Estos patrones de siempre, incluso en los modelos que explicar con precisión la estructura de autocorrelación?

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

Pregunta 2

He montado GLS modelos de la mejor manera posible y plausible de correlación de estructuras con mis datos, pero no son la realidad sustancialmente mejor ajuste que el GLM sin ningún tipo de correlación estructura: sólo una GLS modelo es ligeramente mejor (AIC puntuación = 1.8 inferior), mientras que el resto tienen mayores valores de AIC. Sin embargo, este es sólo el caso cuando todos los modelos están equipados por REML, no ML, donde GLS modelos son claramente mucho mejor, pero entiendo que a partir de las estadísticas de los libros sólo se debe utilizar REML para comparar modelos con diferente correlación de las estructuras y los mismos efectos fijos por razones que no voy a detallar aquí.

Dada la claridad temporalmente autocorrelated la naturaleza de los datos, si no hay ningún modelo es incluso moderadamente mejor que la simple GLM ¿cuál es la forma más adecuada para decidir el modelo a utilizar para la inferencia, suponiendo que yo estoy usando un método adecuado (que finalmente se desea utilizar AIC comparar diferentes combinaciones variables)?

Q1 'simulación' explorando residual de los patrones en los modelos con y sin la adecuada correlación de estructuras

Generar la simulación y la respuesta de la variable con un efecto cíclico de 'tiempo', y un positivo efecto lineal 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 sobre 'el tiempo' con la variación aleatoria:

plot(time,y)

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

plot(x,y)

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

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

El modelo muestra clara patrones cíclicos en los residuos cuando se trazan en contra de 'el tiempo', como sería de esperar:

plot(time, m1$residuals)

Y lo que debería ser un agradable, clara, la falta de cualquier patrón o tendencia en los residuos cuando se trazan contra la 'x':

plot(x, m1$residuals)

Un modelo simple de "y ~ + x" que incluye un proceso autorregresivo de correlación estructura de orden 1 debe ajustarse a los datos mucho mejor que el anterior modelo, ya que de la autocorrelación de la estructura, cuando se evaluó el uso de AIC:

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

Sin embargo, el modelo debe mostrar de manera casi idéntica 'temporalmente' autocorrelated residuales:

plot(time, m2$residuals)

Muchas gracias por cualquier consejo.

Answer 1

Q1

Estás haciendo dos cosas mal aquí. La primera es generalmente una mala cosa; no en general profundizar en el modelo de objetos y de extraer de los componentes. Aprender a usar el extractor de funciones, en este caso, resid(). En este caso, usted está consiguiendo algo útil, pero si había un tipo diferente de modelo de objetos, tales como un GLM de glm(), a continuación, mod$residuals contendría trabajo de los residuos de la última de las NIÑAS de la iteración y son algo que generalmente no quiere!!!

La segunda cosa que usted está haciendo mal es algo que me sorprendió demasiado. Los residuos extraídos (y también han extraído si hubiera usado resid()) son las primas o la respuesta de los residuos. Básicamente esta es la diferencia entre los valores ajustados y los valores observados de la respuesta, teniendo en cuenta los efectos fijos de los términos sólo. Estos valores contienen la misma residual de autocorrelación como la de m1 debido a los efectos fijos (o si se prefiere, el predictor lineal) es el mismo en los dos modelos (~ time + x).

Para obtener los residuos que incluyen la correlación plazo especificado, se necesita la normalizado de los residuos. Te estas haciendo:

resid(m1, type = "normalized")

Este (y otros tipos de residuos disponibles) se describe 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"'.

Por medio de la comparación, aquí están las ACFs de la materia prima (la respuesta) y la cantidad normalizada de residuos

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

Para ver por qué esto está sucediendo, y donde las primas de los residuos no incluyen la correlación plazo, considere el modelo ajustado

$$y = \beta_0 + \beta_1 \mathrm{time} + \beta_2 \mathrm{x} + \varepsilon$$

where

$$ \varepsilon \sim N(0, \sigma^2 \Lambda) $$

and $\Lambda$ is a correlation matrix defined by an AR(1) with parameter $\hat{\rho}$ where the non-diagonal elements of the matrix are filled with values $\rho^{|d|}$, where $d$ is the positive integer separation in time units of pairs of residuals.

The raw residuals, the default returned by resid(m2) are from the linear predictor part only, hence from this bit

$$ \beta_0 + \beta_1 \mathrm{time} + \beta_2 \mathrm{x} $$

and hence they contain none of the information on the correlation term(s) included in $\Lambda$.

Q2

It seems you are trying to fit a non-linear trend with a linear function of time and account for lack of fit to the "trend" with an AR(1) (or other structures). If your data are anything like the example data you give here, I would fit a GAM to allow for a smooth function of the covariates. This model would be

$$y = \beta_0 + f_1(\mathrm{time}) + f_2(\mathrm{x}) + \varepsilon$$

and initially we'll assume the same distribution as for the GLS except that initially we'll assume that $\Lambda = \mathbf{I}$ (an identity matrix, so independent residuals). This model can be fitted using

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

where select = TRUE applies some extra shrinkage to allow the model to remove either of the terms from the model.

This model gives

> 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

and has smooth terms that look like this:

The residuals from this model are also better behaved (raw residuals)

acf(resid(m3))

Now a word of caution; there is an issue with smoothing time series in that the methods that decide how smooth or wiggly the functions are assumes that the data are independent. What this means in practical terms is that the smooth function of time (s(time)) could fit information that is really random autocorrelated error and not only the underlying trend. Hence you should be very careful when fitting smoothers to time series data.

There are a number of ways round this, but one way is to switch to fitting the model via gamm() which calls lme() internally and which allows you to use the correlation argument you used for the gls() model. Here is an example

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))

Note that I have to fix the degrees of freedom for s(time) as there is an identifiability issue with these data. The model could be a wiggly s(time) and no AR(1) ($\rho = 0$) or a linear s(time) (1 degree of freedom) and a strong AR(1) ($\rho >> .5$). Hence I make an educated guess as to the complexity of the underlying trend. (Note I didn't spend much time on this dummy data, but you should look at the data and use your existing knowledge of the variability in time to determine an appropriate number of degrees of freedom for the spline.)

The model with the AR(1) does not represent a significant improvement over the model without the 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

If we look at the estimate for $\hat{\rho}} vemos

> intervals(mm2$lme)
....

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

donde Phi es lo que se llama $\rho$. Por lo tanto, 0 es un valor posible para $\rho$. La estimación es ligeramente mayor que cero, por lo que tendrá un efecto despreciable en el ajuste del modelo y, por tanto, podría dejarlo en el modelo si no hay un fuerte a priori de la razón para asumir residual de autocorrelación.