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

24

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

wprowadź opis zdjęcia tutaj

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ń. wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

     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)

wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

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.

wprowadź opis zdjęcia tutaj

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ąć)]

wprowadź opis zdjęcia tutaj

svannoy
źródło
1
Sformułowanie „jeden z naszych 50 stanów” oznacza, że ​​wszyscy czytelnicy należą do Stanów Zjednoczonych. Widocznie czai się tu także wielu kosmitów.
Nick Cox,
1
Czy to z publicznego zestawu danych? Czy możesz udostępnić dane z tygodnia na tydzień, a nawet z dnia na dzień?
Elvis
1
@Elvis - zamieściłem link do danych dziennej liczby. Dane pochodzą z akt zgonu, które są „publicznym rejestrem”, ale wymagają procesu; jednak zagregowane dane liczbowe nie. PS - sam spróbowałem linku i zadziałało, ale wcześniej nie publikowałem tego w publicznym folderze dropbox, więc daj mi znać, jeśli link nie działa.
svannoy
1
Ponieważ twoje dane są liczone, oczekiwałbym, że wariancja będzie związana ze średnią. Zwykłe modele szeregów czasowych nie uwzględniają tego (możesz jednak spróbować powiedzieć transformację , powiedzmy Freemana-Tukeya ), lub możesz spojrzeć na model szeregu czasowego, który jest przeznaczony do zliczania danych. (Jeśli tego nie zrobisz, może to nie stanowić
większego
1
ytμtVar(yt)=μt

Odpowiedzi:

13

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

sałata(2)πt12)grzech(2)πt12)

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

Elvis
źródło
2
p
2
Całkowicie się zgadzam, właśnie to sugerowałem :) „Z pewnością można powiedzieć, że to sugeruje efekt”; a nie „Jest efekt”! To, co uważam za interesujące, to transformacja trygonometryczna, jest ona bardzo naturalna i nie rozumiem, dlaczego nie jest widywana częściej. Ale to tylko punkt wyjścia ... ładowanie, ocena modelu ... wiele do zrobienia.
Elvis
1
Nie ma problemu! Mój zło wtedy nie byłem w stanie wykryć dowcipu. :)
usεr11852 mówi Przywróć Monic
2
+1. Poisson spotyka Fouriera ... Myślę, że ekonomiści i inni podkreślają zmienne wskaźnikowe, ponieważ sezonowość może być gwałtowna, ale podejście trygonometryczne często pomaga.
Nick Cox,
2
W rzeczy samej. Recenzja samouczka, którą napisałem, jest dostępna na stronie stata-journal.com/article.html?article=st0116
Nick Cox
8

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)

szacowany składnik sezonowy

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.

javlacalle
źródło
Dzięki @javiacalle, zbadam zaproponowane przez ciebie metody. Czy mogę zapytać o pański wniosek z opublikowanego wykresu, czy to fakt, że amplituda wzrasta w miarę upływu lat, jest podstawą waszego komentarza: „widzimy, że szacowany wzór sezonowy nie jest stabilny w czasie”, czy też robi to uwzględniać bardziej subtelną obserwację, że kształt każdego piku jest nieco inny? Zakładam to pierwsze, ale wiemy, dokąd prowadzą nas założenia.
svannoy
2
χ2)
@svannoy Główny wniosek oparty na metodach szeregów czasowych zastosowanych w mojej odpowiedzi jest taki, że w przykładowych danych nie ma wyraźnego wzorca sezonowego. Mimo że cykle sezonowe wyjaśniają część zmienności danych, wzorca sezonowego nie można wiarygodnie zidentyfikować, ponieważ jest on zasłonięty przez wysoki stopień nieregularnych wahań (można to również sprawdzić, wyświetlając funkcję wzmocnienia wybranego modelu ARIMA pokazaną w pytaniu) .
javlacalle
@oDDsKooL Zrobiłem również test chi-kwadrat w dzień tygodnia, sobota / niedziela są nieco poniżej oczekiwań, a poniedziałek / wtorek są nieco powyżej ....
svannoy
6

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.

wprowadź opis zdjęcia tutaj

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.

wprowadź opis zdjęcia tutaj

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 .

wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj

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.

wprowadź opis zdjęcia tutaj

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.

wróżka
źródło
Dzięki za tę odpowiedź i ciekawe komentarze. Nie znam UCM, wydaje się to bardzo interesujące, postaram się o tym pamiętać ...
Elvis
(+1) Interesująca analiza. Nadal byłbym ostrożny, jeśli chodzi o stwierdzenie istotnego deterministycznego wzorca sezonowego, ale twoje wyniki wymagają dokładniejszego przyjrzenia się tej możliwości. Pomocny może być test Canova i Hansena na stabilność sezonową, opisany na przykład tutaj .
javlacalle,
gretlπ/6
1
+1. Wiele interesujących i pomocnych komentarzy. Do twojej listy psychologów, psychiatrów i socjologów można dodać meteorolog / klimatolog. Taki człowiek chciałby dodać, że żadne dwa lata nie są identyczne w cyklach opadów i temperatur. Z grubsza przypuszczałbym, że zimą będzie więcej depresji (krótsze dni itp.), Ale tyle, ile przypuszczam, biorąc pod uwagę pewne dane.
Nick Cox
Dzięki @forecaster, to wiele dodaje do mojej nauki. Jestem psychologiem ze stopniem zdrowia publicznego. Dodałbym epidemiologa do twojej listy. Jak wspomniałem na początku, istnieje wiele mitologii (zwanych też teoretykami) na temat sezonowych trendów i samobójstw. Można argumentować za trendami sezonowymi w dowolnym kierunku, dlatego potrzebujemy analiz ilościowych, aby (nie) potwierdzić. Z punktu widzenia zdrowia publicznego, gdybyśmy znaleźli ostre nieciągłości, moglibyśmy ukierunkować interwencje. Nie widzę tego w tych danych. Z perspektywy teorii samobójstwa potwierdzenie nawet niewielkich trendów może wpłynąć na rozwój teorii.
svannoy
1

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.

wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

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

rnso
źródło
(-1) Istnieją już trzy wysokiej jakości odpowiedzi na to pytanie. Twoja odpowiedź również nie odpowiada na zadane pytanie - możesz opublikować je jako komentarz . Nie udzielasz odpowiedzi, w jaki sposób te dane mogą być analizowane.
Tim
Wcześniej opublikowałem tutaj komentarz (patrz poniżej pytania), ale nie mogę zamieszczać cyfr w komentarzach.
rnso
Chociaż rozumiem ten artykuł wstępny, powiem, że @rnso dostarczył ładny wykres, który ładnie ilustruje potencjalny składnik sezonowy i powinien być częścią mojego oryginalnego postu.
svannoy
Rozumiem to i zgadzam się, ale nadal nie jest to odpowiedź, ale komentarz lub poprawa. @rnso mógł zasugerować ci za pomocą komentarza, że ​​możesz przyjrzeć się lub dołączyć taką fabułę.
Tim
@Glen_b, @ Tim: Dodałem kolejną fabułę, która może być przydatna i której nie mogę komentować.
rnso 18.04.15