Esse é um método apropriado para testar efeitos sazonais nos dados da contagem de suicídios?

24

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

insira a descrição da imagem aqui

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. insira a descrição da imagem aqui

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 stlfunçã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")

insira a descrição da imagem aqui

     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)

insira a descrição da imagem aqui

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

insira a descrição da imagem aqui

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

insira a descrição da imagem aqui

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

daily_suicide_data.csv

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.

insira a descrição da imagem aqui

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.

state.daylight.2002.csv

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

insira a descrição da imagem aqui

svannoy
fonte
1
A expressão "um dos nossos 50 estados" implica que todos os leitores pertencem aos Estados Unidos. Manifestamente, muitos aliens também se escondem aqui.
Nick Cox
1
Isso é de um conjunto de dados público? Você poderia disponibilizar os dados semana a semana ou mesmo dia a dia?
Elvis
1
@ Elvis - publiquei um link para os dados da contagem diária. Os dados provêm de atestados de óbito que são 'registro público', mas requerem um processo para serem obtidos; no entanto, os dados agregados da contagem não. PS - Tentei o link pessoalmente e funcionou, mas não publiquei uma pasta pública da caixa de depósito dessa maneira antes; portanto, informe-me se o link não funcionar.
precisa saber é o seguinte
1
Como seus dados são contados, espero que a variação esteja relacionada à média. Os modelos habituais de séries temporais não explicam isso (no entanto, você pode tentar dizer uma transformação , talvez um Freeman-Tukey , por exemplo) ou pode olhar para um modelo de série temporal projetado para dados de contagem. (Se você não fizer isso, pode não ser um grande problema, já que o número só varia ao longo de um fator de dois ou mais.)
Glen_b -Reinstate Monica
1
ytμtVar(yt)=μt

Respostas:

13

Que tal uma regressão de Poisson?

Criei um quadro de dados contendo seus dados, além de um índice tpara o tempo (em meses) e uma variável monthdayspara o número de dias em cada mês.

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

Então fica assim:

> 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

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:

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

Aqui está o resumo do modelo "pequeno":

> 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

Você pode ver que as duas variáveis ​​têm efeitos marginais bastante significativos. Agora veja o modelo maior:

> 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

Bem, é claro que o monthdaysefeito 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.

> 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

porque(2πt12)pecado(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

Agora compare-o com o modelo nulo:

> 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

Então, pode-se dizer com certeza que isso sugere um efeito sazonal ...

Elvis
fonte
2
p
2
Concordo plenamente, foi o que sugeri :) "Pode-se dizer com certeza que isso sugere um efeito"; e não "Há um efeito"! O que acho interessante é essa transformação trigonométrica, é muito natural e não entendo por que não é vista com mais frequência. Mas este é apenas um ponto de partida ... inicialização, avaliação do modelo ... muito o que fazer.
Elvis
1
Sem probs! Meu mal então, eu era incapaz de detectar a piada. :)
usεr11852 diz Reinstate Monic
2
+1. Poisson conhece Fourier ... Acho que economistas e alguns outros enfatizam as variáveis ​​indicadoras porque a sazonalidade pode ser espetada, mas a abordagem trigonométrica geralmente ajuda.
Nick Cox
2
De fato. Uma revisão tutorial que eu escrevi é acessível em stata-journal.com/article.html?article=st0116
Nick Cox
8

Um teste do qui-quadrado é uma boa abordagem como uma visão preliminar da sua pergunta.

A stldecomposiçã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:

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

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:

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)

componente sazonal estimado

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.

javlacalle
fonte
Obrigado @ javiacalle, estarei investigando os métodos que você sugere. Posso perguntar sobre a sua conclusão do gráfico que você publicou, é o fato de a amplitude aumentar à medida que os anos progridem, que é a base do seu comentário, "podemos ver que o padrão sazonal estimado não é estável ao longo do tempo" ou faz isso incluir a observação mais sutil de que a forma de cada pico é um pouco diferente? Eu assumo o primeiro, mas sabemos aonde as suposições nos levam.
svannoy
2
χ2
@svannoy A principal conclusão baseada nos métodos de séries temporais usados ​​em minha resposta é que não há um padrão sazonal claro nos dados da amostra. Apesar dos ciclos sazonais explicarem parte da variabilidade dos dados, um padrão sazonal não pode ser identificado com segurança, pois é obscurecido por um alto grau de flutuações irregulares (isso também pode ser verificado, exibindo a função de ganho do modelo ARIMA escolhido, mostrado na pergunta) .
Javlacalle
@oDDsKooL Também fiz o teste do qui-quadrado em um dia da semana, sábado / domingo são um pouco abaixo das expectativas e Segunda / Terça, estão logo acima ....
svannoy
6

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:

  • Psicólogo
  • Psiquiatra
  • Sociólogo

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:

  1. Estocástico - muda aleatoriamente e é difícil de prever
  2. Determinista - não muda, perfeitamente previsível. Você pode usar manequim ou trigonometria (sin / cos etc.,) para modelar

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

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

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.

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;

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.

insira a descrição da imagem aqui

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.

insira a descrição da imagem aqui

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 .

insira a descrição da imagem aqui insira a descrição da imagem aqui

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.

insira a descrição da imagem aqui

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.

previsor
fonte
Obrigado por esta resposta e por comentários interessantes. Eu não sei UCM, parece muito interessante, vou tentar manter isso em mente ...
Elvis
(+1) análise interessante. Eu ainda seria cauteloso ao concluir a presença de um padrão sazonal determinístico significativo, mas seus resultados exigem uma análise mais detalhada dessa possibilidade. O teste de Canova e Hansen para estabilidade sazonal pode ser útil, é descrito, por exemplo, aqui .
Javlacalle # 12/15
gretlπ/6
1
+1. Muitos comentários interessantes e úteis. À sua lista de psicólogos, psiquiatras, sociólogos podem ser adicionados meteorologista / climatologista. Essa pessoa gostaria de acrescentar que não há dois anos idênticos nos ciclos de precipitação e temperatura. Eu teria imaginado grosseiramente que havia mais depressão no inverno (dias mais curtos, etc.), mas muito para adivinhar dados alguns dados.
Nick Cox
Obrigado @forecaster, isso adiciona muito ao meu aprendizado. Sou psicóloga, formada em saúde pública. Eu adicionaria epidemiologista à sua lista. Como mencionei no início, há muita mitologia (também conhecida como teorização) sobre tendências sazonais e suicídio. Pode-se argumentar fortemente sobre tendências sazonais, em qualquer direção, por isso precisamos de análises quantitativas para (des) confirmar. Do ponto de vista da saúde pública, se encontrássemos descontinuidades acentuadas, poderíamos direcionar intervenções. Não estou vendo isso nesses dados. Do ponto de vista da teoria do suicídio, confirmar até pequenas tendências poderia informar o desenvolvimento da teoria.
precisa saber é o seguinte
1

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.

insira a descrição da imagem aqui

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:

insira a descrição da imagem aqui

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:

insira a descrição da imagem aqui

Isso mostra um claro aumento nos números em maio e uma queda em outubro.

rnso
fonte
(-1) Já existem três respostas de alta qualidade para esta pergunta. Sua resposta também não responde à pergunta postada - você pode publicá-la como um comentário . Você não fornece respostas sobre como esses dados podem ser analisados.
Tim
Eu já tinha publicado comentários aqui (veja abaixo a pergunta), mas não consigo postar números nos comentários.
rnso
Embora eu compreenda o editorial aqui, direi que o @rnso forneceu um bom gráfico que ilustra bem o potencial componente sazonal e deveria ter sido parte do meu post original.
precisa saber é o seguinte
Entendo isso e concordo, mas ainda assim não é uma resposta, mas um comentário ou melhoria. @rnso poderia ter sugerido a você por meio de comentário que você poderia procurar ou incluir esse enredo.
Tim
@Glen_b, @ Tim: adicionei outro enredo que pode ser útil e que não posso colocar em um comentário.
rnso