Modelos predictivos para la distribución del número de casos diarios de la COVID-19 para los diferentes distritos de Barcelona

Grado en Estadística Título: Modelos predictivos para la distribución del número de casos diarios de la COVID-19 para los diferentes distritos de Barcelona Autor: Pablo Pemán Rodríguez Director: Dr. David Moriña Soler Departamento: Econometria, Estadística y Economía Aplicada Convocatoria: 2n 2020-2021 Resumen La COVID-19 ha afectado a toda la sociedad, es por ello por lo que hemos considerado necesario ajustar modelos que nos permitan explicar y predecir la serie temporal de los casos detectados. Al ser la COVID-19 una enfermedad infecciosa, donde los casos de un día vienen influenciados por los días anteriores, consideramos que los mejores modelos para explicar la serie serán los autorregresivos, en concreto, los AR y los INAR. Analizaremos la serie en función del sexo y el distrito. A nivel de resultados obtenemos diferentes mecanismos en ambos modelos en función del sexo, y también vemos como los barrios con mayor población o densidad poblacional son aquellos con más tendencia a detectar más casos. Comparando los modelos, se observa claramente que el modelo INAR es superior en todos los aspectos al modelo AR, teniendo en cuenta que nuestra serie no es estacionaria, y, por tanto, ninguno es perfecto. Palabras clave: COVID-19, Modelización de enfermedades infecciosas, Análisis estocástico de series temporales, Series temporales de valores enteros Abstract COVID-19 has affected the whole of society, which is why we have considered it necessary to adjust models that allow us to explain and predict the time series of the detected cases. As COVID-19 is an infectious disease, where the cases of one day are influenced by the previous days, we consider that the best models to explain the series will be the autoregressive ones, specifically, the AR and the INAR. We will analyze the series as a function of sex and district. At the level of results, we obtain different mechanisms in both models as a function of sex, and we also see how the districts with higher population or population density are those with a greater tendency to detect more cases. Comparing the models, it is clearly observed that the INAR model is superior in all aspects to the AR model, taking into account that our series is not stationary, and, therefore, neither is perfect. Key words: COVID-19, Infectious disease modeling, Stochastic time series analysis, Integer valued time series. Clasificación AMS: • 37M10 Time series analysis • 62M10 Time series, auto-correlation, regression, etc. • 62M20 Prediction • 62P10 Applications to biology and medical sciences Agradecimientos Primero de todo, agradecer a mi tutor el Dr. David Moriña Soler por haber aceptado ser mi tutor en este trabajo, y por guiarme y ayudarme con todas las dudas que me encontré en la realización de este. Doy las gracias a los profesores del grado de Estadística por hacerme crecer en mi la pasión por el análisis de datos. Finalmente, doy las gracias a mi familia y amigos, por siempre estar a mi lado y apoyarme en los momentos difíciles, sin ellos no estaría donde me encuentro hoy en día. ÍNDICE I. INTRODUCCIÓN ......................................................................................................................... 1 I.1. COVID-19 ................................................................................................................................... 1 I.2. Análisis de series temporales ..................................................................................................... 3 I.3. Objetivos del trabajo .................................................................................................................. 5 I.4. Hipótesis ....................................................................................................................................... 6 II. METODOLOGÍA ..................................................................................................................... 7 II.1. Base de datos .............................................................................................................................. 7 II.2. Modelos AR ................................................................................................................................ 8 II.3. Modelos INAR ......................................................................................................................... 10 II.3. Indicadores para las predicciones .......................................................................................... 12 III. RESULTADOS ........................................................................................................................ 13 III.1. Descriptiva de los Datos ........................................................................................................ 13 III.2. Identificación de los modelos ................................................................................................ 16 III.3. Modelos AR ............................................................................................................................ 19 III.3.1 Estimación de los parámetros ............................................................................................ 19 III.3.2. Validación del modelo ....................................................................................................... 20 III.3.3 Predicciones ....................................................................................................................... 21 III.4. Modelos INAR ........................................................................................................................ 26 III.4.1 Estimación de los parámetros ............................................................................................ 26 III.4.2 Predicciones ....................................................................................................................... 28 III.5. Comparación de los modelos ................................................................................................ 33 IV. CONCLUSIONES ................................................................................................................... 35 V. BIBLIOGRAFIA ......................................................................................................................... 37 VI. ANNEXOS ............................................................................................................................... 38 1 I. INTRODUCCIÓN I.1. COVID-19 La COVID-19 es la enfermedad infecciosa causada por el coronavirus que se ha descubierto más recientemente (SARS-CoV2). Tanto este nuevo virus, como la enfermedad que provoca, eran desconocidos antes de que estallara el brote en Wuhan (China) en diciembre de 2019. Su rápida propagación, su transmisión aérea y la necesidad del mínimo contacto para la infección, añadido a la posibilidad de no manifestarse, o retrasar su manifestación, provocó que en marzo de 2020 la COVID-19 se hubiera propagado por todo el mundo y la Organización Mundial de la Salud (OMS) lo clasificara como pandemia. [1] La COVID-19 ha afectado a todo el mundo, pero sobre todo ha tenido mayor violencia en los países y zonas más desarrolladas. Wuhan es el centro político, cultural, social y económico de China. En Estados Unidos los estados que más sufrieron fueron los de Nueva York, Nueva Jersey, Michigan o California. En España, la Comunidad Autónoma de Madrid y Cataluña sufrieron los mayores ataques. Esto se explica a que todas estas zonas acumulan grandes concentraciones de personas, ya sean residentes o turistas. Todas estas aglomeraciones facilitan la propagación del virus, no solo a nivel regional, sino también mundial. La COVID-19 no solo ha afectado a nivel sanitario, provocando colapsos de los hospitales y millones de muertes alrededor del mundo. También ha afectado a la sociedad y la economía a nivel mundial. Ilustración 1. Tasa interanual del IPI en España A nivel español, podemos comprobar como la producción industrial llego a mínimos históricos, superando incluso los datos correspondientes a la crisis económica que sufrió el país en 2009 (Ilustración 1). 2 Ilustración 2. Variación anual de la compraventa de viviendas en España También podemos comprobar como el número de compraventas de viviendas ha bajado muchísimo, superando también a los valores mínimos que hubo en España durante la crisis económica de 2009 y los posteriores años (Ilustración 2). Ilustración 3. Evolución del paro en la ciudad de Barcelona. Fuente: https://datosmacro.expansion.com/paro/espana/municipios/cataluna/barcelona/barcelona En el gráfico superior podemos ver la afectación social que ha tenido la COVID-19 en la ciudad de Barcelona. Observamos como el paro ha aumentado a partir del primer trimestre de 2020, inicio de la pandemia, y sigue subiendo. (Ilustración 3) 3 En conclusión, la pandemia provocada por la COVID-19 ha afectado negativamente al mundo actual tal y como lo conocemos, y aunque los gobiernos estén intentando paliar los efectos adversos, está claro que la COVID-19 cambiará para siempre la sociedad y la forma de vida actual. Es por ello por lo que creemos que conseguir estimar un modelo capaz de explicar y predecir el comportamiento de la propagación del virus es vital para poder llevar a cabo las medidas preventivas y estrategias necesarias para mitigar al mínimo los daños, tanto sociales, económicos o humanos, que pueda provocar. I.2. Análisis de series temporales Las series temporales se definen como una sucesión de valores de una variable aleatoria ordenados cronológicamente en el tiempo. Su expresión es la siguiente: 𝑥1,𝑥2,𝑥3,…,𝑥t donde t representa la posición temporal que ocupa el valor dentro de la serie. Las series temporales pueden constar de varios componentes, los cuales son importantes para su posterior análisis. Dichos componentes son los siguientes: • Tendencia: Componente de la serie temporal que representa la evolución a largo plazo de la serie, ya sea de forma ascendente o decreciente. • Componente cíclico: Recoge los movimientos oscilatorios por encima y por debajo de la tendencia que se mantienen en periodos superiores al año. • Componente estacional: Recoge los movimientos oscilatorios que se repiten en periodos menores a un año y siempre en el mismo momento temporal. • Componente irregular: Recoge las variaciones esporádicas de la serie, que no están incluidas en ninguno de los componentes estacionales, y, por tanto, aquellas con un carácter residual. 4 Ilustración 4. Número acumulado de casos de COVID-19 en la ciudad de Barcelona El gráfico superior (Ilustración 4) hace referencia a la serie del número de casos de COVID-19 acumulados en la ciudad de Barcelona, y es un muy buen ejemplo para ver el componente de la tendencia en una serie temporal. Además, usaremos este gráfico para introducir un nuevo término, la estacionariedad. Una serie es estacionaria cuando la media y la varianza son constantes a lo largo de toda la serie, y las autocovarianzas entre dos variables aleatorias solo dependen de la distancia temporal que las separa. En nuestro gráfico (Ilustración 4), podemos ver como claramente, al tener una tendencia, esto hace que la media no sea constante a lo largo de la serie. Un componente estacional fuerte podría provocar una serie heterocedástica, es decir, con una varianza no constante. Para paliar estos problemas, normalmente es suficiente con aplicar el logaritmo, o una diferencia estacional, respectivamente. Tenemos dos formas básicas para el análisis de las series temporales dentro del análisis univariante: el análisis determinista y el análisis estocástico. El análisis determinista se basa en gran parte en la descomposición de la serie temporal en todos los componentes explicados anteriormente. En función de la presencia o no de tendencia o estacionalidad, las series se pueden clasificar de la siguiente forma: Tabla 1.1 Tipos de series en función de los componentes presentes Tendencia Si No Estacionalidad Si Tipo IV Tipo II No Tipo III Tipo I 5 La presencia o no de estos componentes puede comprobarse, además de por la observación de la propia serie, por diferentes estadísticos de contraste, en concreto el contraste de Daniel para la tendencia, y el contraste de Kruskal-Wallis para la estacionalidad. Cada tipo de serie tiene diferentes métodos para predecir su comportamiento en función de los componentes entre los que se puede descomponer. Algunos de esto métodos son: método de las medias móviles, método de la tendencia lineal, método ingenuo estacional, método de descomposición …. El problema de este análisis es cuando la descomposición de la serie en alguno de estos componentes no describe con precisión la propia serie. En estas situaciones es cuando el análisis estocástico entra en juego. El análisis estocástico no intenta descomponer la serie, sino que se centra en las autocorrelaciones de los valores. Como bien sabemos, el valor de una serie temporal puede venir influenciado por los valores anteriores de la serie. Por ejemplo, el precio de una acción de bolsa no solo viene influenciado por las acciones que hayan podido ocurrir ese día, sino también por los sucesos de días anteriores. Este tipo de modelos, donde un valor temporal concreto viene influenciado por sus valores pasados se conoce como modelo autorregresivo, AR. Hay otros tipos de modelos, en función de cómo sea la serie y sus autocorrelaciones. Si el valor no viene influenciado por los valores pasados se conoce como modelo de medias móviles MA, y si es mixto se conoce como modelo ARMA. Todos estos modelos comentados son modelos para series continuas, en cambio, para las series que siguen una distribución discreta tenemos los mismos modelos, pero adaptados para utilizar una distribución marginal discreta. Dichos modelos son: modelo autorregresivo para valores enteros no negativos, INAR, modelo de medias móviles para valores enteros no negativos, INMA…. La razón para usar estos modelos, los cuales son la contraparte discreta de los modelos continuos explicados anteriormente, es que dichos modelos utilizan la distribución normal como distribución marginal, y nuestra serie al ser un simple contaje, con valores pequeños, es muy difícil que sea explicada por una distribución marginal normal, y es por ello por lo que también utilizamos los modelos con distribución marginal discreta. La gran diferencia es que los modelos discretos introducen un operador de refinamiento, el cual no está presente en los modelos continuos, y que nos permite asegurar que las estimaciones del modelo serán valores enteros. En concreto, nos centraremos en los modelos AR e INAR, cuya diferencia principal, como hemos comentado anteriormente, es el tipo de series para los que se utilizan. El modelo AR hace referencia a series continuas, mientras que el modelo INAR hace referencia a series discretas. Para una mayor explicación de ambos modelos véase el capítulo II. METODOLOGÍA. I.3. Objetivos del trabajo El análisis temporal de los casos de la COVID-19 en Barcelona, en función del distrito analizado y el sexo, se pretende llevar a cabo mediante la parametrización de modelos capaces de explicar la serie, ajustarse correctamente a la volatilidad de los posibles valores, y predecir con éxito su evolución. Es por ello por lo que los objetivos de este trabajo son: 12 ∑𝑃(𝑋𝑡+𝑘 =𝑗) 𝑙1 𝑗=0 ≈α/2 ,∑𝑃(𝑋𝑡+𝑘 =𝑗) 𝑙2 𝑗=𝑙1≈1 − α En estas ecuaciones vemos cómo vamos sumando las probabilidades de que las predicciones tomen el valor j, hasta que llegamos a la probabilidad α/2. Dicha j, será nuestro límite inferior, a partir del cual calcularemos el límite superior de la misma manera. Sumaremos probabilidades hasta alcanzar la probabilidad 1− α. Con estas fórmulas, también podremos calcular la mediana de las predicciones. La única variación es la de sumar probabilidades hasta llegar a la probabilidad 0.50. ∑𝑃(𝑋𝑡+𝑘 =𝑗) 𝑙1 𝑗=0 ≈0.50 Los apartados correspondientes a la estimación de los parámetros, el ajuste del modelo y las predicciones se podrán consultar a continuación en el capítulo III II.3. Indicadores para las predicciones La eficacia de las predicciones se puede analizar de varias formas. Aquí explicaremos los 3 indicadores que utilizaremos para determinar si las predicciones del modelo se ajustan correctamente a la serie o no. Dichos indicadores serán: los errores de predicción, la cobertura, y el rango de las regiones de predicción. Los errores que estudiaremos son; el error cuadrático medio (EQM), el error absoluto medio (EAM), y el error porcentual absoluto (EPAM). • EQM = 1 𝑘∑(𝑋𝑇 (𝑗)−𝑋𝑡+𝑗)2 𝑘 𝑗=1 • EAM = 1 𝑘∑|𝑋𝑇 (𝑗)−𝑋𝑡+𝑗| 𝑘 𝑗=1 (9) • EPAM = 1 𝑘∑|𝑋𝑇 (𝑗)−𝑋𝑡+𝑗| 𝑋𝑡+𝑗 𝑘 𝑗=1 Donde 𝑋𝑇 (𝑗) hace referencia a la predicción t+j, y k hace referencia al número de predicciones calculadas. En relación con la cobertura, simplemente hemos contado las observaciones entre enero y marzo de 2021 que se encuentran dentro de las regiones de predicción calculadas para dicho tiempo muestral. Para calcular el rango de las regiones de predicción, hemos calculado el rango de cada predicción, restando el límite inferior al límite superior, y una vez hemos calculado el rango para cada predicción hemos calculado la media de todas ellas. El rango lo podemos utilizar para analizar la precisión de los modelos a la hora de predecir. 13 III. RESULTADOS III.1. Descriptiva de los Datos Nuestra base de datos está dividida por los diversos distritos de Barcelona, y en este apartado procederemos a comentar y analizar cada distrito por separado. Primero de todo observaremos de manera gráfica cómo se comportan las series. Ilustración 5. Representación gráfica del número de casos detectados en los distritos de Barcelona En todos los distritos podemos ver como el virus se ha comportado de manera similar, con una ola inicial alrededor de marzo-abril de 2020, inicio de la pandemia. Posteriormente, a finales de julio de 2020 podemos ver como se inicia una segunda ola, que durará hasta inicios de diciembre, y alcanza su pico, y el pico máximo de la serie, en octubre-noviembre. Finalmente, entre mediados y finales de diciembre empieza una tercera ola, muy similar a la segunda e incluso sobrepasándola en algunos distritos, Les Corts o Sarrià-Sant Gervasi, la cual perdura hasta marzo de 2021, fecha final de nuestro estudio. (Para ver las representaciones individuales de cada distrito ir a los Anexos) A nivel de sexos vemos que la COVID-19 se comporta más o menos igual en hombres como en mujeres, con alguna diferencia puntual, por ejemplo, el pico de casos en hombres de Les Corts alrededor de mayo 2020. Esto tiene sentido ya que no se han detectado evidencias para afirmar que la COVID-19 se comporta distinto en función del sexo. Las posibles diferencias que se observen se pueden deber más aun tema social o psicológico, que hace que un sexo respete mejor las medidas que otro, o tenga más propensión a infectarse. Pero aquí no podemos afirmar ni siquiera que estas diferencias sean significativas. Más adelante cuando estimemos los parámetros del modelo podremos ahondar en esta hipótesis. Considero importante explicar, que, aunque la primera ola sea la menor en cuanto a número de casos detectados, fue la más mortal. Este bajo número de casos detectados se puede explicar fácilmente. Al inicio de la pandemia no había métodos necesarios para dicha detección, y, además, los hospitales estuvieron colapsados debido al enorme número de gente infectada, lo 14 que provocó que no se pudiera atender a todo el mundo, en especial, aquellos que no mostraban síntomas claros del virus. Es por ello por lo que el número de casos es inferior a las demás olas, las cuales sucedieron cuando había métodos eficaces de detección, y dichos métodos estaban disponibles para todos aquellos sospechosos de estar infectados. (Para ver una proyección de la mortalidad ir Ilustración Anexos.1) Ahora procederemos a analizar la serie desde un punto de vista estadístico. Como hemos ido comentado anteriormente, para realizar un buen análisis de la serie y poder estimar los modelos que la explican, necesitamos que nuestras series sean estacionarias, para poder analizar de una forma correcta los coeficientes de autocorrelación. La estacionariedad viene dada por una varianza y media constantes. Por desgracia, en nuestros gráficos podemos ver como los casos van aumentado a medida que se detectan más, una pequeña tendencia, y, por tanto, podemos afirmar con claridad que nuestra serie no estacionaria con relación a la media. A nivel de varianza no encontramos muchos problemas. El mayor problema de la varianza es la estacionalidad ya que produciría valores anómalos, lejanos a la media, sin que haya aumentado o disminuido el número de casos. Un ejemplo claro sería una serie temporal con el número de casos por la gripe. Durante todo el año el número de casos sería similar, pero al llegar el invierno aumentará de forma drástica. En nuestro caso vemos como el número de casos se mantiene similar a sus anteriores. Una vez hemos visto la representación gráfica de las series, y hemos determinado que no existe estacionariedad, vamos a analizar más en profundidad cada distrito. Para ello calcularemos el mínimo y máximo de contagios en un día, la media de contagios por día y su desviación estándar, el número de casos totales en el distrito desde el inicio hasta el final del estudio (1/03/2020 – 3/03/2021), y finalmente calcularemos la incidencia de casos por 10.000 habitantes, y relacionaremos todos estos valores con la población total del distrito, su densidad y la renta. En términos generales, podemos observar como el distrito que ha sufrido un mayor impacto por la COVID-19 ha sido el distrito de l’Eixample, superior en todas las categorías a los demás distritos. Hay que tener en cuenta que también es el distrito con mayor población y densidad, y, por tanto, tiene sentido que sea el distrito con mayor afectación. En términos generales Media Desviación estándard Máximo Casos Diarios Número de Casos Totales Casos por 10.000 habitantes Población Densidad Poblacional Renta familiar en miles de € CIUTAT VELLA 11,03 11,20 67 8.116 749,19 108.331 23.908 84.3 EIXAMPLE 26,20 24,38 132 19.282 712,32 270.694 35.556 122.4 GRÀCIA 10,49 10,73 60 7.722 624,50 123.651 29.130 105.3 HORTA-GUINARDÓ 18,33 16,70 90 13.493 771,92 174.799 14.155 78.0. LES CORTS 7,54 9,23 68 5.551 671,75 82.635 13.697 137.3 NOU BARRIS 20,76 19,61 114 15.280 878,10 174.012 20.820 55.0 SANT ANDREU 15,73 14,61 75 11.574 761,57 151.976 22.459 74.6 SANT MARTÍ 22,34 19,94 112 16.443 680,33 241.691 22.129 88.1 SANTS-MONTJUÏC 19,88 18,62 94 14.629 779,86 187.584 8.024 84.6 SARRIÀ-SANT GERVASI 13,35 13,82 71 9.828 650,18 151.157 7.235 182.8 Tabla 3.1 Estadística descriptiva de los casos detectados por distrito 15 podemos ver como el número de casos totales se encuentra en torno a los 8.000-14.000 casos. Dentro de este intervalo encontramos la mayoría de los distritos como: Ciutat Vella, HortaGuinardó, Sant Andreu, Sants-Montjuïc y SarriaSant Gervasi. Sorprendentemente, Les Corts es el distrito, no solo con menor población, sino también con una menor afectación, ya que el total de casos detectados es de 5.551, bastante más bajo que en el siguiente en la lista, Gràcia con 7.772. Por encima de la media no solo encontramos l’Eixample, si no también encontramos Nou Barris y Sant Martí, con 15.280 y 16.443, respectivamente. (Para ver tablas específicas en función del sexo mirar tablas Anexos 3 y 4). A nivel socioeconómico, nuestra hipótesis es que aquellos distritos con una renta familiar superior tendrán menos casos. No solo porque en teoría la población ha de ser menor, ya que no todo el mundo se puede permitir las viviendas en dichos distritos, si no también, porque a medida que aumenta la renta las viviendas tienden a ser más espaciosas, y las familias no se encuentran con la necesidad de compartir residencia con otras burbujas familiares, o incluso alquilar habitaciones. Si observamos los distritos con mayor renta, estos son: Sarrià-Sant Gervasi, Les Corts y L’Eixample. Tanto Les Corts como Sarrià son los distritos con mayores rentas, y los distritos con menos casos, pero también son los distritos con menor población y densidad. Si miramos L’Eixample, aun siendo el tercer distrito con mayor renta, es el distrito con mayor número de casos, población y densidad poblacional. Estas desavenencias en el distrito de l’Eixample se pueden explicar gracias a los barrios que forman l’Eixample. Los barrios de l’Esquerra de l’Eixample, y, sobre todo, el barrio de La Dreta de L’Eixample, donde encontramos Paseo de Gràcia, son barrios de alto nivel adquisitivo, que hacen que la renta del distrito total aumente, aunque en los demás distritos la renta no sea tan alta. Es por ello por lo que, aunque tengamos una renta alta, el número de casos es más elevado de lo que se podía suponer. Aún no podemos afirmar que la renta realmente afecte al número de contagios, porque los distritos que hemos comprobado distaban bastante en cuanto a niveles de población y densidad poblacional. Si miramos otros distritos con características similares, tanto en renta como en población y densidad poblacional, encontramos los distritos de Sant Andreu, Sant Martí y Ciutat Vella. En todos observamos una renta próxima a 80.000€, y una densidad poblacional de 22.000 personas por km2. Observamos como Sant Andreu es el distrito con menos casos totales, y menor renta, pero es el distrito con mayor incidencia de casos por 10.000 habitantes. Esto nos reafirma en la creencia de que, en distritos con menor renta, en las viviendas tienden a vivir más personas, por tanto, parece que se detectan muchos casos en los núcleos familiares y eso hace que la incidencia sea alta, pero después el virus no se contagia hacia otras residencias, reduciendo el número de casos totales, como puede suceder en los otros distritos, donde la incidencia es menor, pero el número de casos totales es mayor. En conclusión, parece ser que la renta, más que afectar al número de casos totales, afecta a la incidencia de casos totales por 10.000 habitantes, ya que las bajas rentas promueven que se creen núcleos residenciales más acinados. Mientras que la población total y, en concreto, su densidad, afectan al número de casos totales detectados, ya que, a más población o densidad poblacional, existe más probabilidad de contagio. Probablemente el distrito más curioso sea el de Sants-Montjuïc, ya es el segundo distrito con más población, pero vemos que es el segundo distrito con una menor densidad. Esto se puede explicar en gran parte al complejo de Montjuïc, el cual no es habitable. Por tanto, podemos intuir que esta densidad es un poco irreal, ya que, aunque el distrito sea muy grande, no todo es zona habitable, y resulta en una densidad poblacional baja. Esto lo que provoca, son núcleos de 16 viviendas concentrados, y provoca que los casos en Sants-Montjuïc no acaben de seguir la lógica establecida por la densidad de población y la población total. Tiene un número elevado de casos por la baja densidad de población, pero si miramos la población total, se observan menos casos que en otros distritos con menos residentes. Otra curiosidad que no se muestra en la tabla, es que en todos los distritos el mínimo de casos detectados en un día fue 0. Esta información no es muy relevante ya que la detección no significa que no hubiera gente nueva contagiada ese día. En definitiva, podemos observar en la tabla como existen diferencias a nivel de contagio en los diferentes distritos, y entre sexos, lo que defiende nuestras hipótesis iniciales de que realmente existen diferencias entre géneros y tipos de distrito. En general, los distritos con mayor población y densidad son los distritos con más casos, mientras que, a nivel de sexo, parece que son las mujeres aquellas con más propensión al contagio (Véanse Tablas Anexos 3 y 4). III.2. Identificación de los modelos Para poder determinar el tipo de modelo y su orden, se procederá a analizar los coeficientes de autocorrelación (ACF) y los coeficientes de autocorrelación parciales (PACF) de cada distrito. Estos gráficos nos indican el grado de correlación que tiene las observaciones de la serie. Para ver si existe correlación tenemos que observar cuantos coeficientes son superiores a 0, dichos coeficientes nos dirán el grado de correlación entre las observaciones. En nuestro caso, serán el número de días anteriores que afectan al número de casos del día analizado. Debido a que la COVID-19 es una enfermedad contagiosa, donde los contagios de un día dependen de los casos detectados en los días anteriores, creemos correcto suponer que los modelos autorregresivos son los indicados para su análisis. Estos modelos se caracterizan por tener un ACF donde el número de coeficientes estadísticamente diferentes de 0, decrecen de forma exponencial hacia 0, y unos PACF con p coeficientes superiores a 0, los cuales nos indican el orden del modelo. Cosa muy importante para destacar que sabemos, y era esperado, es que nuestra serie no es estacionaria, y, por tanto, esperamos que los resultados que obtengamos no serán perfectos al 100%. 17 • ACF Ilustración 6. FAS para las mujeres en los diferentes distritos de Barcelona Ilustración 7. FAS de los hombres en los distintos distritos de Barcelona Si observamos las ilustraciones 6 y 7, podemos ver claramente como los coeficientes de autocorrelación en ningún distrito tienden a 0 a lo largo del tiempo, pero vemos como descienden exponencialmente con el tiempo. También es importante destacar el comportamiento de dichos coeficientes y la forma que dibujan. El primer coeficiente es siempre 1, pero luego observamos que entre los coeficientes 2 y 8, desciende y asciende formando una especie de valle, un comportamiento simétrico. Esta simetría nos indica que puede existir cierta estacionalidad de orden 7, es decir, semanal. Dicha estacionalidad se podrá confirmar cuando analicemos los PACF. El hecho de que los coeficientes no hayan alcanzado el valor 0, se explica por esta posible estacionalidad. Si realmente existe, al haber solo 21 coeficientes solo hemos dado 3 “vueltas”, y, por eso es posible que aún no se puedan considerar 0. Si calculamos más coeficientes seguramente viéramos a los coeficientes alcanzar la región significativa. 18 • PACF Ilustración 8. PACF para las mujeres distintos distritos de Barcelona Ilustración 9. PACF para los hombres en los distintos distritos de Barcelona Al contrario que en el gráfico de la ACF, podemos ver claramente, como en ambos sexos, el número de coeficientes no descienden exponencialmente a 0, si no que únicamente observamos como el primer coeficiente es estadísticamente diferente de 0, y luego los demás coeficientes son estadísticamente equivalentes a 0, exceptuando algún coeficiente específico. Esto nos reafirma definitivamente que nuestros modelos son del tipo AR. (Véase tabla Anexos.2) En la mayoría de los distritos el número de coeficientes superiores a 0, antes de equivaler a 0, son 1. Se pueden apreciar ciertos coeficientes que superan el 0 en momentos puntuales, lo cual nos podría indicar cierta estacionalidad. Es decir, que cada x días hubiera más probabilidad de contagio, por ejemplo, los findes de semana es probable que la gente socialice más, y, por tanto, hay más posibilidades de que los lunes el número de contagios aumente. Esto podría provocar una pequeña modificación en los modelos, ya que el día de la semana afectaría al número de innovaciones que se suceden. 19 Observamos que dichos coeficientes se repiten en los coeficientes 1, 7 y 14, lo que nos indica una estacionalidad de orden 7, tal y como sospechamos anteriormente. Los demás coeficientes que son estadísticamente diferentes de 0, y no ocupan una posición múltiple de 7, se explican mediante el componente aleatorio de la serie, que influye en que algunos coeficientes salgan significativos, pero no tiene mayor importancia. En conclusión, gracias al análisis de los ACF y PACF hemos podido confirmar que realmente los modelos a parametrizar son del tipo autorregresivo, y, además, estos modelos serán de orden 1 con una estacionalidad de orden 7. III.3. Modelos AR Durante todo este apartado, al igual que en los anteriores, hay que recordar que trabajamos con el contaje de casos de la COVID-19 19 en Barcelona, con los datos desagrados por sexo y distrito. III.3.1 Estimación de los parámetros El modelo autorregresivo sigue la fórmula (1) explicada anteriormente. Como podemos observar, al haber determinado que el orden de nuestros modelos es 1, únicamente tendremos que estimar el parámetro ϕ1. Para ello usaremos el método de la máxima verosimilitud. Una vez aplicada el código de R obtenemos los parámetros siguientes: No se observan grandes diferencias en ϕ1 a nivel de sexos, parece que los hombres obtienen valores ligeramente superiores. La mayor separación se encuentra en el distrito de Nou Barris, donde los casos del día anterior afectan mucho menos a las mujeres que a los hombres. Tabla 3.2 Parámetros modelo AR Mujeres Hombres 𝛟𝟏 Desviación estándard Intercept Desviación estándard 𝛟𝟏 Desviación estándard Intercept Desviación estándard CIUTAT VELLA 0,6448 0,0397 10,3741 0,1114 0,6573 0,0391 11,5436 0,1212 EIXAMPLE 0,6660 0,0388 28,3597 0,1108 0,6710 0,0385 23,6629 0,0968 GRÀCIA 0,6417 0,0399 11,6624 0,1065 0,6342 0,0402 9,1789 0,0822 HORTAGUINARDÓ 0,6828 0,0380 20,6341 0,1255 0,6849 0,0379 15,6833 0,0985 LES CORTS 0,5661 0,0428 7,9256 0,1155 0,6412 0,0398 7,0891 0,1168 NOU BARRIS 0,6754 0,0383 22,7039 0,1445 0,7160 0,0363 18,3635 0,1272 SANT ANDREU 0,6634 0,0389 16,7437 0,1182 0,6564 0,0392 14,4677 0,1007 SANT MARTÍ 0,6370 0,0401 24,1985 0,0973 0,6249 0,0406 20,1791 0,0812 SANTSMONTJUÏC 0,5919 0,0420 21,3601 0,1101 0,6015 0,0416 18,1391 0,093 SARRIÀ-SANT GERVASI 0,6717 0,0385 14,6853 0,1163 0,7214 0,036 11,7370 0,1046 20 El parámetro intercept hace referencia a las innovaciones, es decir, los casos nuevos detectados diariamente. Vemos que, en todos los distritos, excepto Ciutat Vella, las mujeres tienden a tener más casos diarios que los hombres. Esto se podría entender como que las mujeres se contagian en mayor abundancia que los hombres. Esta observación se complementa con lo observado en la descriptiva desagregada por sexos (tablas Anexo.3 y Anexos.4), donde las mujeres eran el sexo con más casos detectados, y refuerzan nuestra hipótesis de que existen diferencias entre sexos a la hora de casos detectados. A nivel de significación de los parámetros, hemos calculado su error estándar y hemos comprobado que ϕ1 es significativamente diferente de 0. III.3.2. Validación del modelo Una vez tenemos un posible modelo, hemos de validarlo, comprobar que realmente cumple las condiciones necesarias para poder afirmar que este modelo explica la serie, y, por tanto, poder aplicarlo para intentar predecir los valores futuros de la serie. Hay que destacar que nuestra serie no es estacionaria, y, por tanto, puede fallar alguna validación. La validación consta de diversos pasos: 1. Normalidad de los residuos 2. Autocorrelación residuos III.3.2.1. Normalidad de los residuos Para asegurar la normalidad de los residuos llevaremos a cabo la prueba de Shapiro-Wilk. Tabla 3.3. Valores P-valor prueba de Shapiro_Wilk P-Valor Hombres Mujeres CIUTAT VELLA 1,45E-15 1,40E-15 EIXAMPLE 9,91E-17 1,27E-16 GRÀCIA 2,99E-16 3,88E-18 HORTA-GUINARDÓ 1,46E-14 1,83E-12 LES CORTS 1,75E-20 1,56E-23 NOU BARRIS 7,52E-17 4,88E-16 SANT ANDREU 8,31E-16 9,35E-15 SANT MARTÍ 5,79E-19 7,08E-18 SANTS-MONTJUÏC 8,17E-18 2,76E-16 SARRIÀ-SANT GERVASI 1,39E-16 3,35E-16 Observamos como todos los p-valores de las pruebas son inferiores a 0.05, y, por tanto, tenemos evidencias suficientes para rechazar la hipótesis nula de normalidad, y afirmar que la distribución de la serie no es una normal. Este resultado tiene todo el sentido del mundo ya que nuestra serie es un contaje de valores bajos, y, por tanto, es muy difícil que su 21 distribución marginal siga una normal. Además, sabemos que la serie no es estacionaria y por ello el modelo sabemos que no será perfecto. III.3.2.2. Independencia de los residuos Tenemos que asegurar que los residuos del modelo se comportan como ruido blanco. Para ello llevaremos a cabo la prueba de Box-Pierce. Tabla 3.4. P-valor prueba de Box-Pierce P-Valor Hombres Mujeres CIUTAT VELLA 0,3226 0,8570 EIXAMPLE 0,6071 0,6860 GRÀCIA 0,0502 0,6101 HORTA-GUINARDÓ 0,7243 0,9480 LES CORTS 0,5741 0,3496 NOU BARRIS 0,6489 0,6655 SANT ANDREU 0,1975 0,6889 SANT MARTÍ 0,2015 0,5350 SANTS-MONTJUÏC 0,6671 0,9109 SARRIÀ-SANT GERVASI 0,5314 0,0969 Hemos obtenido todos los p-valores superiores a 0.05, por tanto, con un nivel de significación del 95%, no tenemos suficientes evidencias para rechazar la hipótesis nula de independencia, y, entonces podemos afirmar que los residuos son independientes entre sí. En conclusión, hemos validado correctamente el modelo, a excepción de la normalidad de los residuos, y, a continuación, podemos llevar a cabo las predicciones. III.3.3 Predicciones Para las predicciones, no solo comprobaremos la capacidad predictiva del modelo para valores conocidos (Predicción ex - post) sino que también comprobaremos la eficacia del modelo para valores futuros desconocidos (Predicción ex - ante). Las predicciones seguirán la forma expresada en la fórmula (3). III.3.3.1 Predicción Ex - Post Para la predicción Ex - Post, cogeremos la serie desde el 1/03/2020 hasta el 31/12/2020, y calcularemos las predicciones diarias para los meses de enero y febrero de 2021. Los parámetros de la Tabla 3.2. hacen referencia a los modelos con toda la serie como período 28 algunos tardan un par horas o días, y eso hace que un caso salga detectado, por ejemplo, un jueves cuando en realidad el contagio se produjo el martes. A nivel de sexos, observamos como las mujeres presentan innovaciones mucho más abundantes que los hombres. En la mayoría de los días y distritos, se detectan más nuevos casos diarios de mujeres que de hombres. Esto coincide con lo observado en la descriptiva de la serie (Véanse tablas Anexos.3 y Anexo.4), donde las mujeres eran el sexo con más casos detectados. En conclusión, dichos resultados nos refuerzan la hipótesis que existen diferencias a nivel de contagio entre los sexos. III.4.2 Predicciones Para las predicciones del modelo INAR seguiremos la misma estructura que en los modelos AR. En este caso las predicciones estarán basadas en la fórmula (7). III.4.2.1. Predicciones Ex - Post Para las predicciones Ex - Post usaremos las observaciones desde el 1 de marzo de 2020 hasta el 31 de diciembre de 2020, y predeciremos los valores para los meses de enero y febrero de 2021, y los compararemos con los valores reales observados. Como el caso de los modelos AR, volveremos a parametrizar nuevos modelos, ya que la muestra ha cambiado. (Ver Tabla Anexos.7 y Anexos.8. Parámetros modelo INAR predicciones Ex-Post.) 29 Ilustración 14. Predicciones Ex-Post del modelo INAR para las mujeres Ilustración 15. Predicciones Ex -Post del modelo INAR para los hombres 30 Observamos como el modelo es capaz de adaptarse perfectamente a las tendencias que sigue la serie. Una cosa muy destacable es que podemos observar como al principio de la serie las predicciones son inferiores, pero a medida que pasa el tiempo y los casos detectados descienden las predicciones se ajustan a la perfección. Este evento se podría explicar debido a que el inicio de las predicciones es enero y es donde nos encontramos con la tercera ola. El aumento exponencial durante enero hace que el modelo no sea capaz de explicar a la perfección este período muestral, ya que las observaciones de diciembre no indicaban que se pudiera dar estos valores extremos. Por otro lado, a medida que el número de casos disminuye y toman un comportamiento más lógico, el modelo es capaz de explicar perfectamente la serie y predecirla. Al igual que en los modelos AR, calcularemos los errores de predicción, la cobertura, y, el rango de las regiones de predicción. Los errores de predicción obtenidos son: (Véanse fórmulas 9) Tabla 3.10. Errores de predicción para el modelo INAR EQM EAM EPAM Hombres Mujeres Hombres Mujeres Hombres Mujeres CIUTAT VELLA 61,10 110,19 6,33 7,85 66,69% 56,07% EIXAMPLE 899,92 1181,92 21,31 23,92 46,32% 64,54% GRÀCIA 174,66 209,56 9,27 10,28 48,04% 64,90% HORTA-GUINARDÓ 228,16 372,94 10,87 14,06 41,14% 63,32% LES CORTS 194,51 237,61 9,33 9,97 65,08% 62,54% NOU BARRIS 484,15 680,92 16,66 18,95 62,76% 63,30% SANT ANDREU 340,27 426,98 13,47 14,69 53,88% 59,18% SANT MARTÍ 528,63 705,32 16,66 18,74 54,33% 42,39% SANTS-MONTJUÏC 414,76 496,48 14,56 15,89 58,12% 61,18% SARRIÀ-SANT GERVASI 353,95 544,18 14,69 17,44 66,22% 52,75% Observamos que los EPAM son superiores al 10%, por tanto, nos indican bajas capacidades predictivas de los modelos INAR estimados. Algo destacable es que los EPAM son inferiores para los hombres que para las mujeres. Esto quiere decir que dentro de los modelos INAR, aquellos para los hombres son los más capaces para explicar la serie. (Véase Tabla Anexos.6) A continuación, para aumentar el análisis de la capacidad predictiva de los modelos calcularemos la cobertura de las predicciones, porcentaje de observaciones dentro de las regiones de predicción, y el rango de dichas regiones de predicción. 31 Los resultados obtenidos para la cobertura son: Tabla 3.11. Porcentaje de Cobertura para las predicciones del modelo INAR Cobertura Hombres Mujeres CIUTAT VELLA 61,29% 56,45% EIXAMPLE 40,32% 43,55% GRÀCIA 51,61% 56,45% HORTA-GUINARDÓ 51,61% 43,55% LES CORTS 51,61% 54,84% NOU BARRIS 32,26% 41,94% SANT ANDREU 46,77% 45,16% SANT MARTÍ 43,55% 48,39% SANTS-MONTJUÏC 38,71% 48,39% SARRIÀ-SANT GERVASI 27,42% 32,26% Observamos coberturas aceptables, superiores al 50%, en los distritos de Ciutat Vella, Gràcia, Horta-Guinardó y Les Corts. Pero el problema lo tenemos en los otros distritos, con porcentajes inferiores al 50% alcanzando cotas mínimas del 27% en el distrito de Sarrià-Sant Gervàsi. Dichas coberturas no son especialmente altas, y esto concuerda con lo obtenido en los EPAM, donde se observaba una baja capacidad predictiva. Además, observando las Ilustraciones 14 y 15, sobre todo al principio de la serie, se ve con claridad como las regiones de predicción no incluyen las observaciones, y esto provoca bajas coberturas. A continuación, para poder explicar un poco mejor estas coberturas calcularemos el rango de las regiones de predicción. Los resultados obtenidos son: Tabla 3.12. Rangos de las regiones de predicción modelos INAR Rangos de predicción Hombres Mujeres CIUTAT VELLA 12,26 11,50 EIXAMPLE 16,34 18,55 GRÀCIA 10,08 12,26 HORTA-GUINARDÓ 14,11 16,69 LES CORTS 8,77 9,66 NOU BARRIS 14,44 17,56 SANT ANDREU 13,29 14,63 SANT MARTÍ 15,32 17,56 SANTS-MONTJUÏC 14,71 17,24 SARRIÀ-SANT GERVASI 12,19 13,47 Observamos rangos muy pequeños. El rango más amplio lo encontramos en 18.55 casos para las mujeres en l’Eixample. Estos rangos nos están indicando que nuestros modelos INAR intentan ser muy precisos en sus predicciones, aunque parece que si observamos la cobertura dichos modelos no son capaces de adaptarse bien a los valores de la serie. 32 III.4.2.2. Predicciones Ex - Ante Para la predicción Ex - Ante cogemos la totalidad de nuestra serie, desde el 1 de marzo de 2020 hasta el 3 de marzo de 2021, e intentamos predecir las observaciones que se obtendrían en los próximos 20 días. Hay que tener en cuenta que tenemos estacionalidad de orden 7, y, por tanto, las lambdas se comportaran siguiendo la relación: λ𝑛=λ𝑛+7 Ilustración 16. Predicciones Ante del modelo INAR para las mujeres Ilustración 17. Predicciones Ante del modelo INAR para los hombres Se puede observar cómo claramente los modelos han sido capaces de adaptarse a la tendencia de las últimas observaciones y ofrecer unas predicciones plausibles, que nos pueden permitir llevar estrategias eficientes. También se puede destacar como las mujeres obtienen predicciones 33 superiores a los hombres, lo cual tiene sentido ya que como vimos (Véanse tablas Anexos 3 y 4), las mujeres son aquellas con más casos detectados. En conclusión, parece ser que nuestro modelo tiene una baja capacidad predictiva, tanto, para las predicciones Ex - Post como para las predicciones Ex - Ante, y presenta unas coberturas bajas en algunos distritos, lo cual se debe a que los modelos INAR son muy precisos, con regiones de predicción bastante estrechas, y la tercera ola con valores extremos que encontramos en enero de 2021. III.5. Comparación de los modelos Una vez hemos estimado los dos tipos de modelos: AR e INAR, consideramos oportuno comparar ambos modelos, tanto a nivel de parámetros (Véanse tablas 3.2, 3.8 y 3.9), como a, nivel de predicciones (Véanse tablas 3.5, 3.6, 3.7, 3.10, 3.11, 3.12) para discernir cuál de los dos explica mejor la serie temporal de número de casos diarios en los diferentes distritos de Barcelona. A nivel de parámetros vamos a dividirlos en dos categorías: resultados en función del sexo y en función del distrito. • Sexo: A nivel del parámetro α observamos mayores valores para los hombres, lo que nos empieza a indicar la existencia de diferentes mecanismos para los sexos. Luego a nivel de innovaciones, casos que no depende de los días anteriores, es en las mujeres en las que se detectan más casos, es decir que las mujeres tienen más propensión a contagiarse. Estos resultados concuerdan con los obtenidos en la descriptiva realizada donde observábamos que más de la mitad de los casos detectados eran de mujeres. (Véanse tablas Anexos.3 y Anexos.4) • Distritos: En cuanto a la diferenciación de los parámetros en los diferentes distritos, en ambos modelos obtenemos que el número de casos, tanto los dependientes de días anteriores como las innovaciones aumentan cuanto mayor es la densidad poblacional del distrito. También hemos observado como la renta afecta a los casos a nivel de que, en distritos con renta baja, las agrupaciones familiares suelen ser más extensas que en los distritos de renta más alta, y por ello, cuando un miembro familiar se contagia, rápidamente se contagian más personas que en otros distritos donde las agrupaciones pueden ser inferiores. Renta y densidad poblacional están bastante relacionadas. Distritos con renta alta tienen menor densidad poblacional, ya que poca gente puede permitirse el coste de vida, y, por ende, en los distritos de renta más baja se producen mayores agrupaciones de viviendas, lo que aumenta la densidad poblacional, y facilita la propagación del virus. 34 Ilustración 18. Número de casos de la COVID-19 en relación a los quintiles privación socioeconómica[11] Con la imagen superior (Ilustración 18) no observamos una tendencia ascendente, desde el primer quintil hasta el último, pero si podemos ver como en el último quintil, aquel con mayor privación socioeconómica, y, por tanto, menor renta, alcanzamos el máximo de casos. En conclusión, observamos como la renta influye en el número de casos detectados, a menor renta, mayor número de casos de la COVID-19. A nivel de predicciones, obtenemos resultados malos en ambos modelos. En relación con los errores de predicción, en ambos modelos obtenemos EPAM elevados, inferiores en los modelos INAR, indicando bajas capacidades predictivas. En adición, los modelos AR presentan porcentajes de cobertura muy altos, superiores al 60%, a diferencia de los modelos INAR, donde solo 4 distritos son capaces de superar el 50% y ni siquiera llegan al 60%, pero hay que comentar que dichos porcentajes de los modelos AR son un poco “irreales”, ya que los rangos de las regiones de predicción de los modelos AR son mucho más grandes que los modelos INAR. En los modelos AR el rango varía desde 20 casos diarios en un distrito hasta más de 50 casos en varios distritos, mientras que en los modelos INAR el máximo rango que obtenemos es de 19 casos en un solo distrito, mientras que los demás distritos se mueven en valores en torno a los 10-15 casos de rango. En conclusión, hemos podido ver como ambos modelos se comportan de la misma manera a la hora de estimar de los parámetros, diferenciando correctamente entre sexos y distritos, en la parte regresiva los hombres obtienen valores superiores, mientras que en la parte de las innovaciones es al revés, mientras que en los distritos se observa claramente como los barrios con mayor densidad poblacional obtienen parámetros con valores más altos. A nivel de predicciones obtenemos que ambos modelos tienen una capacidad predictiva muy baja, siendo los modelos INAR los que obtienen un mejor EPAM y más precisos, mientras que los AR son los menos eficaces y más imprecisos. 35 IV. CONCLUSIONES La COVID-19 ha tenido un impacto en todos los niveles de la sociedad, (Véanse Ilustraciones. 1, 2, 3, 4, y Anexos.1) es por ello por lo que hemos considerado necesario parametrizar modelos que puedan ajustarse a la serie y así disponer de instrumentos estadísticos que nos permitan predecir el comportamiento del virus para llevar a cabos estrategias eficientes, y mitigar su impacto. Nuestra hipótesis inicial, es que al ser la COVID-19 una enfermedad infecciosa, donde los casos diarios dependen de los casos detectados los días anteriores, los modelos más adecuados para parametrizar serían los modelos autorregresivos. Para confirmar esta hipótesis hemos observados los coeficientes de autocorrelación (Véanse Ilustraciones 6, 7, 8, y 9). Al analizar los coeficientes también hemos descubierto la existencia de una estacionalidad de orden 7, es decir, semanal. Aunque no fuera una hipótesis inicial, el hecho de que el contagio de la COVID19 se comporte diferente en función del día de la semana tiene todo el sentido del mundo. Los fines de semana, al haber más vida social, la probabilidad de contagio aumenta con respecto a otro día de la semana. En conclusión, utilizaremos los modelos tipo autorregresivos, más concreto el modelo AR, para distribuciones continuas, y el modelo INAR, para distribuciones discretas. Dichos modelos serán de orden 1, y presentarán una estacionalidad de 7 días. Una vez tenemos el tipo de modelos a estimar, hemos de comparar los resultados. Tal y como comentamos al inicio del trabajo, nuestra hipótesis es que tanto el sexo, como el distrito a analizar influyen en el número de casos. En relación con el sexo, dicha hipótesis se defiende al ver como las mujeres son aquellas con más de la mitad de los casos (Véanse Tablas Anexos.3 y 4), y también observando las diferencias en los parámetros del modelo, tanto a nivel regresivo, donde los hombres presentan valores superiores, como a nivel de las innovaciones, donde vemos como las mujeres llevan la delantera (Véanse Tablas 3.2, 3.8, 3.9, Anexos.5, 7, y 8). A nivel de distritos también se observa como la densidad poblacional del distrito es un elemento clave para analizar los contagios. Aquellos distritos con más densidad poblacional son los distritos con mayores casos detectados y parámetros del modelo. También hemos concluido que la renta influye de una forma indirecta, distritos con alta renta son poco poblados, y, por tanto, se dificulta la propagación del virus. (Véase Tabla 3.1) Los modelos AR obtienen resultados optimistas, pero un poco irreales. Obtienen unos EPAM alrededor del 70-100%, lo que nos indica una bajísima capacidad predictiva (Véase Tabla 3.5), pero, por otro lado, observando la representación gráfica se observa como los valores predichos son incluidos en las regiones de predicción (Véanse Ilustraciones 10 y 11). Esto supone una cobertura mínima del 65% en todos los distritos (Véanse tabla 3.6), la cual está influenciada por las enormes regiones de predicción del modelo AR, donde el rango varía desde los 20 casos diarios hasta los 50. (Véase tabla 3.7) En cambio, los modelos INAR presentan mejores EPAM (Véase Tabla 3.10), y por las representaciones gráficas se observa como también se adaptan mucho mejor al comportamiento de la serie (Véanse Ilustraciones 14 y 15), mientras que los modelos AR acaban convergiendo a la media, el modelo INAR es capaz de dibujar las tendencias que sigue la serie. A nivel de cobertura es inferior para los modelos INAR (Véase Tabla 3.11), pero también lo es la región de predicción (Véase Tabla 3.12), por lo que el modelo INAR tiende a ser más preciso que los modelos AR. 36 Tal y como se ha comentado a lo largo del trabajo, esperábamos que la serie no fuera estacionaria y así ha resultado, por tanto, destacamos que ninguno de los dos modelos ha cumplido todas las validaciones y, por ende, ninguna estimación será perfecta. Aun así, recomendamos enormemente el modelo INAR sobre el AR, no solo por los mejores resultados, y el mejor comportamiento observado, sino porque al ser nuestra serie un contaje de valores bajos, un modelo que siga una distribución discreta (INAR), tiene mejor expectativa para explicar la serie que un modelo con distribución marginal continua (AR). 37 V. BIBLIOGRAFIA [1] “Coronavirus: La OMS declara el brote de Covid-19 pandemia.” https://www.redaccionmedica.com/secciones/sanidad-hoy/coronavirus-pandemiabrote-de-covid-19-nivel-mundial-segun-oms-1895 (accessed Mar. 06, 2021). [2] M. Woodward, S. A. E. Peters, and K. Harris, “Social deprivation as a risk factor for COVID-19 mortality among women and men in the UK Biobank: nature of risk and context suggests that social interventions are essential to mitigate the effects of future pandemics,” J. Epidemiol. Community Health, p. jech-2020-215810, Apr. 2021, doi: 10.1136/jech-2020-215810. [3] R Core Team, “R: A Language and Environment for Statistical Computing,” 2019, [Online]. Available: https://www.r-project.org/. [4] H. Wickham., “ggplot2: Elegant Graphics for Data Analysis.,” Springer-Verlag New York, 2016. [5] D. Jin‐Guan and L. Yuan, “THE INTEGER‐VALUED AUTOREGRESSIVE (INAR(p)) MODEL,” J. Time Ser. Anal., vol. 12, no. 2, pp. 129–142, 1991, doi: 10.1111/j.1467-9892.1991.tb00073.x. [6] D. Moriña, P. Puig, J. Ríos, A. Vilella, and A. Trilla, “A statistical model for hospital admissions caused by seasonal diseases,” Stat. Med., vol. 30, no. 26, pp. 3125–3136, Nov. 2011, doi: 10.1002/sim.4336. [7] F. De, C. Matemáticas, and J. A. T. Cazorla, UNIVERSIDAD COMPLUTENSE DE MADRID MODELO “INARMA” PARA SERIES TEMPORALES DE VALORES ENTEROS: ANÁLISIS, PROPIEDADES ASINTOMÁTICAS Y ESTIMACIÓN. MEMORIA PARA OPTAR AL GRADO DE DOCTOR PRESENTADA POR José Nerys Funes Torres Bajo la dirección del doctor. 2010. [8] C. Felipe González López, A. : María, and E. Correal Núñez, “Modelos de Series de Tiempo para Datos Enteros.” [9] D. Moriña, P. Puig, J. Ríos, A. Vilella, and A. Trilla, “A statistical model for hospital admissions caused by seasonal diseases,” Stat. Med., vol. 30, no. 26, pp. 3125–3136, Nov. 2011, doi: 10.1002/sim.4336. [10] R. Bu and B. McCabe, “Model selection, estimation and forecasting in INAR(p) models: A likelihood-based Markov Chain approach,” Int. J. Forecast., vol. 24, no. 1, pp. 151–162, Jan. 2008, doi: 10.1016/j.ijforecast.2007.11.002. [11] “#COVID19aldiaBCN.” https://aspb.shinyapps.io/COVID19_BCN/#Nivell_socioeconòmic_de_l’àrea_de_resid ència (accessed Jun. 17, 2021). 44 casos <- df[df$SexeCodi==s,] casos <- casos$Casos ed<-c(ed,summary(casos),sum(casos),sd(casos),(sum(casos)/poblacio[i])*10000) } } edT <- c() for (i in 1:10){ df <- dadesBCN3[dadesBCN3$SectorSanitariDescripcio==unique(dadesBCN3$SectorSanitariDescripcio)[i],] casos <- df$Casos edT<-c(edT,summary(casos),sum(casos),sd(casos),(sum(casos)/poblacio[i])*10000) } ed<-matrix(ed,nrow =10,byrow=TRUE) edH <- data.frame("Mean"=ed[,4],"Desviación estándard"=ed[,8],"Max"=ed[,6],"Num Casos"=ed[,7],"Població"=poblacio,"Casos por 10.000 habitantes"=ed[,9],"Densidad poblacional"= dens,row.names = unique(dadesBCN3$SectorSanitariDescripcio)) edD <- data.frame("Mean"=ed[,13],"Desviación estándard"=ed[,17],"Max"=ed[,15],"Num Casos"=ed[,16],"Població"=poblacio,"Casos por 10.000 habitantes"=ed[,18],"Densidad poblacional"=dens ,row.names = unique(dadesBCN3$SectorSanitariDescripcio)) edT <- matrix(edT,nrow=10,byrow=TRUE) edT <- data.frame("Media"=edT[,4],"Desviación estándard"=edT[,8],"Máximo"=edT[,6],"Número de Casos Totales"=edT[,7],"Población"=poblacio,"Casos por 10.000 habitantes"=edT[,9],"Densidad poblacional"=dens,"Renta familiar en miles de ???"= renta,row.names = unique(dadesBCN3$SectorSanitariDescripcio)) write.csv(edT,"C:/Users/pablo/OneDrive/Escritorio/TFG/DescriptivaTotal.csv", row.names = TRUE) write.csv(edH,"C:/Users/pablo/OneDrive/Escritorio/TFG/DescriptivaHombres.csv", row.names = TRUE) 45 write.csv(edD,"C:/Users/pablo/OneDrive/Escritorio/TFG/DescriptivaMujeres.csv", row.names = TRUE) #FAS de cada barrio for (s in levels(dadesBCN3$SexeCodi)){ df <- dadesBCN3[dadesBCN3$SexeCodi==s,] par(mfrow=c(2,5)) for(i in 1:10){ dd <- df[df$SectorSanitariDescripcio==unique(df$SectorSanitariDescripcio)[i],] dd <- ts(dd$Casos) acf(dd,ylim=c(-1,1),lag=21, main=c('FAS',s,substring(unique(dadesBCN3$SectorSanitariDescripcio)[i],first=10))) } } #FAP de cada barrio for (s in levels(dadesBCN3$SexeCodi)){ df <- dadesBCN3[dadesBCN3$SexeCodi==s,] par(mfrow=c(2,5)) for(i in 1:10){ dd <- df[df$SectorSanitariDescripcio==unique(df$SectorSanitariDescripcio)[i],] dd <- ts(dd$Casos) pacf(dd,ylim=c(-1,1),lag=21, main=c('FAP',s,substring(unique(dadesBCN3$SectorSanitariDescripcio)[i],first=10))) } } 46 #Modelo AR(1) #Estimación arc<- c() ars <- c() arn <- c() ari <-c() arcv <- c() par(mfrow=c(2,5)) for (s in levels(dadesBCN3$SexeCodi)){ df <- dadesBCN3[dadesBCN3$SexeCodi==s,] for (i in 1:10){ casos <- df[df$SectorSanitariDescripcio==unique(df$SectorSanitariDescripcio)[i],] model<-arima(casos$Casos,c(1,0,0),seasonal=list(order=c(0,0,0),period=7)) arc <- c(arc,model$coef) ars <- c(ars,2*pnorm(c(abs(model$coef)/sqrt(diag(model$var.coef))), mean=0, sd=1, lower.tail=FALSE)) arn<-c(arn,shapiro.test(model$residuals)$p.value) ari<- c(ari,Box.test(model$residuals)$p.value) arcv <- c(arcv,sqrt(diag(model$var.coef))) #Predicciones Ex -Ante serief<-predict(model,n.ahead=20) inf = serief$pred - 2*serief$se inf<- ifelse(inf<0,0,inf) sup = serief$pred + 2*serief$se 47 x<-casos$TipusCasData[349:368] y<-seq.Date(as.Date("2021-03-04"),as.Date("2021-03-23"),1) plot(casos$Casos[349:368],col="cadetblue1", main=c('Prediccion ExAnte',s,substring(unique(dadesBCN3$SectorSanitariDescripcio)[i],first=10)),type="lines",ylab=" Número Casos",xlab="Fecha",xaxt="n",xlim=c(0,40),ylim=c(0,max(sup)+5)) lines(seq(21,40,1),serief$pred, col="red" ) lines(seq(21,40,1),inf, col="darkgreen") lines(seq(21,40,1),sup, col="darkgreen") axis(1,at=seq(1,40,1),labels=c(x,y)) } } #Parámetros significativos y residuos independientes, no normalidad, no hay estacionariedad arc <- matrix(arc,ncol=2,byrow = TRUE) arcH <- data.frame("Alpha"=arc[1:10,1],"Intercept" = arc[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arcD <- data.frame("Alpha"=arc[11:20,1],"Intercept" = arc[11:20,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arcv <- matrix(arcv,ncol=2,byrow = TRUE) arcvH <- data.frame("Alpha"=arcv[1:10,1],"Intercept" = arcv[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arcvD <- data.frame("Alpha"=arcv[11:20,1],"Intercept" = arcv[11:20,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) ars <- matrix(ars,ncol=2,byrow = TRUE) arsH <- data.frame("Alpha"=ars[1:10,1],"Intercept" = ars[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arsD <- data.frame("Alpha"=ars[11:20,1],"Intercept" = ars[11:20,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) 48 arn <- matrix(arn,ncol=1,byrow = TRUE) arnH <- data.frame("Alpha"=arn[1:10,],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arnD <- data.frame("Alpha"=arn[11:20,],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) ari <- matrix(ari,ncol=1,byrow = TRUE) ariH <- data.frame("Alpha"=ari[1:10,],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) ariD <- data.frame("Alpha"=ari[11:20,],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) write.csv(arcH,"C:/Users/pablo/OneDrive/Escritorio/TFG/ParametrosARHombres.csv", row.names = TRUE) write.csv(arcD,"C:/Users/pablo/OneDrive/Escritorio/TFG/ParametrosARMujeres.csv", row.names = TRUE) write.csv(arcvH,"C:/Users/pablo/OneDrive/Escritorio/TFG/ParametrosDEARHombres.csv", row.names = TRUE) write.csv(arcvD,"C:/Users/pablo/OneDrive/Escritorio/TFG/ParametrosDEARMujeres.csv", row.names = TRUE) write.csv(arsH,"C:/Users/pablo/OneDrive/Escritorio/TFG/SignificacionARHombres.csv", row.names = TRUE) write.csv(arsD,"C:/Users/pablo/OneDrive/Escritorio/TFG/SignificacionARMujeres.csv", row.names = TRUE) write.csv(arnH,"C:/Users/pablo/OneDrive/Escritorio/TFG/normalidadARHombres.csv", row.names = TRUE) write.csv(arnD,"C:/Users/pablo/OneDrive/Escritorio/TFG/normalidadARMujeres.csv", row.names = TRUE) write.csv(ariH,"C:/Users/pablo/OneDrive/Escritorio/TFG/independenciaARHombres.csv", row.names = TRUE) write.csv(ariD,"C:/Users/pablo/OneDrive/Escritorio/TFG/independenciaARMujeres.csv", row.names = TRUE) #Predicciones Ex-Post cobertura <- function(li,ls,obs){ kk<-c() 49 k<-0 for(i in 1:length(ls)){ if(li[i] <= obs[i] && obs[i] <= ls[i]){ k<-k+1 } } kk <- c(kk,k) return(kk) } arc<- c() ars <- c() arn <- c() ari <-c() arcv <- c() eqmAR<- c() eamAR <- c() epamAR <- c() cobsAR<-c() riAR<-c() par(mfrow=c(2,5)) for (s in levels(dadesBCN3$SexeCodi)){ df <- dadesBCN3[dadesBCN3$SexeCodi==s,] for (i in 1:10){ casos <- df[df$SectorSanitariDescripcio==unique(df$SectorSanitariDescripcio)[i],] serie<-casos$Casos[1:306] model2<-arima(serie,c(1,0,0),seasonal=list(order=c(0,0,0),period=7)) serie2f<-predict(model2,n.ahead=62,interval="confidence") 50 inf = serie2f$pred - 1.96*serie2f$se inf <-ifelse(inf<0,0,inf) sup = serie2f$pred + 1.96*serie2f$se arc <- c(arc,model2$coef) ars <- c(ars,2*pnorm(c(abs(model2$coef)/sqrt(diag(model2$var.coef))), mean=0, sd=1, lower.tail=FALSE)) arn<-c(arn,shapiro.test(model2$residuals)$p.value) ari<- c(ari,Box.test(model2$residuals)$p.value) arcv <- c(arcv,sqrt(diag(model2$var.coef))) preds <- data.frame(TCasos=serie2f$pred, LL=inf, UL=sup, TipusCasData=unique(dadesBCN4$TipusCasData[dadesBCN4$TipusCasData>"2020-12-31"])) d<- dadesBCN4[dadesBCN4$SectorSanitariDescripcio==names(table(dadesBCN4$SectorSanitariDescripcio)) [i] & dadesBCN4$SexeCodi==s & dadesBCN4$TipusCasData>"2020-12-31", ] assign(paste0("p",s, i), ggplot(d, aes(x=TipusCasData, y=TCasos)) + geom_line()+scale_x_date(date_breaks = "1 week", date_labels = "%m-%d- %Y")+theme(axis.text.x=element_text(angle=60, hjust=1), legend.position = "none")+ geom_ribbon(data=preds,aes(ymin=LL,ymax=UL),alpha=0.3)+geom_line(data=preds,aes(x=TipusCasData , y=TCasos, color = "yellow"))+ xlab("Data")+ylab("Nombre de casos")) #Errores predicción op<-data.frame("Obs"=casos[307:368,]$Casos,"Pred"=serie2f$pred) op <- op[op$Obs!=0,] errors<- op$Obs-op$Pred eqmAR<-c(eqmAR,sum(errors*errors)/length(op$Obs)) eamAR<-c(eamAR,sum(abs(errors))/length(op$Obs)) epamAR<-c(epamAR,sum(abs(errors)/abs(op$Obs))/length(op$Obs)) #Cobertura cobsAR<-c(cobsAR,cobertura(inf,sup,casos$Casos[307:368])) #Rangs Predicció 51 riAR<-c(riAR,mean(sup-inf)) } } ggarrange(pHome1,pHome2,pHome3,pHome4,pHome5,pHome6,pHome7,pHome8,pHome9,pHome10,labels=substr ing(unique(dadesBCN3$SectorSanitariDescripcio),first=10),ncol=5,nrow=2,common.legend = TRUE) ggarrange(pDona1,pDona2,pDona3,pDona4,pDona5,pDona6,pDona7,pDona8,pDona9,pDona10,labels=substr ing(unique(dadesBCN3$SectorSanitariDescripcio),first=10),ncol=5,nrow=2,common.legend = TRUE) #Parametres arc <- matrix(arc,ncol=2,byrow = TRUE) arcH <- data.frame("Alpha"=arc[1:10,1],"Intercept" = arc[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arcD <- data.frame("Alpha"=arc[11:20,1],"Intercept" = arc[11:20,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arcv <- matrix(arcv,ncol=2,byrow = TRUE) arcvH <- data.frame("Alpha"=arcv[1:10,1],"Intercept" = arcv[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arcvD <- data.frame("Alpha"=arcv[11:20,1],"Intercept" = arcv[11:20,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) ars <- matrix(ars,ncol=2,byrow = TRUE) arsH <- data.frame("Alpha"=ars[1:10,1],"Intercept" = ars[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arsD <- data.frame("Alpha"=ars[11:20,1],"Intercept" = ars[11:20,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arn <- matrix(arn,ncol=1,byrow = TRUE) arnH <- data.frame("Alpha"=arn[1:10,],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) arnD <- data.frame("Alpha"=arn[11:20,],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) ari <- matrix(ari,ncol=1,byrow = TRUE) ariH <- data.frame("Alpha"=ari[1:10,],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) ariD <- data.frame("Alpha"=ari[11:20,],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) 52 write.csv(arcH,"C:/Users/pablo/OneDrive/Escritorio/TFG/ParametrosARHombres.csv", row.names = TRUE) write.csv(arcD,"C:/Users/pablo/OneDrive/Escritorio/TFG/ParametrosARMujeres.csv", row.names = TRUE) write.csv(arcvH,"C:/Users/pablo/OneDrive/Escritorio/TFG/ParametrosDEARHombres.csv", row.names = TRUE) write.csv(arcvD,"C:/Users/pablo/OneDrive/Escritorio/TFG/ParametrosDEARMujeres.csv", row.names = TRUE) write.csv(arsH,"C:/Users/pablo/OneDrive/Escritorio/TFG/SignificacionARHombres.csv", row.names = TRUE) write.csv(arsD,"C:/Users/pablo/OneDrive/Escritorio/TFG/SignificacionARMujeres.csv", row.names = TRUE) write.csv(arnH,"C:/Users/pablo/OneDrive/Escritorio/TFG/normalidadARHombres.csv", row.names = TRUE) write.csv(arnD,"C:/Users/pablo/OneDrive/Escritorio/TFG/normalidadARMujeres.csv", row.names = TRUE) write.csv(ariH,"C:/Users/pablo/OneDrive/Escritorio/TFG/independenciaARHombres.csv", row.names = TRUE) write.csv(ariD,"C:/Users/pablo/OneDrive/Escritorio/TFG/independenciaARMujeres.csv", row.names = TRUE) #Errors Predicción eqmAR <- matrix(eqmAR,ncol=2) eqmAR <- data.frame("EQM Hombres"=eqmAR[1:10,1],"EQM Mujeres"=eqmAR[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) write.csv(eqmAR,"C:/Users/pablo/OneDrive/Escritorio/TFG/EQMAR.csv", row.names = TRUE) eamAR <- matrix(eamAR,ncol=2) eamAR <- data.frame("EAM Hombres"=eamAR[1:10,1],"EAM Mujeres"=eamAR[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) write.csv(eamAR,"C:/Users/pablo/OneDrive/Escritorio/TFG/EAMAR.csv", row.names = TRUE) epamAR <- matrix(epamAR,ncol=2) epamAR <- data.frame("EPAM Hombres"=epamAR[1:10,1],"EPAM Mujeres"=epamAR[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) 53 write.csv(epamAR,"C:/Users/pablo/OneDrive/Escritorio/TFG/EPAMARs.csv", row.names = TRUE) #Cobertura cobsAR<- matrix(cobsAR,ncol=2) cobsAR <- data.frame("Cobertura Hombres"=cobsAR[1:10,1],"Cobertura Mujeres"=cobsAR[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) write.csv(cobsAR,"C:/Users/pablo/OneDrive/Escritorio/TFG/CoberturaARs.csv", row.names = TRUE) #Rang Interval Predicció riAR<- matrix(riAR,ncol=2) riAR <- data.frame("Rang Hombres"=riAR[1:10,1],"Rang Mujeres"=riAR[1:10,2],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) write.csv(riAR,"C:/Users/pablo/OneDrive/Escritorio/TFG/RangARs.csv", row.names = TRUE) #Modelo INAR(1) #Función Verosimilitud y estimación inar <- function(param){ lv <- 0 for (t in 2:length(x)){ if(weekdays(d[t])=="lunes"){ xt <- min(x[t-1],x[t]) bin <- 0 for (i in 0:xt ){ bin <- bin + (choose(x[t-1],i) * param[1]^i *(1param[1])^(x[t-1] - i ))* ( exp(- param[2]*pob/10000)* (param[2]*pob/10000)^(x[t]- i ) ) / factorial(x[t]- i ) } lv <- lv + log(bin) } if(weekdays(d[t])=="martes"){ xt <- min(x[t-1],x[t]) bin <- 0 for (i in 0:xt ){ bin <- bin + (choose(x[t-1],i) * param[1]^i *(1param[1])^(x[t-1] - i ))* ( exp(- param[3]*pob/10000)* (param[3]*pob/10000)^(x[t]- i ) ) / factorial(x[t]- i ) } 60 z1<- cbind(z2,deT) z1 <- data.frame("Alpha1"=z2[,1],"Alpha2"= z2[,2],"Lambda"=z2[,3], "SD Alpha1"=z2[,4],"SD Alpha2"=z2[,5],"SD Lambda"=z2[,6],row.names = unique(dadesBCN32$SectorSanitariDescripcio)) write.csv(z1,"C:/Users/pablo/OneDrive/Escritorio/TFG/EstimacionInar.csv", row.names = TRUE) #Predicciones probfunc <- function(xn, prob, p1, lamb, pob){ xt <- min(xn, prob) bin <- 0 for (i in 0:xt){ bin <- bin + (choose(xn,i)*p1^i*(1-p1)^(xn-i))*(exp(-lamb*pob/10000)*(lamb*pob/10000)^(probi))/factorial(prob-i) } return(bin) } pred <- function(xn, p1, lambda, pob,num){ mediana <- vector() li <- vector() ls <- vector() for (npred in 1:num){ sumprob <- 0 k <- 0 while (sumprob < 0.50){ probk <- probfunc(xn, k, p1[npred], lambda[npred], pob) sumprob <- probk + sumprob mediana[npred] <- k k <- k + 1 } sumprob <- 0 k <- 0 while (sumprob < 0.95){ probk <- probfunc(xn, k, p1[npred], lambda[npred], pob) sumprob <- probk + sumprob if (sumprob < 0.025){ li[npred] <- k 61 } k <- k + 1 } k <- 0 sumprob <- 0 while (sumprob < 0.95){ probk <- probfunc(xn, k, p1[npred], lambda[npred], pob) sumprob <- probk + sumprob ls[npred] <- k k <- k + 1 } } li[is.na(li)] <- 0 resultat=list(LowerLim=li, Median=mediana, UpperLim=ls) return(resultat) } #Ex-Ante homes <- vector() for (i in 1:length(table(dadesBCN4$SectorSanitariDescripcio))){ nou_valor <- dadesBCN4$TCasos[dadesBCN4$SexeCodi=="Home" & dadesBCN4$SectorSanitariDescripcio==names(table(dadesBCN4$SectorSanitariDescripcio))[i] & dadesBCN4$TipusCasData=="2021-03-03"] homes <- c(homes, nou_valor) } predH_ante <- list() for (j in 1:length(table(dadesBCN4$SectorSanitariDescripcio))){ poblacio <- unique(dadesBCN4$Poblacio[dadesBCN4$SectorSanitariDescripcio==names(table(dadesBCN4$SectorSani tariDescripcio))[j]]) lambda <- vector() for (i in 1:20){ if (i %in% seq(1, 20, 7)) lambda[i] <- pars[[j]][1, 5] if (i %in% seq(2, 20, 7)) lambda[i] <- pars[[j]][1, 6] if (i %in% seq(3, 20, 7)) lambda[i] <- pars[[j]][1, 7] 62 if (i %in% seq(4, 20, 7)) lambda[i] <- pars[[j]][1, 8] if (i %in% seq(5, 20, 7)) lambda[i] <- pars[[j]][1, 2] if (i %in% seq(6, 20, 7)) lambda[i] <- pars[[j]][1, 3] if (i %in% seq(7, 20, 7)) lambda[i] <- pars[[j]][1, 4] } f <- vector() g <- vector() l <- vector() f[1] <- pars[[j]][1, 1] l[1] <- lambda[1] for (i in 2:20){ f[i] <- f[1]*f[i-1] l[i] <- f[1]*l[i-1] + lambda[i] } predH_ante[[j]] <- pred(homes[j], f, l, poblacio,20) } #Dones dones <- vector() for (i in 1:length(table(dadesBCN4$SectorSanitariDescripcio))){ nou_valor <- dadesBCN4$TCasos[dadesBCN4$SexeCodi=="Dona" & dadesBCN4$SectorSanitariDescripcio==names(table(dadesBCN4$SectorSanitariDescripcio))[i] & dadesBCN4$TipusCasData=="2021-03-03"] dones <- c(dones, nou_valor) } predD_ante <- list() for (j in 1:length(table(dadesBCN4$SectorSanitariDescripcio))){ poblacio <- unique(dadesBCN4$Poblacio[dadesBCN4$SectorSanitariDescripcio==names(table(dadesBCN4$SectorSani tariDescripcio))[j]]) lambda <- vector() s for (i in 1:20){ if (i %in% seq(1, 20, 7)) lambda[i] <- pars[[j]][1, 5] if (i %in% seq(2, 20, 7)) lambda[i] <- pars[[j]][1, 6] 63 if (i %in% seq(3, 20, 7)) lambda[i] <- pars[[j]][1, 7] if (i %in% seq(4, 20, 7)) lambda[i] <- pars[[j]][1, 8] if (i %in% seq(5, 20, 7)) lambda[i] <- pars[[j]][1, 2] if (i %in% seq(6, 20, 7)) lambda[i] <- pars[[j]][1, 3] if (i %in% seq(7, 20, 7)) lambda[i] <- pars[[j]][1, 4] } f <- vector() g <- vector() l <- vector() f[1] <- pars[[j]][1, 1] l[1] <- lambda[1] for (i in 2:20){ f[i] <- f[1]*f[i-1] l[i] <- f[1]*l[i-1] + lambda[i] } predD_ante[[j]] <- pred(dones[j], f, l, poblacio,20) } par(mfrow=c(2,5)) for (s in levels(dadesBCN4$SexeCodi)){ df2<-dadesBCN3[dadesBCN3$SexeCodi==s,] for(i in 1:10){ casos <- df2[df2$SectorSanitariDescripcio==unique(df2$SectorSanitariDescripcio)[i],] x<-casos$TipusCasData[349:368] y<-seq.Date(as.Date("2021-03-04"),as.Date("2021-03-23"),1) if (s=="Home"){ 64 plot(casos$Casos[349:368],col="cadetblue1", main=c('Prediccion ExAnte',s,substring(unique(dadesBCN3$SectorSanitariDescripcio)[i],first=10)),type="lines",ylab=" Número Casos",xlab="Fecha",xaxt="n",xlim=c(0,40),ylim=c(0,max(casos$Casos[349:368]+10))) lines(seq(21,40,1),predH_ante[[i]]$Median, col="red" ) lines(seq(21,40,1),predH_ante[[i]]$LowerLim, col="darkgreen") lines(seq(21,40,1),predH_ante[[i]]$UpperLim, col="darkgreen") axis(1,at=seq(1,40,1),labels=c(x,y)) }else{ plot(casos$Casos[349:368],col="cadetblue1", main=c('Prediccion ExAnte',s,substring(unique(dadesBCN3$SectorSanitariDescripcio)[i],first=10)),type="lines",ylab=" Número Casos",xlab="Fecha",xaxt="n",xlim=c(0,40),ylim=c(0,max(casos$Casos[349:368]+10))) lines(seq(21,40,1),predD_ante[[i]]$Median, col="red" ) lines(seq(21,40,1),predD_ante[[i]]$LowerLim, col="darkgreen") lines(seq(21,40,1),predD_ante[[i]]$UpperLim, col="darkgreen") axis(1,at=seq(1,40,1),labels=c(x,y)) }} } #Ex-Post CV <- dadesBCN3 [dadesBCN3$SectorSanitariDescripcio == "BARCELONA CIUTAT VELLA",] x <- ts(CV[CV$SexeCodi=="Home",]$Casos[1:306]) d <- unique(CV$TipusCasData)[1:306] pob<-poblacio[1] lmvCVH <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deCVH<- sqrt(diag(solve(lmvCVH$hessian))) x <- ts(CV[CV$SexeCodi=="Dona",]$Casos[1:306]) lmvCVD <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deCVD<- sqrt(diag(solve(lmvCVD$hessian))) EIX <- dadesBCN3 [dadesBCN3$SectorSanitariDescripcio == "BARCELONA EIXAMPLE",] x <- ts(EIX[EIX$SexeCodi=="Home",]$Casos[1:306]) d <- unique(EIX$TipusCasData)[1:306] pob <- poblacio[2] lmvEIXH <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deEIXH<- sqrt(diag(inv(lmvEIXH$hessian))) 65 x <- ts(EIX[EIX$SexeCodi=="Dona",]$Casos[1:306]) lmvEIXD <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deEIXD<- sqrt(diag(inv(lmvEIXD$hessian))) GRACIA <- dadesBCN3 [dadesBCN3$SectorSanitariDescripcio == "BARCELONA GRÀCIA",] x <- ts(GRACIA[GRACIA$SexeCodi=="Home",]$Casos[1:306]) d <- unique(GRACIA$TipusCasData)[1:306] pob<-poblacio[3] lmvGRACIAH <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deGRACIAH <- sqrt(diag(inv(lmvGRACIAH$hessian))) x <- ts(GRACIA[GRACIA$SexeCodi=="Dona",]$Casos[1:306]) lmvGRACIAD <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deGRACIAD <- sqrt(diag(inv(lmvGRACIAD$hessian))) HG <- dadesBCN3 [dadesBCN3$SectorSanitariDescripcio == "BARCELONA HORTA-GUINARDÓ",] x <- ts(HG[HG$SexeCodi=="Home",]$Casos[1:306]) d <- unique(HG$TipusCasData)[1:306] pob<-poblacio[4] lmvHGH <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deHGH <- sqrt(diag(inv(lmvHGH$hessian))) x <- ts(HG[HG$SexeCodi=="Dona",]$Casos[1:306]) lmvHGD <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deHGD <- sqrt(diag(inv(lmvHGD$hessian))) LESC <- dadesBCN3 [dadesBCN3$SectorSanitariDescripcio == "BARCELONA LES CORTS",] x <- ts(LESC[LESC$SexeCodi=="Home",]$Casos[1:306]) d <- unique(LESC$TipusCasData)[1:306] pob<-poblacio[5] lmvLESCH <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deLESCH <- sqrt(diag(inv(lmvLESCH$hessian))) x <- ts(LESC[LESC$SexeCodi=="Dona",]$Casos[1:306]) lmvLESCD <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deLESCD <- sqrt(diag(inv(lmvLESCD$hessian))) 66 NB <- dadesBCN3 [dadesBCN3$SectorSanitariDescripcio == "BARCELONA NOU BARRIS",] x <- ts(NB[NB$SexeCodi=="Home",]$Casos[1:306]) d <- unique(NB$TipusCasData)[1:306] pob<-poblacio[6] lmvNBH <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deNBH <- sqrt(diag(inv(lmvNBH$hessian))) x <- ts(NB[NB$SexeCodi=="Dona",]$Casos[1:306]) lmvNBD <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deNBD <- sqrt(diag(inv(lmvNBD$hessian))) SA <- dadesBCN3 [dadesBCN3$SectorSanitariDescripcio == "BARCELONA SANT ANDREU",] x <- ts(SA[SA$SexeCodi=="Home",]$Casos[1:306]) d <- unique(SA$TipusCasData)[1:306] pob<-poblacio[7] lmvSAH <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deSAH <- sqrt(diag(inv(lmvSAH$hessian))) x <- ts(SA[SA$SexeCodi=="Dona",]$Casos[1:306]) lmvSAD <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deSAD <- sqrt(diag(inv(lmvSAD$hessian))) SMRT <- dadesBCN3 [dadesBCN3$SectorSanitariDescripcio == "BARCELONA SANT MARTÍ",] x <- ts(SMRT[SMRT$SexeCodi=="Home",]$Casos[1:306]) d <- unique(SMRT$TipusCasData)[1:306] pob<-poblacio[8] lmvSMRTH <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deSMRTH <- sqrt(diag(inv(lmvSMRTH$hessian))) x <- ts(SMRT[SMRT$SexeCodi=="Dona",]$Casos[1:306]) lmvSMRTD <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deSMRTD <- sqrt(diag(inv(lmvSMRTD$hessian))) SM <- dadesBCN3 [dadesBCN3$SectorSanitariDescripcio == "BARCELONA SANTS-MONTJUÏC",] 67 x <- ts(SM[SM$SexeCodi=="Home",]$Casos[1:306]) d <- unique(SM$TipusCasData)[1:306] pob<-poblacio[9] lmvSMH <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deSMH <-sqrt(diag(solve(lmvSMH$hessian))) x <- ts(SM[SM$SexeCodi=="Dona",]$Casos[1:306]) lmvSMD <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deSMD <-sqrt(diag(solve(lmvSMD$hessian))) SSG <- dadesBCN3 [dadesBCN3$SectorSanitariDescripcio == "BARCELONA SARRIÀ-SANT GERVASI",] x <- ts(SSG[SSG$SexeCodi=="Home",]$Casos[1:306]) d <- unique(SSG$TipusCasData)[1:306] pob<-poblacio[10] lmvSSGH <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deSSGH <- sqrt(diag(inv(lmvSSGH$hessian))) x <- ts(SSG[SSG$SexeCodi=="Dona",]$Casos[1:306]) lmvSSGD <- nlm(inar, p=c(0.5,1,1,1,1,1,1,1),hessian=TRUE) deSSGD <- sqrt(diag(inv(lmvSSGD$hessian))) peHpost <- c(lmvCVH$estimate, lmvEIXH$estimate,lmvGRACIAH$estimate,lmvHGH$estimate,lmvLESCH$estimate,lmvNBH$estimate,lmvSAH$ estimate,lmvSMRTH$estimate,lmvSMH$estimate,lmvSSGH$estimate) peHpost <- matrix(peHpost,ncol=8,byrow = TRUE) peHpost <- data.frame(Alpha=peHpost[,1],Lambda1=peHpost[,2],Lambda2=peHpost[,3],Lambda3=peHpost[,4],Lambd a4=peHpost[,5],Lambda5=peHpost[,6],Lambda6=peHpost[,7],Lambda7=peHpost[,8],row.names= unique(dadesBCN3$SectorSanitariDescripcio)) write.csv(peHpost,"C:/Users/pablo/OneDrive/Escritorio/TFG/ParametresHPostINAR.csv", row.names = TRUE) deTHpost <- c(deCVH,deEIXH,deGRACIAH,deHGH,deLESCH,deNBH,deSAH,deSMRTH,deSMH,deSSGH) deTHpost <- matrix(deTHpost,byrow=TRUE,ncol=8) deTHpost <- data.frame("Alpha error"=deTHpost[,1],"Lambda1 error"= deTHpost[,2],"Lambda2 error"=deTHpost[,3],"Lambda3 error"=deTHpost[,4],"Lambda4 error"=deTH[,5],"Lambda5 error"=deTH[,6],"Lambda6 error"=deTH[,7],"Lambda7 error"=deTH[,8],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) write.csv(deTHpost,"C:/Users/pablo/OneDrive/Escritorio/TFG/DesviacioHPostINAR.csv", row.names = TRUE) 68 peDpost <- c(lmvCVD$estimate, lmvEIXD$estimate,lmvGRACIAD$estimate,lmvHGD$estimate,lmvLESCD$estimate,lmvNBD$estimate,lmvSAD$ estimate,lmvSMRTD$estimate,lmvSMD$estimate,lmvSSGD$estimate) peDpost <- matrix(peDpost,ncol=8,byrow=TRUE) peDpost <- data.frame(Alpha=peDpost[,1],Lambda1=peDpost[,2],Lambda2=peDpost[,3],Lambda3=peDpost[,4],Lambd a4=peDpost[,5],Lambda5=peDpost[,6],Lambda6=peDpost[,7],Lambda7=peDpost[,8],row.names= unique(dadesBCN3$SectorSanitariDescripcio)) write.csv(peDpost,"C:/Users/pablo/OneDrive/Escritorio/TFG/ParametresDPostINAR.csv", row.names = TRUE) deTDpost <- c(deCVD,deEIXD,deGRACIAD,deHGD,deLESCD,deNBD,deSAD,deSMRTD,deSMD,deSSGD) deTDpost <- matrix(deTDpost,byrow=TRUE,ncol=8) deTDpost <- data.frame("Alpha error"=deTDpost[,1],"Lambda1 error"= deTDpost[,2],"Lambda2 error"=deTDpost[,3],"Lambda3 error"=deTDpost[,4],"Lambda4 error"=deTDpost[,5],"Lambda5 error"=deTDpost[,6],"Lambda6 error"=deTDpost[,7],"Lambda7 error"=deTDpost[,8],row.names = unique(dadesBCN3$SectorSanitariDescripcio)) write.csv(deTHpost,"C:/Users/pablo/OneDrive/Escritorio/TFG/desviacioDPostINAR.csv", row.names = TRUE) # Homes homes <- vector() for (i in 1:length(table(dadesBCN4$SectorSanitariDescripcio))){ nou_valor <- dadesBCN4$TCasos[dadesBCN4$SexeCodi=="Home" & dadesBCN4$SectorSanitariDescripcio==names(table(dadesBCN4$SectorSanitariDescripcio))[i] & dadesBCN4$TipusCasData=="2020-12-31"] homes <- c(homes, nou_valor) } predH_post <- list() for (j in 1:length(table(dadesBCN4$SectorSanitariDescripcio))){ poblacio <- unique(dadesBCN4$Poblacio[dadesBCN4$SectorSanitariDescripcio==names(table(dadesBCN4$SectorSani tariDescripcio))[j]]) lambda <- vector() for (i in 1:62){ if (i %in% seq(1, 62, 7)) lambda[i] <- peHpost[j,6] if (i %in% seq(2, 62, 7)) lambda[i] <- peHpost[j,7] if (i %in% seq(3, 62, 7)) lambda[i] <- peHpost[j,8] if (i %in% seq(4, 62, 7)) lambda[i] <- peHpost[j,2] if (i %in% seq(5, 62, 7)) lambda[i] <- peHpost[j,3] 69 if (i %in% seq(6, 62, 7)) lambda[i] <- peHpost[j,4] if (i %in% seq(7, 62, 7)) lambda[i] <- peHpost[j,5] } f <- vector() g <- vector() l <- vector() f[1] <- pars[[j]][1, 1] l[1] <- lambda[1] for (i in 2:62){ f[i] <- f[1]*f[i-1] l[i] <- f[1]*l[i-1] + lambda[i] } predH_post[[j]] <- pred(homes[j], f, l, poblacio,62) } #Dones dones <- vector() for (i in 1:length(table(dadesBCN4$SectorSanitariDescripcio))){ nou_valor <- dadesBCN4$TCasos[dadesBCN4$SexeCodi=="Dona" & dadesBCN4$SectorSanitariDescripcio==names(table(dadesBCN4$SectorSanitariDescripcio))[i] & dadesBCN4$TipusCasData=="2020-12-31"] dones <- c(dones, nou_valor) } predD_post <- list() ### Cada element de la llista és un vector amb les prediccions per cada barri for (j in 1:length(table(dadesBCN4$SectorSanitariDescripcio))) { poblacio <- unique(dadesBCN4$Poblacio[dadesBCN4$SectorSanitariDescripcio==names(table(dadesBCN4$SectorSani tariDescripcio))[j]]) lambda <- vector() # Les prediccions comencen en divendres for (i in 1:62) { if (i %in% seq(1, 62, 7)) lambda[i] <- peDpost[j,6] if (i %in% seq(2, 62, 7)) lambda[i] <- peDpost[j,7]