Ressourcen für die Analyse unterbrochener Zeitreihen in R.

12

Ich bin ziemlich neu in R. Ich habe versucht, mich über Zeitreihenanalysen zu informieren und bin bereits fertig

  1. Shumway und Stoffers Zeitreihenanalyse und ihre Anwendungen 3. Auflage ,
  2. Hyndmans exzellente Prognose: Prinzipien und Praxis
  3. Avril Coghlan verwendet R für die Zeitreihenanalyse
  4. A. Ian McLeod et al. Zeitreihenanalyse mit R.
  5. Angewandte Zeitreihenanalyse von Dr. Marcel Dettling

Bearbeiten: Ich bin nicht sicher, wie ich damit umgehen soll, aber ich habe eine nützliche Ressource außerhalb von Cross Validated gefunden. Ich wollte es hier aufnehmen, falls jemand auf diese Frage stößt.

Segmentierte Regressionsanalyse unterbrochener Zeitreihenstudien in der Medikamentengebrauchsforschung

Ich habe eine univariate Zeitreihe der Anzahl der konsumierten Artikel (Zähldaten), die 7 Jahre lang täglich gemessen wurden. Etwa in der Mitte der Zeitreihe wurde eine Intervention auf die Studienpopulation angewendet. Es wird nicht erwartet, dass dieser Eingriff eine sofortige Wirkung hervorruft, und der Zeitpunkt des Wirkungseintritts ist im Wesentlichen nicht bekannt.

Mit Hyndmans forecastPaket habe ich ein ARIMA-Modell an die Daten vor der Intervention angepasst auto.arima(). Ich bin mir jedoch nicht sicher, wie ich diese Anpassung verwenden soll, um zu beantworten, ob sich der Trend statistisch signifikant geändert hat, und um die Menge zu quantifizieren.

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

Gibt es Ressourcen, die sich speziell mit der Analyse unterbrochener Zeitreihen in R befassen? Ich habe diesen Umgang mit ITS in SPSS gefunden, konnte ihn jedoch nicht in R übersetzen.

dais.johns
quelle
Möchten Sie Rückschlüsse darauf ziehen, ob die Intervention einen statistisch signifikanten Effekt hatte, oder möchten Sie die Intervention modellieren, um bessere Prognosen zu erhalten ? Und könnten Sie möglicherweise die Daten zur Verfügung stellen?
Stephan Kolassa
@StephanKolassa Auf jeden Fall! Mein Ziel ist es, Schlussfolgerungen zu ziehen. Ich werde Dummy-Daten in einer Bearbeitung bereitstellen, um meinen Standpunkt besser zu veranschaulichen.
dais.johns
@StephanKolassa Daten nach bestem Wissen und Gewissen.
dais.johns
Frühere Untersuchungen deuten darauf hin, dass der Interventionseffekt in der Größenordnung von +/- 5% liegt.
dais.johns
@StephanKolassa Bereitstellung der tatsächlich verwendbaren Daten
dais.johns

Antworten:

4

Dies wird als Änderungspunktanalyse bezeichnet. Das R-Paket changepointkann dies für Sie tun: Siehe die Dokumentation hier (einschließlich Verweise auf die Literatur): http://www.lancs.ac.uk/~killick/Pub/KillickEckley2011.pdf

Brent Kerby
quelle
Vielen Dank. Ich untersuche das. Soweit ich das beurteilen kann, würde dies mögliche Änderungspunkte in der Serie berechnen, aber die Trenddifferenz nicht analysieren. Ich entschuldige mich, wenn diese Annahme falsch ist. Ich konnte das Paket nur oberflächlich überprüfen.
dais.johns
Nachdem Sie den Änderungspunkt identifiziert haben, können Sie die Daten in zwei Zeitreihen (vor und nach dem Änderungspunkt) aufteilen und die Parameter der beiden Zeitreihen separat schätzen. Noch ein paar Vorschläge: Da Ihre Daten einen starken saisonalen Trend aufweisen, sollte dieser vor der Änderungspunktanalyse entfernt werden. Wenn Sie ein ARIMA-Modell verwenden möchten, sollte die Differenzierung auch vor der Änderungspunktanalyse durchgeführt werden (oder Sie müssen alternativ ein spezielleres Verfahren verwenden).
Brent Kerby
Vielen Dank für Ihre Vorschläge, die ich umsetzen möchte, und als "beantwortet" markieren, wenn dies das Problem löst.
dais.johns
2

Ich würde ein hierarchisches Modell mit wiederholten Messungen vorschlagen. Diese Methode sollte robuste Ergebnisse liefern, da jeder Einzelne als seine / ihre eigene Kontrolle fungiert. Versuchen Sie, diesen Link von der UCLA zu überprüfen .

kblansit
quelle
0

Für einen Bayes'schen Ansatz können Sie mcpein Poisson- oder Binomial-Modell anpassen (da Sie Zählungen aus Zeiträumen mit festem Intervall haben), wobei die Autoregression auf die Residuen (im Protokollbereich) angewendet wird. Vergleichen Sie dann ein Zwei-Segment-Modell mit einem Ein-Segment-Modell unter Verwendung der Kreuzvalidierung.

Bevor wir beginnen, beachten Sie, dass dieses Modell für diesen Datensatz nicht gut passt und die Kreuzvalidierung instabil aussieht. Daher würde ich in Szenarien mit hohen Einsätzen auf Folgendes verzichten, aber es zeigt einen allgemeinen Ansatz:

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

Für den vorliegenden Datensatz ergibt dies

       elpd_diff se_diff
model2    0.0       0.0 
model1 -459.1      64.3 

Das heißt, ein elpd_diff/se_diffVerhältnis von etwa 7 zugunsten des Nullmodells (keine Änderung). Mögliche Verbesserungen sind:

  • Modellierung des periodischen Trends mit sin()oder cos().
  • Hinzufügen vorheriger Informationen über den wahrscheinlichen Ort der Änderung, z prior = list(cp_1 = dnorm(1999.8, 0.5).

Lesen Sie mehr über das Modellieren der Autoregression, das Durchführen von Modellvergleichen und das Festlegen von Prioritäten für die mcpWebsite . Offenlegung: Ich bin der Entwickler von mcp.

Jonas Lindeløv
quelle