Modelos jerárquicos para predecir el número de infectados por COVID-19 en América del Sur

A finales de 2019 se identificó en Wuhan, China una gran cantidad de personas enfermas con neumonía, estudios de investigadores permitieron determinar que se trataba de un virus totalmente desconocido, un nuevo coronavirus. Coronavirus es una gran familia de virus conocidos por causar enfermedades respiratorias y gastrointestinales que van desde un resfriado común hasta manifestaciones clínicas más severas como las causadas por el síndrome respiratorio de Oriente Medio (MERS) y el síndrome respiratorio agudo grave (SARS). Este coronavirus llamado 2019 novel coronavirus (COVID-19) es de fácil transmisión y la gravedad de sus síntomas ha causado gran incertidumbre y temor debido a que el virus actualmente se ha propagado por más de 180 países.

Por ello, se propone estudiar una base de datos sobre el COVID-19 actualizada día a día en Kaggle nombrada como \(covid\_19\_data.csv\).

Descripción de la base de datos

Esta es una base de datos con 11341 observaciones y 8 variables que son:

• Sno: número de observación.
• ObservationDate: Fecha de observación MM/DD/AAAA.
• Province/State: Provincia o estado de la observación (Podría estar vacía cuando falta).
• Country/Region: País de la observación.
• Last Update: Hora en UTC en la que se actualiza la fila para la provincia o país dado sin estandarizar.
• Confirmed: Número acumulado de casos confirmados hasta esa fecha.
• Deaths: Número acumulado de muertes hasta esa fecha.
• Recovered: Número acumulado de casos recuperados hasta esa fecha.

Objetivo

Modelar la cantidad de infectados por el virus COVID-19 que se pueden detectar en una fecha teniendo en cuenta la cantidad de casos confirmados y su diferencia entre los muertos y recuperados para cada uno de los países de Sur América con el fin de determinar el incremento de infección

Análisis descriptivo

A continuación se presenta un resumen de la base de datos original:

summary(Covid)

##       SNo        Observation_Date      Province_State        Country_Region
##  Min.   :    1   Min.   :2020-01-22           :3595   US            :4796  
##  1st Qu.: 2836   1st Qu.:2020-03-01   Texas   : 267   Mainland China:1920  
##  Median : 5671   Median :2020-03-15   Georgia : 173   Australia     : 331  
##  Mean   : 5671   Mean   :2020-03-08   Virginia: 146   Canada        : 266  
##  3rd Qu.: 8506   3rd Qu.:2020-03-23   Kentucky: 134   France        : 136  
##  Max.   :11341   Max.   :2020-03-23   Missouri: 130   UK            : 102  
##                                       (Other) :6896   (Other)       :3790  
##               Last_Update     Confirmed           Deaths          Recovered    
##  2020-03-23 23:19:34:3179   Min.   :    0.0   Min.   :   0.00   Min.   :    0  
##  2020-03-11T20:00:00: 197   1st Qu.:    1.0   1st Qu.:   0.00   1st Qu.:    0  
##  2020-03-23 23:19:21: 182   Median :    5.0   Median :   0.00   Median :    0  
##  2020-03-16T14:38:45: 118   Mean   :  491.6   Mean   :  17.52   Mean   :  175  
##  2020-03-15T18:20:18:  91   3rd Qu.:   71.0   3rd Qu.:   0.00   3rd Qu.:    2  
##  2020-03-15T18:20:19:  71   Max.   :67800.0   Max.   :6077.00   Max.   :59882  
##  (Other)            :7503

Donde se puede ver la cantidad de observaciones que tiene la base de datos original, junto a la fecha de inicio de los registros con la última registrada. Además, la variable que contiene la provincia o estado son la única variabla con Na’s. Hay un país donde han detectado más de 67000 casos de COVID-19 y otro donde han muerto un poco más de 6000 personas. Lo cual es bastante desconsolador teniendo en cuenta el corto periodo de tiempo entre el 22 de enero hasta el 23 de marzo de 2020.

Con el fin de poder realizar el análisis jerárquico se realizaron algunos arreglos en la base de datos. Se filtró la base para los países de América del Sur únicamente. Además, se eliminaron variables innecesarias como \(SNo\) y \(Last\_Update\) pero se incluyeron dos nuevas que son: \(Infected\) (número de personas infectadas a la fecha) que se estima de la siguiente manera: \(Infected=Confirmed-Deaths-Recovered\) y \(ID\_Country\) que es el identificador del país asignado por orden alfabético. Por último, se hizo que el origen de detección del virus o fecha 0 coincidiera con la fecha 26 de febrero de 2020 puesto que en esta parte del continente se detecto el primer caso ese día. Por ello, para conocer su estructura se presentan las primeras observaciones de la nueva base de datos:

head(covid)

##   Observation_Date ID_Country Country_Region Confirmed Deaths Recovered
## 1                0          1      Argentina         0      0         0
## 2                1          1      Argentina         0      0         0
## 3                2          1      Argentina         0      0         0
## 4                3          1      Argentina         0      0         0
## 5                4          1      Argentina         0      0         0
## 6                5          1      Argentina         0      0         0
##   Infected
## 1        0
## 2        0
## 3        0
## 4        0
## 5        0
## 6        0

De la base de datos anterior, se puede ver también:

El número de muertes por COVID-19 incrementa a medida que el número de casos confirmados aumenta, en otras palabras, el número de contagiados y la cantidad de personas que mueren son cantidades que aumentan proporcionalmente.

En trece de los dieciséis países de América del Sur, por lo menos, una persona posee el virus actualmente. En seis de ellos hay más de mil casos confirmados a la fecha.

En orden decreciente Brazil, Ecuador y Argentina, tienen un número de muertes que ya supero los veinticinco pero países como Chile, Perú y Colombia mantienen tasas de mortalidad bajas aunque la cantidad de infectados sea alta. Comparando Argentina y Colombia que tienen una cantidad de casos confirmados similares se puede notar una gran diferencia entre la mortalidad en cada país.

Así, como ha tenido más muertes, Argentina, ha podido recuperar más pacientes con COVID-19 que ningun otro país y Venezuela que no tiene tantos casos confirmados ha podido recuperar a más personas con la enfermedad.

Veamos ahora, lo que se desea modelar, en la escala normal:

Por las formas exponenciales visibles en los gráficos se proponen dos procedimientos: el primero, tomar la variable \(Infected\) como el logaritmo y el segundo, la variable \(Infected\) como una respuesta \(Poisson\). Siguiendo el primer procedimiento, la gráfica quedaría:

Realizando un test de Shapiro-Wilk, se determinó que la variable \(Infected\) no se distribuye Normal, por lo tanto, se realiza la segunda propuesta de procedimiento.

Modelos a considerar

Modelo 1

\[Infected_{ij} \sim Poisson(\lambda_{ij})\] \[log(\lambda_{ij}) = \beta_0+\beta_1 Observation\_Date_{ij}+ b_{oi}\] \[b_0 \sim N(0, \sigma_{b_0}^2)\] con un vector de parámetros \(\Theta=(\beta_0, \beta_1, \sigma_{b_0})^T\)

Modelo 2

\[Infected_{ij} \sim Poisson(\lambda_{ij})\] \[log(\lambda_{ij}) = \beta_0+\beta_1 Observation\_Date_{ij}+ b_{oi}+b_{1i} Observation\_Date_{ij}\] \[\begin{pmatrix} b_0 \\ b_1 \end{pmatrix} \sim N \left[ \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_{b_0}^2 & \sigma_{b_{01}}^2 \\ \sigma_{b_{01}}^2 & \sigma_{b_1}^2 \end{pmatrix} \right]\] con un vector de parámetros \(\Theta=(\beta_0, \beta_1, \sigma_{b_0}, \sigma_{b_1}, \sigma_{b_{01}})^T\)

Modelo 3

\[Infected_{ij} \sim Poisson(\lambda_{ij})\] \[log(\lambda_{ij}) = \beta_0+\beta_1 Observation\_Date_{ij} + \beta_2Deaths_{ij} + \beta_3Recovered_{ij} + b_{oi}+b_{1i} Observation\_Date_{ij}\] \[\begin{pmatrix} b_0 \\ b_1 \end{pmatrix} \sim N \left[ \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_{b_0}^2 & \sigma_{b_{01}}^2 \\ \sigma_{b_{01}}^2 & \sigma_{b_1}^2 \end{pmatrix} \right]\] con un vector de parámetros \(\Theta=(\beta_0, \beta_1, \beta_2, \beta_3, \sigma_{b_0}, \sigma_{b_1}, \sigma_{b_{01}})^T\)

Construcción de los modelos

Modelo 1

mod1 <- glmer(Infected ~ Observation_Date + (1 | ID_Country), 
              data=covid, 
              family = poisson(link = "log"),
              nAGQ = 1)
summary(mod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Infected ~ Observation_Date + (1 | ID_Country)
##    Data: covid
## 
##      AIC      BIC   logLik deviance df.resid 
##   2493.7   2505.3  -1243.9   2487.7      348 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -6.9732 -0.8589 -0.3160  0.5786  9.1132 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  ID_Country (Intercept) 2.808    1.676   
## Number of obs: 351, groups:  ID_Country, 13
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.842229   0.467902  -6.074 1.24e-09 ***
## Observation_Date  0.298312   0.002138 139.523  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Observtn_Dt -0.106

Modelo 2

mod2 <- glmer(Infected ~ Observation_Date + (Observation_Date | ID_Country), 
              data=covid, 
              family = poisson(link = "log"),
              nAGQ = 1)
summary(mod2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Infected ~ Observation_Date + (Observation_Date | ID_Country)
##    Data: covid
## 
##      AIC      BIC   logLik deviance df.resid 
##   2228.4   2247.7  -1109.2   2218.4      346 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5692 -0.9436 -0.3969  0.4434  8.5034 
## 
## Random effects:
##  Groups     Name             Variance Std.Dev. Corr 
##  ID_Country (Intercept)      2.492418 1.57874       
##             Observation_Date 0.003678 0.06065  -0.36
## Number of obs: 351, groups:  ID_Country, 13
## 
## Fixed effects:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.25797    0.45646  -4.947 7.55e-07 ***
## Observation_Date  0.27171    0.01773  15.323  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Observtn_Dt -0.413

Modelo 3

mod3 <- glmer(Infected ~ Observation_Date + Deaths + Recovered + (Observation_Date | ID_Country), 
              data=covid, 
              family = poisson(link = "log"),
              nAGQ = 1)
summary(mod3)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Infected ~ Observation_Date + Deaths + Recovered + (Observation_Date |  
##     ID_Country)
##    Data: covid
## 
##      AIC      BIC   logLik deviance df.resid 
##   2208.8   2235.8  -1097.4   2194.8      344 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6087 -0.9337 -0.3963  0.4000  8.6509 
## 
## Random effects:
##  Groups     Name             Variance Std.Dev. Corr 
##  ID_Country (Intercept)      2.449785 1.56518       
##             Observation_Date 0.003699 0.06082  -0.34
## Number of obs: 351, groups:  ID_Country, 13
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -2.336606   0.453561  -5.152 2.58e-07 ***
## Observation_Date  0.276224   0.017833  15.490  < 2e-16 ***
## Deaths            0.001654   0.001591   1.040    0.299    
## Recovered        -0.013902   0.002971  -4.679 2.89e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Obsr_D Deaths
## Observtn_Dt -0.404              
## Deaths       0.028 -0.041       
## Recovered    0.044 -0.063  0.002

Análisis de residuales

Modelo 1

Modelo 2

Modelo 3

Comparación de los modelos

Se puede notar en los residuales de los ajustes anteriores que la varianza aún tiene algunos problemas de homogeneidad. Por ello, se propone un último modelo y corregir dicho problema. El resultado del modelo es el siguiente:

mod4 <- glmer(Infected ~ poly(Observation_Date, 4) + (Observation_Date | ID_Country), 
              data=covid, 
              family = poisson(link = "log"),
              nAGQ = 1)
summary(mod4)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Infected ~ poly(Observation_Date, 4) + (Observation_Date | ID_Country)
##    Data: covid
## 
##      AIC      BIC   logLik deviance df.resid 
##   2154.5   2185.4  -1069.2   2138.5      343 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6401 -0.8734 -0.2990  0.3598 12.2654 
## 
## Random effects:
##  Groups     Name             Variance Std.Dev. Corr 
##  ID_Country (Intercept)      2.716115 1.64806       
##             Observation_Date 0.004479 0.06693  -0.44
## Number of obs: 351, groups:  ID_Country, 13
## 
## Fixed effects:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  1.0121     0.4230   2.393 0.016723 *  
## poly(Observation_Date, 4)1  46.4221     3.3563  13.831  < 2e-16 ***
## poly(Observation_Date, 4)2  -5.4227     1.6450  -3.297 0.000979 ***
## poly(Observation_Date, 4)3   3.2578     1.0861   3.000 0.002704 ** 
## poly(Observation_Date, 4)4  -1.9724     0.4286  -4.602 4.19e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) p(O_D,4)1 p(O_D,4)2 p(O_D,4)3
## ply(O_D,4)1 -0.035                              
## ply(O_D,4)2  0.143 -0.516                       
## ply(O_D,4)3 -0.117  0.444    -0.939             
## ply(O_D,4)4  0.088 -0.342     0.771    -0.921

Al comparar los modelos por una prueba de verosimilitud sencilla, se puede determinar que:

anova(mod1, mod2, mod3, mod4)

## Data: covid
## Models:
## mod1: Infected ~ Observation_Date + (1 | ID_Country)
## mod2: Infected ~ Observation_Date + (Observation_Date | ID_Country)
## mod3: Infected ~ Observation_Date + Deaths + Recovered + (Observation_Date | 
## mod3:     ID_Country)
## mod4: Infected ~ poly(Observation_Date, 4) + (Observation_Date | ID_Country)
##      Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)    
## mod1  3 2493.7 2505.3 -1243.9   2487.7                              
## mod2  5 2228.4 2247.7 -1109.2   2218.4 269.340      2  < 2.2e-16 ***
## mod3  7 2208.8 2235.8 -1097.4   2194.8  23.629      2  7.397e-06 ***
## mod4  8 2154.5 2185.4 -1069.2   2138.5  56.265      1  6.333e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Añadir cada una de las variables a cada modelo consecutivamente a llevado a mejorar el modelo, esto se refleja en el análisis de la correlación de Pearson que se muestra a continuación:

# Correlación de Pearson
# Modelo 1
cor(covid$Infected, predict(mod1, type = "response"))

## [1] 0.9919989

# Modelo 2
cor(covid$Infected, predict(mod2, type = "response"))

## [1] 0.9948642

# Modelo 3
cor(covid$Infected, predict(mod3, type = "response"))

## [1] 0.9949198

# Modelo 4
cor(covid$Infected, predict(mod4, type = "response"))

## [1] 0.9955308

Siendo el último resultado, el correspondiente al modelo 4, eso quiere decir que es un mejor modelo que los otros. Sus residuales también muestran gran mejoría mostrando una dispersión más homogénea:

Presentación del mejor modelo

Para terminar se desea visualizar como quedó el modelo. Se presenta aquí:

Viendo las gráficas parecen satisfactorias puesto que el comportamiento de las curvas por grupo sigue la tendencia de los puntos como se debe. Esto podría ser útil para predecir cantidad de infectados con solo darle al modelo la fecha que queremos predecir. Por ejemplo, el último valor observado en la base de datos para la variable \(Observation\_Date\) es \(26\) el registro corresponde al día 23 de marzo de 2020 que realmente pertenece a la información recolectada durante el día 22 de marzo, tomando como \(Country\_Region\) a Colombia por la base se sabe que ese día habían 269 casos de personas infectadas, el modelo predice para el día 23 de marzo 291 infectados aproximadamente. Haciendo lo mismo para el día 25 de marzo, se predicen 384 infectados, según el Ministerio de Salud en Colombia fueron 458 personas con la enfermedad. Esta diferencia se presenta puesto que solo se tienen 27 observaciones por grupo pero sirve para conocer la propagación del virus en distintas poblaciones, también como motivador para continuar con el periodo de cuarentena y seguir las recomendaciones, tanto como buscar nuevas estrategias para lograr que el virus no se siga propagando a la velocidad que lo hace actualmente. Para que el modelo siga trabajando y mejore en sus predicciones, es recomendable actualizar la base de datos diariamente, o en lo posible, cuando se vaya a usar el modelo.