Peras y Manzanas

En esta ocasión seré mucho más breve, utilizaremos la misma metodología que aquí, pero adaptado a un dataset diferente, para este reporte nos basaremos en el reporte técnico diario sobre los casos confirmados de COVID-19 que publica el gobierno mexicano con corte a las 13hrs.

En este dataset reportan la fecha de inicio de los síntomas, pero no la fecha de confirmación; debido a las limitantes de medición, aparece un aplanamiento al final de la curva, por lo que para hacer el modelaje, solo se considerara el 75% de los datos (rango que podemos considerar como curva verdadera).

La diferencia de reportar los datos de esta forma, en contraste con reportar la fecha de confirmación, es que nos da un panorama más realista y formal del comportamiento de la pandemia, mientras reportar la fecha de confirmación es útil para rastrear bajones en la capacidad de medición, con reportar fecha de inicio de los síntomas se corrige este artefacto, a costa de perder el 35% final de los datos.

Para entender mejor:

Ahora, poniendo días del calendario:

La motivación surge de que muchos nos estamos enfocando en los Dashboards que reportan en base a fecha de confirmación, los cuales son muy útiles a gran escala, pero dan la falsa impresión que estos son útiles a escala nacional o para los especialistas en la materia encargados del sistema de salud del país, y hace que algunas personas no especialistas aprovechen la desinformación con una clara tendencia política y hacen fearmongering, minimizando el valor del juicio de los que realmente son especialistas. No hago esto para defender a nadie, ni me autodenomino especialista, al contrario, si los que leen este documento encuentran algo que no cuadra, espero que les sea útil para desmentir cosas.

Lamentablemente el reporte que está disponible al público está en formato PDF, lo que limita la automatización de la recuperación de datos en tiempo real ligeramente, por el momento no tengo un repositorio en GitHub donde pueda vaciar diario las actualizaciones, por lo que voy a dejar el script comentado que tengo para recuperar los datos y hacer la tabla de incidencias.

En el dataset se reportan también las procedencias, para hacer el análisis vamos a usar solo los datos de “contacto”, es decir, aquellos casos que no fueron importados.

Librerías que vamos a usar:

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

Como extraer los datos del PDF en linea, sin descargarlo:

El Script esta comentado porque la sección donde se hace uso del GUI no compila correctamente usando Rmarkdown. Pero funciona perfectamente si lo copias a un script y quitas el # de la parte donde esta el comando.

# El URL cambia a diario
#goburl<-"https://www.gob.mx/cms/uploads/attachment/file/544957/Tabla_casos_positivos_COVID-19_resultado_InDRE_2020.04.02.pdf"

# El rango de paginas del documento tambien cambia a diario
#pgs<-c(1:32)

# Este comando permite seleccionar el area de las tablas del PDF mediante un GUI
# Se recomienda solo seleccionar el area de las tablas, sin nombres de columnas.

#v<-locate_areas(goburl, 
#                pages = pgs)

# Extraemos las tablas del documento PDF con los parametros actualizados,
# y los ponemos en una lista de matrices.
#out <- extract_tables(goburl, 
#                      pages = pgs,
#                      output = "matrix",
#                      method = "stream",
#                      encoding="latin1",
#                      guess = FALSE,
#                      area = v)

#Se recomienda hacer una revision manual para detectar errores.

# Convertimos cada elemento de la lista en un data.frame

#out.1<-lapply(out,function(x) data.frame(x,
#                                         stringsAsFactors = FALSE))

# UDF para reasignar nombres de columnas

#crename<-function(x,y){
#  colnames(x)<-y
#  return(x)
#}

#Nombres de columnas
#cnames<-c("n_caso",
#          "estado",
#          "sexo",
#          "edad",
#          "fecha_i_s",
#          "RT_PCR",
#          "procedencia",
#          "fecha_arr_mex")

# Iterar UDF sobre elementos de la lista
#out.2<-lapply(out.1,function(x) crename(x,
#                                        cnames))

# Generamos una unica tabla uniendo todos los elementos de la lista
# Mediante nombre de columna
#mexcovid<-bind_rows(out.2)

# Filtramos los casos de contacto
#mexcovid<-mexcovid[mexcovid$procedencia=="Contacto",]

# Generamos tabla de incidencias diarias
#y<-as.matrix((table(mexcovid$fecha_i_s)))

#d<-row.names(y)
#incc<-as.numeric(y[,1])

#mexcovid_is<-data.frame(Date=dmy(as.character(d)),
#                        incident_cases=incc,
#                        stringsAsFactors = FALSE)

# Ordenamos por fecha
#mexcovid_is<-mexcovid_is[order(mexcovid_is$Date),]

# Consolidamos fechas
# md<-min(mexcovid_is$Date)
# Md<-max(mexcovid_is$Date)
# t <- 1:as.integer(ymd(Md) - (ymd(md)-days(1)))

# mxcov<-data.frame(Date = ymd(md)-days(1) + days(t))

# mexcovid_is<-merge(mxcov,mexcovid_is,by="Date",all.x = TRUE) %>%
#   mutate(incident_cases = replace_na(incident_cases, 0))


# Generamos incidencias acumuladas
#mexcovid_is$cumulative_cases<-cumsum(mexcovid_is$incident_cases)

# Generamos proporciones acumuladas (para localizar el 0.75 de los datos)
#mexcovid_is$p_dat<-mexcovid_is$cumulative_cases/sum(mexcovid_is$incident_cases)

# Guardamos las tablas en el directorio de trabajo
#write.table(mexcovid,
#            "COVID19MEX_daily_incidence_by_symptoms_date_confirmed.txt",
#            sep = "\t",
#            row.names = FALSE,
#            quote = FALSE)
#write.table(mexcovid_is,
#            "COVID19MEX_cumulative_incidence_by_symptoms_date_confirmed.txt",
#            sep = "\t",
#            row.names = FALSE,
#            quote = FALSE)
# Para consolidar el codigo:
df<-mexcovid_is

# Filtramos datos cuya proporcion acumulada sea menor a 0.76

df<-df[df$p_dat<0.76,]

tail(df)
##          Date incident_cases cumulative_cases     p_dat
## 31 2020-03-18             49              291 0.3277027
## 32 2020-03-19             54              345 0.3885135
## 33 2020-03-20             75              420 0.4729730
## 34 2020-03-21             63              483 0.5439189
## 35 2020-03-22             57              540 0.6081081
## 36 2020-03-23             82              622 0.7004505

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-17"
#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   1   2   2   2   2   3   3   3   4   4   4   4   4   5   6   6   8  10
## [20]  13  16  24  35  43  56  88 117 151 191 242 291 345 420 483 540 622

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.6 y gamma=0.4.

Opt <- optim(c(0.6, 0.4), #(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.5928624 0.4071377

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.456172

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.

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:

Ahora incidencias acumuladas (75% de los datos):

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 se ajusta bastante bien.

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: 5128"
## [1] "Para hoy 2020-04-03 El número de incidencias nuevas ajustadas para COVID-19 en México sera de: 869 nuevas infecciones respecto a ayer: 2020-04-02"

Ahora, vamos a predecir cuanta gente infectada va a haber en 35 días, dada la de duplicación que tenemos hasta hoy, 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:

Disclaimers: