sábado, 30 de diciembre de 2017

Modelando el Suicidio en Yucatán



Tengo la tarea de escribir una entrada para el blog de Foco Rojo donde hablaré un poco de la predicción del suicidio; y me pareció una buena idea contextualizar un poco el tema en nuestro estado. Lo haré ofreciendo información sencilla por medio de Sistemas de Información Geográfico, así como con un pequeño modelamiento estadístico. Aquí les hablaré un poco de como obtuve este modelo.

Paso 1. Obtener los datos apropiados

Afortunadamente, México dispone de unos de los sistemas estadísticos más completos y accesibles a nivel mundial. Por ejemplo, ¿Sabían que pocas organizaciones de estadísticas nacionales disponen también de información georreferenciada? y de las que lo tienen ¿Que son pocos los países que lo ofrecen de forma totalmente gratuita a unos cuantos clicks de distancia? Pues México es uno de esos pocos países.

Gracias a ello me di a la tarea de buscar la información de Suicidios en Yucatán de 1990 hasta el año anterior (2016). Para acceder a ello basta que sigan la ruta "Estadística > Fuente/Proyecto > Registros Administrativos > Estadísticas vitales > Mortalidad > Tabulados > Mortalidad General". Ahí podrán filtrar las estadísticas de mortalidad por estado, municipio, causa, sexo, edad, año, y muchísimas variables más (Solo denle gracias a la Secretaría de Salud [SS] por ser tan organizada con la información que provee al INEGI). En mi caso, descargue en formato de Excel las "Defunciones accidentales y violentas", por estado y por año. Eligí la causa "Suicidio", las cuales corresponden a lo reportado por la SS como "Muertes por lesiones autoinflingidas intencionalmente" que según el CIE-10 corresponden a las categorías X60 a X84.

Y listo, ¿No?... Bueno, no tan fácil. Verán, obviamente de 1990 a 2016 habrá más suicidios, porque bueno, la población crece. Entonces, quizás deberíamos en ves de trabajar con valores absolutos trabajar con datos que estén escalados de alguna forma. Pensé primero trabajar con las proyecciones de población de CONEVAL, pero dado que son proyecciones, mejor preferí trabajar con otro indicador más realista. En mi caso, elegí trabajar con "Mortalidad general", es decir, todas las causas de muerte, y a partir de ello, dividí los suicidios de cada año, entre la mortalidad general de cada año, para obtener una especie de proporción. Dado que los valores eran muy pequeños, decidí reescalarlos por cada 100,000 muertes. Ahora sí ¡Listo!... Ya tenía la tasa de suicidios por cada 100,000 muertes en Yucatán de Enero de 1990 a Diciembre de 2016.

Paso 2. Cargar la serie de tiempo... en R

La verdad es que estoy un poco oxidado en análisis de series de tiempo, y en R más. Así que en parte esta entrada es para no solo comentarles lo que hice, sino también para recordarme en futuras ocasiones que rayos hice.

Dicho lo anterior, lo primero era cargar la base de datos en R (Base = SuicidiosYuc); y ahí vino la primera situación importante. Si cargáramos una base de datos a R, sin importar cual sea, R lo detectará como eso, como un "data frame", y en mi caso, si bien era una base de datos, dado los análisis que quería realizar, tenía que hacer que R lo identificará como lo que es, una "serie de tiempo".

La realidad es que resolver lo anterior fue más fácil de lo que creí. Basta usar el comando ts(base, freq=h start=c(año, mes)), que volverá cualquier data frame en una serie de tiempo. En mi caso, con el comando:
Suicidio <- ts(SuicidioYuc,freq=12,start=c(1990,1))
Creé una nueva base (ahora en formato Serie de Tiempo) llamado Suicidio, a partir de la base original "SuicidioYuc", con una frecuencia de 12 (porque son 12 meses), y cuya temporalidad inicia el mes de Enero (1) del año 1990.

¿Funcionó?... ¡Claro! Para asegurarme hice su respectivo gráfico con el comando:
plot(Suicidio,main="Tasa de Suicidios en Yucatán por cada 100,000 defunciones", xlab="Año", ylab="Tasa de suicidios")
El cual pueden ver justo a continuación:

Paso 3. Analizando los componentes de la Serie de Tiempo


La siguiente tarea, es analizar si podemos modelar los datos en R. Para ello tenemos que conocer sus componentes. Pero déjenme les explico un poco de como funcionan las series de tiempo para que comprendan adecuadamente lo que hice. Pondré una pequeña notación matemática, pero prometo que será sencillo y les aliento a leerlo para que lo entiendan.

Primero, imaginemos que tenemos una serie de datos Yt. Nosotros podríamos modelar la serie de datos Yt si consideráramos el siguiente modelo estadístico:

Yt = Tt + St + et

Donde, Tt, es la tendencia. Es decir, ¿los datos incrementan con el tiempo? ¿Decrecen? ¿O se mantienen sin tendencia? (esto últilmo sería lo que llamamos una serie "estacionaria", porque los datos no cambian, no tienen tendencia). St por otro lado, es la "estacionalidad" (no lo confundan con "estacionaria"). Es decir, ¿La serie tiene ciclos o estaciones? (La S viene de Season, "estación" en inglés). Y e simplemente es un error, un modelo no puede ser perfecto, siempre conlleva a un elemento de error, que en el mejor de los casos deseamos que sea un error completamente azaroso, aleatorio. Ejemplo: imaginen que con todos sus ahorros, abren una tienda, y analizan sus ganacias. La tendencia sería lo que esperaríamos ver en el caso de que su tienda tenga éxito. Es decir, sus ganancias irían aumentando poco a poco conforme la gente vaya conociendo más su tienda. Los ciclos o estaciones serían, por ejemplo, cada quincena, pues resulta que la gente va más esos días y compra más cosas. Entonces, verían que sus ventas tienen un ciclo quincenal. Y el error, bueno, el comportamiento de los clientes seguro tiene algun componente aleatorio, no podemos predecirlo con total exactitud.

Dicho lo anterior, me di a la tarea de analizar los componentes de la serie de tiempo de suicidios en Yucatán que va de 1990 a 2016. Con los siguientes comando pude hacer, primero, la descomposición de la serie en estos 3 elementos; y segundo, graficarlo.
descomp1=decompose(Suicidio)
plot(descomp1)
El resultado, es el siguiente gráfico:



Aquí podemos apreciar, en el gráfico de arriba, la serie original, es decir, la tasa de suicidios por cada 100,000 defunciones en función del tiempo (por eso lleva la etiqueta "original"). Justo debajo, vemos el gráfico que dice "Tendencia". Ahí podemos observar como la tasa de suicidios tiene una tendencia creciente, por lo tanto, cuando lo modelemos, será conveniente incluir esta tendencia. En tercer lugar, tenemos el componente de ciclos, donde parece que sí existen ciclos, y dado que son 27 (como los años que analizamos), seguramente es un ciclo anual. Por lo tanto, habría que considerar esta información al modelarlo. Finalmente, tenemos el componente aleatorio, o el error, donde no observamos algún patrón en especial, por lo que seguramente podremos modelar los datos sin preocuparnos por este componente. Por cierto, este patrón, que no tiene ningún patrón, se le conoce como caminata aleatoria.

Paso 4. Crear una serie de tiempo para modelar y otra para probar

El siguiente paso pues, es crear dos "sub-series" de tiempo. Uno que nos servirá para entrenar el modelo, y otra para probarlo. En mi caso cree la serie de entrenamiento con los datos de 1990 hasta 2015 (Serie: "Sui.train"), y utilice los datos de 2016 para probar el modelo (Serie: "Sui.test"). Lo hice con los siguientes comandos, donde solo tengo que especificar la base serie (Suicidio), y el inicio y final de las sub-series.

Sui.train <- window(Suicidio, start = c(1990,1), end = c(2015,12))
Sui.test <- window(Suicidio, start = c(2016,1), end = c(2016,12))

Paso 5. Modelando la serie de tiempo de forma automática


Hasta ahora no habíamos necesitado algún paquete adicional en R, pero ahora necesitaremos los paquetes "forecast" para hacer predicciones; "tseries" para correr algunas pruebas de hipótesis; y "TSPred" para graficar las predicciones. Así que sí están replicándo lo que hice, bájenlos, y cárguenlos en su memoria de R.

Pensé, tontamente como verán posteriormente, que podría modelarlo fácilmente con el modelizador automático del paquete forecast. Con el comando auto.arima(), simplemente hacemos que R automáticamente prueba varios modelos ARIMA. No voy a adentrarme a lo que es un modelo ARIMA (Modelo Autorregresivo Integrado de Promedio Móviles o AutoRegressive Integrated Moving Average Model), pero básicamente yo le señalo una serie de parámetros para que considere la tendencia y ciclos. De hecho, en realidad al final use un modelo SARIMA, es que lo mismo pero con un componente de ciclos (Seasonal ARIMA). La forma de correr el modelizador automático que seguí fue con el siguiente código:

arima1 <- auto.arima(Sui.train, trace = TRUE, test = "kpss", ic = "aic")
El modelo automático que me dió R, fue un modelo ARIMA (1,1,1). Que básicamente no tomaba en consideración los ciclos. Por lo tanto, no fue extraño que la prueba de Ljung-Box saliera bastante extraña. A grandes rasgos, esta prueba lo que analiza es si los errores se autocorrelacionan. Y si recordamos un poco, los errores son el componente aleatorio de la serie, por lo tanto, no deberían correlacionarse de ninguna forma. Pero ¡Sorpresa!, la prueba de autocorrelación de LJung-Box, tuvo una Ji cuadrada de 63.43, con 36 grados de libertad, y un valor p = 0.003185. Es decir, la autocorrelación fue significativa. Estoy casi seguro que esta autocorrelación se debe a que el modelo automático no está tomando en consideración los ciclos.

Solo para propósitos de seguir ilustrando el proceso (con este modelo mal hecho), le pedí que me hiciera 12 predicciones con el comando:
arima1.forescast <- forecast(arima1, h=12)
Y luego le pedí que me graficara esas predicciones contrastándolas con los datos originales con el paquete y comando siguiente:
library(TSPred)
plotarimapred(Sui.test, arima1, xlim = c(2010,2017), range.percent = 0.05)
Que simplemente le pedí, que con la base de prueba (Sui.test), contrastara el modelo que hicimos  (arima1), que me acote el gráfico de 2010 a 2017 (xlim = c(2010,2017)), y que me incluya los intervalos de confianza del 95% para las predicciones (range.percent = 0.05). El resultado, pueden apreciarlo justo abajo:



El estimado, que es la línea azul, como pueden ver, no es muy bueno, pues no se ajusta del todo a los datos reales de 2016 (la línea negra punteada). Como ya mencioné, seguramente es porque no está considerando los ciclos.

Paso 6. Modelo final

Intenté por varios métodos con el modelizador automático modelar la serie de tiempo, pero en ninguna ocasión consideró los ciclos. Incluso a pesar de que por default corre modelos con ciclos. En fin, la verdad es que me decidí hacer el modelo "a mano", y tras probar con diversos parámetros, el modelo que tuvo un mejor ajuste de todos, fue un modelo SARIMA (1,1,1)(1,0,2)[12]; el cual corrí con el siguiente comando:
arima3 <- arima(Sui.train, order = c(1,1,1), seasonal = c(1,0,2))
Nuevamente probe el supuesto de autocorrelación de los errores con los siguientes comando, el primero para ver el gráfico de autocorrelación, y el segundo para la prueba de Ljung-Box de autocorrelación:
#Gráfico de la función de autocorrelación
acf(arima3$residuals, lag.max = 36, type = "correlation", plot = TRUE)
#Test Ljung-Box para autocorrelación lineal
Box.test(arima3$residuals, lag=36, type = "Ljung-Box")
El resultado podemos apreciarlo en el gráfico de a continuación. Y, en general, ningún coeficiente de correlación fue lo suficientemente grande para pensar que la autocorrelación será significativa.



Lo anterior lo comprobamos con la prueba de Ljun-Box, que me dio un valor de Ji-cuadrada de 48.36(36), con un valor p = 0.08. Lo que significa que no existe autocorrelación... ¡Bravo!

El paso siguiente, fue realizar una predicción utilizando mi nuevo modelo (arima3). Y posteriormente graficarlo contra los datos de prueba, es decir, los verdaderos datos de 2016. Todo ello con los siguientes comandos:
arima3.forescast <- forecast(arima3, h=12)
plotarimapred(Sui.test, arima3, xlim = c(2010,2017), range.percent = 0.05, main = "Predicción de ARIMA(1,1,1)(1,0,2)[12]", xlab = "Año", ylab = "Tasa de Suicidios")
El resultado, puede apreciarse a continuación. Y tal como se mira, este modelo es mucho mejor que el anterior.


Y aunque podemos ver un repunte que sobresale del estimado para 2016, lo cierto es que en general, pareciera que el modelo se ajusta relativamente bien a los datos. En un análisis más sencillo, me tomé la libertad de graficar la distribución de la tasa de suicidios por cada mes por medio de un gráfico de violín en R, así como añadir un gráfico de caja y bigotes a cada distribución. Esto lo hice con el siguiente comando, utilizando el paquete ggplot2:

library(ggplot2)ggplot(data = SUICIDIO, aes(Mes, Tasa, fill = Mes)) + geom_violin(trim=FALSE) +  geom_boxplot(width=0.2, fill="white") +  labs(title="Tasa promedio de suicidio en Yucatán por Mes (1990-2016)",       x="Mes", y = "Tasa") + theme_minimal() + theme(axis.text.x=element_blank())




En este resultado, podemos apreciar de mejor forma los ciclos de la tasa de suicidios. La línea media dentro de cada gráfico de caja y bigotes, nos representa la mediana, y vemos por ejemplo, que en el verano la tasa de suicidios aumenta, y que los meses que más suicidios presentan en promedio son abril y agosto. Además, se puede apreciar como junio y diciembre, son los meses que presentan una mayor y menor variabilidad respectivamente.

Conclusión

En general, observamos una tendencia creciente en los suicidios en Yucatán. Incluso si lo consideramos de forma relativa como lo hice, pues modelé la tasa, no las frecuencias absolutas. Y aunque la tendencia es pequeña, ciertamente hay que considerar ello para la política pública de Yucatán. Segundo, la presencia de ciclos anuales en el modelo indica que podemos analizar estas tendencias y ver sus patrones mensuales. Por ejemplo, en la imagen superior, podemos observar como claramente, los meses donde ocurren más suicidios en promedio es en abril, y agosto (¿temporada vacacional?), y de hecho podemos observar como en general, desde mediados de primavera hasta finales del verano, la tasa de suicidios aumenta. Sin duda es otro elemento más importante para tomar en cuenta en los programas de prevención. ¿Quizás estar más alerta de las personas en riesgo en esos meses?

Aún nos falta mucho por hacer para prevenir el suicidio en nuestro estado, pero un primer paso es conocer mejor este fenómeno. Si bien este fue un ejercicio sencillo y con la intención únicamente de contextualizar la situación del suicidio en nuestro estado, lo cierto es que nos brindó información muy interesante.

Les invito a analizar más datos de forma similar, y si pueden ajustar un mejor modelo que el que yo realicé, ¡Comenten como llegaron a ello!

0 comentarios:

Publicar un comentario