Czy jest to właściwa metoda testowania skutków sezonowych w danych dotyczących liczby samobójstw?

Mam 17 lat (1995–2011) danych dotyczących aktu zgonu związanych ze śmiercią samobójczą dla stanu w USA. Istnieje wiele mitologii na temat samobójstw i miesięcy / pór roku, wiele z nich jest sprzecznych, a literatura I ” Po przejrzeniu recenzji nie rozumiem zastosowanych metod ani nie ufam wynikom.

Dlatego postanowiłem sprawdzić, czy mogę ustalić, czy samobójstwa są bardziej lub mniej prawdopodobne w danym miesiącu w ramach mojego zbioru danych. Wszystkie moje analizy są wykonywane w języku R.

Całkowita liczba samobójstw w danych wynosi 13 909.

Jeśli spojrzysz na rok z najmniejszą liczbą samobójstw, zdarzają się one w ciągu 309/365 dni (85%). Jeśli spojrzysz na rok, w którym popełniono najwięcej samobójstw, zdarzają się one w 339/365 dni (93%).

Każdego roku jest wiele dni bez samobójstw. Jednak po zsumowaniu przez wszystkie 17 lat, dochodzi do samobójstw każdego dnia roku, w tym 29 lutego (chociaż tylko 5, gdy średnia wynosi 38).

Samo zsumowanie liczby samobójstw każdego dnia w roku nie oznacza wyraźnej sezonowości (moim zdaniem).

Zagregowane na poziomie miesięcznym, średnie samobójstwa na miesiąc wahają się od:

(m = 65, sd = 7,4, do m = 72, sd = 11,1)

Moje pierwsze podejście polegało na agregacji zestawu danych według miesięcy dla wszystkich lat i wykonaniu testu chi-kwadrat po obliczeniu oczekiwanych prawdopodobieństw dla hipotezy zerowej, że nie ma systematycznej wariancji w liczbie samobójstw według miesięcy. Obliczałem prawdopodobieństwa dla każdego miesiąca, biorąc pod uwagę liczbę dni (i dostosowując luty dla lat przestępnych).

Wyniki chi-kwadrat nie wykazały istotnych różnic w zależności od miesiąca:

# So does the sample match  expected values?
chisq.test(monthDat$suicideCounts, p=monthlyProb)
# Yes, X-squared = 12.7048, df = 11, p-value = 0.3131

Poniższy obraz pokazuje łączną liczbę miesięcznie. Poziome czerwone linie są ustawione zgodnie z oczekiwanymi wartościami odpowiednio dla lutego, 30 dni i 31 dni. Zgodnie z testem chi-kwadrat, żaden miesiąc nie jest powyżej 95% przedziału ufności dla oczekiwanych zliczeń.

Myślałem, że skończyłem, dopóki nie zacząłem badać danych szeregów czasowych. Jak sobie wyobrażam wiele osób, zacząłem od nieparametrycznej metody sezonowego rozkładu przy użyciu stlfunkcji w pakiecie statystyk.

Aby utworzyć dane szeregów czasowych, zacząłem od zagregowanych danych miesięcznych:

suicideByMonthTs <- ts(suicideByMonth$monthlySuicideCount, start=c(1995, 1), end=c(2011, 12), frequency=12) 

# Plot the monthly suicide count, note the trend, but seasonality?
plot(suicideByMonthTs, xlab="Year",
  ylab="Annual  monthly  suicides")

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1995  62  47  55  74  71  70  67  69  61  76  68  68
1996  64  69  68  53  72  73  62  63  64  72  55  61
1997  71  61  64  63  60  64  67  50  48  49  59  72
1998  67  54  72  69  78  45  59  53  48  65  64  44
1999  69  64  65  58  73  83  70  73  58  75  71  58
2000  60  54  67  59  54  69  62  60  58  61  68  56
2001  67  60  54  57  51  61  67  63  55  70  54  55
2002  65  68  65  72  79  72  64  70  59  66  63  66
2003  69  50  59  67  73  77  64  66  71  68  59  69
2004  68  61  66  62  69  84  73  62  71  64  59  70
2005  67  53  76  65  77  68  65  60  68  71  60  79
2006  65  54  65  68  69  68  81  64  69  71  67  67
2007  77  63  61  78  73  69  92  68  72  61  65  77
2008  67  73  81  73  66  63  96  71  75  74  81  63
2009  80  68  76  65  82  69  74  88  80  86  78  76
2010  80  77  82  80  77  70  81  89  91  82  71  73
2011  93  64  87  75 101  89  87  78 106  84  64  71

A następnie przeprowadził stl()rozkład

# Seasonal decomposition
suicideByMonthFit <- stl(suicideByMonthTs, s.window="periodic")
plot(suicideByMonthFit)

W tym momencie zaniepokoiłem się, ponieważ wydaje mi się, że istnieje zarówno składnik sezonowy, jak i trend. Po wielu badaniach w Internecie postanowiłem postępować zgodnie z instrukcjami Roba Hyndmana i George'a Athanasopoulosa zgodnie z ich tekstem online „Prognozowanie: zasady i praktyka”, w szczególności w celu zastosowania sezonowego modelu ARIMA.

Używałem adf.test()i kpss.test()oceniałem stacjonarność i otrzymałem sprzeczne wyniki. Obaj odrzucili hipotezę zerową (zauważając, że testują hipotezę przeciwną).

adfResults <- adf.test(suicideByMonthTs, alternative = "stationary") # The p < .05 value 
adfResults

    Augmented Dickey-Fuller Test

data:  suicideByMonthTs
Dickey-Fuller = -4.5033, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary

kpssResults <- kpss.test(suicideByMonthTs)
kpssResults

    KPSS Test for Level Stationarity

data:  suicideByMonthTs
KPSS Level = 2.9954, Truncation lag parameter = 3, p-value = 0.01

Następnie użyłem algorytmu z książki, aby sprawdzić, czy potrafię określić stopień różnicowania, który należy wykonać zarówno dla trendu, jak i sezonu. Skończyłem z nd = 1, ns = 0.

Następnie pobiegłem auto.arima, który wybrał model, który miał zarówno trend, jak i składnik sezonowy wraz ze stałą typu „dryfowania”.

# Extract the best model, it takes time as I've turned off the shortcuts (results differ with it on)
bestFit <- auto.arima(suicideByMonthTs, stepwise=FALSE, approximation=FALSE)
plot(theForecast <- forecast(bestFit, h=12))
theForecast

> summary(bestFit)
Series: suicideByMonthFromMonthTs 
ARIMA(0,1,1)(1,0,1)[12] with drift         

Coefficients:
          ma1    sar1     sma1   drift
      -0.9299  0.8930  -0.7728  0.0921
s.e.   0.0278  0.1123   0.1621  0.0700

sigma^2 estimated as 64.95:  log likelihood=-709.55
AIC=1429.1   AICc=1429.4   BIC=1445.67

Training set error measures:
                    ME    RMSE     MAE       MPE     MAPE     MASE       ACF1
Training set 0.2753657 8.01942 6.32144 -1.045278 9.512259 0.707026 0.03813434

Wreszcie spojrzałem na pozostałości po dopasowaniu i jeśli dobrze to rozumiem, ponieważ wszystkie wartości mieszczą się w granicach progu, zachowują się jak biały szum, a zatem model jest dość rozsądny. Przeprowadziłem test Portmanteau, jak opisano w tekście, który miał wartość ap znacznie powyżej 0,05, ale nie jestem pewien, czy mam prawidłowe parametry.

Acf(residuals(bestFit))

Box.test(residuals(bestFit), lag=12, fitdf=4, type="Ljung")

    Box-Ljung test

data:  residuals(bestFit)
X-squared = 7.5201, df = 8, p-value = 0.4817

Po powrocie i ponownym przeczytaniu rozdziału o modelowaniu arima, zdaję sobie sprawę, że auto.arimazdecydowałem się modelować trend i sezon. I zdaję sobie również sprawę, że prognozowanie nie jest specjalnie analizą, którą prawdopodobnie powinienem zrobić. Chcę wiedzieć, czy określony miesiąc (lub bardziej ogólnie pora roku) powinien zostać oznaczony jako miesiąc wysokiego ryzyka. Wydaje się, że narzędzia w literaturze prognostycznej są bardzo trafne, ale może nie najlepsze w przypadku mojego pytania. Wszelkie uwagi są mile widziane.

Publikuję link do pliku csv, który zawiera dzienne liczby. Plik wygląda następująco:

head(suicideByDay)

        date year month day_of_month t count
1 1995-01-01 1995    01           01 1     2
2 1995-01-03 1995    01           03 2     1
3 1995-01-04 1995    01           04 3     3
4 1995-01-05 1995    01           05 4     2
5 1995-01-06 1995    01           06 5     3
6 1995-01-07 1995    01           07 6     2

Daily_suicide_data.csv

Liczba to liczba samobójstw, które miały miejsce tego dnia. „t” to sekwencja liczbowa od 1 do całkowitej liczby dni w tabeli (5533).

Zanotowałem poniższe komentarze i pomyślałem o dwóch rzeczach związanych z modelowaniem samobójstw i pór roku. Po pierwsze, jeśli chodzi o moje pytanie, miesiące są po prostu pełnomocnikami do oznaczania zmiany sezonu, nie jestem zainteresowany, ani żaden konkretny miesiąc nie różni się od innych (to oczywiście interesujące pytanie, ale to nie to, co postanowiłem zbadać). Dlatego myślę, że sensowne jest wyrównanie miesięcy poprzez użycie pierwszych 28 dni wszystkich miesięcy. Kiedy to robisz, stajesz się nieco gorszy, co interpretuję jako więcej dowodów na brak sezonowości. W wynikach poniżej pierwsze dopasowanie jest reprodukcją odpowiedzi poniżej z wykorzystaniem miesięcy z ich prawdziwą liczbą dni, a następnie zestawem danych samobójstwoByShortMonthw których liczba samobójstw została obliczona na podstawie pierwszych 28 dni wszystkich miesięcy. Interesuje mnie, co ludzie sądzą o tym, czy ta zmiana jest dobrym pomysłem, czy nie jest konieczna czy szkodliwa?

> summary(seasonFit)

Call:
glm(formula = count ~ t + days_in_month + cos(2 * pi * t/12) + 
    sin(2 * pi * t/12), family = "poisson", data = suicideByMonth)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4782  -0.7095  -0.0544   0.6471   3.2236  

Coefficients:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)         2.8662459  0.3382020   8.475  < 2e-16 ***
t                   0.0013711  0.0001444   9.493  < 2e-16 ***
days_in_month       0.0397990  0.0110877   3.589 0.000331 ***
cos(2 * pi * t/12) -0.0299170  0.0120295  -2.487 0.012884 *  
sin(2 * pi * t/12)  0.0026999  0.0123930   0.218 0.827541    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 302.67  on 203  degrees of freedom
Residual deviance: 190.37  on 199  degrees of freedom
AIC: 1434.9

Number of Fisher Scoring iterations: 4

> summary(shortSeasonFit)

Call:
glm(formula = shortMonthCount ~ t + cos(2 * pi * t/12) + sin(2 * 
    pi * t/12), family = "poisson", data = suicideByShortMonth)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.2414  -0.7588  -0.0710   0.7170   3.3074  

Coefficients:
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)         4.0022084  0.0182211 219.647   <2e-16 ***
t                   0.0013738  0.0001501   9.153   <2e-16 ***
cos(2 * pi * t/12) -0.0281767  0.0124693  -2.260   0.0238 *  
sin(2 * pi * t/12)  0.0143912  0.0124712   1.154   0.2485    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 295.41  on 203  degrees of freedom
Residual deviance: 205.30  on 200  degrees of freedom
AIC: 1432

Number of Fisher Scoring iterations: 4

Drugą rzeczą, nad którą bardziej się zastanawiałem, jest kwestia używania miesiąca jako zastępcy sezonu. Być może lepszym wskaźnikiem sezonu jest liczba godzin dziennych w danym obszarze. Te dane pochodzą ze stanu północnego, który ma znaczne różnice w świetle dziennym. Poniżej znajduje się wykres światła dziennego z roku 2002.

Kiedy używam tych danych, a nie miesiąca, efekt jest nadal znaczący, ale efekt jest bardzo, bardzo mały. Odchylenie resztkowe jest znacznie większe niż w powyższych modelach. Jeśli godziny dzienne są lepszym modelem dla pór roku, a dopasowanie nie jest tak dobre, czy jest to więcej dowodów na bardzo mały efekt sezonowy?

> summary(daylightFit)

Call:
glm(formula = aggregatedDailyCount ~ t + daylightMinutes, family = "poisson", 
    data = aggregatedDailyNoLeap)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0003  -0.6684  -0.0407   0.5930   3.8269  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      3.545e+00  4.759e-02  74.493   <2e-16 ***
t               -5.230e-05  8.216e-05  -0.637   0.5244    
daylightMinutes  1.418e-04  5.720e-05   2.479   0.0132 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 380.22  on 364  degrees of freedom
Residual deviance: 373.01  on 362  degrees of freedom
AIC: 2375

Number of Fisher Scoring iterations: 4

Wysyłam godziny dzienne na wypadek, gdyby ktoś chciał się tym pobawić. Uwaga: nie jest to rok przestępny, więc jeśli chcesz podać minuty dla lat przestępnych, ekstrapoluj lub odzyskaj dane.

state.daylight.2002.csv

[ Edytuj, aby dodać fabułę z usuniętej odpowiedzi (mam nadzieję, że rnso też nie przeszkadza mi przeniesienie fabuły w usuniętej odpowiedzi tutaj na pytanie. Svannoy, jeśli mimo wszystko nie chcesz dodawać tej odpowiedzi, możesz ją cofnąć)]

Co z regresją Poissona?

Stworzyłem ramkę danych zawierającą twoje dane, a także indeks tczasu (w miesiącach) i zmienną monthdaysdla liczby dni w każdym miesiącu.

T <- read.table("suicide.txt", header=TRUE)
U <- data.frame( year = as.numeric(rep(rownames(T),each=12)), 
         month = rep(colnames(T),nrow(T)), 
         t = seq(0, length = nrow(T)*ncol(T)), 
         suicides = as.vector(t(T)))
U$monthdays <- c(31,28,31,30,31,30,31,31,30,31,30,31)
U$monthdays[ !(U$year %% 4) & U$month == "Feb" ] <- 29

Wygląda to tak:

> head(U,14)
   year month  t suicides monthdays
1  1995   Jan  0       62        31
2  1995   Feb  1       47        28
3  1995   Mar  2       55        31
4  1995   Apr  3       74        30
5  1995   May  4       71        31
6  1995   Jun  5       70        30
7  1995   Jul  6       67        31
8  1995   Aug  7       69        31
9  1995   Sep  8       61        30
10 1995   Oct  9       76        31
11 1995   Nov 10       68        30
12 1995   Dec 11       68        31
13 1996   Jan 12       64        31
14 1996   Feb 13       69        29

Porównajmy teraz model z efektem czasowym i efektem liczby dni z modelem, w którym dodajemy efekt miesiąca:

> a0 <- glm( suicides ~ t + monthdays, family="poisson", data = U )
> a1 <- glm( suicides ~ t + monthdays + month, family="poisson", data = U )

Oto podsumowanie „małego” modelu:

> summary(a0)

Call:
glm(formula = suicides ~ t + monthdays, family = "poisson", data = U)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-2.7163  -0.6865  -0.1161   0.6363   3.2104

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.8060135  0.3259116   8.610  < 2e-16 ***
t           0.0013650  0.0001443   9.461  < 2e-16 ***
monthdays   0.0418509  0.0106874   3.916 9.01e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 302.67  on 203  degrees of freedom
Residual deviance: 196.64  on 201  degrees of freedom
AIC: 1437.2

Number of Fisher Scoring iterations: 4

Widać, że te dwie zmienne mają w znacznym stopniu znaczący efekt krańcowy. Teraz spójrz na większy model:

> summary(a1)

Call:
glm(formula = suicides ~ t + monthdays + month, family = "poisson",
    data = U)

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-2.56164  -0.72363  -0.05581   0.58897   3.09423

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)  1.4559200  2.1586699   0.674    0.500
t            0.0013810  0.0001446   9.550   <2e-16 ***
monthdays    0.0869293  0.0719304   1.209    0.227
monthAug    -0.0845759  0.0832327  -1.016    0.310
monthDec    -0.1094669  0.0833577  -1.313    0.189
monthFeb     0.0657800  0.1331944   0.494    0.621
monthJan    -0.0372652  0.0830087  -0.449    0.653
monthJul    -0.0125179  0.0828694  -0.151    0.880
monthJun     0.0452746  0.0414287   1.093    0.274
monthMar    -0.0638177  0.0831378  -0.768    0.443
monthMay    -0.0146418  0.0828840  -0.177    0.860
monthNov    -0.0381897  0.0422365  -0.904    0.366
monthOct    -0.0463416  0.0830329  -0.558    0.577
monthSep     0.0070567  0.0417829   0.169    0.866
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 302.67  on 203  degrees of freedom
Residual deviance: 182.72  on 190  degrees of freedom
AIC: 1445.3

Number of Fisher Scoring iterations: 4

Oczywiście monthdaysefekt zanika; można to oszacować tylko dzięki latom przestępnym !! Utrzymanie go w modelu (i uwzględnienie lat przestępnych) pozwala na wykorzystanie pozostałych odchyleń do porównania dwóch modeli.

> anova(a0, a1, test="Chisq")
Analysis of Deviance Table

Model 1: suicides ~ t + monthdays
Model 2: suicides ~ t + monthdays + month
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       201     196.65
2       190     182.72 11   13.928    0.237

$\cos\left( {2\pi t \over 12}\right)$$\sin\left( {2\pi t \over 12}\right)$

> a2 <- glm( suicides ~ t + monthdays + cos(2*pi*t/12) + sin(2*pi*t/12),
             family="poisson", data = U )
> summary(a2)

Call:
glm(formula = suicides ~ t + monthdays + cos(2 * pi * t/12) +
    sin(2 * pi * t/12), family = "poisson", data = U)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-2.4782  -0.7095  -0.0544   0.6471   3.2236

Coefficients:
                     Estimate Std. Error z value Pr(>|z|)
(Intercept)         2.8676170  0.3381954   8.479  < 2e-16 ***
t                   0.0013711  0.0001444   9.493  < 2e-16 ***
monthdays           0.0397990  0.0110877   3.589 0.000331 ***
cos(2 * pi * t/12) -0.0245589  0.0122658  -2.002 0.045261 *
sin(2 * pi * t/12)  0.0172967  0.0121591   1.423 0.154874
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 302.67  on 203  degrees of freedom
Residual deviance: 190.37  on 199  degrees of freedom
AIC: 1434.9

Number of Fisher Scoring iterations: 4

Teraz porównaj to z modelem zerowym:

> anova(a0, a2, test="Chisq")
Analysis of Deviance Table

Model 1: suicides ~ t + monthdays
Model 2: suicides ~ t + monthdays + cos(2 * pi * t/12) + sin(2 * pi *
    t/12)
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       201     196.65                   
2       199     190.38  2   6.2698   0.0435 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Można więc z pewnością powiedzieć, że sugeruje to efekt sezonowy ...

Test chi-kwadrat jest dobrym podejściem jako wstępny pogląd na twoje pytanie.

stlRozkładu może być mylące, jako narzędzie do badania obecności sezonowości. Ta procedura umożliwia przywrócenie stabilnego wzorca sezonowego, nawet jeśli jako sygnał wejściowy zostanie przekazany biały szum (losowy sygnał bez struktury). Spróbuj na przykład:

plot(stl(ts(rnorm(144), frequency=12), s.window="periodic"))

Patrzenie na zamówienia wybrane przez automatyczną procedurę wyboru modelu ARIMA jest również nieco ryzykowne, ponieważ sezonowy model ARIMA nie zawsze wiąże się z sezonowością (szczegółowe informacje można znaleźć w tej dyskusji ). W takim przypadku wybrany model generuje cykle sezonowe, a komentarz @RichardHardy jest uzasadniony, jednak potrzebny jest dalszy wgląd, aby stwierdzić, że samobójstwa są motywowane sezonowym wzorcem.

Poniżej podsumowuję niektóre wyniki na podstawie analizy opublikowanej serii miesięcznej. Jest to komponent sezonowy oszacowany na podstawie podstawowego strukturalnego modelu szeregów czasowych:

require(stsm)
m <- stsm.model(model = "llm+seas", y = x)
fit <- stsmFit(m, stsm.method = "maxlik.td.scoring")
plot(tsSmooth(fit)$states[,2], ylab = "")
mtext(text = "seasonal component", side = 3, adj = 0, font = 2)

Podobny komponent został wyodrębniony przy użyciu oprogramowania TRAMO-SEATS z domyślnymi opcjami. Widzimy, że szacowany wzorzec sezonowy nie jest stabilny w czasie, a zatem nie potwierdza hipotezy o powtarzającym się wzorcu liczby samobójstw w miesiącach w okresie próby. Uruchamiając oprogramowanie X-13ARIMA-SEATS z domyślnymi opcjami, połączony test sezonowości stwierdził, że możliwa do zidentyfikowania sezonowość nie występuje.

Edytuj (zobacz tę odpowiedź i mój komentarz poniżej, aby znaleźć możliwą do zidentyfikowania sezonowość ).

Biorąc pod uwagę charakter twoich danych, warto uzupełnić tę analizę opartą na metodach szeregów czasowych o model danych zliczania (np. Model Poissona) i przetestować istotność sezonowości w tym modelu. Dodanie danych licznika tagów do pytania może spowodować zwiększenie liczby wyświetleń i potencjalnych odpowiedzi w tym kierunku.

Jak zauważono w moim komentarzu, jest to bardzo interesujący problem. Wykrywanie sezonowości nie jest zadaniem statystycznym. Rozsądnym podejściem byłoby skonsultowanie teorii i ekspertów, takich jak:

• Psycholog
• Psychiatra
• Socjolog

na temat tego problemu, aby zrozumieć „dlaczego” istniałaby sezonowość w celu uzupełnienia analizy danych. Przechodząc do danych, zastosowałem doskonałą metodę dekompozycji zwaną modelem nieobserwowanych składników (UCM), która jest formą metody przestrzeni stanów. Zobacz także ten bardzo dostępny artykuł Koopmana. Moje podejście jest podobne do @Javlacalle. Zapewnia nie tylko narzędzie do dekompozycji danych szeregów czasowych, ale także obiektywnie ocenia obecność lub brak sezonowości za pomocą testów istotności. Nie jestem wielkim fanem testowania istotności danych nie eksperymentalnych, ale nie znam żadnej innej procedury, w której można by przetestować swoją hipotezę dotyczącą obecności / braku sezonowości na danych szeregów czasowych.

Wielu ignoruje, ale bardzo ważną cechą, którą chcielibyśmy zrozumieć, jest rodzaj sezonowości:

1. Stochastyczny - zmienia się losowo i jest trudny do przewidzenia
2. Deterministyczny - nie zmienia się, doskonale przewidywalny. Możesz modelować za pomocą manekina lub trygonometrii (sin / cos itp.)

W przypadku danych z długich szeregów czasowych, takich jak twoje, możliwe jest, że sezonowość mogła się zmieniać z czasem. Znów UCM jest jedynym znanym mi podejściem, które może wykryć te stochastyczne / deterministyczne zjawiska sezonowości. UCM może rozłożyć problem na następujące „komponenty”:

Time Series Data = level + Slope + Seasonality + Cycle + Causal + Error(Noise).

Możesz również sprawdzić, czy poziom, nachylenie, cykl jest deterministyczny czy stochastyczny. Uwaga level + slope = trend! Poniżej przedstawiam analizę twoich danych za pomocą UCM. Do analizy wykorzystałem SAS.

data input;
format date mmddyy10.;
date = intnx( 'month', '1jan1995'd, _n_-1 );
      input deaths@@;
datalines;
62    47  55  74  71  70  67  69  61  76  68  68
64    69  68  53  72  73  62  63  64  72  55  61
71    61  64  63  60  64  67  50  48  49  59  72
67    54  72  69  78  45  59  53  48  65  64  44
69    64  65  58  73  83  70  73  58  75  71  58
60    54  67  59  54  69  62  60  58  61  68  56
67    60  54  57  51  61  67  63  55  70  54  55
65    68  65  72  79  72  64  70  59  66  63  66
69    50  59  67  73  77  64  66  71  68  59  69
68    61  66  62  69  84  73  62  71  64  59  70
67    53  76  65  77  68  65  60  68  71  60  79
65    54  65  68  69  68  81  64  69  71  67  67
77    63  61  78  73  69  92  68  72  61  65  77
67    73  81  73  66  63  96  71  75  74  81  63
80    68  76  65  82  69  74  88  80  86  78  76
80    77  82  80  77  70  81  89  91  82  71  73
93    64  87  75  101 89  87  78  106 84  64  71
;
run;

ods graphics on;
 proc ucm data = input plots=all; 
      id date interval = month; 
      model deaths ; 
      irregular ; 
      level checkbreak; 
      season length = 12 type=trig var = 0 noest; * Note I have used trigonometry to model seasonality;
   run;

   ods graphics off;

Po kilku iteracjach dotyczących różnych składników i kombinacji, zakończyłem oszczędnym modelem następującej postaci:

Istnieje poziom stochastyczny + deterministyczna sezonowość + pewne wartości odstające, a dane nie mają innych wykrywalnych cech.

Poniżej znajduje się analiza istotności różnych składników. Zauważ, że użyłem trygonometrii (czyli sin / cos w oświadczeniu o sezonowości w PROC UCM) podobnym do @Elvis i @Nick Cox. Możesz także użyć fikcyjnego kodowania w UCM, a kiedy testowałem, oba dały podobne wyniki. W tej dokumentacji znajdziesz różnice między dwoma sposobami modelowania sezonowości w SAS.

Jak pokazano powyżej, masz wartości odstające: dwa impulsy i jedna zmiana poziomu w 2009 r. (Czy bańka ekonomiczna / mieszkaniowa odgrywała rolę po 2009 r.?), Co można wyjaśnić dalszą głęboką analizą nurkowania. Dobrą cechą użycia Proc UCMjest to, że zapewnia doskonałą wydajność graficzną.

Poniżej znajduje się sezonowość oraz połączony wykres trendów i sezonowości. To, co zostało, to hałas .

Ważniejszym testem diagnostycznym, jeśli chcesz użyć wartości p i testowania istotności, jest sprawdzenie, czy twoje resztki są pozbawione wzorów i normalnie rozmieszczone, co jest spełnione w powyższym modelu przy użyciu UCM i jak pokazano poniżej na pozostałych wykresach diagnostycznych, takich jak acf / pacf i inni.

Wniosek : Na podstawie analizy danych z wykorzystaniem UCM i testów istotności dane wydają się mieć sezonowość i obserwujemy dużą liczbę zgonów w miesiącach letnich w maju / czerwcu / lipcu, a najniższą w miesiącach zimowych w grudniu i lutym.

Dodatkowe uwagi : Proszę również wziąć pod uwagę praktyczne znaczenie wielkości wahań sezonowych. Aby zaprzeczyć sprzecznym z prawdą argumentom, skonsultuj się z ekspertami w dziedzinie w celu dalszego uzupełnienia i potwierdzenia twojej hipotezy.

W żadnym wypadku nie mówię, że jest to jedyne podejście do rozwiązania tego problemu. Cechą, która podoba mi się w UCM jest to, że umożliwia jawne modelowanie wszystkich funkcji szeregów czasowych i jest również bardzo wizualna.

Do wstępnej oceny wizualnej można zastosować następujący wykres. Rysując dane miesięczne z krzywą lessa i 95% przedziałem ufności, wydaje się, że wzrost w połowie roku osiąga maksimum w czerwcu. Inne czynniki mogą powodować, że dane mają szeroki rozkład, dlatego trend sezonowy może być maskowany na tym surowym wykresie lessowym. Punkty danych zostały roztrzęsione.

Edycja: Poniższy wykres przedstawia krzywą lessa i przedział ufności dla zmiany liczby obserwacji w stosunku do liczby w poprzednim miesiącu:

Pokazuje to również, że w miesiącach w pierwszej połowie roku liczba przypadków stale rośnie, a spada w drugiej połowie roku. To także sugeruje szczyt w połowie roku. Jednak przedziały ufności są szerokie i wynoszą 0, tzn. Brak zmian w ciągu roku, co wskazuje na brak istotności statystycznej.

Różnicę liczby miesięcy można porównać ze średnimi z poprzednich 3 miesięcy:

Pokazuje to wyraźny wzrost liczby w maju i spadek w październiku.