Tendance des données de séries chronologiques irrégulières

J'ai un ensemble de données de mesures de la température de l'eau prises à partir d'un grand plan d'eau à intervalles irréguliers sur une période de plusieurs décennies. (Galveston Bay, TX si vous êtes intéressé)

Voici la tête des données:

  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

(LA MESURE est la mesure de température qui nous intéresse.)

L'ensemble complet est disponible ici: https://github.com/jscarlton/galvBayData/blob/master/gbtemp.csv

Je voudrais supprimer les effets des variations saisonnières pour observer la tendance (le cas échéant) de la température au fil du temps. La décomposition de séries chronologiques est-elle la meilleure façon de procéder? Comment gérer le fait que les mesures n'ont pas été prises à intervalle régulier? J'espère qu'il existe un package R pour ce type d'analyse, bien que Python ou Stata convienne également. 

(Remarque: pour cette analyse, je choisis d'ignorer la variabilité spatiale dans les mesures. Idéalement, je prendrais cela en compte également, mais je pense que cela serait désespérément complexe.)

Plutôt que d'essayer de décomposer la série chronologique de manière explicite, je suggère plutôt de modéliser les données spatio-temporellement car, comme vous le verrez ci-dessous, la tendance à long terme varie probablement spatialement, la tendance saisonnière varie avec la tendance à long terme et spatialement.

J'ai trouvé que les modèles additifs généralisés (GAM) sont un bon modèle pour ajuster des séries temporelles irrégulières comme vous le décrivez.

Ci-dessous, j'illustre un modèle rapide que j'ai préparé pour les données complètes du formulaire suivant

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

où

• $\alpha$ est l'ordonnée à l'origine du modèle,
• $f_1(\text{ToD}_i)$ est une fonction fluide de l'heure de la journée,
• $f_2(\text{DoY}_i)$ est une fonction régulière du jour de l'année,
• $f_3(\text{Year}_i)$ est une fonction lisse de l'année,
• $f_4(\text{x}_i, \text{y}_i)$ est un lissage 2D de longitude et de latitude,
• $f_5(\text{DoY}_i, \text{Year}_i)$ est un produit tensoriel lisse du jour de l'année et de l'année,
• $f_6(\text{x}_i, \text{y}_i, \text{ToD}_i)$ produit tenseur lisse de l'emplacement et de l'heure de la journée
• $f_7(\text{x}_i, \text{y}_i, \text{DoY}_i)$ produit tensoriel lisse de l'emplacement jour de l'année &
• $f_8(\text{x}_i, \text{y}_i, \text{Year}_i$ produit tenseur lisse d'emplacement et d'année

En effet, les quatre premiers lissages sont les principaux effets de

1. moment de la journée,
2. saison,
3. tendance à long terme,
4. variation spatiale

tandis que les quatre autres produits tensoriels lissent les interactions douces du modèle entre les covariables déclarées, qui modélisent

5. comment le modèle saisonnier de la température varie au fil du temps,
6. comment l’effet de l’heure varie dans l’espace,
7. comment l'effet saisonnier varie spatialement, et
8. comment la tendance à long terme varie spatialement

Les données sont chargées dans R et massées un peu avec le code suivant

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

Le modèle lui-même est ajusté à l'aide de la bam()fonction qui est conçue pour adapter les GAM à de plus grands ensembles de données tels que celui-ci. Vous pouvez également l'utiliser gam()pour ce modèle, mais cela prendra un peu plus de temps pour s'adapter.

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)

Les s()termes sont les effets principaux, tandis que les ti()termes sont des lissages d' interaction des produits tensoriels où les effets principaux des covariables nommées ont été supprimés de la base. Ces ti()lissages sont un moyen d'inclure les interactions des variables déclarées d'une manière numériquement stable.

L' knotsobjet ne fait que définir les points de terminaison du lissage cyclique que j'ai utilisé pour l'effet du jour de l'année - nous voulons que 23:59 le 31 décembre se rejoigne en douceur avec 00:01 le 1er janvier. Cela représente dans une certaine mesure les années bissextiles.

Le résumé du modèle indique que tous ces effets sont importants;

> 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

Une analyse plus approfondie voudrait vérifier si nous avons besoin de toutes ces interactions; certains ti()termes spatiaux n'expliquent que de petites quantités de variation dans les données, comme l'indique$F$statistique; il y a beaucoup de données ici, donc même de petites tailles d'effets peuvent être statistiquement significatives mais sans intérêt.

Cependant, à titre de vérification rapide, la suppression des trois ti()lissages spatiaux ( m.sub) se traduit par un ajustement significativement moins bon, tel qu'évalué par l'AIC:

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

Nous pouvons tracer les effets partiels des cinq premiers lissages en utilisant la plot()méthode - les lissages du produit tensoriel 3D ne peuvent pas être tracés facilement et pas par défaut.

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

L' scale = 0argument ici met toutes les parcelles à leur propre échelle, pour comparer les ampleurs des effets, nous pouvons désactiver cela:

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

On voit maintenant que l'effet saisonnier domine. La tendance à long terme (en moyenne) est indiquée dans le graphique supérieur droit. Cependant, pour vraiment regarder la tendance à long terme, vous devez choisir une station, puis prédire à partir du modèle de cette station, en fixant l'heure du jour et le jour de l'année à certaines valeurs représentatives (midi, pour un jour de l'année en été dire). Là, au début de l'année ou deux de la série, certaines valeurs de températures basses par rapport au reste des enregistrements, sont probablement relevées dans tous les lissages impliqués YEAR. Ces données doivent être examinées de plus près.

Ce n'est pas vraiment l'endroit pour entrer dans cela, mais voici quelques visualisations des ajustements du modèle. Je regarde d'abord le modèle spatial de la température et comment il varie au cours des années de la série. Pour ce faire, je prédis à partir du modèle d'une grille 100x100 sur le domaine spatial, à midi le jour 180 de chaque année:

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)

puis je mets à manquant, NAles valeurs prévues fitpour tous les points de données qui se trouvent à une certaine distance des observations (proportionnel; dist)

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

et joindre les prédictions aux données de prédiction

pred <- cbind(pdata, Fitted = fit)

La définition de valeurs prédites NAcomme celle-ci nous empêche d'extrapoler au-delà du support des données.

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

nous obtenons ce qui suit

Nous pouvons voir la variation des températures d'une année à l'autre avec un peu plus de détails si nous animons plutôt que de facette l'intrigue

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)

Pour examiner plus en détail les tendances à long terme, nous pouvons prévoir des stations particulières. Par exemple, pour STATION_ID13364 et la prévision des jours au cours des quatre trimestres, nous pouvons utiliser les éléments suivants pour préparer les valeurs des covariables que nous voulons prédire à (midi, le jour de l'année 1, 90, 180 et 270, à la station sélectionnée , et évaluer la tendance à long terme à 500 valeurs également espacées)

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

Ensuite, nous prédisons et demandons des erreurs standard, pour former un intervalle de confiance approximativement ponctuel à 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 nous traçons

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

produire

Évidemment, il y a beaucoup plus à modéliser ces données que ce que je montre ici, et nous voudrions vérifier l'autocorrélation résiduelle et le sur-ajustement des splines, mais aborder le problème comme l'un des modèles des caractéristiques des données permet une analyse plus détaillée examen des tendances.

Vous pouvez bien sûr simplement modéliser chacun STATION_IDséparément, mais cela gâcherait les données, et de nombreuses stations ont peu d'observations. Ici, le modèle emprunte à toutes les informations de la station pour combler les lacunes et aider à estimer les tendances d'intérêt.

Quelques notes sur bam()

Le bam()modèle utilise toutes les astuces de mgcv pour estimer rapidement le modèle (plusieurs threads 4 ), la sélection rapide de la fluidité REML ( method = 'fREML') et la discrétisation des covariables. Avec ces options activées, le modèle tient en moins d'une minute sur ma station de travail Xeon double 4 cœurs de l'ère 2013 avec 64 Go de RAM.

Vous pouvez simplement utiliser la fonction decomposequi sépare vos séries chronologiques en trois composantes, tendance, saisonnière et aléatoire. Vous pouvez également extraire les différentes valeurs du résultat et les tracer. Assurez-vous de définir vos données comme une série chronologique. La fonction stlfait essentiellement la même chose mais vous donne plus de flexibilité dans la façon dont vous pouvez décomposer vos données.

Je recommande également le site Web suivant

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

est-ce que cela aide?

Une option serait d'utiliser une stratégie de régression qui traite vos données comme tirées d'une fonction sous-jacente continue (voir la discussion ici: Existe - t-il un étalon-or pour modéliser des séries chronologiques à espacement irrégulier? ). De là, vous pouvez utiliser une méthode comme l'analyse du spectre singulier pour décomposer le signal (pour R: https://cran.r-project.org/web/packages/Rssa/index.html ).