Ressources pour l'analyse des séries chronologiques interrompues dans R

12

Je suis assez nouveau pour R. J'ai tenté de lire l'analyse des séries chronologiques et j'ai déjà terminé

  1. Analyse des séries temporelles de Shumway et Stoffer et ses applications 3e édition ,
  2. L'excellente prévision de Hyndman : principes et pratique
  3. Avril Coghlan utilise R pour l'analyse des séries chronologiques
  4. A. Ian McLeod et al Analyse des séries chronologiques avec R
  5. Analyse des séries chronologiques appliquées du Dr Marcel Dettling

Edit: je ne sais pas comment gérer cela, mais j'ai trouvé une ressource utile en dehors de Cross Validated. Je voulais l'inclure ici au cas où quelqu'un tomberait sur cette question.

Analyse de régression segmentée des études de séries chronologiques interrompues dans la recherche sur l'utilisation des médicaments

J'ai une série chronologique univariée du nombre d'articles consommés (données de comptage) mesurés quotidiennement pendant 7 ans. Une intervention a été appliquée à la population étudiée à peu près au milieu de la série chronologique. Cette intervention ne devrait pas produire d'effet immédiat et le moment du début de l'effet est essentiellement inconnu.

En utilisant le forecastpackage de Hyndman, j'ai adapté un modèle ARIMA aux données de pré-intervention à l'aide auto.arima(). Mais je ne sais pas comment utiliser cet ajustement pour déterminer s'il y a eu un changement statistiquement significatif de tendance et quantifier le montant.

# for simplification I will aggregate to monthly counts
# I can later generalize any teachings the community supplies
count <- c(2464, 2683, 2426, 2258, 1950, 1548, 1108,  991, 1616, 1809, 1688, 2168, 2226, 2379, 2211, 1925, 1998, 1740, 1305,  924, 1487, 1792, 1485, 1701, 1962, 2896, 2862, 2051, 1776, 1358, 1110,  939, 1446, 1550, 1809, 2370, 2401, 2641, 2301, 1902, 2056, 1798, 1198,  994, 1507, 1604, 1761, 2080, 2069, 2279, 2290, 1758, 1850, 1598, 1032,  916, 1428, 1708, 2067, 2626, 2194, 2046, 1905, 1712, 1672, 1473, 1052,  874, 1358, 1694, 1875, 2220, 2141, 2129, 1920, 1595, 1445, 1308, 1039,  828, 1724, 2045, 1715, 1840)
# for explanatory purposes
# month <- rep(month.name, 7)
# year <- 1999:2005
ts <- ts(count, start(1999, 1))
train_month <- window(ts, start=c(1999,1), end = c(2001,1))
require(forecast)
arima_train <- auto.arima(train_month)
fit_month <- Arima(train_month, order = c(2,0,0), seasonal = c(1,1,0), lambda = 0)
plot(forecast(fit_month, 36)); lines(ts, col="red")

Existe-t-il des ressources traitant spécifiquement de l'analyse des séries chronologiques interrompues dans R? J'ai trouvé ce problème avec les STI dans SPSS mais je n'ai pas pu le traduire en R.

dais.johns
la source
Voulez-vous déduire si l'intervention a eu un effet statistiquement significatif, ou voulez-vous modéliser l'intervention pour obtenir de meilleures prévisions ? Et pourriez-vous éventuellement rendre les données disponibles?
Stephan Kolassa
@StephanKolassa Certainement! Mon objectif est de faire de l'inférence. Je fournirai des données factices dans une édition pour mieux illustrer mon propos.
dais.johns
@StephanKolassa Données fournies au mieux de mes capacités.
dais.johns
Des recherches antérieures suggèrent que l'effet de l'intervention soit à l'échelle d'un changement de +/- 5%.
dais.johns
@StephanKolassa Fourni des données utilisables réelles
dais.johns

Réponses:

4

C'est ce qu'on appelle l'analyse des points de changement. Le package R changepointpeut le faire pour vous: voir la documentation ici (y compris les références à la littérature): http://www.lancs.ac.uk/~killick/Pub/KillickEckley2011.pdf

Brent Kerby
la source
Je vous remercie. Je regarde ça. Autant que je sache, cela calculerait les points de changement possibles dans la série, mais n'analyserait pas la différence de tendance. Je m'excuse si cette hypothèse est incorrecte, je n'ai pas pu revoir le paquet autrement que superficiellement.
dais.johns
Après avoir identifié le point de changement, vous pouvez diviser les données en deux séries temporelles (avant et après le point de changement) et estimer les paramètres des deux séries temporelles séparément. Quelques suggestions supplémentaires: comme vos données ont une forte tendance saisonnière, cela devrait être supprimé avant l'analyse du point de changement; et si vous envisagez d'utiliser un modèle ARIMA, la différenciation doit également être effectuée avant l'analyse du point de changement (ou, alternativement, vous devrez utiliser une procédure plus spécialisée).
Brent Kerby
Merci pour vos suggestions que j'essaierai de mettre en œuvre et marquerai comme "répondues" si cela résout le problème.
dais.johns
2

Je proposerais un modèle hiérarchique de mesures répétées. Cette méthode devrait fournir des résultats solides car chaque individu agira comme son propre contrôle. Essayez de vérifier ce lien depuis UCLA.

kblansit
la source
0

Pour une approche bayésienne, vous pouvez utiliser mcppour ajuster un modèle de Poisson ou binomial (parce que vous avez des comptes à partir de périodes à intervalle fixe) avec une autorégression appliquée aux résidus (dans l'espace logarithmique). Comparez ensuite un modèle à deux segments à un modèle à un segment à l'aide de la validation croisée.

Avant de commencer, notez que pour ce jeu de données, ce modèle ne convient pas bien et la validation croisée semble instable. Je m'abstiendrai donc d'utiliser les éléments suivants dans des scénarios à enjeux élevés, mais cela illustre une approche générale:

# Fit the change point model
library(mcp)
model_full = list(
  count ~ 1 + ar(1),  # intercept and AR(1)
  ~ 1  # New intercept
)
fit_full = mcp(model_full, data = df, family = poisson(), par_x = "year")


# Fit the null model
model_null = list(
  count ~ 1 + ar(1)  # just a stable AR(1)
)
fit_null = mcp(model_null, data = df, family = poisson(), par_x = "year")

# Compare predictive performance using LOO cross-validation
fit_full$loo = loo(fit_full)
fit_null$loo = loo(fit_null)
loo::loo_compare(fit_full$loo, fit_null$loo)

Pour le présent ensemble de données, cela se traduit par

       elpd_diff se_diff
model2    0.0       0.0 
model1 -459.1      64.3 

C'est-à-dire un elpd_diff/se_diffrapport d'environ 7 en faveur du modèle nul (pas de changement). Les améliorations possibles incluent:

  • modéliser la tendance périodique en utilisant sin()ou cos().
  • l'ajout d'informations préalables sur l'emplacement probable du changement, par exemple prior = list(cp_1 = dnorm(1999.8, 0.5).

En savoir plus sur la modélisation de l'autorégression, la comparaison de modèles et la définition des priorités du mcpsite Web . Divulgation: je suis le développeur de mcp.

Jonas Lindeløv
la source