¿Es este un método apropiado para evaluar los efectos estacionales en los datos de recuento de suicidios?

Tengo 17 años (1995 a 2011) de datos de certificados de defunción relacionados con muertes por suicidio para un estado en los EE. UU. Hay mucha mitología sobre suicidios y los meses / estaciones, en gran parte contradictorios, y de la literatura I ' He revisado, no tengo una idea clara de los métodos utilizados o la confianza en los resultados.

Así que me dispuse a ver si puedo determinar si los suicidios tienen más o menos probabilidades de ocurrir en un mes determinado dentro de mi conjunto de datos. Todos mis análisis se realizan en R.

El número total de suicidios en los datos es de 13,909.

Si observa el año con la menor cantidad de suicidios, ocurren en 309/365 días (85%). Si observa el año con la mayor cantidad de suicidios, ocurren en 339/365 días (93%).

Entonces, hay una buena cantidad de días al año sin suicidios. Sin embargo, cuando se agregan en los 17 años, hay suicidios todos los días del año, incluido el 29 de febrero (aunque solo 5 cuando el promedio es 38).

Simplemente sumar el número de suicidios en cada día del año no indica una clara estacionalidad (a mi entender).

Agregados a nivel mensual, el promedio de suicidios por mes varía de:

(m = 65, sd = 7.4, a m = 72, sd = 11.1)

Mi primer enfoque fue agregar el conjunto de datos por mes para todos los años y hacer una prueba de chi-cuadrado después de calcular las probabilidades esperadas para la hipótesis nula, de que no había una variación sistemática en los recuentos de suicidios por mes. Calculé las probabilidades para cada mes teniendo en cuenta la cantidad de días (y ajustando febrero por años bisiestos).

Los resultados de chi-cuadrado no indicaron variación significativa por mes:

# So does the sample match  expected values?
chisq.test(monthDat$suicideCounts, p=monthlyProb)
# Yes, X-squared = 12.7048, df = 11, p-value = 0.3131

La siguiente imagen indica recuentos totales por mes. Las líneas rojas horizontales se colocan en los valores esperados para febrero, 30 días y 31 días, respectivamente. De acuerdo con la prueba de chi-cuadrado, ningún mes está fuera del intervalo de confianza del 95% para los recuentos esperados.

Pensé que había terminado hasta que comencé a investigar datos de series de tiempo. Como imagino que mucha gente lo hace, comencé con el método de descomposición estacional no paramétrico usando la stlfunción en el paquete de estadísticas.

Para crear los datos de series temporales, comencé con los datos mensuales 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

Y luego realizó la stl()descomposición

# Seasonal decomposition
suicideByMonthFit <- stl(suicideByMonthTs, s.window="periodic")
plot(suicideByMonthFit)

En este punto me preocupé porque me parece que hay un componente estacional y una tendencia. Después de mucha investigación en Internet, decidí seguir las instrucciones de Rob Hyndman y George Athanasopoulos tal como se exponen en su texto en línea "Pronósticos: principios y práctica", específicamente para aplicar un modelo ARIMA estacional.

Solía adf.test()y kpss.test()para evaluar la estacionariedad y obtuve resultados contradictorios. Ambos rechazaron la hipótesis nula (señalando que prueban la hipótesis opuesta).

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

Luego usé el algoritmo en el libro para ver si podía determinar la cantidad de diferenciación que se necesitaba hacer tanto para la tendencia como para la temporada. Terminé con nd = 1, ns = 0.

Luego corrí auto.arima, que eligió un modelo que tenía una tendencia y un componente estacional junto con una constante de 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, miré los residuos del ajuste y si entiendo esto correctamente, ya que todos los valores están dentro de los límites del umbral, se comportan como ruido blanco y, por lo tanto, el modelo es bastante razonable. Ejecuté una prueba de portmanteau como se describe en el texto, que tenía un valor p muy por encima de 0.05, pero no estoy seguro de tener los parámetros correctos.

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

Después de volver a leer el capítulo sobre modelado de arima nuevamente, ahora me doy cuenta de que auto.arimaelegí modelar la tendencia y la temporada. Y también me estoy dando cuenta de que el pronóstico no es específicamente el análisis que probablemente debería estar haciendo. Quiero saber si un mes específico (o más generalmente la época del año) debe marcarse como un mes de alto riesgo. Parece que las herramientas en la literatura de pronósticos son muy pertinentes, pero quizás no sean las mejores para mi pregunta. Cualquier y toda entrada es muy apreciada.

Estoy publicando un enlace a un archivo csv que contiene los recuentos diarios. El archivo se ve así:

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

El recuento es la cantidad de suicidios que ocurrieron ese día. "t" es una secuencia numérica del 1 al número total de días en la tabla (5533).

Tomé nota de los comentarios a continuación y pensé en dos cosas relacionadas con el modelado del suicidio y las estaciones. Primero, con respecto a mi pregunta, los meses son simplemente representantes para marcar el cambio de estación, no estoy interesado en si un mes en particular es diferente o no (eso es una pregunta interesante, pero no es lo que me propuse investigar). Por lo tanto, creo que tiene sentido igualar los meses simplemente usando los primeros 28 días de todos los meses. Cuando haces esto, obtienes un ajuste ligeramente peor, lo que estoy interpretando como más evidencia de la falta de estacionalidad. En el resultado a continuación, el primer ajuste es una reproducción de una respuesta a continuación utilizando meses con su número real de días, seguido de un conjunto de datos suicideByShortMonthen el que se calcularon los recuentos de suicidios a partir de los primeros 28 días de todos los meses. Me interesa lo que la gente piense acerca de si este ajuste es una buena idea, no es necesario o dañino.

> 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

La segunda cosa que he investigado más es la cuestión de usar el mes como proxy de la temporada. Quizás un mejor indicador de la temporada es la cantidad de horas de luz que recibe un área. Estos datos provienen de un estado del norte que tiene una variación sustancial en la luz del día. A continuación se muestra un gráfico de la luz del día del año 2002.

Cuando uso estos datos en lugar del mes del año, el efecto sigue siendo significativo, pero el efecto es muy, muy pequeño. La desviación residual es mucho mayor que los modelos anteriores. Si las horas del día son un mejor modelo para las estaciones y el ajuste no es tan bueno, ¿es esto más evidencia de un efecto estacional muy pequeño?

> 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

Estoy publicando las horas del día en caso de que alguien quiera jugar con esto. Tenga en cuenta que este no es un año bisiesto, por lo que si desea incluir los minutos para los años bisiestos, extrapole o recupere los datos.

state.daylight.2002.csv

[ Editar para agregar la trama de la respuesta eliminada (espero que a rnso no le importe que mueva la trama de la respuesta eliminada aquí arriba a la pregunta. Svannoy, si no quieres que esto se agregue después de todo, puedes revertirlo)]

¿Qué pasa con una regresión de Poisson?

Creé un marco de datos que contiene sus datos, más un índice tpara el tiempo (en meses) y una variable monthdayspara la cantidad de días en cada mes.

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

Entonces se ve así:

> 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

Ahora comparemos un modelo con un efecto de tiempo y un efecto de varios días con un modelo en el que agreguemos un efecto de mes:

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

Aquí está el resumen del modelo "pequeño":

> 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

Puede ver que las dos variables tienen efectos marginales en gran medida significativos. Ahora mira el modelo más grande:

> 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

Bueno, por supuesto, el monthdaysefecto se desvanece; ¡solo se puede estimar gracias a los años bisiestos! Mantenerlo en el modelo (y tener en cuenta los años bisiestos) permite utilizar las desviaciones residuales para comparar los dos 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

$\cos\left( {2\pi t \over 12}\right)$$\sin\left( {2\pi t \over 12}\right)$

> 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

Ahora compárelo con el 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

Entonces, seguramente se puede decir que esto sugiere un efecto estacional ...

Una prueba de chi-cuadrado es un buen enfoque como vista preliminar de su pregunta.

La stldescomposición puede ser engañosa como herramienta para evaluar la presencia de estacionalidad. Este procedimiento logra devolver un patrón estacional estable incluso si se pasa como entrada un ruido blanco (señal aleatoria sin estructura). Prueba por ejemplo:

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

Mirar las órdenes elegidas por un procedimiento de selección automática de modelo ARIMA también es un poco arriesgado, ya que un modelo ARIMA estacional no siempre involucra estacionalidad (para más detalles, vea esta discusión ). En este caso, el modelo elegido genera ciclos estacionales y el comentario de @RichardHardy es razonable, sin embargo, se requiere una mayor comprensión para concluir que los suicidios son impulsados ​​por un patrón estacional.

A continuación, resumo algunos resultados basados ​​en el análisis de las series mensuales que publicaste. Este es el componente estacional estimado sobre el modelo de serie temporal estructural básico:

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)

Se extrajo un componente similar utilizando el software TRAMO-SEATS con opciones predeterminadas. Podemos ver que el patrón estacional estimado no es estable en el tiempo y, por lo tanto, no respalda la hipótesis de un patrón recurrente en el número de suicidios a través de los meses durante el período de la muestra. Al ejecutar el software X-13ARIMA-SEATS con opciones predeterminadas, la prueba combinada de estacionalidad concluyó que la estacionalidad identificable no está presente.

Editar (vea esta respuesta y mi comentario a continuación para obtener una definición de estacionalidad identificable ).

Dada la naturaleza de sus datos, valdría la pena complementar este análisis basado en métodos de series de tiempo con un modelo de datos de conteo (por ejemplo, modelo de Poisson) y probar la importancia de la estacionalidad en ese modelo. Agregar los datos de conteo de etiquetas a su pregunta puede generar más vistas y posibles respuestas en esta dirección.

Como se señaló en mi comentario, este es un problema muy interesante. Detectar la estacionalidad no es un ejercicio estadístico solo. Un enfoque razonable sería consultar a la teoría y a expertos como:

• Psicólogo
• Psiquiatra
• Sociólogo

sobre este problema para entender "por qué" habría estacionalidad para complementar el análisis de datos. En cuanto a los datos, utilicé un excelente método de descomposición llamado modelo de componentes no observados (UCM), que es una forma de método de espacio de estado. Vea también este artículo muy accesible de Koopman. Mi enfoque es similar a @Javlacalle. No solo proporciona una herramienta para descomponer datos de series temporales, sino que también evalúa objetivamente la presencia o ausencia de estacionalidad a través de pruebas de significación. No soy un gran admirador de las pruebas de significación en datos no experimentales, pero no conozco ningún otro procedimiento en el que pueda probar su hipótesis sobre la presencia / ausencia de estacionalidad en datos de series de tiempo.

Muchos ignoran pero una característica muy importante que uno quisiera entender es el tipo de estacionalidad:

1. Estocástico: cambios aleatorios y difíciles de predecir
2. Determinista: no cambia, perfectamente predecible. Podrías usar un ficticio o trigonometría (sin / cos, etc.) para modelar

Para datos de series temporales largas como la suya, es posible que la estacionalidad haya cambiado con el tiempo. Una vez más, la UCM es el único enfoque que sé que puede detectar esta estacionalidad estocástica / determinista. UCM puede descomponer su problema en los siguientes "componentes":

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

También puede probar si el nivel, la pendiente, el ciclo es determinista o estocástico. Tenga en cuenta que level + slope = trend. A continuación, presento un análisis de sus datos utilizando UCM. Usé SAS para hacer el análisis.

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;

Después de varias iteraciones considerando diferentes componentes y combinaciones, terminé con un modelo parsimonioso de la siguiente forma:

Hay un nivel estocástico + estacionalidad determinista + algunos valores atípicos y los datos no tienen otras características detectables.

A continuación se muestra el análisis de significancia de varios componentes. Tenga en cuenta que usé trigonometría (que es sin / cos en la declaración de estacionalidad en PROC UCM) similar a @Elvis y @Nick Cox. También podría usar la codificación ficticia en UCM y cuando lo probé, ambos dieron resultados similares. Consulte esta documentación para conocer las diferencias entre las dos formas de modelar la estacionalidad en SAS.

Como se muestra arriba, tiene valores atípicos: dos pulsos y un cambio de nivel en 2009 (¿La economía / la burbuja de la vivienda jugaron un papel después de 2009?), Lo que podría explicarse mediante un análisis más profundo. Una buena característica del uso Proc UCMes que proporciona una excelente salida gráfica.

A continuación se muestra la estacionalidad y una trama combinada de tendencia y estacionalidad. Lo que queda es ruido .

Una prueba de diagnóstico más importante si desea utilizar valores de p y pruebas de significación es verificar si sus residuos no tienen patrón y están distribuidos normalmente, lo que se cumple en el modelo anterior con UCM y como se muestra a continuación en las gráficas de diagnóstico de residuos como acf / pacf y otros.

Conclusión : Basado en el análisis de datos usando UCM y pruebas de significación, los datos parecen tener estacionalidad y vemos un alto número de muertes en los meses de verano de mayo / junio / julio y el más bajo en los meses de invierno de diciembre y febrero.

Consideraciones adicionales : considere también la importancia práctica de la magnitud de la variación estacional. Para negar argumentos contrafácticos, consulte a expertos en dominios para complementar y validar aún más su hipótesis.

De ninguna manera estoy diciendo que este es el único enfoque para resolver este problema. La característica que me gusta de UCM es que le permite modelar explícitamente todas las características de series temporales y también es muy visual.

Para la estimación visual inicial, se puede usar el siguiente gráfico. Al trazar los datos mensuales con la curva de loess y su intervalo de confianza del 95%, parece que hay un aumento a mediados de año que alcanza su punto máximo en junio. Otros factores pueden estar causando que los datos tengan una distribución amplia, por lo tanto, la tendencia estacional puede estar enmascarada en este gráfico de loess de datos sin procesar. Los puntos de datos han sido alterados.

Editar: La siguiente gráfica muestra la curva de loess y el intervalo de confianza para el cambio en el número de casos del número en el mes anterior:

Esto también muestra que durante los meses en la primera mitad del año, el número de casos sigue aumentando, mientras que disminuyen en la segunda mitad del año. Esto también sugiere un pico a mediados de año. Sin embargo, los intervalos de confianza son amplios y abarcan 0, es decir, sin cambios, durante todo el año, lo que indica una falta de significación estadística.

La diferencia del número de un mes se puede comparar con el promedio de los valores de los 3 meses anteriores:

Esto muestra un claro aumento en los números en mayo y una caída en octubre.