Tengo 17 años (1995 a 2011) de datos de certificados de defunción relacionados con muertes por suicidio para un estado de EE.UU. Hay mucha mitología por ahí sobre los suicidios y los meses/temporadas, mucha de ella contradictoria, y de la literatura que he revisado, no tengo una idea clara de los métodos utilizados o la confianza en los resultados.
Así que me he propuesto ver si puedo determinar si es más o menos probable que se produzcan suicidios 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 nos fijamos en el año con menos suicidios, éstos se producen en 309/365 días (85%). Si nos fijamos en el año con más suicidios, se producen en 339/365 días (93%).
Así que hay un buen número de días al año sin suicidios. Sin embargo, cuando se suman los 17 años, hay suicidios todos los días del año, incluido el 29 de febrero (aunque sólo 5 cuando la media es de 38).
La simple suma del número de suicidios en cada día del año no indica una estacionalidad clara (a mi entender).
Agregados a nivel mensual, la media de suicidios por mes oscila entre:
(m=65, sd=7,4, a m=72, sd=11,1)
Mi primera aproximación consistió en agregar el conjunto de datos por meses para todos los años y hacer una prueba de ji-cuadrado tras calcular las probabilidades esperadas para la hipótesis nula, que no había varianza sistemática en los recuentos de suicidios por mes. Calculé las probabilidades para cada mes teniendo en cuenta el número de días (y ajustando febrero para los años bisiestos).
Los resultados del chi-cuadrado no indicaron ninguna variación significativa en función del mes:
# So does the sample match expected values?
chisq.test(monthDat$suicideCounts, p=monthlyProb)
# Yes, X-squared = 12.7048, df = 11, p-value = 0.3131La imagen siguiente indica los recuentos totales por mes. Las líneas rojas horizontales se sitúan en los valores previstos para febrero, los meses de 30 días y los meses de 31 días, respectivamente. De acuerdo con la prueba chi-cuadrado, ningún mes queda fuera del intervalo de confianza del 95% para los recuentos previstos. 
Pensaba que había terminado hasta que empecé a investigar datos de series temporales. Como imagino que hace mucha gente, empecé con el método de descomposición estacional no paramétrica utilizando el stl del paquete stats.
Para crear los datos de las series temporales, empecé 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 71Y 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 tanto un componente estacional como una tendencia. Después de mucho buscar en Internet, decidí seguir las instrucciones de Rob Hyndman y George Athanasopouls que figuran en su texto en línea "Forecasting: principles and practice", concretamente aplicar un modelo ARIMA estacional.
Utilicé adf.test() y kpss.test() para evaluar estacionariedad y obtuve resultados contradictorios. Ambos rechazaron la hipótesis nula (señalando que ponen a prueba la hipótesis contraria).
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.01A continuación, utilicé el algoritmo del libro para ver si podía determinar la cantidad de diferenciación que había que hacer tanto para la tendencia como para la estación. Terminé con nd = 1, ns = 0.
Luego corrí auto.arima que eligió un modelo con un componente tendencial y otro 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.03813434Por último, he mirado los residuos del ajuste y, si lo he entendido bien, como todos los valores están dentro de los límites del umbral, se comportan como ruido blanco y, por tanto, el modelo es bastante razonable. He ejecutado un 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.4817Después de volver a leer el capítulo sobre el modelado arima, ahora me doy cuenta de que auto.arima optó por modelar la tendencia y la temporada. Y también me estoy dando cuenta de que la previsión no es específicamente el análisis que probablemente debería estar haciendo. Quiero saber si un mes concreto (o, más en general, una época del año) debería marcarse como mes de alto riesgo. Parece que las herramientas de la literatura de previsión son muy pertinentes, pero quizás no las mejores para mi pregunta. Agradeceré cualquier aportación.
Estoy publicando un enlace a un archivo csv que contiene los recuentos diarios. El archivo tiene este aspecto:
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 2El recuento es el número de suicidios ocurridos ese día. "t" es una secuencia numérica de 1 al número total de días de la tabla (5533).
He tomado nota de los comentarios de abajo y he pensado en dos cosas relacionadas con el suicidio de modelos y las estaciones. En primer lugar, con respecto a mi pregunta, los meses son simplemente aproximaciones para marcar el cambio de estación, no estoy interesado en si un mes en particular es o no diferente de los demás (por supuesto, es una pregunta interesante, pero no es lo que me propuse investigar). Por lo tanto, creo que tiene sentido igualar los meses utilizando simplemente los 28 primeros días de todos los meses. Al hacer esto, se obtiene un ajuste ligeramente peor, lo que interpreto como una prueba más de la falta de estacionalidad. En la salida de abajo, el primer ajuste es una reproducción de una respuesta de abajo utilizando los meses con su verdadero número de días, seguido de un conjunto de datos suicideByShortMonth en el que los recuentos de suicidios se calcularon a partir de los 28 primeros días de todos los meses. Me interesa saber qué opina la gente sobre si este ajuste es una buena idea, no es necesario o es perjudicial.
> 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: 4La segunda cosa que he investigado más es la cuestión de utilizar el mes como sustituto de la temporada. Quizá un mejor indicador de la temporada sea el número de horas de luz que recibe una zona. Estos datos proceden de un estado septentrional que tiene una variación sustancial de la luz diurna. A continuación se muestra un gráfico de la luz del día a partir del año 2002.
Cuando utilizo 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 en los modelos anteriores. Si las horas de luz son un modelo mejor para las estaciones, y el ajuste no es tan bueno, ¿es esto una prueba más de que el efecto estacional es 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: 4Pongo las horas de luz por si alguien quiere jugar con esto. Ojo, este no es un año bisiesto, así que si quieres poner los minutos de los años bisiestos, extrapola o recupera los datos.
[ Editar para añadir 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 añada después de todo, puedes revertirlo)].
