Tenho 17 anos (1995 a 2011) de dados de atestados de óbito relacionados a mortes por suicídio em um estado nos EUA. Existe muita mitologia por aí sobre suicídios e os meses / estações do ano, muitas das quais contraditórias e a literatura. Ao analisar, não entendo claramente os métodos usados ou a confiança nos resultados.
Então, decidi ver se consigo determinar se há mais ou menos probabilidade de suicídio em um determinado mês dentro do meu conjunto de dados. Todas as minhas análises são feitas em R.
O número total de suicídios nos dados é 13.909.
Se você olhar o ano com menos suicídios, eles ocorrerão em 309/365 dias (85%). Se você observar o ano com mais suicídios, eles ocorrerão em 339/365 dias (93%).
Portanto, há um número razoável de dias por ano sem suicídios. No entanto, quando agregados em todos os 17 anos, há suicídios todos os dias do ano, incluindo 29 de fevereiro (embora apenas 5 quando a média é 38).
Simplesmente somar o número de suicídios em cada dia do ano não indica uma sazonalidade clara (para os meus olhos).
Agregados no nível mensal, a média de suicídios por mês varia de:
(m = 65, sd = 7,4, m = 72, sd = 11,1)
Minha primeira abordagem foi agregar os dados estabelecidos por mês para todos os anos e fazer um teste qui-quadrado após calcular as probabilidades esperadas para a hipótese nula, de que não havia variação sistemática nas contagens de suicídio por mês. Calculei as probabilidades de cada mês, levando em consideração o número de dias (e ajustando fevereiro por anos bissextos).
Os resultados do qui-quadrado não indicaram variação significativa por mês:
# So does the sample match expected values?
chisq.test(monthDat$suicideCounts, p=monthlyProb)
# Yes, X-squared = 12.7048, df = 11, p-value = 0.3131
A imagem abaixo indica a contagem total por mês. As linhas vermelhas horizontais são posicionadas nos valores esperados para fevereiro, 30 dias e 31 dias, respectivamente. Consistente com o teste do qui-quadrado, nenhum mês está fora do intervalo de confiança de 95% para as contagens esperadas.
Eu pensei que tinha terminado até começar a investigar dados de séries temporais. Como imagino que muitas pessoas fazem, comecei com o método de decomposição sazonal não paramétrico usando a stl
função no pacote de estatísticas.
Para criar os dados da série temporal, comecei com os dados mensais agregados:
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
E então realizou a stl()
decomposição
# Seasonal decomposition
suicideByMonthFit <- stl(suicideByMonthTs, s.window="periodic")
plot(suicideByMonthFit)
Nesse ponto, fiquei preocupado porque me parece que há um componente sazonal e uma tendência. Depois de muita pesquisa na Internet, decidi seguir as instruções de Rob Hyndman e George Athanasopoulos, conforme dispostas em seu texto on-line "Previsão: princípios e prática", especificamente para aplicar um modelo ARIMA sazonal.
Eu usei adf.test()
e kpss.test()
avaliei a estacionariedade e obtive resultados conflitantes. Ambos rejeitaram a hipótese nula (observando que testam a hipótese oposta).
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
Em seguida, usei o algoritmo do livro para verificar se era possível determinar a quantidade de diferenciação necessária para a tendência e a estação. Eu terminei com nd = 1, ns = 0.
Eu então corri auto.arima
, que escolheu um modelo que possuía uma tendência e um componente sazonal, juntamente com uma constante do tipo "deriva".
# 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
Finalmente, observei os resíduos do ajuste e, se entendi corretamente, como todos os valores estão dentro dos limites, eles se comportam como ruído branco e, portanto, o modelo é bastante razoável. Fiz um teste do portmanteau, conforme descrito no texto, que tinha um valor de p bem acima de 0,05, mas não tenho certeza se tenho os parâmetros corretos.
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
Depois de voltar e ler o capítulo sobre modelagem de arima, percebo agora que auto.arima
optou por modelar tendência e estação. E também estou percebendo que a previsão não é especificamente a análise que eu provavelmente deveria estar fazendo. Quero saber se um mês específico (ou, geralmente, época do ano) deve ser sinalizado como um mês de alto risco. Parece que as ferramentas na literatura de previsão são altamente pertinentes, mas talvez não sejam as melhores para minha pergunta. Toda e qualquer informação é muito apreciada.
Estou postando um link para um arquivo CSV que contém as contagens diárias. O arquivo fica assim:
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
Contagem é o número de suicídios que aconteceram naquele dia. "t" é uma sequência numérica de 1 ao número total de dias na tabela (5533).
Tomei nota dos comentários abaixo e pensei em duas coisas relacionadas à modelagem de suicídio e temporadas. Primeiro, no que diz respeito à minha pergunta, os meses são simplesmente proxies para marcar a mudança de estação, não estou interessado em saber se um mês em particular é diferente dos outros (é claro que é uma pergunta interessante, mas não é o que me propus a investigar). Portanto, acho que faz sentido equalizar os meses simplesmente usando os primeiros 28 dias de todos os meses. Quando você faz isso, fica um pouco pior, o que estou interpretando como mais uma evidência da falta de sazonalidade. Na saída abaixo, o primeiro ajuste é uma reprodução de uma resposta abaixo usando meses com seu número real de dias, seguido por um conjunto de dados suicideByShortMonthem que as contagens de suicídio foram calculadas desde os primeiros 28 dias de todos os meses. Estou interessado no que as pessoas pensam sobre se esse ajuste é uma boa ideia, não é necessária ou prejudicial?
> 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
A segunda coisa que eu examinei mais é a questão de usar o mês como proxy para a temporada. Talvez um melhor indicador de estação seja o número de horas de luz do dia que uma área recebe. Esses dados são provenientes de um estado do norte que possui uma variação substancial na luz do dia. Abaixo está um gráfico da luz do dia de 2002.
Quando uso esses dados em vez do mês do ano, o efeito ainda é significativo, mas o efeito é muito, muito pequeno. O desvio residual é muito maior que os modelos acima. Se o horário de verão é um modelo melhor para as estações do ano e o ajuste não é tão bom, isso é mais uma evidência de um efeito sazonal muito pequeno?
> 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
Estou postando o horário de verão, caso alguém queira brincar com isso. Observe que este não é um ano bissexto; portanto, se você deseja colocar minutos para os anos bissextos, extrapole ou recupere os dados.
[ Editar para adicionar um gráfico da resposta excluída (espero que rnso não me importe de mover o gráfico na resposta excluída aqui para a pergunta. Svannoy, se você não quiser que isso seja adicionado, afinal, você pode revertê-lo)]
fonte
Respostas:
Que tal uma regressão de Poisson?
Criei um quadro de dados contendo seus dados, além de um índice
t
para o tempo (em meses) e uma variávelmonthdays
para o número de dias em cada mês.Então fica assim:
Agora vamos comparar um modelo com um efeito de tempo e um número de dias com um modelo no qual adicionamos um efeito de mês:
Aqui está o resumo do modelo "pequeno":
Você pode ver que as duas variáveis têm efeitos marginais bastante significativos. Agora veja o modelo maior:
Bem, é claro que o
monthdays
efeito desaparece; só pode ser estimado graças a anos bissextos !! Mantê-lo no modelo (e levando em consideração os anos bissextos) permite usar os desvios residuais para comparar os dois modelos.Agora compare-o com o modelo nulo:
Então, pode-se dizer com certeza que isso sugere um efeito sazonal ...
fonte
Um teste do qui-quadrado é uma boa abordagem como uma visão preliminar da sua pergunta.
A
stl
decomposição pode ser enganosa como uma ferramenta para testar a presença de sazonalidade. Este procedimento consegue retornar um padrão sazonal estável, mesmo que um ruído branco (sinal aleatório sem estrutura) seja passado como entrada. Tente por exemplo:Observar as ordens escolhidas por um procedimento de seleção automática de modelo ARIMA também é um pouco arriscado, pois um modelo ARIMA sazonal nem sempre envolve sazonalidade (para detalhes, consulte esta discussão ). Nesse caso, o modelo escolhido gera ciclos sazonais e o comentário de @RichardHardy é razoável; no entanto, é necessário um insight adicional para concluir que os suicídios são motivados por um padrão sazonal.
Abaixo, sumario alguns resultados com base na análise das séries mensais que você publicou. Esse é o componente sazonal estimado no modelo estrutural básico de série temporal:
Um componente semelhante foi extraído usando o software TRAMO-SEATS com opções padrão. Podemos ver que o padrão sazonal estimado não é estável ao longo do tempo e, portanto, não suporta a hipótese de um padrão recorrente no número de suicídios ao longo de meses durante o período da amostra. Executando o software X-13ARIMA-SEATS com opções padrão, o teste combinado de sazonalidade concluiu que a sazonalidade identificável não está presente.
Editar (veja esta resposta e meu comentário abaixo para obter uma definição de sazonalidade identificável ).
Dada a natureza de seus dados, vale a pena complementar esta análise com base em métodos de séries temporais com um modelo para dados de contagem (por exemplo, modelo de Poisson) e testar a significância da sazonalidade nesse modelo. Adicionar os dados de contagem de tags à sua pergunta pode resultar em mais visualizações e possíveis respostas nessa direção.
fonte
Como observado no meu comentário, este é um problema muito interessante. Detectar a sazonalidade não é apenas um exercício estatístico. Uma abordagem razoável seria consultar a teoria e especialistas como:
sobre esse problema para entender "por que" haveria sazonalidade para complementar a análise de dados. Chegando aos dados, usei um excelente método de decomposição chamado modelo de componentes não observados (UCM), que é uma forma de método de espaço de estado. Veja também este artigo muito acessível de Koopman. Minha abordagem é semelhante à @Javlacalle. Ele não apenas fornece uma ferramenta para decompor dados de séries temporais, mas também avalia objetivamente a presença ou ausência de sazonalidade por meio de testes de significância. Não sou muito fã de testes de significância em dados não experimentais, mas não conheço outro procedimento que possa testar sua hipótese sobre presença / ausência de sazonalidade em dados de séries temporais.
Muitos ignoram, mas uma característica muito importante que alguém gostaria de entender é o tipo de sazonalidade:
Para dados demorados de séries temporais como a sua, é possível que a sazonalidade tenha mudado ao longo do tempo. Novamente, a UCM é a única abordagem que eu sei que pode detectar essa sazonalidade estocástica / determinística. O UCM pode decompor seu problema nos seguintes "componentes":
Você também pode testar se o nível, inclinação, ciclo é determinístico ou estocástico. Por favor note isso
level + slope = trend
. A seguir, apresento algumas análises sobre seus dados usando o UCM. Eu usei o SAS para fazer a análise.Após várias iterações considerando diferentes componentes e combinações, terminei com um modelo parcimonioso da seguinte forma:
Há um nível estocástico + sazonalidade determinística + alguns valores discrepantes e os dados não possuem outros recursos detectáveis.
Abaixo está a análise de significância de vários componentes. Observe que eu usei trigonometria (que é sin / cos na declaração de sazonalidade no PROC UCM) semelhante a @Elvis e @Nick Cox. Você também pode usar codificação fictícia no UCM e, quando eu testei, ambos deram resultados semelhantes. Consulte esta documentação para obter diferenças entre as duas maneiras de modelar a sazonalidade no SAS.
Como mostrado acima, você tem discrepâncias: dois pulsos e uma mudança de nível em 2009 (a bolha da economia / habitação teve um papel após 2009 ??), o que pode ser explicado por uma análise mais aprofundada. Uma boa característica do uso
Proc UCM
é que ele fornece excelente saída gráfica.Abaixo está a sazonalidade e um gráfico combinado de tendência e sazonalidade. O que sobra é barulho .
Um teste de diagnóstico mais importante, se você deseja usar os valores de p e o teste de significância, é verificar se seus resíduos são sem padrão e normalmente distribuídos, o que é satisfeito no modelo acima usando o UCM e como mostrado abaixo nos gráficos de diagnóstico residual como acf / pacf e outros.
Conclusão : Com base na análise de dados usando UCM e testes de significância, os dados parecem ter sazonalidade e observamos um alto número de mortes nos meses de verão de maio / junho / julho e menores nos meses de inverno de dezembro e fevereiro.
Considerações adicionais : Considere também o significado prático da magnitude da variação sazonal. Para negar argumentos contrafactuais, consulte especialistas em domínio para complementar e validar ainda mais sua hipótese.
Não estou dizendo que essa é a única abordagem para resolver esse problema. O recurso que eu gosto no UCM é que ele permite modelar explicitamente todos os recursos da série temporal e também é altamente visual.
fonte
gretl
Para estimativa visual inicial, o gráfico a seguir pode ser usado. Traçando os dados mensais com curva de loess e seu intervalo de confiança de 95%, parece que há um aumento no meio do ano em junho. Outros fatores podem estar causando uma ampla distribuição dos dados; portanto, a tendência sazonal pode estar sendo mascarada nesse gráfico de perda de dados brutos. Os pontos de dados foram tremidos.
Editar: O gráfico a seguir mostra a curva de loess e o intervalo de confiança para alteração no número de casos em relação ao número do mês anterior:
Isso também mostra que, durante os meses da primeira metade do ano, o número de casos continua aumentando, enquanto eles caem na segunda metade do ano. Isso também sugere um pico no meio do ano. No entanto, os intervalos de confiança são amplos e ultrapassam 0, ou seja, nenhuma alteração ao longo do ano, indicando uma falta de significância estatística.
A diferença do número de um mês pode ser comparada com a média dos valores dos três meses anteriores:
Isso mostra um claro aumento nos números em maio e uma queda em outubro.
fonte