¿Es posible especificar un término de autocorrelación anidado al trabajar con una estructura jerárquica (GAM)?

Question

Estoy modelando la ocurrencia de una especie en 5 sitios diferentes de forma horaria (presencia/ausencia), basado en una serie de predictores temporales (por ejemplo, momento del año, ciclo día/noche, mareas...). Las covariables están indicadas por x1, x2, ... en el código a continuación. Para información, tengo ~70 000 puntos de datos.

Estoy utilizando una estructura de HGAM, como se introdujo en Pedersen et al. 2019. Para cada covariable, primero investigué diferentes opciones de especificación (suavizador global o no, penalización compartida o no...), y seleccioné la mejor basada en AIC y mi pregunta de investigación.

Cuando pongo todos los términos juntos en el modelo, termino con una estructura como esta:

model <- bam(response ~ offset(log(offset)) + s(year, bs="re") + 
      Site + s(x1, m=2, bs="cc", k=8) + 
      s(x1, Site, bs = "fs", xt=list(bs="cc"), m=2) + 
      s(x2, bs = "cc", by=Site, m=2, k=8) + 
      s(x3, m=2, bs="cc", k=10) + 
      s(x3, by = Site, bs = "cc", m=1, k=10) + 
      s(x4, bs = "cc", by=Site, m=2, k=8) + 
      s(x5, bs = "tp", by=Site, m=2, k=8), family = "binomial", 
      data = data.all, method="REML", cluster=cl, select=TRUE)

La devianza explicada del modelo es solo alrededor del 13%. Dada la alta resolución temporal (horaria), espero cierta autocorrelación temporal en los residuos. Leí en la documentación de bam(), y en otros posts (aquí y allí), que esto se podría especificar usando el argumento rho en bam() siguiendo el orden del conjunto de datos.

• Sin embargo, no estoy seguro de si esto se aplica cuando el modelo tiene una estructura jerárquica. Sé que es factible en gamm(), pero el problema entonces es que no puedo especificar múltiples términos de interacción suavizador-factor como está escrito actualmente en el modelo...
• ¿Es algo que se podría abordar en brms? ¿O de alguna otra manera?

Para información, aquí está la salida de las funciones de gráficos acf y pcaf:

Answer 1

1. La estructura jerárquica es principalmente para los efectos "fijos", pero debes considerarla cuando se trata de una estructura de autocorrelación. Con ambos correlation en gamm() y rho en bam(), a menos que se indique lo contrario, el modelo asumirá una serie temporal larga única.

   Este comportamiento puede ser el que deseas, pero es posible que, por ejemplo en tu caso, quieras anidar la estructura de correlación dentro del día en lugar de simplemente el orden temporal de las muestras. Dicha anidación indicaría que la estructura de autocorrelación opera en el nivel de agrupamiento más bajo (dentro del día) — no puedes optar por anidarla dentro del año cuando tienes datos a nivel de subdiario a menos que crees una variable que ordene las muestras a nivel de subdiario dentro del año.

   Con correlation utilizarías corAR1(form =~ t | year) para indicar que deseas un AR(1) anidado dentro del año, donde t ordena las observaciones dentro del year. Para obtener esto anidado dentro del site solo necesitas agregar el lado derecho de la fórmula. Una forma de hacer esto es crear una variable year_site = interaction(year, site, drop = TRUE) en tus datos y luego modificar la fórmula para que sea corAR1(form =~ t | year_site). Puede haber formas más directas de hacer esto dentro de la fórmula, pero todos los ejemplos que he visto utilizan un solo factor de agrupamiento. Con rho debes crear un vector lógico que sea FALSE en todas partes, excepto en la primera observación de cada "serie" en la que deseas que la estructura de correlación esté anidada.

   Es importante tener en cuenta que estás estimando el mismo parámetro único $\rho$ independientemente de cómo especifiques el AR(1), que opera dentro de cada nivel del factor de agrupamiento (o de manera equivalente para otros términos de ARMA). En ese sentido, los parámetros son globales en la forma en que describen una estructura común de autocorrelación dentro de cada nivel de la jerarquía.

2. brms funciona esencialmente de la misma manera. Consulta https://paul-buerkner.github.io/brms/reference/autocor-terms.html y luego las funciones de correlación individuales vinculadas desde esa página que documentan argumentos como time (para ordenar las observaciones) y gr para el factor de agrupamiento.