R0 en el tiempo en rangos y con dos valores.
DataIntelligence
05-09-2020
Aplicaremos el calculo del R0 no a los rangos de series de tiempo sino a dos momentos, a dos dias. Demostraremos que el resultado es idéntico al de si hiciéramos el cálculo sobre rangos de tiempo.
library(knitr)
library(kableExtra)
library(janitor)
library(dplyr)
library(readxl)
library(xlsx)
# la siguiente linea impide el despliegue de advertencias en forma local:
options(warn = -1)
library(tidyverse)
library(gghighlight)
library(magrittr)
library(readr)
library(ape)
library(ggdendro)
library(rmarkdown)
library(ggplot2)
library(gganimate)
library(lubridate)
library(plotly)
theme_set(theme_bw()) covid19 <-
read_xlsx("C:/Users/usuario/Desktop/GitHub-trabajo/Algoritmos_R_y_Python/R0/R0_R/salida.xlsx")
head(covid19)## # A tibble: 6 x 26
## Region `Codigo region` Comuna `Codigo comuna` Poblacion Fecha Acumulado Nuevos
## <chr> <dbl> <chr> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 Arica~ 15 Arica 15101 247552 2020~ 6 6
## 2 Arica~ 15 Arica 15101 247552 2020~ 6 0
## 3 Arica~ 15 Arica 15101 247552 2020~ 12 6
## 4 Arica~ 15 Arica 15101 247552 2020~ 41 29
## 5 Arica~ 15 Arica 15101 247552 2020~ 63 22
## 6 Arica~ 15 Arica 15101 247552 2020~ 87 24
## # ... with 18 more variables: Fallecido <dbl>, `Nuevos Fallecidos` <dbl>,
## # Activo <dbl>, `Nuevos Activo` <dbl>, Recuperados <dbl>, `Nuevos
## # Recuperados` <dbl>, `% Fallecido` <dbl>, `% Activo` <dbl>, `%
## # Recuperados` <dbl>, `% Nuevo Casos` <dbl>, `% Nuevos Fallecido` <dbl>, `%
## # Nuevos Activo` <dbl>, `% Nuevos Recuperados` <dbl>, `Casos/1MM hab` <dbl>,
## # `Fallecidos/1MM hab` <dbl>, `Activos/1MM hab` <dbl>, `Recuperados/1MM
## # hab` <dbl>, r0 <dbl>covid19 <- covid19 %>% filter(Comuna == 'Punta Arenas')
covid19 <- subset(covid19, Comuna != "Chillán")
covid19 <- subset(covid19, Fecha> "2020-05-22")
#covid19 <- subset(covid19, Fecha > "43934")
# covid19 <- covid19[5,]
# covid19
R0_v <- c()
comuna_v <- c()
codigo_comuna_v <- c()
a2 <- 0
s2 <- 0
s_v2 <- c()
y_bueno2 <- c()
R0_v2 <- c()
comuna_v2 <- c()
codigo_comuna_v2 <- c()
a <- 0
s <- 0
s_v <- c()
y_bueno <- c()
Fecha <- as.Date(covid19$Fecha, origin = "1899-12-30")ux <- unique(covid19$"Codigo comuna")
for(t in ux)
{
# susceptibles.por.dia2 <- 0
if (t != "No Informado")
{
comuna_en_cuestion = covid19[covid19$"Codigo comuna" == t, ]
longitud = nrow(comuna_en_cuestion)
comuna <- comuna_en_cuestion$Comuna
codigo_comuna <- comuna_en_cuestion$"Codigo comuna"
# dia <- comuna_en_cuestion$"Día"
fecha <- comuna_en_cuestion$Fecha
if (nrow(comuna_en_cuestion) != 0 && (comuna_en_cuestion$"Poblacion" != 0))
{
population = as.numeric(comuna_en_cuestion$"Poblacion")
infectados.por.dia <<- cumsum(as.numeric(comuna_en_cuestion$"Nuevos"))
fallecidos.por.dia <<- cumsum(as.numeric(comuna_en_cuestion$"Nuevos Fallecidos"))
infectados.por.dia2 = as.numeric(infectados.por.dia + fallecidos.por.dia)
recuperados.por.dia = cumsum(as.numeric(comuna_en_cuestion$"Nuevos Recuperados"))
# susceptibles.por.dia debe ir acumulando
susceptibles.por.dia = population - infectados.por.dia2 - recuperados.por.dia
# Creamos la tabla SIR:
tabla.comuna =
data.frame(
unique(comuna_en_cuestion$Fecha),
susceptibles.por.dia,
infectados.por.dia2,
# infec.acum,
recuperados.por.dia,
population
)
names(tabla.comuna) =
c("Fecha",
"Susceptibles",
"Infectados",
# "Infectados acumulados",
"Recuperados",
"Poblacion")
x = tabla.comuna$Recuperados
# se generan singularidades
y = population * log(tabla.comuna$Susceptibles)
}
}
print(tabla.comuna)
# # susceptibles.por.dia debe ir acumulando la resta de infectados
# for (t in 1:nrow(comuna_en_cuestion))
# {
# s <<- s + tabla.comuna$Infectados[t] + tabla.comuna$Recuperados[t]
# s_v[t] <- tabla.comuna$Poblacion - s
# y_bueno[t] = population * log(s_v[t])
# }
#
# dataf = data.frame(s_v, tabla.comuna$Infectados, y_bueno, tabla.comuna$Recuperados)
# print(dataf)
m <- seq(2, longitud)
for (i in m)
{
xx <- tabla.comuna$Recuperados[1:i]
# xx <- tabla.comuna$Recuperados[i-1:i]
yy <- y[1:i]
# yy <- y[i-1:i]
if(!is.null(yy)){
estimacion.R0 = -summary(lm(yy ~ xx))$coefficients[2]
R0_v[i] <- estimacion.R0
comuna_v <- comuna
codigo_comuna_v <- codigo_comuna
fecha_v <- fecha
# dia_v <- dia
}
}
for (i in m)
{
#xx <- tabla.comuna$Recuperados[1:i]
xxx <- tabla.comuna$Recuperados[i-1:i]
#yy <- y[1:i]
yyy <- y[i-1:i]
if(!is.null(yyy)){
estimacion.R02 = -summary(lm(yyy ~ xxx))$coefficients[2]
R0_v2[i] <- estimacion.R02
comuna_v2 <- comuna
codigo_comuna_v2 <- codigo_comuna
fecha_v2 <- fecha
# dia_v <- dia
}
}
eee <- data.frame(comuna_v, R0_v , fecha_v)
write.table(
eee,
'R0.csv',
col.names = F,
append = T,
sep = ","
)
eee2 <- data.frame(comuna_v2, R0_v2 , fecha_v)
write.table(
eee2,
'R02.csv',
col.names = F,
append = T,
sep = ","
)
R0_v <- c()
comuna_v <- c()
codigo_comuna_v <- c()
a <- 0
s <- 0
s_v <- c()
y <- c()
R0_v2 <- c()
comuna_v2 <- c()
codigo_comuna_v2 <- c()
a2 <- 0
s2 <- 0
s_v2 <- c()
y2 <- c()
}## Fecha Susceptibles Infectados Recuperados Poblacion
## 1 2020-05-25 141939 23 22 141984
## 2 2020-05-29 141879 54 51 141984
## 3 2020-06-01 141844 64 76 141984
## 4 2020-06-05 141800 79 105 141984
## 5 2020-06-08 141769 94 121 141984
## 6 2020-06-12 141736 129 119 141984
## 7 2020-06-15 141647 207 130 141984
## 8 2020-06-19 141549 274 161 141984
## 9 2020-06-23 141464 328 192 141984
## 10 2020-06-28 141315 381 288 141984
## 11 2020-07-01 141267 404 313 141984
## 12 2020-07-05 141182 432 370 141984
## 13 2020-07-10 141088 479 417 141984
## 14 2020-07-13 141037 512 435 141984
## 15 2020-07-17 140987 546 451 141984
## 16 2020-07-20 140937 568 479 141984
## 17 2020-07-24 140870 589 525 141984
## 18 2020-07-27 140820 613 551 141984
## 19 2020-07-31 140671 716 597 141984
## 20 2020-08-03 140583 784 617 141984
## 21 2020-08-07 140404 906 674 141984
## 22 2020-08-10 140209 1037 738 141984
## 23 2020-08-14 139999 1159 826 141984
## 24 2020-08-17 139785 1287 912 141984
## 25 2020-08-21 139400 1530 1054 141984
## 26 2020-08-24 138981 1838 1165 141984
## 27 2020-08-28 138424 2188 1372 141984
## 28 2020-08-31 137952 2507 1525 141984head(eee, 50)## comuna_v R0_v fecha_v
## 1 Punta Arenas NA 2020-05-25
## 2 Punta Arenas 2.070059 2020-05-29
## 3 Punta Arenas 1.768603 2020-06-01
## 4 Punta Arenas 1.653190 2020-06-05
## 5 Punta Arenas 1.664298 2020-06-08
## 6 Punta Arenas 1.834052 2020-06-12
## 7 Punta Arenas 2.218345 2020-06-15
## 8 Punta Arenas 2.627003 2020-06-19
## 9 Punta Arenas 2.853279 2020-06-23
## 10 Punta Arenas 2.571856 2020-06-28
## 11 Punta Arenas 2.467455 2020-07-01
## 12 Punta Arenas 2.333386 2020-07-05
## 13 Punta Arenas 2.253917 2020-07-10
## 14 Punta Arenas 2.225834 2020-07-13
## 15 Punta Arenas 2.221018 2020-07-17
## 16 Punta Arenas 2.210058 2020-07-20
## 17 Punta Arenas 2.181527 2020-07-24
## 18 Punta Arenas 2.158974 2020-07-27
## 19 Punta Arenas 2.167281 2020-07-31
## 20 Punta Arenas 2.192280 2020-08-03
## 21 Punta Arenas 2.234368 2020-08-07
## 22 Punta Arenas 2.287342 2020-08-10
## 23 Punta Arenas 2.332252 2020-08-14
## 24 Punta Arenas 2.370421 2020-08-17
## 25 Punta Arenas 2.416511 2020-08-21
## 26 Punta Arenas 2.492524 2020-08-24
## 27 Punta Arenas 2.560274 2020-08-28
## 28 Punta Arenas 2.623609 2020-08-31head(eee2, 50)## comuna_v2 R0_v2 fecha_v
## 1 Punta Arenas NA 2020-05-25
## 2 Punta Arenas NaN 2020-05-29
## 3 Punta Arenas 2.070059 2020-06-01
## 4 Punta Arenas 1.768603 2020-06-05
## 5 Punta Arenas 1.653190 2020-06-08
## 6 Punta Arenas 1.664298 2020-06-12
## 7 Punta Arenas 1.834052 2020-06-15
## 8 Punta Arenas 2.218345 2020-06-19
## 9 Punta Arenas 2.627003 2020-06-23
## 10 Punta Arenas 2.853279 2020-06-28
## 11 Punta Arenas 2.571856 2020-07-01
## 12 Punta Arenas 2.467455 2020-07-05
## 13 Punta Arenas 2.333386 2020-07-10
## 14 Punta Arenas 2.253917 2020-07-13
## 15 Punta Arenas 2.225834 2020-07-17
## 16 Punta Arenas 2.221018 2020-07-20
## 17 Punta Arenas 2.210058 2020-07-24
## 18 Punta Arenas 2.181527 2020-07-27
## 19 Punta Arenas 2.158974 2020-07-31
## 20 Punta Arenas 2.167281 2020-08-03
## 21 Punta Arenas 2.192280 2020-08-07
## 22 Punta Arenas 2.234368 2020-08-10
## 23 Punta Arenas 2.287342 2020-08-14
## 24 Punta Arenas 2.332252 2020-08-17
## 25 Punta Arenas 2.370421 2020-08-21
## 26 Punta Arenas 2.416511 2020-08-24
## 27 Punta Arenas 2.492524 2020-08-28
## 28 Punta Arenas 2.560274 2020-08-31Y y X de la regresion lineal
p <- ggplot(eee, aes(x=fecha_v, y=R0_v)) + geom_point(aes(y = R0_v), colour = "red")
p
o <- ggplot(eee2, aes(x=fecha_v, y=R0_v2)) + geom_point(aes(y = R0_v2), colour = "blue")
o
#income.graph <- income.graph + geom_smooth(method="loess", col="red", size = 0.5)
# income.graph