Preparación de datos

Paquetes

Para poder generar este modelo requeriremos los siguientes paquetes

library(tidyverse)
library(sf)
library(lubridate)

Datos

Para la procesar el modelo requerimos de las siguiente base de datos:

Regiones

Regiones <- read_rds("Regiones2.rds")

Esta base de datos es un sf que tiene la población de los tres grupos de edad (menores de 25, personas entre 25 y 65 y mayores de 65), número de infectados al 17 marzo del 2020, los polígonos de cada región para su visualización. Podemos ver los datos en la siguiente tabla:

Datos por región
Region Adult Over_65 Under_25 Poblacion Infectados Prevalencia
Región de Arica y Parinacota 116055 24577 85436 226068 1 4.40e-06
Región de Tarapacá 176871 24971 128716 330558 0 0.00e+00
Región de Antofagasta 339740 45441 222353 607534 2 3.30e-06
Región de Magallanes y Antártica Chilena 93100 19380 54053 166533 2 1.20e-05
Región de Aysén del Gral.Ibañez del Campo 57415 9290 36453 103158 1 9.70e-06
Región de Atacama 150896 28110 107162 286168 1 3.50e-06
Región de Coquimbo 391067 89543 276976 757586 2 2.60e-06
Región de Valparaíso 947125 247113 621664 1815902 4 2.20e-06
Región Metropolitana de Santiago 3879060 767377 2466371 7112808 247 3.47e-05
Región de Los Lagos 446193 92888 289627 828708 16 1.93e-05
Región de Los Ríos 200546 48428 135863 384837 1 2.60e-06
Región de La Araucanía 493193 120381 343650 957224 7 7.30e-06
Región del Bío-Bío 819734 183145 553926 1556805 14 9.00e-06
Región de Ñuble 255047 65116 160446 480609 28 5.83e-05
Región del Maule 555289 128570 361091 1044950 14 1.34e-05
Región del Libertador Bernardo O’Higgins 493147 108926 312482 914555 2 2.20e-06

El tener esta base de datos nos permite visualizar fácilmente los datos como un mapa, como se ve en la siguiente figura:

Viajes entre regiones:

En el archivo Viajes_Regiones.rds tenemos una tabla (tipo edgelist) de los números promedio de viajes entre regiones. Esta matriz poblacional fue transformada a probabilidades de movimineto para ser trabajada a partir del modelo realizado por Arenas et al. 2020

Viajes_Regiones <- read_rds("Viajes_Regiones.rds")

En la siguiente tabla se ven las 20 primeras lecturas:

Datos de viajes entre regiones, promedio de personas al día
origen destino n_personas
Región de Valparaíso Región Metropolitana de Santiago 27176.095
Región Metropolitana de Santiago Región de Valparaíso 26876.810
Región Metropolitana de Santiago Región del Libertador Bernardo O’Higgins 20039.000
Región del Libertador Bernardo O’Higgins Región Metropolitana de Santiago 19877.190
Región Metropolitana de Santiago Región de Atacama 13219.238
Región de Atacama Región Metropolitana de Santiago 11941.952
Región del Libertador Bernardo O’Higgins Región del Maule 10899.857
Región de Ñuble Región del Bío-Bío 10309.143
Región del Maule Región del Libertador Bernardo O’Higgins 10300.381
Región del Bío-Bío Región de Ñuble 9997.762
Región del Libertador Bernardo O’Higgins Región de Atacama 7452.905
Región de Atacama Región del Libertador Bernardo O’Higgins 7330.762
Región del Maule Región de Ñuble 6617.571
Región de Antofagasta Región de Atacama 6193.429
Región de Atacama Región de Antofagasta 6190.333
Región de Ñuble Región del Maule 6186.143
Región de Los Lagos Región de Los Ríos 5268.048
Región de Los Ríos Región de Los Lagos 5208.000
Región de La Araucanía Región del Bío-Bío 4908.905
Región del Bío-Bío Región de La Araucanía 4763.048

De acuerdo a lo expresado por Arenas et al 2020, estos datos deben ser transformados a una matriz \(R_{i,j}\). Esta matriz de probabilidades (estocástica), donde \(R_{i,j}\) es la probabilidad que si una persona de \(i\) se mueve (lo cual hace con probabilidad \(p\)) tenga como destino \(j\). La suma en \(j\) de \(R_{i,j}\) es igual a 1.

Como ejemplo para calcular esto, lo vemos para la región de Valparaíso como \(i\) (origen), hacia todo \(j\), destino. Para esto tomamos la base de datos Viajes_Regiones y tomamos todos los viajes que salen de Valparaíso, y dividimos los viajes por el total de viajes de salida, esto nos da:

Prob_Valpo <- Viajes_Regiones %>% dplyr::filter(origen ==   "Región de Valparaíso") %>% mutate(p = n_personas/sum(n_personas)) %>% dplyr::select(destino, p)

Lo cual dá el siguiente data.frame

Probabilidades de salidas de Valparaíso
destino p
Región Metropolitana de Santiago 0.713
Región de Coquimbo 0.088
Región de Atacama 0.085
Región de Arica y Parinacota 0.068
Región del Libertador Bernardo O’Higgins 0.046
Región de Antofagasta 0.000
Región del Maule 0.000
Región de Tarapacá 0.000
Región del Bío-Bío 0.000
Región de Ñuble 0.000
Región de Los Lagos 0.000
Región de La Araucanía 0.000
Región de Magallanes y Antártica Chilena 0.000
Región de Los Ríos 0.000
Región de Aysén del Gral.Ibañez del Campo 0.000

De la misma forma, podemos hacerlo para todas las regiones:

#obtenemos todos los nombres de las regiones
Nombres <- Regiones$Region


#Generamos una lista de Probabilidades
Probs <-list()

#Hacemos el proceso que teníamos arriba para todas las regiones


for(i in 1:length(Nombres)){
  Probs[[i]] <- Viajes_Regiones %>% dplyr::filter(origen == Nombres[i]) %>% mutate(p = n_personas/sum(n_personas)) %>% dplyr::select(destino, p)
  colnames(Probs[[i]])[2] <- Nombres[i]
}

#Luego unimos todos estos dataframes con un full_join y remplazamos los na por 0 en caso de = origen destino
Probs <- Probs %>% reduce(full_join) %>% mutate_if(is.numeric, list(~replace_na(., 0)))

Lo cual entrega el siguiente data frame de probabilidades

Probabilidades de movimiento entre regiones
destino Región de Arica y Parinacota Región de Tarapacá Región de Antofagasta Región de Magallanes y Antártica Chilena Región de Aysén del Gral.Ibañez del Campo Región de Atacama Región de Coquimbo Región de Valparaíso Región Metropolitana de Santiago Región de Los Lagos Región de Los Ríos Región de La Araucanía Región del Bío-Bío Región de Ñuble Región del Maule Región del Libertador Bernardo O’Higgins
Región de Valparaíso 0.343 0.003 0.006 0.002 0.000 0.071 0.346 0.000 0.424 0.000 0.000 0.000 0.000 0.000 0.001 0.034
Región de Coquimbo 0.288 0.001 0.014 0.000 0.000 0.083 0.000 0.088 0.007 0.000 0.000 0.000 0.000 0.000 0.000 0.002
Región de Atacama 0.283 0.097 0.765 0.629 0.030 0.000 0.365 0.085 0.209 0.119 0.048 0.131 0.082 0.034 0.052 0.173
Región de Tarapacá 0.074 0.000 0.117 0.000 0.000 0.012 0.001 0.000 0.004 0.000 0.000 0.000 0.000 0.000 0.009 0.072
Región Metropolitana de Santiago 0.008 0.041 0.092 0.162 0.008 0.300 0.048 0.713 0.000 0.030 0.006 0.026 0.026 0.007 0.046 0.461
Región de Antofagasta 0.004 0.162 0.000 0.003 0.000 0.155 0.014 0.000 0.014 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Región del Libertador Bernardo O’Higgins 0.000 0.553 0.001 0.002 0.000 0.184 0.005 0.046 0.316 0.000 0.001 0.002 0.003 0.003 0.534 0.000
Región del Bío-Bío 0.000 0.000 0.001 0.003 0.000 0.030 0.000 0.000 0.006 0.001 0.005 0.465 0.000 0.588 0.009 0.001
Región del Maule 0.000 0.031 0.000 0.001 0.000 0.021 0.000 0.000 0.005 0.000 0.001 0.004 0.032 0.353 0.000 0.253
Región de Magallanes y Antártica Chilena 0.000 0.000 0.000 0.000 0.001 0.007 0.000 0.000 0.001 0.009 0.000 0.000 0.000 0.000 0.000 0.000
Región de La Araucanía 0.000 0.000 0.000 0.001 0.000 0.030 0.000 0.000 0.004 0.005 0.395 0.000 0.275 0.012 0.004 0.001
Región de Los Ríos 0.000 0.000 0.000 0.001 0.000 0.010 0.000 0.000 0.001 0.461 0.000 0.356 0.002 0.002 0.001 0.000
Región de Ñuble 0.000 0.000 0.000 0.001 0.000 0.012 0.000 0.000 0.001 0.000 0.001 0.006 0.578 0.000 0.343 0.001
Región de Arica y Parinacota 0.000 0.112 0.005 0.000 0.000 0.052 0.220 0.068 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Región de Los Lagos 0.000 0.000 0.000 0.188 0.960 0.029 0.000 0.000 0.006 0.000 0.543 0.009 0.001 0.001 0.001 0.000
Región de Aysén del Gral.Ibañez del Campo 0.000 0.000 0.000 0.008 0.000 0.003 0.000 0.000 0.001 0.373 0.000 0.000 0.000 0.000 0.000 0.000

Generación de condiciones iniciales

Para generar un data frame de las condiciones iniciales, hacemos lo siguente

df_out <- Regiones %>% as.data.frame() %>% 
  #Transformamos el sf en un data.frame lo hace más rápido
           select(Region, Poblacion, Adult, Over_65, Under_25, Infectados) %>%  
  #Seleccionamos las columnas necesarias
          mutate(Infectados_Adult = round((Adult/Poblacion)*Infectados), Infectados_Over_65 = round((Over_65/Poblacion)*Infectados),Infectados_Under_25 = round((Under_25/Poblacion)*Infectados), Infectados_Adult = Infectados_Adult + (Infectados - (Infectados_Adult + Infectados_Under_25 + Infectados_Over_65))) %>%  
  #Calculamos los infectados por grupo etareo
    mutate(Suceptibles_Adult = Adult - Infectados_Adult, Suceptibles_Under_25 = Under_25 - Infectados_Under_25, Suceptibles_Over_65 = Over_65 - Infectados_Over_65, Expuestos_Adult = 0, Expuestos_Under_25 = 0, Expuestos_Over_65 = 0, Asintomaticos_Adult = 0, Asintomaticos_Under_25 = 0, Asintomaticos_Over_65 = 0, UCI_Adult = 0, UCI_Under_25 = 0, UCI_Over_65 = 0, Muerto_Adult = 0, Muerto_Under_25 = 0, Muerto_Over_65 = 0, Recuperados_Adult = 0, Recuperados_Under_25 = 0, Recuperados_Over_65 = 0) %>%  
  #Agregamos el resto de las etapas
  mutate(Time = 1) # Todo esto es el tiempo 1

Esto hace que las condiciones inciales sean las siguientes:

Condiciones inciales por región
Region Poblacion Adult Over_65 Under_25 Infectados Infectados_Adult Infectados_Over_65 Infectados_Under_25 Suceptibles_Adult Suceptibles_Under_25 Suceptibles_Over_65 Expuestos_Adult Expuestos_Under_25 Expuestos_Over_65 Asintomaticos_Adult Asintomaticos_Under_25 Asintomaticos_Over_65 UCI_Adult UCI_Under_25 UCI_Over_65 Muerto_Adult Muerto_Under_25 Muerto_Over_65 Recuperados_Adult Recuperados_Under_25 Recuperados_Over_65 Time
Región de Arica y Parinacota 226068 116055 24577 85436 1 1 0 0 116054 85436 24577 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de Tarapacá 330558 176871 24971 128716 0 0 0 0 176871 128716 24971 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de Antofagasta 607534 339740 45441 222353 2 1 0 1 339739 222352 45441 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de Magallanes y Antártica Chilena 166533 93100 19380 54053 2 1 0 1 93099 54052 19380 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de Aysén del Gral.Ibañez del Campo 103158 57415 9290 36453 1 1 0 0 57414 36453 9290 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de Atacama 286168 150896 28110 107162 1 1 0 0 150895 107162 28110 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de Coquimbo 757586 391067 89543 276976 2 1 0 1 391066 276975 89543 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de Valparaíso 1815902 947125 247113 621664 4 2 1 1 947123 621663 247112 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región Metropolitana de Santiago 7112808 3879060 767377 2466371 247 134 27 86 3878926 2466285 767350 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de Los Lagos 828708 446193 92888 289627 16 8 2 6 446185 289621 92886 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de Los Ríos 384837 200546 48428 135863 1 1 0 0 200545 135863 48428 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de La Araucanía 957224 493193 120381 343650 7 3 1 3 493190 343647 120380 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región del Bío-Bío 1556805 819734 183145 553926 14 7 2 5 819727 553921 183143 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región de Ñuble 480609 255047 65116 160446 28 15 4 9 255032 160437 65112 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región del Maule 1044950 555289 128570 361091 14 7 2 5 555282 361086 128568 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Región del Libertador Bernardo O’Higgins 914555 493147 108926 312482 2 1 0 1 493146 312481 108926 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

Para facilidad programática, dividiremos df_out en una lista con 3 data frames, cada uno con las condiciones inciales para los menores de 25, las personas entre 25 y 65 y las presonas mayores de 65 (estas son las llamadas generaciones y se denotará con el supraíndice g o h dependiendo de la ecuación).

Under_25 <- df_out %>% dplyr::select(Region, Poblacion, contains("Under_25"), Time) %>% rename(Generacion = Under_25)

Adult <- df_out %>% dplyr::select(Region, Poblacion, contains("Adult"), Time) %>% rename(Generacion = Adult)

Over_65 <- df_out %>% dplyr::select(Region, Poblacion, contains("Over_65"), Time) %>% rename(Generacion = Over_65)

#Cambiamos los nombres de las columnas para que todos tengas cosas iguales

colnames(Under_25) <- str_remove_all(string = colnames(Under_25), pattern = "_Under_25")

colnames(Adult) <- str_remove_all(string = colnames(Adult), pattern = "_Adult")

colnames(Over_65) <- str_remove_all(string = colnames(Over_65), pattern = "_Over_65")

### Lo agregamos todo a una lista llamada df_out con los nombres de las edades


df_out <- list(Under_25, Adult, Over_65)

names(df_out) <- c("Under_25", "Adult", "Over_65")

Aquí vemos el grupo entre 25 y 65 como ejemplo

Condiciones inciales por región para las personas entre 25 y 65
Region Poblacion Generacion Infectados Suceptibles Expuestos Asintomaticos UCI Muerto Recuperados Time
Región de Arica y Parinacota 226068 116055 1 116054 0 0 0 0 0 1
Región de Tarapacá 330558 176871 0 176871 0 0 0 0 0 1
Región de Antofagasta 607534 339740 1 339739 0 0 0 0 0 1
Región de Magallanes y Antártica Chilena 166533 93100 1 93099 0 0 0 0 0 1
Región de Aysén del Gral.Ibañez del Campo 103158 57415 1 57414 0 0 0 0 0 1
Región de Atacama 286168 150896 1 150895 0 0 0 0 0 1
Región de Coquimbo 757586 391067 1 391066 0 0 0 0 0 1
Región de Valparaíso 1815902 947125 2 947123 0 0 0 0 0 1
Región Metropolitana de Santiago 7112808 3879060 134 3878926 0 0 0 0 0 1
Región de Los Lagos 828708 446193 8 446185 0 0 0 0 0 1
Región de Los Ríos 384837 200546 1 200545 0 0 0 0 0 1
Región de La Araucanía 957224 493193 3 493190 0 0 0 0 0 1
Región del Bío-Bío 1556805 819734 7 819727 0 0 0 0 0 1
Región de Ñuble 480609 255047 15 255032 0 0 0 0 0 1
Región del Maule 1044950 555289 7 555282 0 0 0 0 0 1
Región del Libertador Bernardo O’Higgins 914555 493147 1 493146 0 0 0 0 0 1

Parámetros y modelo

Usamos los parámetros del modelo de Arenas

#Infectividad de asintomaticos
betaA = 0.06

#Infectividad de Infectados
betaI = 0.06

#numero de contactos promedio
K_g = c(11.8, 13.3, 6.6)

#Tasa latente
Eta = 1/2.34

#Tasa de infeccion asintomática
Alpha_g = c(1/5.06, 1/2.86, 1/2.86)

#Tasa de escape
Mu_g = c(1, 1/3.2, 1/3.2)

#Fracción de casos que requieren UCI
Gamma_g = c(0.002,0.05,0.36)

#Fatalidad en la UCI
Omega_g = 0.42

#Tasa de mortalidad

Psi_g = 1/7

# Tasa de salida de la UCI

Chi_g = 0.1

# Factor de densidad

Epsilon = 0.01

#Factor de movilidad

p_G = c(0, 1, 0)

###Tamaño promedio de hogar

Sigma = 2.5

## Tasa de confinamiento

K_0 = 0


## Matriz de contactos entre generaciones


C_G_H <- matrix(data = c(0.5980, 0.3849, 0.0171, 0.2440, 0.7210, 0.0350, 0.1919, 0.5705, 0.2376), nrow = 3, ncol = 3, byrow = T)

Cálculo de proporción del grupo Susceptibles

Para desarrollar esta ecuación partimos desde la ecuación 14 de Arenas et al. donde calculamos \(n_{j,i}^{m,j}(t)\) que es el número esperado de individuos, del grupo etareo \(h\) en estado infeccioso \(m\) (m pudiendo ser del grupo Infectados o Asintomáticos) que se mueven de la región \(j\) a la \(i\),

Primero lo calcularemos para \(m = I\) para cada región y cada generación

\[n_{j,i}^{I,h}(t) = n_j^h\times\rho_j^{I,h}(t)\times[(1 - p^h)\delta_{i,j} + p^h\times R_{j,i}^h]\]

Ejemplo movimiento adultos de la región Metropolitana a la región de Valparaíso

Esto se calcula con el siguiente código, donde \(i\) es el destino (Valparaíso) y \(j\) el origen (Metropolitana) y \(h\) es el grupo etáreo, en este caso los adultos.

Mat  <- matrix(rep(0,(length(Nombres)*length(Nombres))), nrow = length(Nombres), ncol = length(Nombres))

colnames(Mat) <- Nombres
rownames(Mat) <- Nombres

N_I_h_j_i <- list(Mat, Mat, Mat)
names(N_I_h_j_i) <- c("Under_25", "Adult", "Over_65")

N_A_h_j_i <- list(Mat, Mat, Mat)
names(N_A_h_j_i) <- c("Under_25", "Adult", "Over_65")


for(x in 1:length(N_I_h_j_i)){
  for(i in 1:length(Nombres)){
    for(j in 1:length(Nombres)){
      N_I_h_j_i[[x]][j,i] <- (df_out[[x]] %>% dplyr::filter(Region == Nombres[j]) %>% pull(Generacion))*((df_out[[x]] %>% dplyr::filter(Region == Nombres[j]) %>% pull(Infectados))/(df_out[[x]] %>% dplyr::filter(Region == Nombres[j]) %>% pull(Generacion)))*((1-p_G[x])*0 + p_G[x]*Probs %>% dplyr::filter(destino == Nombres[i]) %>% pull(Nombres[j]))
    }
  }
}

Como ejemplo mostraremos el \(n_{j,i}^{I,h}(t)\), para los adultos esta matríz es

para adultos
Región de Arica y Parinacota Región de Tarapacá Región de Antofagasta Región de Magallanes y Antártica Chilena Región de Aysén del Gral.Ibañez del Campo Región de Atacama Región de Coquimbo Región de Valparaíso Región Metropolitana de Santiago Región de Los Lagos Región de Los Ríos Región de La Araucanía Región del Bío-Bío Región de Ñuble Región del Maule Región del Libertador Bernardo O’Higgins
Región de Arica y Parinacota 0.000 0.074 0.004 0.000 0.000 0.283 0.288 0.343 0.008 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Región de Tarapacá 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Región de Antofagasta 0.005 0.117 0.000 0.000 0.000 0.765 0.014 0.006 0.092 0.000 0.000 0.000 0.001 0.000 0.000 0.001
Región de Magallanes y Antártica Chilena 0.000 0.000 0.003 0.000 0.008 0.629 0.000 0.002 0.162 0.188 0.001 0.001 0.003 0.001 0.001 0.002
Región de Aysén del Gral.Ibañez del Campo 0.000 0.000 0.000 0.001 0.000 0.030 0.000 0.000 0.008 0.960 0.000 0.000 0.000 0.000 0.000 0.000
Región de Atacama 0.052 0.012 0.155 0.007 0.003 0.000 0.083 0.071 0.300 0.029 0.010 0.030 0.030 0.012 0.021 0.184
Región de Coquimbo 0.220 0.001 0.014 0.000 0.000 0.365 0.000 0.346 0.048 0.000 0.000 0.000 0.000 0.000 0.000 0.005
Región de Valparaíso 0.136 0.000 0.001 0.000 0.000 0.170 0.175 0.000 1.425 0.000 0.000 0.000 0.000 0.000 0.001 0.092
Región Metropolitana de Santiago 0.156 0.561 1.920 0.193 0.082 27.968 0.888 56.863 0.000 0.817 0.082 0.475 0.749 0.131 0.717 42.397
Región de Los Lagos 0.000 0.000 0.001 0.075 2.982 0.954 0.000 0.002 0.237 0.000 3.691 0.039 0.010 0.003 0.004 0.002
Región de Los Ríos 0.000 0.000 0.000 0.000 0.000 0.048 0.000 0.000 0.006 0.543 0.000 0.395 0.005 0.001 0.001 0.001
Región de La Araucanía 0.000 0.000 0.000 0.000 0.000 0.392 0.000 0.001 0.078 0.028 1.069 0.000 1.394 0.017 0.013 0.007
Región del Bío-Bío 0.000 0.000 0.002 0.001 0.000 0.572 0.001 0.003 0.180 0.010 0.012 1.928 0.000 4.048 0.224 0.019
Región de Ñuble 0.000 0.000 0.001 0.000 0.000 0.513 0.001 0.005 0.100 0.014 0.024 0.180 8.824 0.000 5.295 0.044
Región del Maule 0.000 0.064 0.000 0.000 0.000 0.361 0.001 0.008 0.324 0.006 0.008 0.029 0.061 2.401 0.000 3.737
Región del Libertador Bernardo O’Higgins 0.000 0.072 0.000 0.000 0.000 0.173 0.002 0.034 0.461 0.000 0.000 0.001 0.001 0.001 0.253 0.000 
for(x in 1:length(N_A_h_j_i)){
  for(i in 1:length(Nombres)){
    for(j in 1:length(Nombres)){
      N_A_h_j_i[[x]][j,i] <- (df_out[[x]] %>% dplyr::filter(Region == Nombres[j]) %>% pull(Generacion))*((df_out[[x]] %>% dplyr::filter(Region == Nombres[j]) %>% pull(Asintomaticos))/(df_out[[x]] %>% dplyr::filter(Region == Nombres[j]) %>% pull(Generacion)))*((1-p_G[x])*0 + p_G[x]*Probs %>% dplyr::filter(destino == Nombres[i]) %>% pull(Nombres[j]))
    }
  }
}

Seguimos en el mismo ejemplo con la ecuación 12

\[(n_i^g)^{eff} = \sum_j[(1 - p^g)\delta_{i,j} + p^g\times R_{j,i}^g]\times n_j^g\]

Esto es, para cada region y generación:

#Para cada generacion
for(x in 1:3){
  #creamos la columna a llenar 
  df_out[[x]]$n_i_g_eff <- NA
  #Para cada region
  for(R in 1:nrow(df_out[[x]])){
   n_i_g_eff <- vector()
    for(i in 1:length(Nombres)){
      n_i_g_eff[i] <- ((1 - p_G[x])*ifelse(Nombres[i] == df_out[[x]][R,]$Region, 1, 0) + p_G[x]*Probs %>% dplyr::filter(destino == df_out[[x]][R,]$Region) %>% pull(Nombres[i]))*(df_out[[x]] %>% dplyr::filter(Region == Nombres[i]) %>% pull(Generacion))
    }
    df_out[[x]][R,]$n_i_g_eff <- sum(n_i_g_eff) 
  }
}

Esto da un valor para \((n_i^g)^{eff}\), donde \(i\) es la región de Valparaíso y \(g\) es la población de entre 25 y 65 años de (la cual es la población de los adultos que potencialmente se podrían mover a valparaíso desde todas las otras regiones regiones).

Condiciones inciales por región para las personas entre 25 y 65, agregando n_i_g_eff
Region Poblacion Generacion Infectados Suceptibles Expuestos Asintomaticos UCI Muerto Recuperados Time n_i_g_eff
Región de Arica y Parinacota 226068 116055 1 116054 0 0 0 0 0 1 183742.34
Región de Tarapacá 330558 176871 0 176871 0 0 0 0 0 1 107419.96
Región de Antofagasta 607534 339740 1 339739 0 0 0 0 0 1 114852.87
Región de Magallanes y Antártica Chilena 166533 93100 1 93099 0 0 0 0 0 1 11129.46
Región de Aysén del Gral.Ibañez del Campo 103158 57415 1 57414 0 0 0 0 0 1 170022.35
Región de Atacama 286168 150896 1 150895 0 0 0 0 0 1 1720074.20
Región de Coquimbo 757586 391067 1 391066 0 0 0 0 0 1 161037.71
Región de Valparaíso 1815902 947125 2 947123 0 0 0 0 0 1 1852565.55
Región Metropolitana de Santiago 7112808 3879060 134 3878926 0 0 0 0 0 1 1097512.11
Región de Los Lagos 828708 446193 8 446185 0 0 0 0 0 1 216349.94
Región de Los Ríos 384837 200546 1 200545 0 0 0 0 0 1 388007.95
Región de La Araucanía 957224 493193 3 493190 0 0 0 0 0 1 331529.37
Región del Bío-Bío 1556805 819734 7 819727 0 0 0 0 0 1 413186.05
Región de Ñuble 480609 255047 15 255032 0 0 0 0 0 1 674020.49
Región del Maule 1044950 555289 7 555282 0 0 0 0 0 1 273393.78
Región del Libertador Bernardo O’Higgins 914555 493147 1 493146 0 0 0 0 0 1 1699633.88

Seguimos con la ecuación 11:

\[n_i^{eff} = \sum_{g=1}^{N_G}(n_i^g)^{eff}\]

Esto se calcula así:

n_i_eff <- data.frame(Region = df_out$Adult$Region, n_i_eff = df_out[[1]]$n_i_g_eff + df_out[[2]]$n_i_g_eff+ df_out[[3]]$n_i_g_eff)

Esto hace que la población efectiva de las 3 generaciones desde distintas generaciones (\(n_i^{eff}\)), lo cual guardamos en un data.frame, el que se ve así

n_i_eff para cada region
Region n_i_eff
Región de Arica y Parinacota 293755.34
Región de Tarapacá 261106.96
Región de Antofagasta 382646.87
Región de Magallanes y Antártica Chilena 84562.46
Región de Aysén del Gral.Ibañez del Campo 215765.35
Región de Atacama 1855346.20
Región de Coquimbo 527556.71
Región de Valparaíso 2721342.55
Región Metropolitana de Santiago 4331260.11
Región de Los Lagos 598864.94
Región de Los Ríos 572298.94
Región de La Araucanía 795560.37
Región del Bío-Bío 1150257.05
Región de Ñuble 899582.49
Región del Maule 763054.78
Región del Libertador Bernardo O’Higgins 2121041.88

Para poder ocupar la función de la ecuación 13, generaremos una función en R

\[f(x) = 1 + (1 - e^{-\epsilon x})\]

Func <- function(x, epsilon = 0.01){
  y <- 1 + (1- exp(-epsilon*x))
  return(y)
}

Ahora seguimos con la ecución 10

\[z^g = \frac{N^g}{\sum_{i = 1}^N f(\frac{n_i^{eff}}{S_i})(n_i^g)^{eff}}\]

Entonces \(z^g\) para la generación de entre 25 y 65 es

z_g <- vector()

for(x in 1:length(df_out)){
  z_g[x] <- (df_out[[x]] %>% summarise(N = sum(Generacion)) %>% pull(N))/sum(Func(n_i_eff$n_i_eff/(df_out[[1]]$Suceptibles + df_out[[2]]$Suceptibles + df_out[[3]]$Suceptibles))*df_out[[x]]$n_i_g_eff)
}

El valor de \(z^g\) es un valor único por cada generación \(g\), esto se almacenó en el vector z_g, donde sus valores son 0.99, 0.977, 0.99 para cada generación respectivamente.

Con eso, podemos ya calcular la ecuación 9, la cual es:

\[P_i^g(t) = 1 - \prod_{h = 1}^{N_G}\prod_{j = 1}^{N}(1 - \beta_A)^{z^g(k^g)f(\frac{n_i^{eff}}{S_i})*C^{g,h}\frac{n_{j,i}^{A,h}(t)}{(n_i^h)^eff}}\times(1 - \beta_I)^{z^g(k^g)f(\frac{n_i^{eff}}{S_i})*C^{g,h}\frac{n_{j,i}^{I,h}(t)}{(n_i^h)^eff}}\]

for (g in 1:3) {
    df_out[[g]]$P_G <- NA
    Temp2 <- list()
    for (h in 1:3) {
        temp <- list()
        for (j in 1:16) {
            temp[[j]] <- (1 - betaA)^(z_g[g] * K_g[g] * Func(x = n_i_eff$n_i_eff/(df_out[[1]]$Suceptibles + 
                df_out[[2]]$Suceptibles + df_out[[3]]$Suceptibles)) * C_G_H[g, 
                h] * (N_A_h_j_i[[h]][j, ]/df_out[[h]]$n_i_g_eff)) * (1 - betaI)^(z_g[g] * 
                K_g[g] * Func(x = n_i_eff$n_i_eff/(df_out[[1]]$Suceptibles + 
                df_out[[2]]$Suceptibles + df_out[[3]]$Suceptibles)) * C_G_H[g, 
                h] * (N_I_h_j_i[[h]][j, ]/df_out[[h]]$n_i_g_eff))
        }
        Temp2[[h]] <- reduce(temp, `*`)
    }
    
    df_out[[g]]$P_G <- 1 - reduce(Temp2, `*`)
}

\(P_i^g(t)\) Es la probabilidad de que un agente de la generación \(G\) se infecte del patógeno dentro de la región \(i\), el resultado para cada generación queda almacenado en la columna P_G de cada dataframe de la lista df_out

Esto nos lleva a la ecuación 8 la cual es

\[\Pi_i^g(t) = (1-p_g)\times P_i^g(t) + p_g \times \sum_{j = 1}^N R_{i,j}^g \times P_j^g(t)\]

Esto lo calculamos así

for(g in 1:3){
  df_out[[g]]$PI = NA
  for(i in 1:nrow(Regiones)){
    df_out[[g]]$PI[i] <- p_G[g]*(df_out[[g]] %>% dplyr::filter(Region == Nombres[i]) %>% pull(P_G)) + p_G[g]*sum((Probs %>% dplyr::filter(destino == Nombres[i]) %>% select_if(is.numeric) %>% reduce(c))*df_out[[g]]$P_G)
  }
}

Calculo proporciones t + 1

Ahora podemos al fin ir a las primeras ecuaciones donde calculamos las proporciones de los principales grupos según el modelo: susceptibles, expuestos, asintomáticos, infectados, hospitalizados en UCI, muertos y recuperados. Tal como lo hicimos anteriormente, acá se presentan las proporciones y números totales por grupo, pero solo para la generación de adultos (personas entre 25 y 60 años)

Susceptibles

\[\rho_i^{S,g}(t+1) = \rho_i^{S,g}(t)*(1-\Pi_i^g(t))\]

Para \(t+1\) entonces:

p_susceptibles <- (df_out[[2]]$Suceptibles/df_out[[2]]$Generacion)*(1-df_out[[2]]$PI)
n_susceptibles <- p_susceptibles*df_out[[2]]$Generacion

Expuestos

\[\rho_i^{E,g}(t+1) = \rho_i^{S,g}(t)*\Pi_i^g(t) + (1-\eta^g)*\rho_i^{E,g}(t)\]

p_expuestos <- ((df_out[[2]]$Suceptibles/df_out[[2]]$Generacion)*df_out[[2]]$PI) + (1 - Eta)*(df_out[[2]]$Expuestos/df_out[[2]]$Generacion)
n_expuestos <- p_expuestos*df_out[[2]]$Generacion

Asintomáticos

\[\rho_i^{A,g}(t+1) = \eta^g\rho_i^{E,g}(t)+(1-\alpha^g)\times \rho_i^{A,g}(t)\]

p_asintomaticos <- (Eta*(df_out[[2]]$Expuestos/df_out[[2]]$Generacion)) + (1 - Alpha_g[2])*(df_out[[2]]$Asintomaticos/df_out[[2]]$Generacion)

n_asintomaticos <- p_asintomaticos*df_out[[2]]$Generacion

Infectados

\[\rho_i^{I,g}(t+1) = \alpha^g\rho_i^{A,g}(t)+(1-\mu^g)\times \rho_i^{I,g}(t)\]

p_infectados <- (Alpha_g[2]*(df_out[[2]]$Asintomaticos/df_out[[2]]$Generacion)) + (1 - Mu_g[2])*(df_out[[2]]$Infectados/df_out[[2]]$Generacion)

n_infectados <- p_infectados*df_out[[2]]$Generacion

Hospitalizados

\[\rho_i^{H,g}(t+1) = \mu^g\times \gamma^g \times \rho_i^{I,g}(t) + \omega^g(1 - \psi^g)\times \rho_i^{H,g}(t) + (1- \omega^g)(1-\chi^g)\rho_i^{H,g}(t)\]

p_hospital <- Mu_g[2]*Gamma_g[2]*(df_out[[2]]$Infectados/df_out[[2]]$Generacion) + Omega_g*(1 - Psi_g)*(df_out[[2]]$UCI/df_out[[2]]$Generacion) + (1 - Omega_g)*(1 - Chi_g)*(df_out[[2]]$UCI/df_out[[2]]$Generacion)

n_hospital <- p_hospital*df_out[[2]]$Generacion

Muertos

\[\rho_i^{D,g}(t+1) = \omega^g \times \psi^g \times \rho_i^{H,g}(t) + \rho_i^{D,g}(t)\]

p_muertos <- Omega_g*Psi_g*(df_out[[2]]$UCI/df_out[[2]]$Generacion) + (df_out[[2]]$Muerto/df_out[[2]]$Generacion)
n_muertos <- p_muertos*df_out[[2]]$Generacion

Recuperados

\[\rho_i^{R,g}(t+1) = \mu^g (1 -\gamma^g)\times \rho_i^{I,g}(t) + (1 - \omega^g) \chi^g \times \rho_i^{H,g}(t) + \rho_i^{R,g}(t)\]

p_recuperados <- Mu_g[2]*(1 - Gamma_g[2])*(df_out[[2]]$Infectados/df_out[[2]]$Generacion) + (1 - Omega_g)*Chi_g*(df_out[[2]]$UCI/df_out[[2]]$Generacion) + (df_out[[2]]$Recuperados/df_out[[2]]$Generacion)

n_recuperados <- p_recuperados*df_out[[2]]$Generacion