El modelo SIR tiene por objetivo predecir si una epidemia se va a expandir o va a disminuir. Para ello establece tres funciones dependientes del tiempo que representan:
\(S\left( t \right)\): Variacion del numero de personas susceptibles de contraer la enfermedad en el transcurso del tiempo.
\(I\left( t \right)\): Variacion de la cantidad de infectados en el transcurso del tiempo.
\(R\left( t \right)\): Variacion de la cantidad de recuperados en el transcurso del tiempo.
Y tres parametros que representan:
\(\beta\): El numero medio de contactos por unidad de tiempo.
\(\gamma\): El numero de recuperados por unidad de tiempo dividido por el total de infectados.
\(N\): La poblacion total
El modelo supone que los susceptibles se infectan a cierta tasa y se recuperan a otra (tasas de transicion). Asi:
\[ S\left( t \right) \Rightarrow TDT 1 \Rightarrow I\left( t \right) \Rightarrow TDT 2 \Rightarrow R \left( t \right)\]
La tasa de transicion 1 (TDT1) se define como:\(\frac { \beta I\left(t\right)}{N}\)
La tasa de transicion 2 (TDT2) se define como:\(\gamma\)
El sistema de ecuaciones diferenciales no lineales que representan lo anterior y que definen al modelo son:
\[\frac{dS\left( t \right)}{dt} = - \frac {\beta I\left( t \right) S\left( t \right)}{N} \] \[\frac{dI\left( t \right)}{dt} = \frac {\beta I\left( t \right) S\left( t \right)}{N} - \gamma I\left( t \right)\] \[\frac{dR\left( t \right)}{dt} = \gamma I\left( t \right)\] Observemos que en realidad es un sistema de dos ecuaciones diferenciales. Despejemos \(\gamma I\left( t \right)\) de la segunda ecuacion y reemplacemos el valor en la tercera:
\[\frac{dR\left( t \right)}{dt} = \gamma I\left( t \right) = \frac {\beta I\left( t \right) S\left( t \right)}{N} - \frac{dI\left( t \right)}{dt}\] \[\frac{dR\left( t \right)}{dt} = \frac {\beta I\left( t \right) S\left( t \right)}{N} - \frac{dI\left( t \right)}{dt}\] Reemplacemos el valor del primer termino del lado derecho por la primera ecuacion:
\[\frac{dR\left( t \right)}{dt} = - \frac{dS\left( t \right)}{dt} - \frac{dI\left( t \right)}{dt}\]
Multipliquemos por \(dt\) la expresion e integremos en forma indefinida:
\[ \int dR\left( t \right) = - \int dS\left( t \right) - \int dI\left( t \right)\] \[ R\left( t \right) = - S \left( t \right) - I\left( t \right) +C\] C en igual a N pues: \[ N = R\left( t \right)+ S \left( t \right)+ I\left( t \right) \] Entonces: \[ R\left( t \right) = - S \left( t \right) - I\left( t \right) +N\] Volvamos al sistema de ecuaciones diferenciales originales, despejemos el valor de \(I\left( t \right)\) de la tercera ecuacion y dividamos la primera por el:
\[\frac{dS\left( t \right)}{dt} = - \frac {\beta I\left( t \right) S\left( t \right)}{N} \]
\[\frac{dI\left( t \right)}{dt} = \frac {\beta I\left( t \right) S\left( t \right)}{N} - \gamma I\left( t \right)\]
\[\frac{dR\left( t \right)}{dt} = \gamma I\left( t \right)\]
\[ I\left( t \right) = \frac{dR\left( t \right)}{ \gamma dt } \]
\[ \frac{ \frac{dS\left( t \right)}{dt}} { \frac{dR\left( t \right)}{ \gamma dt }}= \frac{ - \frac {\beta I\left( t \right) S\left( t \right)}{N}} {I\left( t \right)} \]
\[ \frac{dS\left( t \right)}{dt} \frac{ \gamma dt } {dR\left( t \right)} = - \frac {\beta I\left( t \right) S\left( t \right)}{N I\left( t \right)} \] Al lado izquierdo se van los \(dt\) y al lado derecho se van los \(I(t)\).
\[ \frac{ dS\left( t \right) \gamma } {dR\left( t \right)} = - \frac {\beta S\left( t \right)}{N } \]
Reagrupamos.
\[ \frac{dS\left( t \right)}{S\left( t \right)} {} = - \frac {\beta dR\left( t \right)}{\gamma N } \] Integramos definidamente sobre el tiempo entre \(0\) y \(t\):
\[ \int_0^t \frac{dS\left( t \right) }{S\left( t \right)} {} = - \frac {\beta }{\gamma N } \int_0^t dR\left( t \right) \]
\[\ln S(t)\Big|_0^t= - \frac {\beta R\left( t \right) }{\gamma N }\Big|_0^t\] \[\ln S(t)- \ln S(0)= - \frac {\beta }{\gamma N } (R(t)-R(0))\] Lo que es igual a:
\[\ln \frac { S(t)}{ S(0)}= - \frac {\beta }{\gamma N } (R(t)-R(0))\] elevamos a \(e\) y llevamos el denominador del lado izquierdo al derecho:
\[e^{\ln \frac { S(t)}{ S(0)}} = e^{- \frac {\beta }{\gamma N } (R(t)-R(0))}\] \[ S(t) = S(0) e^{- \frac {\beta }{\gamma N } (R(t)-R(0))}\] Por definicion, en el momento \(t=0\):
\(S(0) = N\), y
\(R(0) = 0\), por lo que la ecuacion queda:
\[ S(t) = N e^{- \frac {\beta }{\gamma N } (R(t))}\] La expresion \(\frac {\beta}{\gamma }\) se denomina numero basico de reproduccion \(R_0\) y es fundamental como lo veremos a continuacion:
Reemplacemos el valor de \(R_0\) en la segunda ecuacion diferencial del sistema:
\[\frac{dI\left( t \right)}{dt} = \frac {\beta I\left( t \right) S\left( t \right)}{N} - \gamma I\left( t \right)\]
\[\frac{dI\left( t \right)}{dt} = \frac {\gamma R_0 I\left( t \right) S\left( t \right)}{N} - \gamma I\left( t \right)\]
\[\frac{dI\left( t \right)}{dt} = [\frac { R_0 S\left( t \right)}{N}-1] \gamma I\left( t \right)\] La tasa de propagacion de la infeccion queda definida entonces por la expresion:
\[\frac { R_0 S\left( t \right)}{N}-1\] Si:
\[\frac { R_0 S\left( t \right)}{N}-1 >0 \]
o bien \[\frac { R_0 S\left( t \right)}{N} >1 \]
\[\frac{dI\left( t \right)}{dt} > 0\]
y la tasa de propagacion de la enfermedad crecera.
Por el contrario, si \[\frac { R_0 S\left( t \right)}{N}<1\] decrecera.
Supongamos que: \[\frac { R_0 S\left( t \right)}{N}<1\]
Por lo tanto:
\[ R_0 < \frac {N}{ S\left( t \right)}\]
ya obtuvimos segun (26) que:\[ S(t) = N e^{- \frac {\beta }{\gamma N } (R(t))}\]
entonces:
\[ { R_0 } < e^{ \frac {\beta }{\gamma N } (R(t))} \]
\[ { R_0 } < e^{ \frac {R_0 R(t)} {N}}\] lo que indica que en el momento t la epidemia tenbdera a disminuir y viceversa.
La estimacion de \(R_0\) la obtenemos de (37):
\[ S(t) = N e^{- \frac {R_0 R(t)}{ N } }\] \[ \ln \frac {S(t)}{N} = \ln e^{- \frac {R_0 R(t)}{ N } }\] \[ \ln S(t) = \ln N - \frac {R_0 R(t)}{ N } \] \[ N \ln S(t) = N \ln N - {R_0 R(t)} \]
y tenemos una funcion lineal sobre la cual podemos aplicar una regresion. El numero basico de reproduccion (cambiandole el signo) es la pendiente de la recta de regresion.
El paso final es compara los valores de \(R_0\) con \(e^{\frac {R_0R(t)}{ N }}\), donde t es el momento final.
#Cargamos la base de datos:
covid19 <- read_csv2("C:/Users/usuario/Documents/GitHub/Analisis_covid19_con_R/analisis_13_mayo_SIR/covid19cc.csv")
str(covid19)## tibble [11,484 x 18] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Region : chr [1:11484] "Maule" "Maule" "Metropolitana" "Metropolitana" ...
## $ Comuna : chr [1:11484] "Talca" "Talca" "Vitacura" "Providencia" ...
## $ Codigo_comuna : num [1:11484] 7101 7101 13132 13123 7101 ...
## $ Fec : num [1:11484] 1 2 2 3 3 3 4 4 4 5 ...
## $ Fecha : num [1:11484] 43893 43894 43894 43895 43895 ...
## $ Clave : num [1:11484] 7.10e+08 7.10e+08 1.31e+09 1.31e+09 7.10e+08 ...
## $ Casos_Diarios : num [1:11484] 2 0 1 1 0 0 0 0 1 0 ...
## $ Casos_Acumulados : num [1:11484] 2 2 1 1 2 1 1 2 2 1 ...
## $ acumulados_corregidos: num [1:11484] 0.845 0.845 1.033 0.634 0.845 ...
## $ Muertes : num [1:11484] 0 0 0 0 0 0 0 0 0 0 ...
## $ Casos_Activos : num [1:11484] 2 2 1 1 2 1 1 2 2 1 ...
## $ activos_corregidos : num [1:11484] 0.845 0.845 1.033 0.634 0.845 ...
## $ Recuperados : num [1:11484] 0 0 0 0 0 0 0 0 0 0 ...
## $ Recuperados_Diarios : num [1:11484] 0 0 0 0 0 0 0 0 0 0 ...
## $ Dia : num [1:11484] 1 2 1 1 3 2 2 4 3 3 ...
## $ Fallecidos_Diarios : num [1:11484] 0 0 0 0 0 0 0 0 0 0 ...
## $ on-off : chr [1:11484] "START" "ON" "START" "START" ...
## $ Max_Dia : num [1:11484] 68 68 67 66 68 67 66 68 67 66 ...
## - attr(*, "spec")=
## .. cols(
## .. Region = col_character(),
## .. Comuna = col_character(),
## .. Codigo_comuna = col_double(),
## .. Fec = col_double(),
## .. Fecha = col_double(),
## .. Clave = col_double(),
## .. Casos_Diarios = col_double(),
## .. Casos_Acumulados = col_double(),
## .. acumulados_corregidos = col_double(),
## .. Muertes = col_double(),
## .. Casos_Activos = col_double(),
## .. activos_corregidos = col_double(),
## .. Recuperados = col_double(),
## .. Recuperados_Diarios = col_double(),
## .. Dia = col_double(),
## .. Fallecidos_Diarios = col_double(),
## .. `on-off` = col_character(),
## .. Max_Dia = col_double()
## .. )#Seleccionamos la comuna de Santiago:
covid19.metro = covid19[covid19$Comuna=="Santiago",]
tabla.metro = data.frame(covid19.metro)
head(tabla.metro,10)## Region Comuna Codigo_comuna Fec Fecha Clave Casos_Diarios
## 1 Metropolitana Santiago 13101 12 43904 1310143904 1
## 2 Metropolitana Santiago 13101 13 43905 1310143905 0
## 3 Metropolitana Santiago 13101 14 43906 1310143906 5
## 4 Metropolitana Santiago 13101 15 43907 1310143907 2
## 5 Metropolitana Santiago 13101 16 43908 1310143908 5
## 6 Metropolitana Santiago 13101 17 43909 1310143909 0
## 7 Metropolitana Santiago 13101 18 43910 1310143910 0
## 8 Metropolitana Santiago 13101 19 43911 1310143911 3
## 9 Metropolitana Santiago 13101 20 43912 1310143912 0
## 10 Metropolitana Santiago 13101 21 43913 1310143913 0
## Casos_Acumulados acumulados_corregidos Muertes Casos_Activos
## 1 1 0.1987491 0 1
## 2 1 0.1987491 0 1
## 3 6 1.1924944 0 6
## 4 8 1.5899926 0 8
## 5 13 2.5837380 0 13
## 6 13 2.5837380 0 13
## 7 13 2.5837380 0 13
## 8 16 3.1799852 0 16
## 9 16 3.1799852 0 16
## 10 16 3.1799852 0 16
## activos_corregidos Recuperados Recuperados_Diarios Dia Fallecidos_Diarios
## 1 0.1987491 0 0 1 0
## 2 0.1987491 0 0 2 0
## 3 1.1924944 0 0 3 0
## 4 1.5899926 0 0 4 0
## 5 2.5837380 0 0 5 0
## 6 2.5837380 0 0 6 0
## 7 2.5837380 0 0 7 0
## 8 3.1799852 0 0 8 0
## 9 3.1799852 0 0 9 0
## 10 3.1799852 0 0 10 0
## on.off Max_Dia
## 1 START 57
## 2 ON 57
## 3 ON 57
## 4 ON 57
## 5 ON 57
## 6 ON 57
## 7 ON 57
## 8 ON 57
## 9 ON 57
## 10 ON 57