Seré breve, me dio curiosidad saber cuál sería la tendencia de infección del COVID-19 en México de hoy a algunos días (o meses), teniendo como escenarios de contagio:

Estoy consciente que ya existen modelos matemáticos que intentan predecir estas tendencias, y uno que otro perdido que intenta comparar con España a ojo de buen cubero, sin hacer pruebas de bondad de ajuste (quizás en algún momento haga un ejercicio de machine learning para ver cómo utilizar adecuadamente los datos de otros países), pero creí que sería un ejercicio interesante de exponer.

La motivación proviene de que este método no se encuentra explicado en español, ni en un lenguaje lo suficientemente simple para seguir el rollo.

Este no es un curso, pero si ya le sabes a R puedes usar este script para aplicarlo a tus análisis, no asumo nada de los que vayan a leer esto.

Omitiré casi todos los detalles estadísticos del proceso, ya que es un ejercicio de divulgación, no uno técnico, para los que les dé curiosidad los detalles matemáticos y programáticos:

Librerías que vamos a usar:

library(lubridate)
library(tidyverse)
library(deSolve)
library(plyr)
library(tidyr)
library(magrittr)
library(dplyr)
library(readr)
library(ggplot2)

Obtener datos en tiempo real:

Queremos ahorrarnos la molestia de descargar los datos cada vez que queramos analizarlos. Mejor automatizamos el proceso:

Los datos provienen del repositorio en github de la Johns Hopkins University.

#URL de los datos (github de la JHU)
jhu_url <- paste("https://raw.githubusercontent.com/CSSEGISandData/", 
                 "COVID-19/master/csse_covid_19_data/", "csse_covid_19_time_series/", 
                 "time_series_covid19_confirmed_global.csv", sep = "")

#Leemos los datos directamente de la fuente para tenerlos al dia

# Nos libramos de malos formatos de nombres de columnas
df <- read_csv(jhu_url) %>% 
  dplyr::rename("province" = "Province/State",
         "country_region" = "Country/Region") %>%
 
   #Transformar datos de "time series" a formato largo
  pivot_longer(-c(province,
                  country_region, 
                  Lat, 
                  Long), 
               names_to = "Date", 
               values_to = "cumulative_cases") %>% 
  
  # Ajustar formato de fechas
  mutate(Date = mdy(Date)) %>% 
  
  #Filtrar datos de México (funciona para cualquier otro país).
  filter(country_region =="Mexico") %>% 
  arrange(country_region, Date) %>% 
  group_by(country_region) %>% 
  
  #Calcular datos de incidencia a partir de los datos de incidencia acumulativa
  mutate(incident_cases = c(0, 
                            diff(cumulative_cases))) %>% 
  ungroup() %>% 
  
  # Quitar datos de provincia, latitud y longitud.
  # Estos son utiles si quieres hacer mapas,
  # O analizar los datos de provincia para los paises
  # que tienen disponible esta granularidad.
  select(-c(province, Lat, Long))

#Y asi quedan
tail(df)
## # A tibble: 6 x 4
##   country_region Date       cumulative_cases incident_cases
##   <chr>          <date>                <dbl>          <dbl>
## 1 Mexico         2020-03-29              848            131
## 2 Mexico         2020-03-30              993            145
## 3 Mexico         2020-03-31             1094            101
## 4 Mexico         2020-04-01             1215            121
## 5 Mexico         2020-04-02             1378            163
## 6 Mexico         2020-04-03             1510            132

Ahora la tabla de Muertos

#URL de los datos (github de la JHU)
jhu_url <- paste("https://raw.githubusercontent.com/CSSEGISandData/", 
                 "COVID-19/master/csse_covid_19_data/", "csse_covid_19_time_series/", 
                 "time_series_covid19_deaths_global.csv", sep = "")

#Leemos los datos directamente de la fuente para tenerlos al dia

# Nos libramos de malos formatos de nombres de columnas
df.d <- read_csv(jhu_url) %>% 
  dplyr::rename("province" = "Province/State",
         "country_region" = "Country/Region") %>%
 
   #Transformar datos de "time series" a formato largo
  pivot_longer(-c(province,
                  country_region, 
                  Lat, 
                  Long), 
               names_to = "Date", 
               values_to = "cumulative_cases") %>% 
  
  # Ajustar formato de fechas
  mutate(Date = mdy(Date)) %>% 
  
  #Filtrar datos de México (funciona para cualquier otro país).
  filter(country_region =="Mexico") %>% 
  arrange(country_region, Date) %>% 
  group_by(country_region) %>% 
  
  #Calcular datos de incidencia a partir de los datos de incidencia acumulativa
  mutate(incident_cases = c(0, 
                            diff(cumulative_cases))) %>% 
  ungroup() %>% 
  
  # Quitar datos de provincia, latitud y longitud.
  # Estos son utiles si quieres hacer mapas,
  # O analizar los datos de provincia para los paises
  # que tienen disponible esta granularidad.
  select(-c(province, Lat, Long))

#Y asi quedan
tail(df.d)
## # A tibble: 6 x 4
##   country_region Date       cumulative_cases incident_cases
##   <chr>          <date>                <dbl>          <dbl>
## 1 Mexico         2020-03-29               16              4
## 2 Mexico         2020-03-30               20              4
## 3 Mexico         2020-03-31               28              8
## 4 Mexico         2020-04-01               29              1
## 5 Mexico         2020-04-02               37              8
## 6 Mexico         2020-04-03               50             13

Si quieres saber que es incidencia cumulativa sigue este link.

Modelo SIR

Construimos función para nuestro modelo Susceptibles-Infectados-Recuperados (SIR). Este modelo nos va a ayudar a predecir la tendencia de la incidencia, mediante un ajuste a la misma. Este modelo es muy usado en epidemiología. Para ver detalles de este y otros modelos usados en epidemiología sigue este link. O los que listé al inicio del documento.

SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta * I * S/N
    dI <- beta * I * S/N - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

Sacamos la incidencia acumulada desde el día en el que se dio el primer caso, hasta el dato más actualizado.

#Fecha de inicio
sir_start_date <- "2020-02-28"
#Fecha actualizada
sir_end_date<- max(df$Date)
#Sacamos los datos de incidencia acumulada
Infected <- df %>% 
  filter(Date >= ymd(sir_start_date), Date <= ymd(sir_end_date)) %>% 
  pull(cumulative_cases)
#Y asi queda
Infected
##  [1]    1    4    5    5    5    5    5    6    6    7    7    7    8   12   12
## [16]   26   41   53   82   93  118  164  203  251  316  367  405  475  585  717
## [31]  848  993 1094 1215 1378 1510

Creamos un vector de incrementos de día de la misma longitud que nuestro vector de incidencias acumuladas, lo vamos a usar más adelante

Day <- 1:(length(Infected))

Vamos a especificar los valores iniciales para S, I y R:

# En este caso, vamos a suponer que la población suceptible va a ser el 70%
# de la población urbana de México.

# N= ((población total en México)*proporción de población urbana)*proporción suceptible
N<-((128932753)*0.74)*0.7


# Los guardamos en un vector de inicialización:
init <- c(S = N - Infected[1], I = Infected[1], R = 0)

Suma de cuadrados residuales (RSS)

Definimos una función para calcular la suma de cuadrados residuales (RSS), al cual pasamos los parámetros beta y gamma que serán optimizados para ajustarse a los datos de incidencia acumulada.

RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[, 3]
  sum((Infected - fit)^2)
}

Ahora, vamos a encontrar los valores de beta y gamma que den como resultado el RSS más pequeño mediante un proceso de optimización, el cual representa el mejor ajuste a nuestros datos. Se sugiere empezar con valores de 0.5 para cada parámetro, pero si es necesario probar con otros, la suma de ambos debe dar 1, este proceso se hace hasta que exista convergencia. En este caso el RSS más pequeño se encontró cuando los valores iniciales de beta=0.5 y gamma=0.5.

Opt <- optim(c(0.5, 0.5), #(beta,gamma)
             RSS, #Función a optimizar
             method = "L-BFGS-B", 
             lower = c(0,0), #Limite minimo
             upper = c(1, 1)) #Limite maximo

# Hay que ver si hay convergencia
Opt$message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
#GUardamos los parametros optimos en Opt_par
Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par
##      beta     gamma 
## 0.6072464 0.3927536

Tasa de contagio (R0)

Estos valores por si mismos no nos dicen nada, pero el cosiente es R0, lo que es la tasa de contagio.

# R0= beta/gamma
RO<-Opt_par[1]/Opt_par[2]

#R0=
as.numeric(RO)
## [1] 1.546126

De momento, la OMS estima que la tasa de contagio del virus es de 1.4 a 2.5, así que por ahora la estimación cae dentro del rango esperado. No encontré ninguna fuente que corroborara la tasa de contagio para COVID-19 en México.

Ajuste de modelo SIR, predicciones y visualizaciones

Ahora vamos a ver como se ajusta el modelo, veremos las predicciones del modelo hasta el día de hoy:

# Tiempo en dias para el ajuste del modelo
t <- 1:as.integer(today() - (ymd(sir_start_date)-days(1)))

# Obtener valores ajustados mediante nuestro modelo SIR con tiempo t y parametros optimizados.
fitted_cumulative_incidence <- data.frame(ode(y = init, 
                                              times = t, 
                                              func = SIR, 
                                              parms = Opt_par))
# Vamos a unir la tabla del modelo ajustado de incidencias y
# la tabla de incidencias observadas mediante la columna de fechas
fitted_cumulative_incidence<-fitted_cumulative_incidence %>% 
  mutate(Date = ymd(sir_start_date)-days(1) + days(t))

df.cc<-select(df,Date, cumulative_cases)

fitted_cumulative_incidence<-merge(fitted_cumulative_incidence,
                                   df.cc,
                                   by="Date",
                                   all.x = TRUE)

Visualizamos Gráficos:

Primero el de incidencia de casos por día:

Del primero al 7 de marzo ocurrieron los primeros casos de importación.

Se observa que a mediados de Marzo se empieza a notar el inicio de la fase exponencial, característica de el contagio comunitario, es decir, que ya se dan casos que no se contagiaron fuera de México, si no dentro.

Ahora veamos el gráfico de Incidencias ajustadas con el modelo SIR vs. observadas

La línea punteada azul representa la fecha cuando entramos a Fase 2.

Parece ser que el modelo ya no se ajusta bien, quizas es necesario otro modelo de predicción.

Para hoy, la cantidad de infectados usando este modelo se estima que será:

## [1] "Para hoy 2020-04-03 El número de incidencias acumulativas ajustadas para COVID-19 en México sera de: 1821"
## [1] "Para hoy 2020-04-03 El número de incidencias nuevas ajustadas para COVID-19 en México sera de: 351 nuevas infecciones respecto a ayer: 2020-04-02"

Ahora, vamos a predecir cuanta gente infectada va a haber en 35 días, al mismo tiempo que veremos cuando es probable que se abrume al ya de por si abrumado sistema de salud pública de México, donde:

La línea punteada roja marca la fecha donde entramos a Fase 2.

Como se observa, es probable que el sistema de salud este limitado para finales de abril.

Por eso mencionan que la la primera quincena de abril es la más importante, porque para la segunda ya vamos a estar en un punto incomodo.

Ahora veremos cómo es que se comporta la curva en los diferentes escenarios de susceptibilidad, donde:

Ahora podemos ver como reducir la suceptibilidad ayuda a que la curva sea menos pronunciada.

Espero sea informativo.

Disclaimers: