Analisis_COVID19.Rmd

---
title: "Analisis COVID19"
author: "Ricardo Alberich, Juan Gabriel Gomila y Arnau Mir"
date: "4/6/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(reticulate)
library(knitr)
#library(kableExtra)
library(tidyverse) # %>%
library(magrittr)  # %<>%
library(lubridate)
library(rnaturalearth)
library(car)
library(plotly)
library(xts)
library(dygraphs)

use_python("/Users/juangabriel/opt/anaconda3/bin/python3")
```

# Carga y limpieza preliminar de los datos

Los datos que se van a analizar en este documento, proceden de la compilación hecha por usuarios de [Kaggle](https://www.kaggle.com/imdevskp/corona-virus-report). La fecha del análisis empieza el 6 de Abril de 2020, utilizando la versión número 73 recopilada en la web anterior.

```{python}
import pandas as pd
datos = pd.read_csv("covid_19_clean_complete.csv")
datos.head(10)
```

```{r}
pd <- import("pandas")
datos <- pd$read_csv("covid_19_clean_complete.csv")
kable(head(datos, 10))
```

```{r}
datos <- read.csv("covid_19_clean_complete.csv", stringsAsFactors = FALSE)
datos %>% head(10) %>% kable()
```

## Estructura de los datos

```{r}
str(datos)
colnames(datos) = c("Provincia_Estado",
                    "Pais_Region",
                    "Latitud", # N+ o S-
                    "Longitud", # E+ o W-
                    "Fecha",
                    "Casos_Confirmados",
                    "Casos_Muertos",
                    "Casos_Recuperados"
                    )
datos %>% head() %>% kable() # %>% kable_styling()
```

* Cualitativas se convierten con `factor` o bien `as.factor`.
* Ordinales se convierten con `ordered`.
* Cuantitativos se convierten con `as.numeric`.

```{r}
datos$Provincia_Estado %<>% factor()
datos$Pais_Region %<>% factor()
#datos$Fecha %<>% as.Date(format="%m/%d/%y")
datos$Fecha %<>% mdy()
str(datos)
```


$$Confirmados = Muertos + Recuperados + Enfermos$$

```{r}
datos %<>%
  mutate(Casos_Enfermos = Casos_Confirmados - Casos_Muertos - Casos_Recuperados)

datos %>%
  filter(Casos_Confirmados > 10000) %>%
  head(10) %>%
  kable()

datos %>% 
  filter(Casos_Enfermos < 0) %>%
  arrange(Provincia_Estado, Fecha) %>%
  kable()

datos %>%
  filter(Provincia_Estado == "Hainan") %>%
  kable()

datos %>%
  filter(Provincia_Estado == "Hainan", Casos_Enfermos < 0) %>%
  mutate(Casos_Recuperados = Casos_Recuperados + Casos_Enfermos,
         Casos_Enfermos = 0) %>%
  kable()
```

# Análisis geográfico

```{r}
#datos_europa = datos[datos$Latitud > 38 & datos$Longitud > -25 & datos$Longitud < 30 , ]

datos_europa = datos %>%
  filter(Latitud > 38, between(Longitud, -25, 30))

nrow(datos_europa)

table(datos_europa$Pais_Region) %>%
  as.data.frame() %>%
  filter(Freq > 0) %>%
  kable()


datos_europa %>%
  filter(Fecha == ymd("2020-03-15")) %>%
  kable()
```

$$d(x, y) = \sqrt{(x_{Lat}-y_{Lat})^2 + (x_{Long}-y_{Long})^2}$$
```{r}
distancia_grados = function(x, y){
  sqrt((x[1]-y[1])^2 + (x[2]-y[2])^2)
}

distancia_grados_potsdam = function(x){
  potsdam = c(52.366956, 13.906734)
  distancia_grados(x, potsdam)
}

dist_potsdam = apply(cbind(datos_europa$Latitud, datos_europa$Longitud),
                     MARGIN = 1, 
                     FUN = distancia_grados_potsdam)

datos_europa %<>%
  mutate(dist_potsdam = dist_potsdam)

datos_europa %>%
  filter(between(Fecha, dmy("2-3-2020"), dmy("7-3-2020")),
         dist_potsdam < 4) %>%
  kable()
```


```{r}
world <- ne_countries(scale = "medium", returnclass = "sf")

datos$Pais_Region = factor(datos$Pais_Region, levels = c(levels(datos$Pais_Region), "United States"))

datos[datos$Pais_Region=="US",]$Pais_Region = "United States"

world %>%
  inner_join(datos, by = c("name" = "Pais_Region")) %>%
  filter(Fecha == dmy("30-03-2020")) %>%
  ggplot() +
  geom_sf(color = "black", aes(fill = Casos_Confirmados)) +
#  coord_sf(crs="+proj=laea +lat_0=50 +lon_0=10 +units=m +ellps=GRS80") +
  scale_fill_viridis_c(option="plasma", trans = "sqrt") +
  xlab("Longitud") + ylab("Latitud") +
  ggtitle("Mapa del mundo ", subtitle = "COVID 19") -> g

ggplotly(g)
```

```{r}
datos %>%
  filter(Fecha == dmy("30-03-2020")) %>%
  ggplot(aes(Longitud, Latitud)) +
  geom_point(aes(size = log(Casos_Confirmados+1), colour = log(Casos_Muertos+1))) +
  coord_fixed() +
  theme(legend.position = "bottom") -> g

ggplotly(g)
```

```{r}

thr = 1000

datos %>%
  filter(Fecha == ymd("2020-04-05"),
         Casos_Confirmados > thr) %>%
  mutate(Prop_Muertos = Casos_Muertos / Casos_Confirmados, 
         Ranking = dense_rank(desc(Prop_Muertos))) %>%
  arrange(Ranking) %>%
  head(20) %>%
  kable()
```


```{r}
datos$lat_class = cut(datos$Latitud, 
                      breaks = seq(from = -90, to = 90, by = 10))
datos$long_class = cut(datos$Longitud, 
                       breaks = seq(from = -180, to = 180, by = 10))
tt = table(datos$lat_class, datos$long_class)
tt = tt[nrow(tt):1, ]
mosaicplot(t(tt), shade = TRUE)
```


## Análisis de datos temporal

```{r}
datos_por_fecha = aggregate(
  cbind(Casos_Confirmados, Casos_Muertos, Casos_Recuperados) ~ Fecha,
  data = datos, 
  FUN = sum
)
datos_por_fecha$Casos_Enfermos = datos_por_fecha$Casos_Confirmados - datos_por_fecha$Casos_Muertos - datos_por_fecha$Casos_Recuperados
head(datos_por_fecha)
tail(datos_por_fecha)

barplot(Casos_Confirmados ~ Fecha, data = datos_por_fecha)

plot(Casos_Confirmados ~ Fecha, data = datos_por_fecha, col = "blue", type = "l", main = "Casos documentados por día en todo el mundo", xlab = "Fecha", ylab = "Número de personas", log = "y")
lines(Casos_Muertos ~ Fecha, data = datos_por_fecha, col = "red")
lines(Casos_Recuperados ~ Fecha, data = datos_por_fecha, col = "green")

legend("topleft", c("Confirmados", "Muertos", "Recuperados"), 
       col = c("blue", "red", "green"), pch = 1, lwd = 2)

datos_por_fecha_ts <- xts(x = datos_por_fecha[, 2:5],
                          order.by = datos_por_fecha$Fecha)
dygraph(datos_por_fecha_ts) %>%
  dyOptions(labelsUTC = TRUE, labelsKMB = TRUE,
            fillGraph = TRUE, fillAlpha = 0.05, 
            drawGrid = FALSE, colors = "#D9AE55") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE) %>%
  dyRoller(rollPeriod = 2)

```

```{r}
datos_spain = datos %>% 
  filter(Pais_Region == "Spain") %>%
  select(Fecha, starts_with("Casos_"))

plot(x = datos_spain$Fecha, y = datos_spain$Casos_Confirmados,
     main = "Casos confirmados en España", type = "s", 
     col = "blue", lwd = 2)

datos_por_fecha_ts <- xts(x = datos_spain[, 2:5],
                          order.by = datos_spain$Fecha)
dygraph(datos_por_fecha_ts) %>%
  dyOptions(labelsUTC = TRUE, labelsKMB = TRUE,
            fillGraph = TRUE, fillAlpha = 0.05, 
            drawGrid = FALSE, colors = "#D9AE55") %>%
  dyRangeSelector() %>%
  dyCrosshair(direction = "vertical") %>%
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE) %>%
  dyRoller(rollPeriod = 2)

barplot(as.matrix(t(datos_spain[, 3:5])),
        names = datos_spain$Fecha, 
        col = c("red", "green", "yellow"),
        main = "Estudio de casos por tipo en España", 
        xlab ="Fecha", ylab = "Número de personas")
legend("topleft", c("Muertos", "Recuperados", "Enfermos"),
       col = c("red", "green", "yellow"), lwd = 2, pch = 1
       )
```


```{r}
datos_spain %<>%
  mutate(Nuevos_Casos_Confirmados = Casos_Confirmados - lag(Casos_Confirmados, n = 1),
         Nuevos_Casos_Muertos = Casos_Muertos - lag(Casos_Muertos, n = 1),
         Nuevos_Casos_Recuperados = Casos_Recuperados - lag(Casos_Recuperados, n = 1)
         )

plot(Nuevos_Casos_Confirmados ~ Fecha, data = datos_spain,
     type = "l", col ="blue", 
     xlab = "Fecha", ylab = "Nuevos casos", 
     main = "Nuevos registros en España")
lines(Nuevos_Casos_Muertos ~ Fecha, data = datos_spain,
      type = "l", col = "red")
lines(Nuevos_Casos_Recuperados ~ Fecha, data = datos_spain,
      type = "l", col = "green")

legend("topleft", c("Confirmados", "Muertos", "Recuperados"), 
       col = c("blue", "red", "green"), 
       lwd = 2, pch = 1)
```

## Análisis por Cohortes

```{r}
primer_contagio = datos %>%
  group_by(Pais_Region) %>%
  filter(Casos_Confirmados > 0) %>%
  summarise(Primer_Contagio = min(Fecha)-1)


data_first = datos %>%
  inner_join(primer_contagio, by = "Pais_Region") %>%
  mutate(Dias_Desde_PC = as.numeric(Fecha - Primer_Contagio)) %>%
  filter(Dias_Desde_PC >= 0) %>%
  group_by(Dias_Desde_PC, Pais_Region) %>%
  summarise(Casos_Confirmados = sum(Casos_Confirmados),
            Casos_Muertos = sum(Casos_Muertos),
            Casos_Recuperados = sum(Casos_Recuperados),
            Casos_Enfermos = sum(Casos_Enfermos))


data_first %>%
  #filter(Pais_Region %in% c("Spain", "Italy", "China", "United States", "Germany")) %>%
  ggplot(aes(x = Dias_Desde_PC, y = Casos_Confirmados)) +
  geom_line(aes(col = Pais_Region)) +  
  xlab("Días desde el primer contagio") +
  ylab("Número de personas contagiadas") + 
  ggtitle("Análisis por Cohortes") +
  theme(legend.position = "none") -> g

ggplotly(g)
```













## Modelos de regresión simple

* $x$: Variable Independiente: número de días desde el origen de la pandemia
* $y$: Variable Dependiente: número de casos confirmados. 

$$y = f(x)$$

```{r}
datos_spain$Dias = as.numeric(datos_spain$Fecha - dmy("22/01/2020"))
```

### Regresión Lineal

$$y = ax+b, a,b\in \mathbb R$$

$$min_{a,b\in\mathbb R} \sum_{i=1}^n (y_i-(ax_i+b))^2$$

```{r}
mod1 <- lm(Casos_Confirmados ~ Dias, data = datos_spain)
summary(mod1)
```

$$Casos\ Confirmados = `r mod1$coefficients[2]` Dias + `r mod1$coefficients[1]`$$

```{r}
plot(datos_spain$Dias, datos_spain$Casos_Confirmados)
abline(mod1, col = "red")

plot(mod1$residuals~mod1$fitted.values, xlab = "Valores Ajustados", ylab = "Residuos del modelo")

residuos = mod1$residuals

qqPlot(residuos, distribution = "norm", 
       mean = mean(residuos), sd = sd(residuos))
```

### Regresión exponencial
$$log(y) = ax+b, a,b \in \mathbb R$$
$$y = e^{ax+b} = m e^{ax}$$


```{r}
mod2 <- lm(log(Casos_Confirmados) ~ Dias, data = datos_spain[datos_spain$Casos_Confirmados>0, ])
summary(mod2)
```

$$Casos\ Confirmados = `r exp(mod2$coefficients[1])` \cdot e^{`r mod2$coefficients[2]`\cdot x}$$

```{r}
plot(datos_spain$Dias, datos_spain$Casos_Confirmados)
lines(exp(mod2$coefficients[1])*exp(mod2$coefficients[2]*datos_spain$Dias), col = "blue")

plot(mod2$residuals ~ mod2$fitted.values, xlab = "Valores Ajustados", ylab = "Residuos del modelo")
residuos = mod2$residuals
qqPlot(residuos, distribution = "norm", 
       mean = mean(residuos), sd = sd(residuos))
```

### Modelo potencial
$$log(y) = a\cdot log(x)+b, a,b\in \mathbb R$$
$$y = e^{a\cdot log(x)+b} = e^b\cdot e^{log(x)^a} = m\cdot x^a$$

```{r}
mod3 <- lm(log(Casos_Confirmados) ~ log(Dias),
           data = datos_spain[datos_spain$Casos_Confirmados > 0, ])
summary(mod3)
```

$$Casos\ Confirmados = `r exp(mod3$coefficients[1])`\cdot Dias^{`r mod2$coefficients[2]`}$$

```{r}
plot(datos_spain$Dias, datos_spain$Casos_Confirmados)
lines(exp(mod3$coefficients[1])*datos_spain$Dias^mod3$coefficients[2], col = "green")

plot(mod3$residuals~mod3$fitted.values, 
     xlab = "Valores Ajustados", ylab = "Residuos del modelo")
residuos = mod3$residuals
qqPlot(residuos, distribution = "norm", mean = mean(residuos), sd = sd(residuos))
```


```{r}
mi_model <- lm(log(Casos_Confirmados) ~ Dias + log(Dias) + I(Dias^2) + I(Dias^3) + sqrt(Dias),
               data = datos_spain[datos_spain$Casos_Confirmados> 0, ] )
summary(mi_model)
```

```{r}
start_date = ymd('2020-01-22')
end_date = ymd('2020-04-30')

dates = seq(start_date+1, end_date, by = "1 day")
days_since_start = as.numeric(dates - start_date)
new_data = data.frame(Dias = days_since_start)

pred1 = predict(mod1, newdata = new_data)
pred2 = exp(predict(mod2, newdata = new_data))
pred3 = exp(predict(mod3, newdata = new_data))
pred4 = exp(predict(mi_model, newdata = new_data))

datos_por_fecha_ts = xts(x = data.frame(Real = c(datos_spain$Casos_Confirmados, rep(NA,length(pred1)- length(datos_spain$Casos_Confirmados))),
                                        Mod_Lin = pred1,
                                        #Mod_Exp = pred2, 
                                        Mod_Pot = pred3, 
                                        Mod_Mixt = pred4),
                         order.by = dates)


dygraph(datos_por_fecha_ts)
```