Tendência em dados irregulares de séries temporais

Eu tenho um conjunto de dados de medições de temperatura da água tiradas de um grande corpo de água em intervalos irregulares durante um período de décadas. (Galveston Bay, TX, se você estiver interessado)

Aqui está o chefe dos dados:

  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

(MEDIÇÃO é a medição de temperatura de interesse.)

O conjunto completo está disponível aqui: https://github.com/jscarlton/galvBayData/blob/master/gbtemp.csv

Gostaria de remover os efeitos da variação sazonal para observar a tendência (se houver) da temperatura ao longo do tempo. A decomposição de séries temporais é a melhor maneira de fazer isso? Como lidar com o fato de que as medições não foram realizadas em intervalos regulares? Espero que exista um pacote R para esse tipo de análise, embora Python ou Stata também fiquem bem. 

(Nota: para esta análise, estou optando por ignorar a variabilidade espacial nas medições. Idealmente, eu também consideraria isso, mas acho que isso seria irremediavelmente complexo.)

Em vez de tentar decompor explicitamente a série cronológica, sugiro que você modele os dados espaço-temporalmente porque, como você verá abaixo, a tendência de longo prazo provavelmente varia espacialmente, a tendência sazonal varia de acordo com a tendência de longo prazo e espacialmente.

Eu descobri que os modelos aditivos generalizados (GAMs) são um bom modelo para ajustar séries temporais irregulares, como você descreve.

Abaixo ilustro um modelo rápido que preparei para os dados completos do seguinte formulário

$$\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}$$

Onde

• $\alpha$ é o modelo de interceptação,
• $f_1(\text{ToD}_i)$ é uma função suave da hora do dia,
• $f_2(\text{DoY}_i)$ é uma função suave do dia do ano,
• $f_3(\text{Year}_i)$ é uma função suave do ano,
• $f_4(\text{x}_i, \text{y}_i)$ é um 2D suave de longitude e latitude,
• $f_5(\text{DoY}_i, \text{Year}_i)$ é um produto tensorial suave de dia do ano e ano,
• $f_6(\text{x}_i, \text{y}_i, \text{ToD}_i)$ produto tensorial suave de localização e hora do dia
• $f_7(\text{x}_i, \text{y}_i, \text{DoY}_i)$ produto tensorial suave do local dia do ano &
• $f_8(\text{x}_i, \text{y}_i, \text{Year}_i$ produto tensorial suave de localização e ano

Efetivamente, os quatro primeiros smooths são os principais efeitos de

1. hora do dia,
2. estação,
3. tendência de longo prazo,
4. variação espacial

enquanto os quatro produtos tensores restantes suavizam o modelo de interações suaves entre as covariáveis ​​indicadas, que modelam

5. como o padrão sazonal de temperatura varia ao longo do tempo,
6. como o efeito da hora do dia varia espacialmente,
7. como o efeito sazonal varia espacialmente e
8. como a tendência de longo prazo varia espacialmente

Os dados são carregados no R e massageados um pouco com o seguinte 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))

O próprio modelo é ajustado usando a bam()função projetada para ajustar os GAMs a conjuntos de dados maiores como esse. Você também pode usar gam()para esse modelo, mas levará um pouco mais de tempo para ajustar.

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)

Os s()termos são os principais efeitos, enquanto os ti()termos são interação do produto tensorial, onde os principais efeitos das covariáveis ​​nomeadas foram removidos da base. Esses ti()suaves são uma maneira de incluir interações das variáveis ​​declaradas de maneira numericamente estável.

O knotsobjeto está apenas definindo os pontos de extremidade da suavização cíclica que usei para o efeito do dia do ano - queremos que as 23:59 de 31 de dezembro se unam suavemente às 00:01 de 1º de janeiro. Isso explica, em certa medida, os anos bissextos.

O resumo do modelo indica que todos esses efeitos são 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

Uma análise mais cuidadosa gostaria de verificar se precisamos de todas essas interações; alguns dos ti()termos espaciais explicam apenas pequenas quantidades de variação nos dados, conforme indicado pela estatística ; há muitos dados aqui; mesmo os tamanhos de efeitos pequenos podem ser estatisticamente significativos, mas desinteressantes.$F$

Como uma verificação rápida, no entanto, a remoção dos três ti()suavizações espaciais ( m.sub) resulta em um ajuste significativamente pior conforme avaliado pela AIC:

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

Podemos plotar os efeitos parciais dos cinco primeiros suaves usando o plot()método - o produto do tensor 3D não pode ser plotado facilmente e não por padrão.

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

O scale = 0argumento lá coloca todos os gráficos em sua própria escala, para comparar as magnitudes dos efeitos, podemos desativar isso:

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

Agora podemos ver que o efeito sazonal domina. A tendência de longo prazo (em média) é mostrada no gráfico superior direito. No entanto, para realmente observar a tendência de longo prazo, você precisa escolher uma estação e prever a partir do modelo dessa estação, fixando a hora do dia e o dia do ano com alguns valores representativos (meio-dia, para um dia do ano no verão) dizer). No início do ano ou dois da série, existem alguns valores de temperatura baixa em relação ao restante dos registros, o que provavelmente está sendo detectado em todos os processos de transição YEAR. Esses dados devem ser analisados ​​mais de perto.

Este não é realmente o lugar para entrar nisso, mas aqui estão algumas visualizações do modelo. Primeiro, analiso o padrão espacial da temperatura e como ela varia ao longo dos anos da série. Para fazer isso, prevejo a partir do modelo de uma grade 100x100 no domínio espacial, no meio do dia no dia 180 de cada ano:

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)

em seguida, defino como ausentes NAos valores previstos fitpara todos os pontos de dados que estão a alguma distância das observações (proporcional; dist)

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

e associe as previsões aos dados de previsão

pred <- cbind(pdata, Fitted = fit)

Definir valores previstos NAcomo esse nos impede de extrapolar além do suporte dos dados.

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

obtemos o seguinte

Podemos ver a variação ano a ano nas temperaturas com mais detalhes se animarmos, em vez de facetá-lo,

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 analisar as tendências de longo prazo com mais detalhes, podemos prever estações específicas. Por exemplo, para STATION_ID13364 e prevendo dias nos quatro trimestres, podemos usar o seguinte para preparar valores das covariáveis ​​que desejamos prever (meio-dia, no dia do ano 1, 90, 180 e 270, na estação selecionada e avaliar a tendência de longo prazo em 500 valores igualmente espaçados)

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

Prevemos e solicitamos erros padrão, para formar um intervalo de confiança aproximado de 95%

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 traçamos

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

produzindo

Obviamente, há muito mais para modelar esses dados do que o que mostro aqui, e gostaríamos de verificar se há autocorrelação residual e sobreajuste dos splines, mas abordar o problema como uma das modelagens dos recursos dos dados permite uma análise mais detalhada exame das tendências.

É claro que você poderia apenas modelar cada um STATION_IDseparadamente, mas isso descartaria os dados, e muitas estações têm poucas observações. Aqui, o modelo empresta todas as informações da estação para preencher as lacunas e ajudar na estimativa das tendências de interesse.

Algumas notas sobre bam()

O bam()modelo está usando todos os truques de mgcv para estimar o modelo rapidamente (vários threads 4 ), seleção rápida de suavidade REML ( method = 'fREML') e discretização de covariáveis. Com essas opções ativadas, o modelo se encaixa em menos de um minuto na minha estação de trabalho Xeon de quatro núcleos e era de 2013 com 64 GB de RAM.

Você pode usar a função decomposeque separa suas séries temporais em três componentes: tendência, sazonal e aleatória. Você também pode extrair os diferentes valores do resultado e plotá-los. Certifique-se de definir seus dados como uma série temporal. A função stlestá basicamente fazendo o mesmo, mas oferece mais flexibilidade na forma de decompor seus dados.

Eu também recomendo o seguinte site

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

Isso ajuda?

Uma opção seria usar uma estratégia de regressão que trata seus dados como extraídos de uma função subjacente contínua (consulte a discussão aqui: existe algum padrão-ouro para modelar séries temporais com espaçamento irregular? ). A partir daí, você pode usar um método como o Singular Spectrum Analysis para decompor o sinal (para R: https://cran.r-project.org/web/packages/Rssa/index.html ).