Tendencia en datos de series temporales irregulares

Tengo un conjunto de datos de mediciones de temperatura del agua tomadas de un gran cuerpo de agua a intervalos irregulares durante un período de décadas. (Galveston Bay, TX si estás interesado)

Aquí está la cabeza de los datos:

  STATION_ID     DATE  TIME LATITUDE LONGITUDE YEAR MONTH DAY SEASON MEASUREMENT
1      13296  6/20/91 11:04 29.50889 -94.75806 1991     6  20 Summer        28.0
2      13296  3/17/92  9:30 29.50889 -94.75806 1992     3  17 Spring        20.1
3      13296  9/23/91 11:24 29.50889 -94.75806 1991     9  23   Fall        26.0
4      13296  9/23/91 11:24 29.50889 -94.75806 1991     9  23   Fall        26.0
5      13296  6/20/91 11:04 29.50889 -94.75806 1991     6  20 Summer        28.0
6      13296 12/17/91 10:15 29.50889 -94.75806 1991    12  17 Winter        13.0

(MEDICIÓN es la medición de temperatura de interés).

El conjunto completo está disponible aquí: https://github.com/jscarlton/galvBayData/blob/master/gbtemp.csv

Me gustaría eliminar los efectos de la variación estacional para observar la tendencia (si la hay) en la temperatura a lo largo del tiempo. ¿Es una descomposición de series de tiempo la mejor manera de hacer esto? ¿Cómo manejo el hecho de que las mediciones no se tomaron a intervalos regulares? Espero que haya un paquete R para este tipo de análisis, aunque Python o Stata también estarían bien. 

(Nota: para este análisis, elijo ignorar la variabilidad espacial en las mediciones. Idealmente, también lo explicaría, pero creo que hacerlo sería irremediablemente complejo).

En lugar de tratar de descomponer las series de tiempo explícitamente, sugeriría que modele los datos espacio-temporalmente porque, como verá a continuación, la tendencia a largo plazo probablemente varía espacialmente, la tendencia estacional varía con la tendencia a largo plazo y espacialmente

He descubierto que los modelos aditivos generalizados (GAM) son un buen modelo para ajustar series temporales irregulares como las que usted describe.

A continuación, ilustro un modelo rápido que preparé para los datos completos del siguiente formulario

$$\begin{align} \begin{split} \mathrm{E}(y_i) & = \alpha + f_1(\text{ToD}_i) + f_2(\text{DoY}_i) + f_3(\text{Year}_i) + f_4(\text{x}_i, \text{y}_i) + \\ & \quad f_5(\text{DoY}_i, \text{Year}_i) + f_6(\text{x}_i, \text{y}_i, \text{ToD}_i) + \\ & \quad f_7(\text{x}_i, \text{y}_i, \text{DoY}_i) + f_8(\text{x}_i, \text{y}_i, \text{Year}_i) \end{split} \end{align}$$

dónde

• $\alpha$ es la intercepción del modelo,
• $f_1(\text{ToD}_i)$ es una función suave de la hora del día,
• $f_2(\text{DoY}_i)$ es una función fluida del día del año,
• $f_3(\text{Year}_i)$ es una función fluida del año,
• $f_4(\text{x}_i, \text{y}_i)$ es una suavidad 2D de longitud y latitud,
• $f_5(\text{DoY}_i, \text{Year}_i)$ es un producto tensor uniforme del día del año y año,
• $f_6(\text{x}_i, \text{y}_i, \text{ToD}_i)$ producto tensor sin problemas de ubicación y hora del día
• $f_7(\text{x}_i, \text{y}_i, \text{DoY}_i)$ producto tensor sin problemas de ubicación día del año &
• $f_8(\text{x}_i, \text{y}_i, \text{Year}_i$ producto tensor sin problemas de ubicación y año

Efectivamente, los primeros cuatro suavizados son los principales efectos de

1. hora del día
2. temporada,
3. tendencia a largo plazo,
4. variación espacial

mientras que los cuatro productos tensoriales restantes suavizan las interacciones suaves entre las covariables establecidas, qué modelo

5. cómo el patrón estacional de temperatura varía con el tiempo,
6. cómo el efecto de la hora del día varía espacialmente,
7. cómo varía espacialmente el efecto estacional, y
8. cómo la tendencia a largo plazo varía espacialmente

Los datos se cargan en R y se masajean un poco con el siguiente código

library('mgcv')
library('ggplot2')
library('viridis')
theme_set(theme_bw())
library('gganimate')

galveston <- read.csv('gbtemp.csv')
galveston <- transform(galveston,
                       datetime = as.POSIXct(paste(DATE, TIME),
                                             format = '%m/%d/%y %H:%M', tz = "CDT"))
galveston <- transform(galveston,
                       STATION_ID = factor(STATION_ID),
                       DoY = as.numeric(format(datetime, format = '%j')),
                       ToD = as.numeric(format(datetime, format = '%H')) +
                            (as.numeric(format(datetime, format = '%M')) / 60))

El modelo en sí se ajusta utilizando la bam()función diseñada para ajustar los GAM a conjuntos de datos más grandes como este. También puede usarlo gam()para este modelo, pero tardará un poco más en adaptarse.

knots <- list(DoY = c(0.5, 366.5))
M <- list(c(1, 0.5), NA)
m <- bam(MEASUREMENT ~
         s(ToD, k = 10) +
         s(DoY, k = 30, bs = 'cc') +
         s(YEAR, k = 30) +
         s(LONGITUDE, LATITUDE, k = 100, bs = 'ds', m = c(1, 0.5)) +
         ti(DoY, YEAR, bs = c('cc', 'tp'), k = c(15, 15)) +
         ti(LONGITUDE, LATITUDE, ToD, d = c(2,1), bs = c('ds','tp'), 
            m = M, k = c(20, 10)) +
         ti(LONGITUDE, LATITUDE, DoY, d = c(2,1), bs = c('ds','cc'),
            m = M, k = c(25, 15)) +
         ti(LONGITUDE, LATITUDE, YEAR, d = c(2,1), bs = c('ds','tp'),
            m = M), k = c(25, 15)),
         data = galveston, method = 'fREML', knots = knots,
         nthreads = 4, discrete = TRUE)

Los s()términos son los efectos principales, mientras que los ti()términos son tensor interacción del producto suaviza donde los efectos principales de las covariables nombradas se han eliminado de la base. Estos ti()suavizados son una forma de incluir interacciones de las variables indicadas de forma numéricamente estable.

El knotsobjeto es simplemente establecer los puntos finales de la suavidad cíclica que utilicé para el efecto del día del año: queremos que las 23:59 del 31 de diciembre se unan sin problemas con el 01 de enero a las 00:01. Esto explica en cierta medida los años bisiestos.

El resumen del modelo indica que todos estos efectos son significativos;

> summary(m)

Family: gaussian 
Link function: identity 

Formula:
MEASUREMENT ~ s(ToD, k = 10) + s(DoY, k = 12, bs = "cc") + s(YEAR, 
    k = 30) + s(LONGITUDE, LATITUDE, k = 100, bs = "ds", m = c(1, 
    0.5)) + ti(DoY, YEAR, bs = c("cc", "tp"), k = c(12, 15)) + 
    ti(LONGITUDE, LATITUDE, ToD, d = c(2, 1), bs = c("ds", "tp"), 
        m = list(c(1, 0.5), NA), k = c(20, 10)) + ti(LONGITUDE, 
    LATITUDE, DoY, d = c(2, 1), bs = c("ds", "cc"), m = list(c(1, 
        0.5), NA), k = c(25, 12)) + ti(LONGITUDE, LATITUDE, YEAR, 
    d = c(2, 1), bs = c("ds", "tp"), m = list(c(1, 0.5), NA), 
    k = c(25, 15))

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 21.75561    0.07508   289.8   <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(ToD)                        3.036   3.696    5.956 0.000189 ***
s(DoY)                        9.580  10.000 3520.098  < 2e-16 ***
s(YEAR)                      27.979  28.736   59.282  < 2e-16 ***
s(LONGITUDE,LATITUDE)        54.555  99.000    4.765  < 2e-16 ***
ti(DoY,YEAR)                131.317 140.000   34.592  < 2e-16 ***
ti(ToD,LONGITUDE,LATITUDE)   42.805 171.000    0.880  < 2e-16 ***
ti(DoY,LONGITUDE,LATITUDE)   83.277 240.000    1.225  < 2e-16 ***
ti(YEAR,LONGITUDE,LATITUDE)  84.862 329.000    1.101  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.94   Deviance explained = 94.2%
fREML =  29807  Scale est. = 2.6318    n = 15276

Un análisis más cuidadoso querría verificar si necesitamos todas estas interacciones; algunos de los ti()términos espaciales explican solo pequeñas cantidades de variación en los datos, como lo indica la estadística ; hay una gran cantidad de datos aquí, por lo que incluso los tamaños de efectos pequeños pueden ser estadísticamente significativos pero poco interesantes.$F$

Sin embargo, como una verificación rápida, al eliminar los tres ti()suavizados espaciales ( m.sub), se obtiene un ajuste significativamente peor según lo evaluado por AIC:

> AIC(m, m.sub)
            df      AIC
m     447.5680 58583.81
m.sub 239.7336 59197.05

Podemos trazar los efectos parciales de los primeros cinco suavizados utilizando el plot()método: los suavizadores del producto tensor 3D no se pueden graficar fácilmente y no de forma predeterminada.

plot(m, pages = 1, scheme = 2, shade = TRUE, scale = 0)

El scale = 0argumento allí pone todas las tramas en su propia escala, para comparar las magnitudes de los efectos, podemos desactivar esto:

plot(m, pages = 1, scheme = 2, shade = TRUE)

Ahora podemos ver que domina el efecto estacional. La tendencia a largo plazo (en promedio) se muestra en la gráfica superior derecha. Sin embargo, para ver realmente la tendencia a largo plazo, debe elegir una estación y luego predecir del modelo para esa estación, fijando la hora del día y el día del año en algunos valores representativos (mediodía, para un día del año en verano decir). A principios de año o dos de la serie tiene algunos valores de baja temperatura en relación con el resto de los registros, que probablemente se está recogiendo en todos los suavizados que involucran YEAR. Estos datos deberían analizarse más de cerca.

Este no es realmente el lugar para entrar en eso, pero aquí hay un par de visualizaciones de los ajustes del modelo. Primero miro el patrón espacial de temperatura y cómo varía a lo largo de los años de la serie. Para hacer eso, predigo a partir del modelo para una cuadrícula de 100x100 sobre el dominio espacial, al mediodía del día 180 de cada año:

pdata <- with(galveston,
              expand.grid(ToD = 12,
                          DoY = 180,
                          YEAR = seq(min(YEAR), max(YEAR), by = 1),
                          LONGITUDE = seq(min(LONGITUDE), max(LONGITUDE), length = 100),
                          LATITUDE  = seq(min(LATITUDE), max(LATITUDE), length = 100)))
fit <- predict(m, pdata)

luego puse en falta, NAlos valores predichos fitpara todos los puntos de datos que se encuentran a cierta distancia de las observaciones (proporcional; dist)

ind <- exclude.too.far(pdata$LONGITUDE, pdata$LATITUDE,
                       galveston$LONGITUDE, galveston$LATITUDE, dist = 0.1)
fit[ind] <- NA

y unir las predicciones a los datos de predicción

pred <- cbind(pdata, Fitted = fit)

Establecer valores pronosticados NAcomo este nos impide extrapolar más allá del soporte de los datos.

Usando ggplot2

ggplot(pred, aes(x = LONGITUDE, y = LATITUDE)) +
    geom_raster(aes(fill = Fitted)) + facet_wrap(~ YEAR, ncol = 12) +
    scale_fill_viridis(name = expression(degree*C), option = 'plasma',
                       na.value = 'transparent') +
    coord_quickmap() +
    theme(legend.position = 'top', legend.key.width = unit(2, 'cm'))

obtenemos lo siguiente

Podemos ver la variación de las temperaturas de un año a otro con un poco más de detalle si animamos en lugar de facetar la trama

p <- ggplot(pred, aes(x = LONGITUDE, y = LATITUDE, frame = YEAR)) +
    geom_raster(aes(fill = Fitted)) +
    scale_fill_viridis(name = expression(degree*C), option = 'plasma',
                       na.value = 'transparent') +
    coord_quickmap() +
    theme(legend.position = 'top', legend.key.width = unit(2, 'cm'))+
    labs(x = 'Longitude', y = 'Latitude')

gganimate(p, 'galveston.gif', interval = .2, ani.width = 500, ani.height = 800)

Para ver las tendencias a largo plazo con más detalle, podemos predecir para estaciones particulares. Por ejemplo, para STATION_ID13364 y ​​predecir días en los cuatro trimestres, podríamos usar lo siguiente para preparar los valores de las covariables que queremos predecir (mediodía, el día del año 1, 90, 180 y 270, en la estación seleccionada y evaluar la tendencia a largo plazo en 500 valores igualmente espaciados)

pdata <- with(galveston,
              expand.grid(ToD = 12,
                          DoY = c(1, 90, 180, 270),
                          YEAR = seq(min(YEAR), max(YEAR), length = 500),
                          LONGITUDE = -94.8751,
                          LATITUDE  = 29.50866))

Luego predecimos y pedimos errores estándar, para formar un intervalo de confianza aproximado del 95% puntual

fit <- data.frame(predict(m, newdata = pdata, se.fit = TRUE))
fit <- transform(fit, upper = fit + (2 * se.fit), lower = fit - (2 * se.fit))
pred <- cbind(pdata, fit)

que tramamos

ggplot(pred, aes(x = YEAR, y = fit, group = factor(DoY))) +
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = 'grey', alpha = 0.5) +
    geom_line() + facet_wrap(~ DoY, scales = 'free_y') +
    labs(x = NULL, y = expression(Temperature ~ (degree * C)))

productor

Obviamente, hay mucho más para modelar estos datos que lo que muestro aquí, y nos gustaría verificar la autocorrelación residual y el sobreajuste de las splines, pero abordar el problema como uno de modelar las características de los datos permite una información más detallada examen de las tendencias.

Por supuesto, podría modelar cada uno por STATION_IDseparado, pero eso arrojaría datos, y muchas estaciones tienen pocas observaciones. Aquí el modelo toma prestada toda la información de la estación para llenar los vacíos y ayudar a estimar las tendencias de interés.

Algunas notas sobre bam()

El bam()modelo está utilizando todos los trucos de mgcv para estimar el modelo rápidamente (múltiples hilos 4 ), selección rápida de suavidad REML ( method = 'fREML') y discretización de covariables. Con estas opciones activadas, el modelo cabe en menos de un minuto en mi estación de trabajo Xeon dual de 4 núcleos de la era 2013 con 64 Gb de RAM.

Simplemente puede usar la función decomposeque separa sus series de tiempo en tres componentes, tendencia, estacional y aleatorio. También puede extraer los diferentes valores del resultado y trazarlos. Asegúrese de definir sus datos como una serie de tiempo. La función stlbásicamente hace lo mismo, pero le brinda más flexibilidad sobre cómo puede descomponer sus datos.

También recomiendo el siguiente sitio web

https://www.otexts.org/fpp/6

¿Esto ayuda?

Una opción sería usar una estrategia de regresión que trate sus datos como extraídos de una función subyacente continua (vea la discusión aquí: ¿Hay algún estándar de oro para modelar series de tiempo espaciadas irregularmente? ). Desde allí, puede usar un método como Singular Spectrum Analysis para descomponer la señal (para R: https://cran.r-project.org/web/packages/Rssa/index.html ).