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 stl
funkcji 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.arima
zdecydował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
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.
[ 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ąć)]
źródło
Odpowiedzi:
Co z regresją Poissona?
Stworzyłem ramkę danych zawierającą twoje dane, a także indeks
t
czasu (w miesiącach) i zmiennąmonthdays
dla liczby dni w każdym miesiącu.Wygląda to tak:
Porównajmy teraz model z efektem czasowym i efektem liczby dni z modelem, w którym dodajemy efekt miesiąca:
Oto podsumowanie „małego” modelu:
Widać, że te dwie zmienne mają w znacznym stopniu znaczący efekt krańcowy. Teraz spójrz na większy model:
Oczywiście
monthdays
efekt 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.Teraz porównaj to z modelem zerowym:
Można więc z pewnością powiedzieć, że sugeruje to efekt sezonowy ...
źródło
Test chi-kwadrat jest dobrym podejściem jako wstępny pogląd na twoje pytanie.
stl
Rozkł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: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:
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.
źródło
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:
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:
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”:
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.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 UCM
jest 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.
źródło
gretl
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.
źródło