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. La fecha del análisis empieza el 30 de abril de 2020, utilizando la versión 83 recopilada en la web anterior

Aqui se pueden combinar R y Python

#import pandas as pd
#datos = pd.read_csv("covid_19_clean_complete.csv")
#datos.head(10)#

#pd <- import("pandas")
#datos <- pd$read_csv("covid_19_clean_complete.csv")
#kable(head(datos))

Vamos a hacerlo con el metodo tradicional

path <- '~/Desktop/DATA ANALYST, SCIENCE y CERTIFICADOS/DS4B VIRTUAL CAMP + ESTADÍSTICA + CODING/Predicción y cuadro de mando CORONAVIRUS/CURSO UDEMY EN R CORONAVIRUS/COVID19/covid_19_clean_complete.csv'
read_lines(path,n_max = 5)

## [1] "Province/State,Country/Region,Lat,Long,Date,Confirmed,Deaths,Recovered"
## [2] ",Afghanistan,33.0,65.0,1/22/20,0,0,0"                                  
## [3] ",Albania,41.1533,20.1683,1/22/20,0,0,0"                                
## [4] ",Algeria,28.0339,1.6596,1/22/20,0,0,0"                                 
## [5] ",Andorra,42.5063,1.5218,1/22/20,0,0,0"

datos <- read.csv("covid_19_clean_complete.csv")

glimpse(datos)

## Observations: 25,938
## Variables: 8
## $ Province.State <fct> , , , , , , , , Australian Capital Territory, New Sout…
## $ Country.Region <fct> Afghanistan, Albania, Algeria, Andorra, Angola, Antigu…
## $ Lat            <dbl> 33.0000, 41.1533, 28.0339, 42.5063, -11.2027, 17.0608,…
## $ Long           <dbl> 65.0000, 20.1683, 1.6596, 1.5218, 17.8739, -61.7964, -…
## $ Date           <fct> 1/22/20, 1/22/20, 1/22/20, 1/22/20, 1/22/20, 1/22/20, …
## $ Confirmed      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Deaths         <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Recovered      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …

view(datos)
datos %>% head(10) %>% kable()

Province.State	Country.Region	Lat	Long	Date
	Afghanistan	33.0000	65.0000	1/22/20
	Albania	41.1533	20.1683	1/22/20
	Algeria	28.0339	1.6596	1/22/20
	Andorra	42.5063	1.5218	1/22/20
	Angola	-11.2027	17.8739	1/22/20
	Antigua and Barbuda	17.0608	-61.7964	1/22/20
	Argentina	-38.4161	-63.6167	1/22/20
	Armenia	40.0691	45.0382	1/22/20
Australian Capital Territory	Australia	-35.4735	149.0124	1/22/20
New South Wales	Australia	-33.8688	151.2093	1/22/20

Estructura de los datos

str(datos)

## 'data.frame':    25938 obs. of  8 variables:
##  $ Province.State: Factor w/ 81 levels "","Alberta","Anguilla",..: 1 1 1 1 1 1 1 1 6 49 ...
##  $ Country.Region: Factor w/ 185 levels "Afghanistan",..: 1 2 3 4 5 6 7 8 9 9 ...
##  $ Lat           : num  33 41.2 28 42.5 -11.2 ...
##  $ Long          : num  65 20.17 1.66 1.52 17.87 ...
##  $ Date          : Factor w/ 99 levels "1/22/20","1/23/20",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Confirmed     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Deaths        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Recovered     : int  0 0 0 0 0 0 0 0 0 0 ...

colnames(datos) = c("Provincia_Estado",
                    "Pais_Region",
                    "Latitud",
                    "Longitud",
                    "Fecha",
                    "Casos_Confirmados",
                    "Casos_Muertos",
                    "Casos_Recuperados")
datos %>% head() %>% kable() %>% kable_styling()

Pais_Region	Latitud	Longitud	Fecha
Afghanistan	33.0000	65.0000	1/22/20
Albania	41.1533	20.1683	1/22/20
Algeria	28.0339	1.6596	1/22/20
Andorra	42.5063	1.5218	1/22/20
Angola	-11.2027	17.8739	1/22/20
Antigua and Barbuda	17.0608	-61.7964	1/22/20

Las variables cualitativas se convierten con factor o as.factor`
Las variables ordinales se convierten con ordered
Las variables cuantitavias se convierten con as.numeric

** En este caso todos los datos se han cargado con el tipo correcto porque hemos utilizado el read.csv, el cual tiene por defecto el paramentro StringAsFactor como TRUE, esto hace que todos los campos que vienen cn formato de texto o entre comillas sean directamente factores. Los unicos casos donde esto no nos interesa es cuando el String identifica al objeto; Por ejemplo un numero de matricula.**

Si ponemos StringAsFactors = FALSE se nos cargaran como character seguramente, por lo que podríamos modificar esos campos de la siguiente manera:

#datos$Provincia_Estado = as.factor(datos$Provincia_Estado)
#datos$Pais_Region = as.factor(datos$Pais_Region)

#str(datos)

Con el doble PIPE se puede hacer la misma accion que arriba, pero mas facil. Esto hace que la accion vaya hacia adelante y despues vuelva.

datos$Provincia_Estado %<>% as.factor()
datos$Pais_Region %<>% as.factor()
#datos$Fecha %<>% as.Date(format="%m/%d/%y") # "%m/%d/%y %H:%M:%S" si tenemos hora
datos$Fecha %<>% mdy()
str(datos)

## 'data.frame':    25938 obs. of  8 variables:
##  $ Provincia_Estado : Factor w/ 81 levels "","Alberta","Anguilla",..: 1 1 1 1 1 1 1 1 6 49 ...
##  $ Pais_Region      : Factor w/ 185 levels "Afghanistan",..: 1 2 3 4 5 6 7 8 9 9 ...
##  $ Latitud          : num  33 41.2 28 42.5 -11.2 ...
##  $ Longitud         : num  65 20.17 1.66 1.52 17.87 ...
##  $ Fecha            : Date, format: "2020-01-22" "2020-01-22" ...
##  $ Casos_Confirmados: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Casos_Muertos    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Casos_Recuperados: int  0 0 0 0 0 0 0 0 0 0 ...

Podemos hacer operacion con las variables fechas de la libreria lubridate.

d2 <- dmy("28/04/20")
d1 <- dmy("21/01/20")
d2-d1 # Expresa la diferencia de dias

## Time difference of 98 days

is.difftime(d2-d1)

## [1] TRUE

days(d2-d1)

## [1] "98d 0H 0M 0S"

\[ CasosConfirmados = Muertos + Recuperados + Enfermos \] Vamos a crear una nueva columna de Casos enfermos y despues ver si tenemos datos anomalos

datos %<>%
  mutate(Casos_Enfermos = Casos_Confirmados - Casos_Muertos - Casos_Recuperados)

datos %>%
  filter(Casos_Confirmados > 10000) %>% 
  tail() %>%
  kable() %>%
  kable_styling()

	Pais_Region	Latitud	Longitud	Fecha	Casos_Confirmados	Casos_Muertos	Casos_Recuperados	Casos_Enfermos
868	Sweden	63.0000	16.0000	2020-04-29	20302	2462	1005	16835
869	Switzerland	46.8182	8.2275	2020-04-29	29407	1716	22600	5091
870	Turkey	38.9637	35.2433	2020-04-29	117589	3081	44040	70468
871	United Arab Emirates	24.0000	54.0000	2020-04-29	11929	98	2329	9502
872	United Kingdom	55.3781	-3.4360	2020-04-29	165221	26097	0	139124
873	US	37.0902	-95.7129	2020-04-29	1039909	60967	120720	858222

datos %>%
  filter(Casos_Enfermos < 0) %>% # Miramos si hay casos anomalos
  arrange(Provincia_Estado, Fecha) %>%
  kable() %>%
  kable_styling()

Provincia_Estado	Pais_Region	Latitud	Longitud	Fecha	Casos_Confirmados	Casos_Muertos	Casos_Recuperados	Casos_Enfermos
Diamond Princess	Canada	0.0000	0.0000	2020-03-22	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-03-23	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-03-24	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-03-25	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-03-26	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-03-27	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-03-28	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-03-29	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-03-30	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-03-31	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-01	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-02	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-03	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-04	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-05	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-06	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-07	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-08	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-09	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-10	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-11	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-12	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-13	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-14	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-15	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-16	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-17	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-18	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-19	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-20	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-21	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-22	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-23	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-24	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-25	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-26	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-27	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-28	0	1	0	-1
Diamond Princess	Canada	0.0000	0.0000	2020-04-29	0	1	0	-1
Hainan	China	19.1959	109.7453	2020-03-24	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-25	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-26	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-27	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-28	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-29	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-30	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-31	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-04-01	168	6	168	-6

datos %>%
  filter(Provincia_Estado == 'Hainan') %>%
  kable() %>%
  kable_styling()

Provincia_Estado	Pais_Region	Latitud	Longitud	Fecha	Casos_Confirmados	Casos_Muertos	Casos_Recuperados	Casos_Enfermos
Hainan	China	19.1959	109.7453	2020-01-22	4	0	0	4
Hainan	China	19.1959	109.7453	2020-01-23	5	0	0	5
Hainan	China	19.1959	109.7453	2020-01-24	8	0	0	8
Hainan	China	19.1959	109.7453	2020-01-25	19	0	0	19
Hainan	China	19.1959	109.7453	2020-01-26	22	0	0	22
Hainan	China	19.1959	109.7453	2020-01-27	33	1	0	32
Hainan	China	19.1959	109.7453	2020-01-28	40	1	0	39
Hainan	China	19.1959	109.7453	2020-01-29	43	1	0	42
Hainan	China	19.1959	109.7453	2020-01-30	46	1	1	44
Hainan	China	19.1959	109.7453	2020-01-31	52	1	1	50
Hainan	China	19.1959	109.7453	2020-02-01	62	1	1	60
Hainan	China	19.1959	109.7453	2020-02-02	64	1	4	59
Hainan	China	19.1959	109.7453	2020-02-03	72	1	4	67
Hainan	China	19.1959	109.7453	2020-02-04	80	1	5	74
Hainan	China	19.1959	109.7453	2020-02-05	99	1	5	93
Hainan	China	19.1959	109.7453	2020-02-06	106	1	8	97
Hainan	China	19.1959	109.7453	2020-02-07	117	2	10	105
Hainan	China	19.1959	109.7453	2020-02-08	124	2	14	108
Hainan	China	19.1959	109.7453	2020-02-09	131	3	19	109
Hainan	China	19.1959	109.7453	2020-02-10	138	3	19	116
Hainan	China	19.1959	109.7453	2020-02-11	144	3	20	121
Hainan	China	19.1959	109.7453	2020-02-12	157	4	27	126
Hainan	China	19.1959	109.7453	2020-02-13	157	4	30	123
Hainan	China	19.1959	109.7453	2020-02-14	159	4	43	112
Hainan	China	19.1959	109.7453	2020-02-15	162	4	39	119
Hainan	China	19.1959	109.7453	2020-02-16	162	4	52	106
Hainan	China	19.1959	109.7453	2020-02-17	163	4	59	100
Hainan	China	19.1959	109.7453	2020-02-18	163	4	79	80
Hainan	China	19.1959	109.7453	2020-02-19	168	4	84	80
Hainan	China	19.1959	109.7453	2020-02-20	168	4	86	78
Hainan	China	19.1959	109.7453	2020-02-21	168	4	95	69
Hainan	China	19.1959	109.7453	2020-02-22	168	4	104	60
Hainan	China	19.1959	109.7453	2020-02-23	168	5	106	57
Hainan	China	19.1959	109.7453	2020-02-24	168	5	116	47
Hainan	China	19.1959	109.7453	2020-02-25	168	5	124	39
Hainan	China	19.1959	109.7453	2020-02-26	168	5	129	34
Hainan	China	19.1959	109.7453	2020-02-27	168	5	131	32
Hainan	China	19.1959	109.7453	2020-02-28	168	5	133	30
Hainan	China	19.1959	109.7453	2020-02-29	168	5	148	15
Hainan	China	19.1959	109.7453	2020-03-01	168	5	149	14
Hainan	China	19.1959	109.7453	2020-03-02	168	5	151	12
Hainan	China	19.1959	109.7453	2020-03-03	168	5	155	8
Hainan	China	19.1959	109.7453	2020-03-04	168	5	158	5
Hainan	China	19.1959	109.7453	2020-03-05	168	6	158	4
Hainan	China	19.1959	109.7453	2020-03-06	168	6	158	4
Hainan	China	19.1959	109.7453	2020-03-07	168	6	158	4
Hainan	China	19.1959	109.7453	2020-03-08	168	6	159	3
Hainan	China	19.1959	109.7453	2020-03-09	168	6	159	3
Hainan	China	19.1959	109.7453	2020-03-10	168	6	159	3
Hainan	China	19.1959	109.7453	2020-03-11	168	6	159	3
Hainan	China	19.1959	109.7453	2020-03-12	168	6	160	2
Hainan	China	19.1959	109.7453	2020-03-13	168	6	160	2
Hainan	China	19.1959	109.7453	2020-03-14	168	6	160	2
Hainan	China	19.1959	109.7453	2020-03-15	168	6	160	2
Hainan	China	19.1959	109.7453	2020-03-16	168	6	161	1
Hainan	China	19.1959	109.7453	2020-03-17	168	6	161	1
Hainan	China	19.1959	109.7453	2020-03-18	168	6	161	1
Hainan	China	19.1959	109.7453	2020-03-19	168	6	161	1
Hainan	China	19.1959	109.7453	2020-03-20	168	6	161	1
Hainan	China	19.1959	109.7453	2020-03-21	168	6	161	1
Hainan	China	19.1959	109.7453	2020-03-22	168	6	161	1
Hainan	China	19.1959	109.7453	2020-03-23	168	6	161	1
Hainan	China	19.1959	109.7453	2020-03-24	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-25	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-26	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-27	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-28	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-29	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-30	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-03-31	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-04-01	168	6	168	-6
Hainan	China	19.1959	109.7453	2020-04-02	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-03	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-04	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-05	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-06	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-07	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-08	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-09	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-10	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-11	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-12	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-13	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-14	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-15	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-16	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-17	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-18	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-19	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-20	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-21	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-22	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-23	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-24	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-25	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-26	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-27	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-28	168	6	162	0
Hainan	China	19.1959	109.7453	2020-04-29	168	6	162	0

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

##   Provincia_Estado Pais_Region Latitud Longitud      Fecha Casos_Confirmados
## 1           Hainan       China 19.1959 109.7453 2020-03-24               168
## 2           Hainan       China 19.1959 109.7453 2020-03-25               168
## 3           Hainan       China 19.1959 109.7453 2020-03-26               168
## 4           Hainan       China 19.1959 109.7453 2020-03-27               168
## 5           Hainan       China 19.1959 109.7453 2020-03-28               168
## 6           Hainan       China 19.1959 109.7453 2020-03-29               168
## 7           Hainan       China 19.1959 109.7453 2020-03-30               168
## 8           Hainan       China 19.1959 109.7453 2020-03-31               168
## 9           Hainan       China 19.1959 109.7453 2020-04-01               168
##   Casos_Muertos Casos_Recuperados Casos_Enfermos
## 1             6               168              0
## 2             6               168              0
## 3             6               168              0
## 4             6               168              0
## 5             6               168              0
## 6             6               168              0
## 7             6               168              0
## 8             6               168              0
## 9             6               168              0
##   Casos_Recuperados == Casos_Recuperados + Casos_Enfermos
## 1                                                   FALSE
## 2                                                   FALSE
## 3                                                   FALSE
## 4                                                   FALSE
## 5                                                   FALSE
## 6                                                   FALSE
## 7                                                   FALSE
## 8                                                   FALSE
## 9                                                   FALSE

Como los casos de Hainan ya hemos descubierto que se debe a un error, lo hemos modificado para no borrarlo.

Análisis geográfico de Europa

Como yo quiero todas las filas que pertenecen a Europa, tengo que filtrar en [X,] ya que el seguno concepto son las columnas

# datos_europa <- datos[datos$Latitud > 38 & datos$Longitud > -25 & datos$Longitud < 30, ] # Forma normal de filtrarlo

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

nrow(datos_europa)

## [1] 4455

table(datos_europa$Pais_Region) %>%
  as.data.frame() %>%   # Transformamos a data frame porque el filter solo se aplica a df
  filter(Freq > 0) %>%
  kable() %>%
  kable_styling()

Var1	Freq
Albania	99
Andorra	99
Austria	99
Belarus	99
Belgium	99
Bosnia and Herzegovina	99
Bulgaria	99
Croatia	99
Czechia	99
Denmark	198
Estonia	99
Finland	99
France	99
Germany	99
Greece	99
Holy See	99
Hungary	99
Iceland	99
Ireland	99
Italy	99
Kosovo	99
Latvia	99
Liechtenstein	99
Lithuania	99
Luxembourg	99
Moldova	99
Monaco	99
Montenegro	99
Netherlands	99
North Macedonia	99
Norway	99
Poland	99
Portugal	99
Romania	99
San Marino	99
Serbia	99
Slovakia	99
Slovenia	99
Spain	99
Sweden	99
Switzerland	99
United Kingdom	297

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

Provincia_Estado	Pais_Region	Latitud	Longitud	Fecha	Casos_Confirmados	Casos_Muertos	Casos_Recuperados	Casos_Enfermos
	Albania	41.15330	20.16830	2020-03-15	42	1	0	41
	Andorra	42.50630	1.52180	2020-03-15	1	0	1	0
	Austria	47.51620	14.55010	2020-03-15	860	1	6	853
	Belarus	53.70980	27.95340	2020-03-15	27	0	3	24
	Belgium	50.83330	4.00000	2020-03-15	886	4	1	881
	Bosnia and Herzegovina	43.91590	17.67910	2020-03-15	24	0	0	24
	Bulgaria	42.73390	25.48580	2020-03-15	51	2	0	49
	Croatia	45.10000	15.20000	2020-03-15	49	0	1	48
	Czechia	49.81750	15.47300	2020-03-15	253	0	0	253
Faroe Islands	Denmark	61.89260	-6.91180	2020-03-15	11	0	0	11
	Denmark	56.26390	9.50180	2020-03-15	864	2	1	861
	Estonia	58.59530	25.01360	2020-03-15	171	0	1	170
	Finland	64.00000	26.00000	2020-03-15	244	0	10	234
	France	46.22760	2.21370	2020-03-15	4499	91	12	4396
	Germany	51.00000	9.00000	2020-03-15	5795	11	46	5738
	Greece	39.07420	21.82430	2020-03-15	331	4	8	319
	Holy See	41.90290	12.45340	2020-03-15	1	0	0	1
	Hungary	47.16250	19.50330	2020-03-15	32	1	1	30
	Iceland	64.96310	-19.02080	2020-03-15	171	5	8	158
	Ireland	53.14240	-7.69210	2020-03-15	129	2	0	127
	Italy	43.00000	12.00000	2020-03-15	24747	1809	2335	20603
	Latvia	56.87960	24.60320	2020-03-15	30	0	1	29
	Liechtenstein	47.14000	9.55000	2020-03-15	4	0	0	4
	Lithuania	55.16940	23.88130	2020-03-15	12	0	1	11
	Luxembourg	49.81530	6.12960	2020-03-15	59	1	0	58
	Moldova	47.41160	28.36990	2020-03-15	23	0	0	23
	Monaco	43.73330	7.41670	2020-03-15	2	0	0	2
	Montenegro	42.50000	19.30000	2020-03-15	0	0	0	0
	Netherlands	52.13260	5.29130	2020-03-15	1135	20	2	1113
	North Macedonia	41.60860	21.74530	2020-03-15	14	0	1	13
	Norway	60.47200	8.46890	2020-03-15	1221	3	1	1217
	Poland	51.91940	19.14510	2020-03-15	119	3	0	116
	Portugal	39.39990	-8.22450	2020-03-15	245	0	2	243
	Romania	45.94320	24.96680	2020-03-15	131	0	9	122
	San Marino	43.94240	12.45780	2020-03-15	101	5	4	92
	Serbia	44.01650	21.00590	2020-03-15	48	0	0	48
	Slovakia	48.66900	19.69900	2020-03-15	54	0	0	54
	Slovenia	46.15120	14.99550	2020-03-15	219	1	0	218
	Spain	40.00000	-4.00000	2020-03-15	7798	289	517	6992
	Sweden	63.00000	16.00000	2020-03-15	1022	3	1	1018
	Switzerland	46.81820	8.22750	2020-03-15	2200	14	4	2182
Channel Islands	United Kingdom	49.37230	-2.36440	2020-03-15	3	0	0	3
Isle of Man	United Kingdom	54.23610	-4.54810	2020-03-15	0	0	0	0
	United Kingdom	55.37810	-3.43600	2020-03-15	1140	21	18	1101
	Kosovo	42.60264	20.90298	2020-03-15	0	0	0	0

\[ d(x, y) = /sqrt\] Aqui se hace una funcion como distancia euclidea para ver a que distancia tenia contagiados desde postdam

distancia_grados = function(x, y){
  sqrt((x[1]-y[1])^2 + (x[2]- y[2])^2)
}

distancia_grados_postdam = function(x){
  postdam = c(52.366956, 13.906734)
  distancia_grados(x, postdam)
}

dist_postdam = apply(cbind(datos_europa$Latitud, datos_europa$Longitud),
                     MARGIN = 1,
                     FUN = distancia_grados_postdam) # CBIND junta las columnas

datos_europa %<>%
  mutate(dist_postdam = dist_postdam)

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

Pais_Region	Latitud	Longitud	Fecha	Casos_Confirmados	Casos_Enfermos	dist_postdam
Czechia	49.8175	15.473	2020-03-02	3	3	2.992142
Czechia	49.8175	15.473	2020-03-03	5	5	2.992142
Czechia	49.8175	15.473	2020-03-04	8	8	2.992142
Czechia	49.8175	15.473	2020-03-05	12	12	2.992142
Czechia	49.8175	15.473	2020-03-06	18	18	2.992142
Czechia	49.8175	15.473	2020-03-07	19	19	2.992142

Vamos a hacer un mapa con rnaturalearth

#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")) %>% # Aqui cruzamos la tabla world con la de datos para que coincidan el name de la tabla world con el nombre de mi tabla datos
  #ggplot() +
  #geom_sf(color = "black", aes(fill = Casos_Confirmados)) +
  #coord_sf(crs = "+proj=laea +lat_0=50 +lan_0=10[ +units=m +ellps=GRS80") +
  #scale_fill_viridis_c(option = "plasma", trans = "sqrt") +
  #xlab("Longitud") + ylab("Latitud") +
  #ggtitle("Mapa del mundo", subtitle = "COVID19")

Si lo que quiero es representar un mapa directamente con los datos que yo tengo y de forma minimalista, procedemos de la siguiente manera:

datos %>%
  filter(Fecha == dmy("30-03-2020")) %>%
  ggplot(aes(Longitud, Latitud)) +
  geom_point(aes(size = Casos_Confirmados, color = Casos_Muertos)) +
  coord_fixed() +
  theme(legend.position = "bottom")

# Para que la estetica de los puntos sea directamente proporcional a los casos confirmados. El coord fixed estira el mapa

# Se puede pasar a logaritmo los casos confirmados y muertos para tener el mapa con puntos mas grandes: size = log(Casos_Confirmados +1), color = log(Casos_Muertos+1))

MAPAS INTERACTIVOS con GGPLOTLY

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")) %>% # Aqui cruzamos la tabla world con la de datos para que coincidan el name de la tabla world con el nombre de mi tabla datos
  ggplot() +
  geom_sf(color = "black", aes(fill = Casos_Confirmados)) +
  # coord_sf(crs = "+proj=laea +lat_0=50 +lan_0=10[ +units=m +ellps=GRS80") +
  scale_fill_viridis_c(option = "plasma", trans = "sqrt") +
  xlab("Longitud") + ylab("Latitud") +
  ggtitle("Mapa del mundo", subtitle = "COVID19") -> g

## Warning: Column `name`/`Pais_Region` joining character vector and factor,
## coercing into character vector

ggplotly(g)

Ranking sobre la proporcion de muertos en un dia determinado, top 10/20

datos %>%
  filter(Fecha == ymd("2020-03-25"),
         Casos_Confirmados > 0) %>%
  mutate(Prop_Muertos = Casos_Muertos / Casos_Confirmados *100,
         Ranking = dense_rank(desc(Prop_Muertos))) %>%
  arrange(Ranking) %>% 
  head(20) %>% # Anade columnas al dataset, no se guardan
  kable() %>% kable_styling()

Provincia_Estado	Pais_Region	Latitud	Longitud	Fecha	Casos_Confirmados	Casos_Muertos	Casos_Recuperados	Casos_Enfermos	Prop_Muertos	Ranking
	Gambia	13.4432	-15.3101	2020-03-25	3	1	0	2	33.333333	1
	Sudan	12.8628	30.2176	2020-03-25	3	1	0	2	33.333333	1
	Zimbabwe	-20.0000	30.0000	2020-03-25	3	1	0	2	33.333333	1
	Cabo Verde	16.5388	-23.0418	2020-03-25	4	1	0	3	25.000000	2
	Guyana	5.0000	-58.7500	2020-03-25	5	1	0	4	20.000000	3
	Gabon	-0.8037	11.6094	2020-03-25	6	1	0	5	16.666667	4
Curacao	Netherlands	12.1696	-68.9900	2020-03-25	6	1	0	5	16.666667	4
	Niger	17.6078	8.0817	2020-03-25	7	1	0	6	14.285714	5
	Bangladesh	23.6850	90.3563	2020-03-25	39	5	7	27	12.820513	6
Cayman Islands	United Kingdom	19.3133	-81.2546	2020-03-25	8	1	0	7	12.500000	7
	San Marino	43.9424	12.4578	2020-03-25	208	21	4	183	10.096154	8
	Italy	43.0000	12.0000	2020-03-25	74386	7503	9362	57521	10.086575	9
	Iraq	33.0000	44.0000	2020-03-25	346	29	103	214	8.381503	10
	Paraguay	-23.4425	-58.4438	2020-03-25	37	3	0	34	8.108108	11
	Iran	32.0000	53.0000	2020-03-25	27017	2077	9625	15315	7.687752	12
	Spain	40.0000	-4.0000	2020-03-25	49515	3647	5367	40501	7.365445	13
	Indonesia	-0.7893	113.9213	2020-03-25	790	58	31	701	7.341772	14
	Algeria	28.0339	1.6596	2020-03-25	302	21	65	216	6.953642	15
	Philippines	13.0000	122.0000	2020-03-25	636	38	26	572	5.974843	16
	Netherlands	52.1326	5.2913	2020-03-25	6412	356	3	6053	5.552090	17

Como agrupar datos cuantitativos: MOSAIC PLOT. Cuando necesitamos cortar valores continuos se suele utilizar la funcion CUT

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,] # Crear una sequencia 
mosaicplot(t(tt), shade = TRUE)

Análisis de datos temporal

DATOS A NIVEL MUNDIAL

datos_por_fecha <- aggregate(
  cbind(Casos_Confirmados, Casos_Muertos, Casos_Recuperados, Casos_Enfermos) ~ 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)

##        Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos
## 1 2020-01-22               555            17                28            510
## 2 2020-01-23               654            18                30            606
## 3 2020-01-24               941            26                35            880
## 4 2020-01-25              1434            42                38           1354
## 5 2020-01-26              2118            56                51           2011
## 6 2020-01-27              2927            82                58           2787

tail(datos_por_fecha)

##         Fecha Casos_Confirmados Casos_Muertos Casos_Recuperados Casos_Enfermos
## 94 2020-04-24           2811598        196718            769361        1845519
## 95 2020-04-25           2897619        202868            796131        1898620
## 96 2020-04-26           2972358        206568            842917        1922873
## 97 2020-04-27           3041759        211167            869418        1961174
## 98 2020-04-28           3116393        217153            902901        1996339
## 99 2020-04-29           3193881        227638            945671        2020572

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") # Lo que hacemos con LOG = Y es pasar a escala logarítimica el eje Y para ver que hay diferentes velocidades de contagiados a partir de determinadas fechas.
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) # Aqui le indico la leyenda. Pch = Point character y Line wide para ver mas gordito o menos la linea

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.1, # Si se quita este ultimo son solo lineas en vez de areas
            drawGrid = FALSE) %>% # Esto es para quitar la cuadricula
  dyRangeSelector() %>% # Estos es para hacer zoom o no abajo
  dyCrosshair(direction = "vertical") %>% # Linea orientativa
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE) %>% # Estos son para el tamaño del circulo y que se queden los puntos aunque no tenga el raton encima
  dyRoller(rollPeriod = 1)

MAPA INTERACTIVO CON DATOS SOBRE ESPAÑA

datos_spain <- datos %>%
  filter(Pais_Region == "Spain") %>% # Con el filter seleccionamos filas
  select(Fecha, starts_with("Casos_")) # Sirve para seleccionar algunas variables del dataset (por columnas que me interesen)


plot(x = datos_spain$Fecha, y = datos_spain$Casos_Confirmados,
     main = "Casos confirmados en España", type = "s", col = "blue", lwd = 2) # la "s" es de STEP, sube cada dia un escalón

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.1, # Si se quita este ultimo son solo lineas en vez de areas
            drawGrid = FALSE) %>% # Esto es para quitar la cuadricula
  dyRangeSelector() %>% # Estos es para hacer zoom o no abajo
  dyCrosshair(direction = "vertical") %>% # Linea orientativa
  dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2,
              hideOnMouseOut = FALSE) %>% # Estos son para el tamaño del circulo y que se queden los puntos aunque no tenga el raton encima
  dyRoller(rollPeriod = 1)

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") # Normalmente se hace de una matriz, pero claro necesitamos transponer los datos con la t y decirle todas las filas pero solo de la columna 2 a la 5

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

EL PROBLEMA DE ESTE ULTIMO BARPLOT APILADO ES QUE NO SABES SI REALMENTE EL NUMERO DE ENFERMOS HA CRECIDO O ES QUE LA BARRA DE RECUPERADOS EMPUJA HACIA ARRIBA LA DE ENFERMOS

Gracias a la funcion LEAD y LAG del paquete DPLYR, aplicados a un vector de datos como pueden ser los casos confirmados, nos hace retroceder en el tiempo

lag(datos_spain$Casos_Confirmados, n = 1) # Nos ha desplazado un dia hacia atras

##  [1]     NA      0      0      0      0      0      0      0      0      0
## [11]      0      1      1      1      1      1      1      1      1      2
## [21]      2      2      2      2      2      2      2      2      2      2
## [31]      2      2      2      2      2      6     13     15     32     45
## [41]     84    120    165    222    259    400    500    673   1073   1695
## [51]   2277   2277   5232   6391   7798   9942  11748  13910  17963  20410
## [61]  25374  28768  35136  39885  49515  57786  65719  73235  80110  87956
## [71]  95923 104118 112065 119199 126168 131646 136675 141942 148220 153222
## [81] 158273 163027 166831 170099 172541 177644 184948 190839 191726 198674
## [91] 200210 204178 208389 213024 219764 223759 226629 229422 232128

lead(datos_spain$Casos_Confirmados, n = 1) # Nos adelanta el dia hacia adelante

##  [1]      0      0      0      0      0      0      0      0      0      1
## [11]      1      1      1      1      1      1      1      2      2      2
## [21]      2      2      2      2      2      2      2      2      2      2
## [31]      2      2      2      6     13     15     32     45     84    120
## [41]    165    222    259    400    500    673   1073   1695   2277   2277
## [51]   5232   6391   7798   9942  11748  13910  17963  20410  25374  28768
## [61]  35136  39885  49515  57786  65719  73235  80110  87956  95923 104118
## [71] 112065 119199 126168 131646 136675 141942 148220 153222 158273 163027
## [81] 166831 170099 172541 177644 184948 190839 191726 198674 200210 204178
## [91] 208389 213024 219764 223759 226629 229422 232128 236899     NA

datos_spain_1 <- 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),
         Nuevos_Casos_Enfermos = Casos_Enfermos - lag(Casos_Enfermos, n = 1))
  
datos_spain_1 %>%
  kable() %>%
  kable_styling()

Fecha	Casos_Confirmados	Casos_Muertos	Casos_Recuperados	Casos_Enfermos	Nuevos_Casos_Confirmados	Nuevos_Casos_Muertos	Nuevos_Casos_Recuperados	Nuevos_Casos_Enfermos
2020-01-22	0	0	0	0	NA	NA	NA	NA
2020-01-23	0	0	0	0	0	0	0	0
2020-01-24	0	0	0	0	0	0	0	0
2020-01-25	0	0	0	0	0	0	0	0
2020-01-26	0	0	0	0	0	0	0	0
2020-01-27	0	0	0	0	0	0	0	0
2020-01-28	0	0	0	0	0	0	0	0
2020-01-29	0	0	0	0	0	0	0	0
2020-01-30	0	0	0	0	0	0	0	0
2020-01-31	0	0	0	0	0	0	0	0
2020-02-01	1	0	0	1	1	0	0	1
2020-02-02	1	0	0	1	0	0	0	0
2020-02-03	1	0	0	1	0	0	0	0
2020-02-04	1	0	0	1	0	0	0	0
2020-02-05	1	0	0	1	0	0	0	0
2020-02-06	1	0	0	1	0	0	0	0
2020-02-07	1	0	0	1	0	0	0	0
2020-02-08	1	0	0	1	0	0	0	0
2020-02-09	2	0	0	2	1	0	0	1
2020-02-10	2	0	0	2	0	0	0	0
2020-02-11	2	0	0	2	0	0	0	0
2020-02-12	2	0	0	2	0	0	0	0
2020-02-13	2	0	0	2	0	0	0	0
2020-02-14	2	0	0	2	0	0	0	0
2020-02-15	2	0	2	0	0	0	2	-2
2020-02-16	2	0	2	0	0	0	0	0
2020-02-17	2	0	2	0	0	0	0	0
2020-02-18	2	0	2	0	0	0	0	0
2020-02-19	2	0	2	0	0	0	0	0
2020-02-20	2	0	2	0	0	0	0	0
2020-02-21	2	0	2	0	0	0	0	0
2020-02-22	2	0	2	0	0	0	0	0
2020-02-23	2	0	2	0	0	0	0	0
2020-02-24	2	0	2	0	0	0	0	0
2020-02-25	6	0	2	4	4	0	0	4
2020-02-26	13	0	2	11	7	0	0	7
2020-02-27	15	0	2	13	2	0	0	2
2020-02-28	32	0	2	30	17	0	0	17
2020-02-29	45	0	2	43	13	0	0	13
2020-03-01	84	0	2	82	39	0	0	39
2020-03-02	120	0	2	118	36	0	0	36
2020-03-03	165	1	2	162	45	1	0	44
2020-03-04	222	2	2	218	57	1	0	56
2020-03-05	259	3	2	254	37	1	0	36
2020-03-06	400	5	2	393	141	2	0	139
2020-03-07	500	10	30	460	100	5	28	67
2020-03-08	673	17	30	626	173	7	0	166
2020-03-09	1073	28	32	1013	400	11	2	387
2020-03-10	1695	35	32	1628	622	7	0	615
2020-03-11	2277	54	183	2040	582	19	151	412
2020-03-12	2277	55	183	2039	0	1	0	-1
2020-03-13	5232	133	193	4906	2955	78	10	2867
2020-03-14	6391	195	517	5679	1159	62	324	773
2020-03-15	7798	289	517	6992	1407	94	0	1313
2020-03-16	9942	342	530	9070	2144	53	13	2078
2020-03-17	11748	533	1028	10187	1806	191	498	1117
2020-03-18	13910	623	1081	12206	2162	90	53	2019
2020-03-19	17963	830	1107	16026	4053	207	26	3820
2020-03-20	20410	1043	1588	17779	2447	213	481	1753
2020-03-21	25374	1375	2125	21874	4964	332	537	4095
2020-03-22	28768	1772	2575	24421	3394	397	450	2547
2020-03-23	35136	2311	2575	30250	6368	539	0	5829
2020-03-24	39885	2808	3794	33283	4749	497	1219	3033
2020-03-25	49515	3647	5367	40501	9630	839	1573	7218
2020-03-26	57786	4365	7015	46406	8271	718	1648	5905
2020-03-27	65719	5138	9357	51224	7933	773	2342	4818
2020-03-28	73235	5982	12285	54968	7516	844	2928	3744
2020-03-29	80110	6803	14709	58598	6875	821	2424	3630
2020-03-30	87956	7716	16780	63460	7846	913	2071	4862
2020-03-31	95923	8464	19259	68200	7967	748	2479	4740
2020-04-01	104118	9387	22647	72084	8195	923	3388	3884
2020-04-02	112065	10348	26743	74974	7947	961	4096	2890
2020-04-03	119199	11198	30513	77488	7134	850	3770	2514
2020-04-04	126168	11947	34219	80002	6969	749	3706	2514
2020-04-05	131646	12641	38080	80925	5478	694	3861	923
2020-04-06	136675	13341	40437	82897	5029	700	2357	1972
2020-04-07	141942	14045	43208	84689	5267	704	2771	1792
2020-04-08	148220	14792	48021	85407	6278	747	4813	718
2020-04-09	153222	15447	52165	85610	5002	655	4144	203
2020-04-10	158273	16081	55668	86524	5051	634	3503	914
2020-04-11	163027	16606	59109	87312	4754	525	3441	788
2020-04-12	166831	17209	62391	87231	3804	603	3282	-81
2020-04-13	170099	17756	64727	87616	3268	547	2336	385
2020-04-14	172541	18056	67504	86981	2442	300	2777	-635
2020-04-15	177644	18708	70853	88083	5103	652	3349	1102
2020-04-16	184948	19315	74797	90836	7304	607	3944	2753
2020-04-17	190839	20002	74797	96040	5891	687	0	5204
2020-04-18	191726	20043	74797	96886	887	41	0	846
2020-04-19	198674	20453	77357	100864	6948	410	2560	3978
2020-04-20	200210	20852	80587	98771	1536	399	3230	-2093
2020-04-21	204178	21282	82514	100382	3968	430	1927	1611
2020-04-22	208389	21717	85915	100757	4211	435	3401	375
2020-04-23	213024	22157	89250	101617	4635	440	3335	860
2020-04-24	219764	22524	92355	104885	6740	367	3105	3268
2020-04-25	223759	22902	95708	105149	3995	378	3353	264
2020-04-26	226629	23190	117727	85712	2870	288	22019	-19437
2020-04-27	229422	23521	120832	85069	2793	331	3105	-643
2020-04-28	232128	23822	123903	84403	2706	301	3071	-666
2020-04-29	236899	24275	132929	79695	4771	453	9026	-4708

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

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

Assignment 1: Tasa de Variación Media

Se define la tasa de variación media como el incremento de nuevos casos de contagio con respecto al día anterior. Con esta definición en mente vamos a hacer los siguientes ejercicios

Questions for this assignment Teniendo en cuenta que la TVM del número de contagios para un día t se podría definir como:

TVM(t) = (Casos en el día t - Casos en el día t-1)/Casos en el día t

utiliza las funciones lead o lag para definir la TVM del número de nuevos casos de contagio en España (o en tu propio país)

Representa utilizando el gráfico adecuado la TVM y discute cómo ha ido evolucionando a lo largo del tiempo.

¿Se observa algún cambio en la tendencia que parezca indicar que el contagio se está frenando y que las medidas de contención son efectivas? Justifica la respuesta.

CALCULO TASA DE VARIACIÓN MEDIA

#TVM <- datos_spain_1$Casos_Confirmados - lag(datos_spain_1$Casos_Confirmados, n = 1) / datos_spain_1$Casos_Confirmados
options(scipen=999)
TVM <- datos_spain_1 %>%
  mutate(TVM = (Casos_Confirmados - lag(Casos_Confirmados, n = 1)) / Casos_Confirmados) 

plot(TVM ~ Fecha, data = TVM,
     type = "l", col = "blue",
     xlab = "Fecha", ylab = "% Casos Confirmados",
     main = "Tasa de variación media en España")

ANÁLISIS POR COHORTES - SEGMENTAR DEPENDIENDO DE DIFERENTES VARIABLES

TENEMOS QUE SACAR PRIMERO LA FECHA INICIAL DE CADA PAIS CUANDO SE ORIGINÓ LA PANDEMIA

primer_contagio <- datos %>%
  group_by(Pais_Region) %>%
  filter(Casos_Confirmados > 0) %>%
  summarise(primer_contagio = min(Fecha)-1) # Aqui saco por paises el primer dia de caso contagiado pero le resto -1 para obtener el dia 0.

primer_contagio %>%
  kable() %>%
  kable_styling()# El ultimo de normalidad antes de que se produjera el primer contagio

Pais_Region	primer_contagio
Afghanistan	2020-02-23
Albania	2020-03-08
Algeria	2020-02-24
Andorra	2020-03-01
Angola	2020-03-19
Antigua and Barbuda	2020-03-12
Argentina	2020-03-02
Armenia	2020-02-29
Australia	2020-01-25
Austria	2020-02-24
Azerbaijan	2020-02-29
Bahamas	2020-03-15
Bahrain	2020-02-23
Bangladesh	2020-03-07
Barbados	2020-03-16
Belarus	2020-02-27
Belgium	2020-02-03
Belize	2020-03-22
Benin	2020-03-15
Bhutan	2020-03-05
Bolivia	2020-03-10
Bosnia and Herzegovina	2020-03-04
Botswana	2020-03-29
Brazil	2020-02-25
Brunei	2020-03-08
Bulgaria	2020-03-07
Burkina Faso	2020-03-09
Burma	2020-03-26
Burundi	2020-03-30
Cabo Verde	2020-03-19
Cambodia	2020-01-26
Cameroon	2020-03-05
Canada	2020-01-25
Central African Republic	2020-03-14
Chad	2020-03-18
Chile	2020-03-02
China	2020-01-21
Colombia	2020-03-05
Congo (Brazzaville)	2020-03-14
Congo (Kinshasa)	2020-03-10
Costa Rica	2020-03-05
Cote d’Ivoire	2020-03-10
Croatia	2020-02-24
Cuba	2020-03-11
Cyprus	2020-03-08
Czechia	2020-02-29
Denmark	2020-02-26
Diamond Princess	2020-02-06
Djibouti	2020-03-17
Dominica	2020-03-21
Dominican Republic	2020-02-29
Ecuador	2020-02-29
Egypt	2020-02-13
El Salvador	2020-03-18
Equatorial Guinea	2020-03-14
Eritrea	2020-03-20
Estonia	2020-02-26
Eswatini	2020-03-13
Ethiopia	2020-03-12
Fiji	2020-03-18
Finland	2020-01-28
France	2020-01-23
Gabon	2020-03-13
Gambia	2020-03-16
Georgia	2020-02-25
Germany	2020-01-26
Ghana	2020-03-13
Greece	2020-02-25
Grenada	2020-03-21
Guatemala	2020-03-13
Guinea	2020-03-12
Guinea-Bissau	2020-03-24
Guyana	2020-03-11
Haiti	2020-03-19
Holy See	2020-03-05
Honduras	2020-03-10
Hungary	2020-03-03
Iceland	2020-02-27
India	2020-01-29
Indonesia	2020-03-01
Iran	2020-02-18
Iraq	2020-02-23
Ireland	2020-02-28
Israel	2020-02-20
Italy	2020-01-30
Jamaica	2020-03-10
Japan	2020-01-21
Jordan	2020-03-02
Kazakhstan	2020-03-12
Kenya	2020-03-12
Kosovo	2020-03-25
Kuwait	2020-02-23
Kyrgyzstan	2020-03-17
Laos	2020-03-23
Latvia	2020-03-01
Lebanon	2020-02-20
Liberia	2020-03-15
Libya	2020-03-23
Liechtenstein	2020-03-03
Lithuania	2020-02-27
Luxembourg	2020-02-28
Madagascar	2020-03-19
Malawi	2020-04-01
Malaysia	2020-01-24
Maldives	2020-03-07
Mali	2020-03-24
Malta	2020-03-06
Mauritania	2020-03-13
Mauritius	2020-03-17
Mexico	2020-02-27
Moldova	2020-03-07
Monaco	2020-02-28
Mongolia	2020-03-09
Montenegro	2020-03-16
Morocco	2020-03-01
Mozambique	2020-03-21
MS Zaandam	2020-03-27
Namibia	2020-03-13
Nepal	2020-01-24
Netherlands	2020-02-26
New Zealand	2020-02-27
Nicaragua	2020-03-18
Niger	2020-03-19
Nigeria	2020-02-27
North Macedonia	2020-02-25
Norway	2020-02-25
Oman	2020-02-23
Pakistan	2020-02-25
Panama	2020-03-09
Papua New Guinea	2020-03-19
Paraguay	2020-03-07
Peru	2020-03-05
Philippines	2020-01-29
Poland	2020-03-03
Portugal	2020-03-01
Qatar	2020-02-28
Romania	2020-02-25
Russia	2020-01-30
Rwanda	2020-03-13
Saint Kitts and Nevis	2020-03-24
Saint Lucia	2020-03-13
Saint Vincent and the Grenadines	2020-03-13
San Marino	2020-02-26
Sao Tome and Principe	2020-04-05
Saudi Arabia	2020-03-01
Senegal	2020-03-01
Serbia	2020-03-05
Seychelles	2020-03-13
Sierra Leone	2020-03-30
Singapore	2020-01-22
Slovakia	2020-03-05
Slovenia	2020-03-04
Somalia	2020-03-15
South Africa	2020-03-04
South Korea	2020-01-21
South Sudan	2020-04-04
Spain	2020-01-31
Sri Lanka	2020-01-26
Sudan	2020-03-12
Suriname	2020-03-13
Sweden	2020-01-30
Switzerland	2020-02-24
Syria	2020-03-21
Taiwan*	2020-01-21
Tanzania	2020-03-15
Thailand	2020-01-21
Timor-Leste	2020-03-21
Togo	2020-03-05
Trinidad and Tobago	2020-03-13
Tunisia	2020-03-03
Turkey	2020-03-10
Uganda	2020-03-20
Ukraine	2020-03-02
United Arab Emirates	2020-01-28
United Kingdom	2020-01-30
Uruguay	2020-03-13
Uzbekistan	2020-03-14
Venezuela	2020-03-13
Vietnam	2020-01-22
West Bank and Gaza	2020-03-04
Western Sahara	2020-04-04
Yemen	2020-04-09
Zambia	2020-03-17
Zimbabwe	2020-03-19
United States	2020-01-21

data_first =  datos %>%
  inner_join(primer_contagio, by = "Pais_Region") %>%
  mutate(dias_desde_PrimerC = as.numeric(Fecha - primer_contagio)) %>%
  filter(dias_desde_PrimerC >=0) %>%
  group_by(dias_desde_PrimerC, 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")) %>%
  ggplot(aes(x=dias_desde_PrimerC, y=Casos_Confirmados)) +
  geom_line(aes(col = Pais_Region)) +
  xlab("Dias desde el primer contagio") +
  ylab("Número de personas contagiadas") +
  ggtitle("Análisi por cohortes") -> g#+
  #theme(legend.position = "none")
ggplotly(g)

Modelos de regresión simple

\(x\): Variable Independiente: Número de dias desde el origen de la pandemia
\(Y\): Variable Dependiente: Número de casos confirmados

\[y = f(x)\]

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

Regresion Lineal

mod1 <- lm(formula = Casos_Confirmados ~Dias, data = datos_spain)
summary(mod1)

## 
## Call:
## lm(formula = Casos_Confirmados ~ Dias, data = datos_spain)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -62873 -35811   4958  35507  61926 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) -61925.6     7929.1   -7.81     0.00000000000679 ***
## Dias          2476.9      139.8   17.72 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39750 on 97 degrees of freedom
## Multiple R-squared:  0.764,  Adjusted R-squared:  0.7616 
## F-statistic:   314 on 1 and 97 DF,  p-value: < 0.00000000000000022

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") # Al verse dos patrones ya sabemos que es un modelo pesimo

residuos = mod1$residuals
residuos

##             1             2             3             4             5 
##  61925.589091  59448.658701  56971.728312  54494.797922  52017.867532 
##             6             7             8             9            10 
##  49540.937143  47064.006753  44587.076364  42110.145974  39633.215584 
##            11            12            13            14            15 
##  37157.285195  34680.354805  32203.424416  29726.494026  27249.563636 
##            16            17            18            19            20 
##  24772.633247  22295.702857  19818.772468  17342.842078  14865.911688 
##            21            22            23            24            25 
##  12388.981299   9912.050909   7435.120519   4958.190130   2481.259740 
##            26            27            28            29            30 
##      4.329351  -2472.601039  -4949.531429  -7426.461818  -9903.392208 
##            31            32            33            34            35 
## -12380.322597 -14857.252987 -17334.183377 -19811.113766 -22284.044156 
##            36            37            38            39            40 
## -24753.974545 -27228.904935 -29688.835325 -32152.765714 -34590.696104 
##            41            42            43            44            45 
## -37031.626494 -39463.556883 -41883.487273 -44323.417662 -46659.348052 
##            46            47            48            49            50 
## -49036.278442 -51340.208831 -53417.139221 -55272.069610 -57167.000000 
##            51            52            53            54            55 
## -59643.930390 -59165.860779 -60483.791169 -61553.721558 -61886.651948 
##            56            57            58            59            60 
## -62557.582338 -62872.512727 -61296.443117 -61326.373506 -58839.303896 
##            61            62            63            64            65 
## -57922.234286 -54031.164675 -51759.095065 -44606.025455 -38811.955844 
##            66            67            68            69            70 
## -33355.886234 -28316.816623 -23918.747013 -18549.677403 -13059.607792 
##            71            72            73            74            75 
##  -7341.538182  -1871.468571   2785.601039   7277.670649  10278.740260 
##            76            77            78            79            80 
##  12830.809870  15620.879481  19421.949091  21947.018701  24521.088312 
##            81            82            83            84            85 
##  26798.157922  28125.227532  28916.297143  28881.366753  31507.436364 
##            86            87            88            89            90 
##  36334.505974  39748.575584  38158.645195  42629.714805  41688.784416 
##            91            92            93            94            95 
##  43179.854026  44913.923636  47071.993247  51335.062857  52853.132468 
##            96            97            98            99 
##  53246.202078  53562.271688  53791.341299  56085.410909

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

## [1] 57 56

No es una buena prediccion ya que los residuos no siguen una normal, aunque la significatividad y el R2 era aceptable.

REGRESIÓN EXPONENCIAL: Buscamos estimar con una regresion lineal el logaritmo de la variable dependiente en funcion de la independiente.

\[LOG(Y) = ax+b, a,b \in \mathbb R \]

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

## 
## Call:
## lm(formula = log(Casos_Confirmados) ~ Dias, data = datos_spain[datos_spain$Casos_Confirmados > 
##     0, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.77186 -1.14690  0.08563  1.29283  2.03263 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) -2.501856   0.359002  -6.969       0.000000000583 ***
## Dias         0.180093   0.006003  29.998 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.455 on 87 degrees of freedom
## Multiple R-squared:  0.9118, Adjusted R-squared:  0.9108 
## F-statistic: 899.9 on 1 and 87 DF,  p-value: < 0.00000000000000022

EL MODELO ANTERIOR DE REGRESION SIMPLE ME APORTABA UN R CUADRADO DEL 0.6, EN ESTE TENEMOS UN MULTIPLE R-SQUARED DE: 0.91 LO QUE SIGNIFICA QUE PARCE MUCHO MEJOR.

Representamos la predicción

plot(datos_spain$Dias, datos_spain$Casos_Confirmados) # X = Dias Y = Confirm
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))

## 99 34 
## 89 24

Modelo potencial

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

## 
## Call:
## lm(formula = log(Casos_Confirmados) ~ log(Dias), data = datos_spain[datos_spain$Casos_Confirmados > 
##     0, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8814 -0.7533  0.2524  1.0041  4.6383 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept) -22.4061     1.0891  -20.57 <0.0000000000000002 ***
## log(Dias)     7.7165     0.2803   27.53 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.573 on 87 degrees of freedom
## Multiple R-squared:  0.897,  Adjusted R-squared:  0.8958 
## F-statistic: 757.8 on 1 and 87 DF,  p-value: < 0.00000000000000022

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))

## 11 12 
##  1  2

Comparar los modelos visualmente

start_date = ymd("2020-01-22")
end_date = ymd("2020-05-5")

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)

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, ])

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_Pon = pred3,
                                        Mod_Mix = pred4),
                                        order.by = dates) #Aglutinaos por columnas las cuatro pedicciones con los datos reales.

dygraph(datos_por_fecha_ts)

Estudio de la homogenidad de la infeccion

library(wbstats)
pop_data <- wb(indicator = "SP.POP.TOTL", startdate = 2018, enddate = 2019)

covid19=read.csv("covid_19_clean_complete.csv")

pop_data[pop_data$country == "Spain",]

##     iso3c date    value indicatorID         indicator iso2c country
## 453   ESP 2019 47076781 SP.POP.TOTL Population, total    ES   Spain
## 454   ESP 2018 46797754 SP.POP.TOTL Population, total    ES   Spain

paises=unique(covid19$Country.Region)
covid19.limpio = c()
for (i in 1:length(paises)){
  if(length(which(paises[i] %in% pop_data$country))>0){
    covid19.limpio = rbind(covid19.limpio, covid19[covid19$Country.Region==paises[1],])
  }
}

covid19.limpio$Date = as.Date(as.character(covid19.limpio$Date), "%m/%d/%Y")

infectados.totales.por.dia = aggregate(covid19.limpio$Confirmed ~covid19.limpio$Date, FUN = sum)

fallecidos.totales.por.dia= aggregate(covid19.limpio$Deaths ~covid19.limpio$Date, FUN = sum)

recuperados.totales.por.dia= aggregate(covid19.limpio$Recovered ~covid19.limpio$Date, FUN = sum)

tabla.totales= data.frame(infectados.totales.por.dia[,1], infectados.totales.por.dia[,2],fallecidos.totales.por.dia[,2],recuperados.totales.por.dia[,2])

names(tabla.totales) = c("Fecha", "Infectados", "Fallecidos", "Recuperados")

x=tabla.totales[,1]
ggplot(tabla.totales, aes(x)) +
  geom_line(aes(y=tabla.totales$Infectados, colour="Infectados")) +
  geom_line(aes(y=tabla.totales$Fallecidos, colour="Fallecidos")) +
  geom_line(aes(y=tabla.totales$Recuperados, colour="Recuperados")) +
  xlab("Fecha") + ylab("Frecuencias") + 
  scale_color_manual(values = c("red", "blue", "green")) -> g

ggplotly(g)

## Warning: Use of `tabla.totales$Infectados` is discouraged. Use `Infectados`
## instead.

## Warning: Use of `tabla.totales$Fallecidos` is discouraged. Use `Fallecidos`
## instead.

## Warning: Use of `tabla.totales$Recuperados` is discouraged. Use `Recuperados`
## instead.

Numero de casos confirmados por pais

fecha = "0020-03-30"

confirmados.por.pais = aggregate(covid19.limpio$Confirmed[covid19.limpio$Date==fecha] ~ covid19.limpio$Country.Region[covid19.limpio$Date==fecha], FUN = sum)
names(confirmados.por.pais)=c("Pais","Confirmados")

tail(confirmados.por.pais)

##          Pais Confirmados
## 1 Afghanistan       26860

Modelado matmático SIR

covid19$Date=as.Date(as.character(covid19$Date), "%m/%d/%Y")
str(covid19)

## 'data.frame':    25938 obs. of  8 variables:
##  $ Province.State: Factor w/ 81 levels "","Alberta","Anguilla",..: 1 1 1 1 1 1 1 1 6 49 ...
##  $ Country.Region: Factor w/ 185 levels "Afghanistan",..: 1 2 3 4 5 6 7 8 9 9 ...
##  $ Lat           : num  33 41.2 28 42.5 -11.2 ...
##  $ Long          : num  65 20.17 1.66 1.52 17.87 ...
##  $ Date          : Date, format: "0020-01-22" "0020-01-22" ...
##  $ Confirmed     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Deaths        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Recovered     : int  0 0 0 0 0 0 0 0 0 0 ...

library(wbstats)

pop_data2 <- wb(indicator = "SP.POP.TOTL", startdate = 2018, enddate = 2019)

pop_data[pop_data$country=="Spain",]

##     iso3c date    value indicatorID         indicator iso2c country
## 453   ESP 2019 47076781 SP.POP.TOTL Population, total    ES   Spain
## 454   ESP 2018 46797754 SP.POP.TOTL Population, total    ES   Spain

Creamos una tabla solo para España

covid19.España = covid19[covid19$Country.Region=="Spain",]

Ahora calculamos los infectados y los recuperados por días:

infectados.por.dia = aggregate(covid19.España$Confirmed ~covid19.España$Date, FUN = sum)
fallecidos.por.dia = aggregate(covid19.España$Deaths ~covid19.España$Date, FUN = sum)
infectados.por.dia2 = infectados.por.dia[,2] + fallecidos.por.dia[,2] # Se consideran a los fallecidos como infectados
recuperados.por.dia = aggregate(covid19.España$Recovered ~covid19.España$Date, FUN = sum)

Ahora calculamos los habitantes totales

habitantes=pop_data$value[pop_data$country=="Spain"]
Susceptibles.por.dia = habitantes-infectados.por.dia2-recuperados.por.dia[,2]

## Warning in habitantes - infectados.por.dia2: longitud de objeto mayor no es
## múltiplo de la longitud de uno menor

Creamos la tabla tabla.españa con toda la info

tabla.España = data.frame(unique(covid19.España$Date), Susceptibles.por.dia,infectados.por.dia2,recuperados.por.dia[,2])
names(tabla.España) = c("Fecha","Susceptibles","Infectados","Recuperados")
head(tabla.España,10)

##         Fecha Susceptibles Infectados Recuperados
## 1  0020-01-22     47076781          0           0
## 2  0020-01-23     46797754          0           0
## 3  0020-01-24     47076781          0           0
## 4  0020-01-25     46797754          0           0
## 5  0020-01-26     47076781          0           0
## 6  0020-01-27     46797754          0           0
## 7  0020-01-28     47076781          0           0
## 8  0020-01-29     46797754          0           0
## 9  0020-01-30     47076781          0           0
## 10 0020-01-31     46797754          0           0

Para estimar el valor R0 de España, calculando la pendiente de la recta de regresión de la variable Y = N In S(t) como función de la variable X = R(t)

x=tabla.España$Recuperados
y=habitantes*log(tabla.España$Susceptibles)

## Warning in habitantes * log(tabla.España$Susceptibles): longitud de objeto mayor
## no es múltiplo de la longitud de uno menor

summary(lm( y ~ x))

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2690768 -2602175  2515603  2605359  2607213 
## 
## Coefficients:
##                  Estimate    Std. Error  t value            Pr(>|t|)    
## (Intercept) 829111952.514    310262.829 2672.289 <0.0000000000000002 ***
## x                  -2.661         7.501   -0.355               0.724    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2631000 on 97 degrees of freedom
## Multiple R-squared:  0.001295,   Adjusted R-squared:  -0.009001 
## F-statistic: 0.1258 on 1 and 97 DF,  p-value: 0.7236

X ESTIMATE (Coefficients) = -3.46 es el R0 y vemos que la regresion es bastante buena porque el R cuadrado es de 0.96 (Multiple R-Squared)

El valor R0 estimado será:

(estimacion.R0= -summary(lm(y~x))$coefficients[2])

## [1] 2.660555

Pera ver si la epidemia se expandirá o no, hemos de comparar el valor anterior con: ….

dia.ultimo=length(tabla.España[1,]) # Ultimo dia
exp(estimacion.R0*tabla.España$Recuperados[dia.ultimo]/habitantes)

## [1] 1 1

Desgraciadamente vemos que el R0 es mayor que 1 por lo que el virus seguira expandiendose

Serie temporal

Una serie temporal es una sucesión de observaciones ordenadas a lo largo del tiempo.

Normalmente se considera una serie temporal para estudiar la evolución de un fenómeno o una variable a lo largo del tiempo.

El objetivo de una serie es tratar de hallar patrones de comportamiento del fenómeno o de la variable a estudiar y de realizar predicciones futuras. Puede considerarse un proceso estocastico en el que se estudian la evolucion en el tiempo de una sucesión de variables aleatorias. Dichas variables normalmente estan correlacionadas o existe una dependencia lineal entre las mismas

En una serie temporal Y(t) se intenta identificar cuatro componentes.

Tendencia(t) como evolucionará la serie
Estacionalidad(t) Oscilaciones que se producen y se repiten
Cíclica(t) son las variaciones cíclicas que se producen a largo plazo
Aleatoria(t) nos da la aportacion aleatoria

Hay dos modelos básicos para modelizar una serie temporal

El modelo aditivo: En este modelo suponemos que la serie temporal Y(t) es la suma de las cuatro componentes.

Y = T+E+C+R

El modelo multiplicativo: Suponemos que la serie temporal Y(t) es el producto de las tres componentes no aleatorias mas la componente aleatoria

Y = TCE+R

En el modelo aditivo la componente estacional mantiene su amplitud mientras que en modelo multiplicativo, dicha componente tiene una amplitud que se va amplificando en el tiempo

Serie temporal ajustada

Serie original de la que hemos eliminado la componente estacional

Modelo Holt_Winters

Doble suavizado temporal. Suponemos que tenemos una Y sin estacionalidad. El doble suavizado consiste en calcular dos series:

S(t): la serie suavizada y
B(t): la estimacion de la tendencia de la serie

S(1) = Y(1), B(1) = Y(1) - Y (0)

Calculo de la serie temporal de infectados

infectados = ts(tabla.España$Infectados, frequency = 7, start = c(1.3)) # Frequency seria por semana y el start quiere decir que el primer dato es de la semana 1 el dia 3
infectados

## Time Series:
## Start = 1.3 
## End = 15.3 
## Frequency = 7 
##  [1]      0      0      0      0      0      0      0      0      0      0
## [11]      1      1      1      1      1      1      1      1      2      2
## [21]      2      2      2      2      2      2      2      2      2      2
## [31]      2      2      2      2      6     13     15     32     45     84
## [41]    120    166    224    262    405    510    690   1101   1730   2331
## [51]   2332   5365   6586   8087  10284  12281  14533  18793  21453  26749
## [61]  30540  37447  42693  53162  62151  70857  79217  86913  95672 104387
## [71] 113505 122413 130397 138115 144287 150016 155987 163012 168669 174354
## [81] 179633 184040 187855 190597 196352 204263 210841 211769 219127 221062
## [91] 225460 230106 235181 242288 246661 249819 252943 255950 261174

library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

autoplot(infectados)

Calculo de las componentes de la serie

Nuestra serie tempoal es estacional con variacion semanal. Calculemos a continuación sus componentes suponiendo que el modelo es aditivo.

components = decompose(infectados, type = "additive")
components$seasonal

## Time Series:
## Start = 1.3 
## End = 15.3 
## Frequency = 7 
##  [1] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
##  [8] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [15] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [22] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [29] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [36] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [43] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [50] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [57] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [64] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [71] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [78] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [85] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [92] -259.5426  136.8090  581.2926  311.6324   96.8059 -246.3998 -620.5976
## [99] -259.5426

Como comenzamos un miercoles, se observa que los jueves, viernes y sabado hay un incremento de infectados mientras que los martes por ejemplo hay un descenso

autoplot(components)

La tendencia es creciente y la aleatoriedad del modelo crece en las semanas 8,9 y 10

Ajuste de la serie temporal

infectados.ajustados = infectados-components$seasonal # Restamos las variaciones estacionales
autoplot(infectados.ajustados)

Predicciones

(prediccion.infectados = HoltWinters(infectados.ajustados, gamma = FALSE))

## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = infectados.ajustados, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.7934889
##  beta : 0.6958784
##  gamma: FALSE
## 
## Coefficients:
##       [,1]
## a 261133.1
## b   4125.8

Hemos usado un doble suavizado exponencial ya que el suavizado simple no funciona bien cuando hay una tendencia en la serie temporal como es nuestro caso.

El valor de alpha es alto lo que nos dice que las predicciones usan el nivel actual de la serie de infectados con un peso elevado, es decir, para predecir el valor, de un dia en particular, usa el valor de los infectados del dia actual con un peso de 0.79.

Como el valor de beta = 1 tenemos que la pendiente de la tendencia de la prediccion en un dia particular vale 1.

El valor cometido en la preidccion se puede ver con SSE

prediccion.infectados$SSE

## [1] 138061092

sqrt(prediccion.infectados$SSE)

## [1] 11749.94

Nos estamos equivocando en unos 11749 infectados.

plot(prediccion.infectados)

En el gráfico vemos que se ajusta bastante bien.

Si queremos hacer predicciones mas alla del ultimo dia tenemos que:

(prediccion.infectados.semama = forecast:::forecast.HoltWinters(prediccion.infectados,h=7)) # H = 7 es la prediccion de la siguiente semana.

##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 15.44286       265258.9 263725.9 266791.9 262914.4 267603.5
## 15.58571       269384.7 266814.6 271954.9 265454.0 273315.5
## 15.72857       273510.5 269628.5 277392.6 267573.4 279447.6
## 15.87143       277636.3 272234.7 283037.9 269375.3 285897.3
## 16.01429       281762.1 274665.7 288858.5 270909.1 292615.1
## 16.15714       285887.9 276940.9 294834.9 272204.7 299571.2
## 16.30000       290013.7 279073.9 300953.6 273282.7 306744.8

La primera columna son los dias que hacemos la prediccion, la segunda la predccion, la tercera y la cuarta son los extremos del intervalo de confianza del valor predicho al 80%. La quinta y la sexta son los extremos del intervalo de confianza de valor predicho al 95%

autoplot(prediccion.infectados.semama)

Comprobación del modelo con jung-box.

ggAcf(prediccion.infectados.semama$residuals[3:75], lag.max = 7)

Tenemos problemas porque la autocorrelacion cuando tenemos un lag de 2 dias de los residuos es significativa.

Box.test(prediccion.infectados.semama$residuals, lag = 7, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  prediccion.infectados.semama$residuals
## X-squared = 3.9914, df = 7, p-value = 0.7808

shapiro.test(prediccion.infectados.semama$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  prediccion.infectados.semama$residuals
## W = 0.89231, p-value = 0.0000008632

No sigue una normalidad porque nos da un p-valor muy pequeño. Por lo tanto el modelo no es adecuado en este caso.

Analisis estadístico COVID19

Víctor López

30/4/2020