Zasoby do analizy przerwanych szeregów czasowych w R.


12

Jestem dość nowy w R. Próbowałem przeczytać analizę szeregów czasowych i już skończyłem

  1. Analiza szeregów czasowych Shumwaya i Stoffera i jej zastosowania 3. edycja ,
  2. Doskonałe prognozy Hyndmana : zasady i praktyka
  3. Avril Coghlan używa R do analizy szeregów czasowych
  4. A. Ian McLeod i in. Analiza szeregów czasowych z R.
  5. Analiza zastosowanych szeregów czasowych dr Marcela Dettlinga

Edycja: Nie jestem pewien, jak sobie z tym poradzić, ale znalazłem przydatny zasób poza Cross Validated. Chciałem to tutaj zamieścić, na wypadek, gdyby ktoś natknął się na to pytanie.

Analiza regresji segmentowej przerwanych badań szeregów czasowych w badaniach nad stosowaniem leków

Mam jednoznaczny szereg czasowy liczby zużytych produktów (dane zliczeń) mierzonych codziennie przez 7 lat. Interwencję zastosowano w badanej populacji mniej więcej w połowie szeregów czasowych. Nie należy oczekiwać, że ta interwencja przyniesie natychmiastowy efekt, a czas pojawienia się efektu jest zasadniczo niepoznawalny.

Korzystając z forecastpakietu Hyndmana, dopasowałem model ARIMA do danych sprzed interwencji przy użyciu auto.arima(). Ale nie jestem pewien, jak wykorzystać to dopasowanie, aby odpowiedzieć, czy nastąpiła statystycznie istotna zmiana trendu i określić kwotę.

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

Czy są jakieś zasoby dotyczące konkretnie analizy przerwanych szeregów czasowych w R? Znalazłem to zajmowanie się ITS w SPSS, ale nie byłem w stanie przetłumaczyć tego na R.


Czy chcesz wnioskować, czy interwencja miała statystycznie znaczący wpływ, czy też chcesz modelować interwencję, aby uzyskać lepsze prognozy ? Czy możesz udostępnić dane?
Stephan Kolassa,

@StephanKolassa Na pewno! Moim celem jest wnioskowanie. Dostarczę pozorowane dane w Edycji, aby lepiej zilustrować mój punkt widzenia.
dais.johns,

@StephanKolassa Dane dostarczone zgodnie z moimi najlepszymi umiejętnościami.
dais.johns,

Wcześniejsze badania sugerują, że wpływ interwencji jest na skalę +/- 5% zmiany.
dais.johns,

@StephanKolassa Podane rzeczywiste dane użyteczne
dais.johns

Odpowiedzi:


4

Jest to znane jako analiza punktu zmiany. Pakiet R changepointmoże to dla Ciebie zrobić: zapoznaj się z dokumentacją tutaj (w tym odniesienia do literatury): http://www.lancs.ac.uk/~killick/Pub/KillickEckley2011.pdf


Dziękuję Ci. Patrzę na to. O ile mogę powiedzieć, obliczy to możliwe punkty zmiany w szeregu, ale nie przeanalizuje różnicy trendów. Przepraszam, jeśli to założenie jest nieprawidłowe, nie byłem w stanie przejrzeć pakietu inaczej niż powierzchownie.
dais.johns,

Po zidentyfikowaniu punktu zmiany można podzielić dane na dwa szeregi czasowe (przed i po punkcie zmiany) i osobno oszacować parametry dwóch szeregów czasowych. Kilka dodatkowych sugestii: ponieważ Twoje dane wykazują silną tendencję sezonową, należy to usunąć przed analizą punktu zmiany; a jeśli zamierzasz użyć modelu ARIMA, różnicowanie należy również wykonać przed analizą punktu zmiany (lub alternatywnie będziesz musiał zastosować bardziej wyspecjalizowaną procedurę).
Brent Kerby,

Dziękuję za sugestie, które spróbuję wdrożyć i oznaczę jako „odpowiedz”, jeśli to rozwiąże problem.
dais.johns,

2

Sugerowałbym model hierarchiczny z powtarzanymi pomiarami. Ta metoda powinna zapewnić solidne wyniki, ponieważ każda osoba będzie działać jako własna kontrola. Spróbuj sprawdzić ten link w UCLA.


0

W przypadku podejścia bayesowskiego można użyć, mcpaby dopasować model Poissona lub dwumianowy (ponieważ masz liczby z okresów o ustalonym interwale) z autoregresją zastosowaną do reszt (w przestrzeni dziennika). Następnie porównaj model dwusegmentowy z modelem jednosegmentowym za pomocą weryfikacji krzyżowej.

Zanim zaczniemy, zauważ, że dla tego zestawu danych ten model nie pasuje dobrze, a sprawdzanie poprawności wygląda niestabilnie. Powstrzymałbym się więc od stosowania następujących scenariuszy w scenariuszach z dużymi stawkami, ale ilustruje to ogólne podejście:

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

W przypadku obecnego zestawu danych powoduje to:

       elpd_diff se_diff
model2    0.0       0.0 
model1 -459.1      64.3 

Tj. elpd_diff/se_diffStosunek około 7 na korzyść modelu zerowego (bez zmian). Możliwe ulepszenia obejmują:

  • modelowanie trendu okresowego za pomocą sin()lub cos().
  • dodanie wcześniejszych informacji o prawdopodobnej lokalizacji zmiany, np prior = list(cp_1 = dnorm(1999.8, 0.5).

Czytaj więcej na temat modelowania autoregresji, robiąc porównanie modelu i ustawienie priors na mcpstronie internetowej . Ujawnienie: Jestem deweloperem mcp.

Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.