У меня есть 17 лет (с 1995 по 2011) данных свидетельств о смерти, связанных со смертями от самоубийств для штата в США. Существует много мифологий о самоубийствах и месяцах / сезонах, большая часть которых противоречива, и литературы, которую я ' После проверки я не получил четкого представления о применяемых методах или уверенности в результатах.
Поэтому я решил посмотреть, смогу ли я определить, могут ли самоубийства быть более или менее вероятными в каком-либо конкретном месяце в моем наборе данных. Все мои анализы сделаны в R.
Общее количество самоубийств в данных 13,909.
Если вы посмотрите на год с наименьшим количеством самоубийств, они происходят в 309/365 дней (85%). Если вы посмотрите на год с наибольшим количеством самоубийств, они происходят в течение 339/365 дней (93%).
Таким образом, каждый год существует достаточное количество дней без самоубийств. Однако при агрегировании за все 17 лет самоубийства происходят каждый день в году, включая 29 февраля (хотя только 5, когда в среднем 38).
Простое суммирование количества самоубийств в каждый день года не указывает на явную сезонность (на мой взгляд).
Среднемесячные самоубийства в месяц в совокупности варьируются от:
(m = 65, sd = 7,4, до m = 72, sd = 11,1)
Мой первый подход состоял в том, чтобы агрегировать набор данных по месяцам за все годы и провести тест хи-квадрат после вычисления ожидаемых вероятностей для нулевой гипотезы о том, что не было систематической дисперсии числа самоубийств по месяцам. Я рассчитал вероятности для каждого месяца с учетом количества дней (и поправив февраль на високосные годы).
Результаты хи-квадрат не показали значительных изменений по месяцам:
# So does the sample match expected values?
chisq.test(monthDat$suicideCounts, p=monthlyProb)
# Yes, X-squared = 12.7048, df = 11, p-value = 0.3131
Изображение ниже показывает общее количество в месяц. Горизонтальные красные линии расположены в ожидаемых значениях для февраля, 30-дневного месяца и 31-дневного месяца соответственно. В соответствии с критерием хи-квадрат ни один месяц не выходит за пределы 95% доверительного интервала для ожидаемых значений.
Я думал, что закончил, пока не начал исследовать данные временных рядов. Как и многие другие, я начал с непараметрического метода сезонной декомпозиции, используя stl
функцию из пакета статистики.
Чтобы создать данные временных рядов, я начал с агрегированных ежемесячных данных:
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
А затем выполнил stl()
разложение
# Seasonal decomposition
suicideByMonthFit <- stl(suicideByMonthTs, s.window="periodic")
plot(suicideByMonthFit)
В этот момент я забеспокоился, потому что мне кажется, что есть и сезонная составляющая, и тенденция. После долгих интернет-исследований я решил следовать указаниям Роба Хиндмана и Джорджа Атанасопулоса, изложенным в их тексте онлайн «Прогнозирование: принципы и практика», в частности, чтобы применить сезонную модель ARIMA.
Я использовал adf.test()
и kpss.test()
для оценки стационарности и получил противоречивые результаты. Они оба отвергли нулевую гипотезу (отметив, что они проверяют противоположную гипотезу).
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
Затем я использовал алгоритм из книги, чтобы посмотреть, смогу ли я определить количество различий, которое необходимо сделать как для тренда, так и для сезона. Я закончил с nd = 1, ns = 0.
Затем я побежал auto.arima
, выбрав модель, в которой были и тренд, и сезонный компонент, и константа типа «дрейф».
# 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
Наконец, я посмотрел на остатки от подгонки и, если я правильно понимаю, так как все значения находятся в пределах пороговых значений, они ведут себя как белый шум, и, таким образом, модель является довольно разумной. Я запустил тест portmanteau, как описано в тексте, у которого значение ap было намного выше 0,05, но я не уверен, что у меня правильные параметры.
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
Вернувшись назад и снова прочитав главу о моделировании аримы, я понимаю, что все- auto.arima
таки решил смоделировать тренд и сезон. И я также понимаю, что прогнозирование - это не тот анализ, который мне, вероятно, следует делать. Я хочу знать, должен ли определенный месяц (или, более обычно, время года) быть помечен как месяц высокого риска. Кажется, что инструменты в литературе по прогнозированию очень уместны, но, возможно, не самый лучший для моего вопроса. Любой вклад очень ценится.
Я публикую ссылку на CSV-файл, который содержит ежедневные счета. Файл выглядит так:
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
Подсчет - это число самоубийств, которые произошли в этот день. «t» - числовая последовательность от 1 до общего количества дней в таблице (5533).
Я принял к сведению комментарии ниже и подумал о двух вещах, связанных с моделированием самоубийства и времен года. Во-первых, что касается моего вопроса, месяцы - это просто прокси для обозначения смены сезона, я не заинтересован в том, отличается ли месяц от других (это, конечно, интересный вопрос, но это не то, что я намеревался задать) исследования). Следовательно, я думаю, что имеет смысл выравнивать месяцы, просто используя первые 28 дней всех месяцев. Когда вы делаете это, вы становитесь немного хуже, что я интерпретирую как еще одно свидетельство отсутствия сезонности. В приведенных ниже выходных данных первое совпадение представляет собой воспроизведение из ответа ниже с использованием месяцев с их истинным числом дней, за которым следует набор данных suicideByShortMonthв котором число самоубийств рассчитывалось с первых 28 дней всех месяцев. Меня интересует, что люди думают о том, действительно ли это регулирование является хорошей идеей, не нужно или вредно?
> 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
Второе, на что я обратил внимание, это вопрос использования месяца в качестве прокси для сезона. Возможно, лучшим показателем времени года является количество дневных часов, которые получает район. Эти данные поступают из северного штата, в котором наблюдается значительное изменение дневного света. Ниже приведен график дневного света 2002 года.
Когда я использую эти данные, а не месяц года, эффект все еще значителен, но эффект очень и очень мал. Остаточное отклонение намного больше, чем у моделей выше. Если дневные часы - лучшая модель для сезонов, и подгонка не так хороша, является ли это большим свидетельством очень небольшого сезонного эффекта?
> 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
Я публикую дневные часы на случай, если кто-нибудь захочет поиграть с этим. Обратите внимание, что это не високосный год, поэтому, если вы хотите указать минуты високосных лет, либо экстраполируйте, либо извлеките данные.
[ Изменить, чтобы добавить заговор из удаленного ответа (надеюсь, rnso не возражает, чтобы я переместил заговор в удаленном ответе здесь на вопрос. Сванной, если вы не хотите, чтобы это добавлялось в конце концов, вы можете отменить его)]
источник
Ответы:
Как насчет регрессии Пуассона?
Я создал фрейм данных, содержащий ваши данные, а также индекс
t
для времени (в месяцах) и переменнуюmonthdays
для количества дней в каждом месяце.Так это выглядит так:
Теперь давайте сравним модель с эффектом времени и числом дней с моделью, в которую мы добавили эффект месяца:
Вот краткое изложение «маленькой» модели:
Вы можете видеть, что две переменные имеют в значительной степени предельные эффекты. Теперь посмотрим на большую модель:
Ну, конечно,
monthdays
эффект исчезает; это можно оценить только благодаря високосным годам !! Сохранение его в модели (и с учетом високосных лет) позволяет использовать остаточные отклонения для сравнения двух моделей.Теперь сравните это с нулевой моделью:
Так что можно с уверенностью сказать, что это говорит о сезонном эффекте ...
источник
Тест хи-квадрат - это хороший подход в качестве предварительного взгляда на ваш вопрос.
stl
Разложение может ввести в заблуждение , в качестве инструмента для проверки на наличие сезонности. Этой процедуре удается вернуть стабильную сезонную модель, даже если в качестве входного сигнала передается белый шум (случайный сигнал без структуры). Попробуйте например:Просмотр заказов, выбранных с помощью процедуры автоматического выбора модели ARIMA, также немного рискован, поскольку сезонная модель ARIMA не всегда предполагает сезонность (подробности см. В этом обсуждении ). В этом случае выбранная модель генерирует сезонные циклы, и комментарий @RichardHardy является разумным, однако требуется дополнительное понимание, чтобы сделать вывод, что самоубийства обусловлены сезонным характером.
Ниже я суммирую некоторые результаты, основанные на анализе ежемесячных серий, которые вы опубликовали. Это сезонная составляющая, рассчитанная по базовой модели структурных временных рядов:
Аналогичный компонент был извлечен с использованием программного обеспечения TRAMO-SEATS с опциями по умолчанию. Мы можем видеть, что предполагаемая сезонная картина нестабильна во времени и, следовательно, не поддерживает гипотезу повторяющейся картины количества самоубийств за месяцы в течение периода выборки. Запустив программное обеспечение X-13ARIMA-SEATS с опциями по умолчанию, комбинированный тест на сезонность показал, что идентифицируемой сезонности нет.
Изменить (см. Этот ответ и мой комментарий ниже для определения идентифицируемой сезонности ).
Учитывая характер ваших данных, стоило бы дополнить этот анализ, основанный на методах временных рядов, моделью для подсчета данных (например, моделью Пуассона) и проверить значимость сезонности в этой модели. Добавление данных подсчета тегов к вашему вопросу может привести к увеличению числа просмотров и потенциальных ответов в этом направлении.
источник
Как отмечается в моем комментарии, это очень интересная проблема. Выявление сезонности - это не только статистическое упражнение. Разумным подходом было бы консультироваться с теорией и экспертами, такими как:
по этой проблеме, чтобы понять «почему» будет сезонность, чтобы дополнить анализ данных. Что касается данных, я использовал отличный метод декомпозиции, называемый моделью ненаблюдаемых компонентов (UCM), который является формой метода пространства состояний. Смотрите также эту очень доступную статью Koopman. Мой подход похож на @Javlacalle. Он не только предоставляет инструмент для разложения данных временных рядов, но также объективно оценивает наличие или отсутствие сезонности посредством проверки значимости. я не большой поклонник тестирования значимости неэкспериментальных данных, но я не знаю ни одной другой процедуры, чтобы вы могли проверить свою гипотезу о наличии / отсутствии сезонности на данных временных рядов.
Многие игнорируют, но очень важная особенность, которую хотелось бы понять, это тип сезонности:
Для длинных данных временного ряда, таких как ваш, возможно, что сезонность могла бы измениться с течением времени. Опять же, UCM - единственный известный мне подход, который может обнаружить эти стохастические / детерминированные сезонности. UCM может разложить вашу проблему на следующие «компоненты»:
Вы также можете проверить, является ли уровень, наклон, цикл детерминированным или стохастическим. Пожалуйста, обратите внимание, что
level + slope = trend
. Ниже я приведу некоторый анализ ваших данных с использованием UCM. Я использовал SAS для анализа.После нескольких итераций с учетом различных компонентов и комбинаций я закончил экономной моделью следующего вида:
Существует стохастический уровень + детерминированная сезонность + некоторые выбросы, и у данных нет никаких других обнаруживаемых признаков.
Ниже приведен анализ значимости различных компонентов. Обратите внимание, что я использовал тригонометрию (это sin / cos в выражении сезонности в PROC UCM), аналогично @Elvis и @Nick Cox. Вы также можете использовать фиктивное кодирование в UCM, и когда я тестировал, оба дали схожие результаты. См. Эту документацию для различий между двумя способами моделирования сезонности в SAS.
Как показано выше, у вас есть выбросы: два импульса и одно смещение уровня в 2009 году (сыграла ли роль экономика / жилищный пузырь после 2009 года ??), что можно объяснить дальнейшим глубоким анализом погружений. Хорошая особенность использования
Proc UCM
является то, что он обеспечивает превосходный графический вывод.Ниже представлена сезонность и комбинированный тренд и сезонность сюжета. Все, что осталось, это шум .
Более важным диагностическим тестом, если вы хотите использовать значения p и тестирование значимости, является проверка того, являются ли ваши остатки без шаблонов и нормально ли распределены, что удовлетворяется в приведенной выше модели с использованием UCM и как показано ниже на графиках остаточных значений, таких как acf / pacf и другие.
Вывод : на основе анализа данных с использованием UCM и тестирования значимости данные, как представляется, имеют сезонность, и мы видим большое количество смертей в летние месяцы май / июнь / июль и самое низкое в зимние месяцы декабрь и февраль.
Дополнительные соображения : Просьба также учитывать практическую значимость величины сезонных колебаний. Чтобы опровергнуть противоречивые аргументы, проконсультируйтесь с экспертами в предметной области, чтобы дополнительно дополнить и подтвердить свою гипотезу.
Я ни в коем случае не говорю, что это единственный подход к решению этой проблемы. Особенность, которая мне нравится в UCM, заключается в том, что она позволяет вам явно моделировать все функции временных рядов, а также очень наглядна.
источник
gretl
Для первоначальной визуальной оценки можно использовать следующий график. Построение ежемесячных данных с кривой лёсса и 95% -ным доверительным интервалом показывает, что в июне наблюдается пик роста в середине года. Другие факторы могут быть причиной того, что данные будут иметь широкое распространение, поэтому сезонный тренд может маскироваться на этом графике лёсса необработанных данных. Точки данных были повреждены.
Изменить: На следующем графике показаны кривая лёсса и доверительный интервал для изменения количества случаев от числа в предыдущем месяце:
Это также показывает, что в течение месяцев в первой половине года число случаев продолжает расти, в то время как во второй половине года они снижаются. Это также говорит о пике в середине года. Однако доверительные интервалы широки и пересекают 0, т.е. не меняются в течение года, что указывает на отсутствие статистической значимости.
Разницу в количестве месяцев можно сравнить со средними значениями за предыдущие 3 месяца:
Это показывает явный рост числа в мае и падение в октябре.
источник