Deaths Predicted for next 10 days in Italy

Kết quả của sử dụng Auto ARIMA (Hyndman & Khandakar, 2008) để dự báo số người chết vì Coronavirus trong 10 ngày tới ở Table 2. Chú ý rằng cho đến hiện tại (ngày 18 - 03) thì dữ liệu mới chỉ có đến ngày 17 - 03 (tức trễ 1 ngày). Dưới đây là R codes:

# Load some packages: 
library(tidyverse)
library(lubridate)

# Current time: 
TODAY <- today()

# Data links (source: https://github.com/CSSEGISandData/COVID-19): 
link1 <- "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Recovered.csv"
link2 <- "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv"
link3 <- "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Deaths.csv"

# Load and clean data: 

link2 %>% 
  read_csv() %>% 
  rename(Province = `Province/State`, Country = `Country/Region`) %>% 
  gather(Date, Confirmed, -Province, -Country, -Lat, -Long) %>% 
  mutate(Date = mdy(Date)) -> df_confirmed


link3 %>% 
  read_csv() %>% 
  rename(Province = `Province/State`, Country = `Country/Region`) %>% 
  gather(Date, Deaths, -Province, -Country, -Lat, -Long) %>% 
  mutate(Date = mdy(Date)) -> df_deaths


# Function for combine data sets by country selected: 

combineData <- function(countryName) {
  
  df_confirmed %>% 
    filter(Country == countryName) %>% 
    group_by(Country, Date) %>% 
    summarise(Confirmed = sum(Confirmed)) %>% 
    ungroup() -> Confirmed
  
  df_deaths %>% 
    filter(Country == countryName) %>% 
    group_by(Country, Date) %>% 
    summarise(Deaths = sum(Deaths)) %>% 
    ungroup() -> Deaths
  
  full_join(Confirmed, Deaths %>% select(-Country), by = "Date") -> df_total
  
  return(df_total)
}


# Use the function: 

lapply(c("China", "Italy"), combineData) -> data_countries
do.call("bind_rows", data_countries) -> df_countries

italy_deaths <- df_countries %>% filter(Country == "Italy")


library(fpp2)


train <- italy_deaths$Deaths[1:50] %>% ts()
test <- italy_deaths$Deaths[51:56] %>% ts()
h <- length(test)


my_arima <- auto.arima(train)

# Use the model for forecasting: 
predicted_arima <- forecast(my_arima, h = h)$mean %>% as.vector()


# Data Frame for comparing: 
df_arima <- data_frame(Date = italy_deaths$Date[51:56], 
                       Actual = test %>% as.vector(), 
                       Predicted = predicted_arima %>% round(0), 
                       Error = Predicted - Actual,
                       Error_Percent = round(Error / Actual, 2))

# Model Performance: 
df_arima %>% knitr::kable(caption = "Table 1: Deaths by Coronavirus Predicted vs. Actuals")

Table 1: Deaths by Coronavirus Predicted vs. Actuals
Date	Actual	Predicted	Error	Error_Percent
2020-03-12	827	1059	232	0.28
2020-03-13	1266	1304	38	0.03
2020-03-14	1441	1567	126	0.09
2020-03-15	1809	1836	27	0.01
2020-03-16	2158	2114	-44	-0.02
2020-03-17	2503	2394	-109	-0.04

train_full <- italy_deaths$Deaths%>% ts()
my_arima_full <- auto.arima(train_full)

# Use the model for forecasting for next 10 days from 2020-03-17: 
predicted_arima <- forecast(my_arima_full, h = 10)$mean %>% as.vector()


df_arima_next10 <- data_frame(Date = italy_deaths$Date %>% max() + 1:10, 
                              Predicted = predicted_arima %>% round(0))


df_arima_next10 %>% knitr::kable(caption = "Table 2: Deaths by Coronavirus Predicted for Next 10 Days")

Table 2: Deaths by Coronavirus Predicted for Next 10 Days
Date	Predicted
2020-03-18	2896
2020-03-19	3283
2020-03-20	3669
2020-03-21	4085
2020-03-22	4482
2020-03-23	4895
2020-03-24	5309
2020-03-25	5720
2020-03-26	6139
2020-03-27	6555

df_predicted <- df_arima_next10 %>% 
  rename(Deaths = Predicted) %>% 
  mutate(Type = "Predicted")


italy_deaths %>% 
  mutate(Type = "Actual") -> df_actual


library(hrbrthemes)
my_font <- "Ubuntu Condensed"

ggplot() + 
  geom_point(data = df_actual, aes(Date, Deaths, color = Type), size = 2) + 
  geom_point(data = df_predicted, aes(Date, Deaths, color = Type), size = 2) + 
  theme_modern_rc(base_family = my_font) + 
  scale_y_continuous(limits = c(0, 7000), breaks = seq(0, 7000, 1000)) + 
  scale_color_manual(values = c("cyan", "red")) + 
  labs(title = "Deaths by Coronavirus Predicted for next 10 days in Italy") + 
  theme(legend.title = element_blank()) + 
  theme(plot.margin = unit(c(1.2, 1.2, 1.2, 1.2), "cm")) + 
  theme(legend.position = "top") + 
  theme(axis.text.y = element_text(color = "white", size = 12)) + 
  theme(axis.text.x = element_text(color = "white", size = 12)) + 
  theme(axis.title.x = element_blank()) + 
  theme(axis.title.y = element_blank()) + 
  theme(legend.text = element_text(size = 11, face = "bold", color = "white", family = my_font)) -> f1 


plotly::ggplotly(f1)

References

Kourentzes, N., Barrow, D. K., & Crone, S. F. (2014). Neural network ensemble operators for time series forecasting. Expert Systems with Applications, 41(9), 4235-4244.
Hyndman, R. J., and Y. Khandakar (2008):“Automatic time series forecasting: The forecast package for R,”. Journal of Statistical Software, 26(3).

LS0tDQp0aXRsZTogJ0RlYXRocyBieSBDb3JvbmF2aXJ1cyBQcmVkaWN0ZWQgZm9yIG5leHQgMTAgZGF5cyBpbiBJdGFseScNCmF1dGhvcjogJ0F1dGhvcjogTmd1eWVuIENoaSBEdW5nJw0Kc3VidGl0bGU6ICJSIE1hY2hpbmUgTGVhcm5pbmcgU2VyaWVzIg0Kb3V0cHV0Og0KICBodG1sX2RvY3VtZW50OiANCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlDQogICAgY29kZV9mb2xkaW5nOiBoaWRlDQogICAgaGlnaGxpZ2h0OiB6ZW5idXJuDQogICAgIyBudW1iZXJfc2VjdGlvbnM6IHllcw0KICAgIHRoZW1lOiAiZmxhdGx5Ig0KICAgIHRvYzogVFJVRQ0KICAgIHRvY19mbG9hdDogVFJVRQ0KLS0tDQoNCmBgYHtyIHNldHVwLGluY2x1ZGU9RkFMU0V9DQprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUsIHdhcm5pbmcgPSBGQUxTRSwgbWVzc2FnZSA9IEZBTFNFLCBmaWcud2lkdGggPSAxMCwgZmlnLmhlaWdodCA9IDYpDQpgYGANCg0KDQojIERlYXRocyBQcmVkaWN0ZWQgZm9yIG5leHQgMTAgZGF5cyBpbiBJdGFseQ0KDQpL4bq/dCBxdeG6oyBj4bunYSBz4butIGThu6VuZyBBdXRvIEFSSU1BIChIeW5kbWFuICYgS2hhbmRha2FyLCAyMDA4KSDEkeG7gyBk4buxIGLDoW8gc+G7kSBuZ8aw4budaSBjaOG6v3QgdsOsIENvcm9uYXZpcnVzIHRyb25nIDEwIG5nw6B5IHThu5tpIOG7nyBUYWJsZSAyLiBDaMO6IMO9IHLhurFuZyBjaG8gxJHhur9uIGhp4buHbiB04bqhaSAobmfDoHkgMTggLSAwMykgdGjDrCBk4buvIGxp4buHdSBt4bubaSBjaOG7iSBjw7MgxJHhur9uIG5nw6B5IDE3IC0gMDMgKHThu6ljIHRy4buFIDEgbmfDoHkpLiBExrDhu5tpIMSRw6J5IGzDoCBSIGNvZGVzOiANCg0KDQpgYGB7cn0NCg0KIyBMb2FkIHNvbWUgcGFja2FnZXM6IA0KbGlicmFyeSh0aWR5dmVyc2UpDQpsaWJyYXJ5KGx1YnJpZGF0ZSkNCg0KIyBDdXJyZW50IHRpbWU6IA0KVE9EQVkgPC0gdG9kYXkoKQ0KDQojIERhdGEgbGlua3MgKHNvdXJjZTogaHR0cHM6Ly9naXRodWIuY29tL0NTU0VHSVNhbmREYXRhL0NPVklELTE5KTogDQpsaW5rMSA8LSAiaHR0cHM6Ly9yYXcuZ2l0aHVidXNlcmNvbnRlbnQuY29tL0NTU0VHSVNhbmREYXRhL0NPVklELTE5L21hc3Rlci9jc3NlX2NvdmlkXzE5X2RhdGEvY3NzZV9jb3ZpZF8xOV90aW1lX3Nlcmllcy90aW1lX3Nlcmllc18xOS1jb3ZpZC1SZWNvdmVyZWQuY3N2Ig0KbGluazIgPC0gImh0dHBzOi8vcmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbS9DU1NFR0lTYW5kRGF0YS9DT1ZJRC0xOS9tYXN0ZXIvY3NzZV9jb3ZpZF8xOV9kYXRhL2Nzc2VfY292aWRfMTlfdGltZV9zZXJpZXMvdGltZV9zZXJpZXNfMTktY292aWQtQ29uZmlybWVkLmNzdiINCmxpbmszIDwtICJodHRwczovL3Jhdy5naXRodWJ1c2VyY29udGVudC5jb20vQ1NTRUdJU2FuZERhdGEvQ09WSUQtMTkvbWFzdGVyL2Nzc2VfY292aWRfMTlfZGF0YS9jc3NlX2NvdmlkXzE5X3RpbWVfc2VyaWVzL3RpbWVfc2VyaWVzXzE5LWNvdmlkLURlYXRocy5jc3YiDQoNCiMgTG9hZCBhbmQgY2xlYW4gZGF0YTogDQoNCmxpbmsyICU+JSANCiAgcmVhZF9jc3YoKSAlPiUgDQogIHJlbmFtZShQcm92aW5jZSA9IGBQcm92aW5jZS9TdGF0ZWAsIENvdW50cnkgPSBgQ291bnRyeS9SZWdpb25gKSAlPiUgDQogIGdhdGhlcihEYXRlLCBDb25maXJtZWQsIC1Qcm92aW5jZSwgLUNvdW50cnksIC1MYXQsIC1Mb25nKSAlPiUgDQogIG11dGF0ZShEYXRlID0gbWR5KERhdGUpKSAtPiBkZl9jb25maXJtZWQNCg0KDQpsaW5rMyAlPiUgDQogIHJlYWRfY3N2KCkgJT4lIA0KICByZW5hbWUoUHJvdmluY2UgPSBgUHJvdmluY2UvU3RhdGVgLCBDb3VudHJ5ID0gYENvdW50cnkvUmVnaW9uYCkgJT4lIA0KICBnYXRoZXIoRGF0ZSwgRGVhdGhzLCAtUHJvdmluY2UsIC1Db3VudHJ5LCAtTGF0LCAtTG9uZykgJT4lIA0KICBtdXRhdGUoRGF0ZSA9IG1keShEYXRlKSkgLT4gZGZfZGVhdGhzDQoNCg0KIyBGdW5jdGlvbiBmb3IgY29tYmluZSBkYXRhIHNldHMgYnkgY291bnRyeSBzZWxlY3RlZDogDQoNCmNvbWJpbmVEYXRhIDwtIGZ1bmN0aW9uKGNvdW50cnlOYW1lKSB7DQogIA0KICBkZl9jb25maXJtZWQgJT4lIA0KICAgIGZpbHRlcihDb3VudHJ5ID09IGNvdW50cnlOYW1lKSAlPiUgDQogICAgZ3JvdXBfYnkoQ291bnRyeSwgRGF0ZSkgJT4lIA0KICAgIHN1bW1hcmlzZShDb25maXJtZWQgPSBzdW0oQ29uZmlybWVkKSkgJT4lIA0KICAgIHVuZ3JvdXAoKSAtPiBDb25maXJtZWQNCiAgDQogIGRmX2RlYXRocyAlPiUgDQogICAgZmlsdGVyKENvdW50cnkgPT0gY291bnRyeU5hbWUpICU+JSANCiAgICBncm91cF9ieShDb3VudHJ5LCBEYXRlKSAlPiUgDQogICAgc3VtbWFyaXNlKERlYXRocyA9IHN1bShEZWF0aHMpKSAlPiUgDQogICAgdW5ncm91cCgpIC0+IERlYXRocw0KICANCiAgZnVsbF9qb2luKENvbmZpcm1lZCwgRGVhdGhzICU+JSBzZWxlY3QoLUNvdW50cnkpLCBieSA9ICJEYXRlIikgLT4gZGZfdG90YWwNCiAgDQogIHJldHVybihkZl90b3RhbCkNCn0NCg0KDQojIFVzZSB0aGUgZnVuY3Rpb246IA0KDQpsYXBwbHkoYygiQ2hpbmEiLCAiSXRhbHkiKSwgY29tYmluZURhdGEpIC0+IGRhdGFfY291bnRyaWVzDQpkby5jYWxsKCJiaW5kX3Jvd3MiLCBkYXRhX2NvdW50cmllcykgLT4gZGZfY291bnRyaWVzDQoNCml0YWx5X2RlYXRocyA8LSBkZl9jb3VudHJpZXMgJT4lIGZpbHRlcihDb3VudHJ5ID09ICJJdGFseSIpDQoNCg0KbGlicmFyeShmcHAyKQ0KDQoNCnRyYWluIDwtIGl0YWx5X2RlYXRocyREZWF0aHNbMTo1MF0gJT4lIHRzKCkNCnRlc3QgPC0gaXRhbHlfZGVhdGhzJERlYXRoc1s1MTo1Nl0gJT4lIHRzKCkNCmggPC0gbGVuZ3RoKHRlc3QpDQoNCg0KbXlfYXJpbWEgPC0gYXV0by5hcmltYSh0cmFpbikNCg0KIyBVc2UgdGhlIG1vZGVsIGZvciBmb3JlY2FzdGluZzogDQpwcmVkaWN0ZWRfYXJpbWEgPC0gZm9yZWNhc3QobXlfYXJpbWEsIGggPSBoKSRtZWFuICU+JSBhcy52ZWN0b3IoKQ0KDQoNCiMgRGF0YSBGcmFtZSBmb3IgY29tcGFyaW5nOiANCmRmX2FyaW1hIDwtIGRhdGFfZnJhbWUoRGF0ZSA9IGl0YWx5X2RlYXRocyREYXRlWzUxOjU2XSwgDQogICAgICAgICAgICAgICAgICAgICAgIEFjdHVhbCA9IHRlc3QgJT4lIGFzLnZlY3RvcigpLCANCiAgICAgICAgICAgICAgICAgICAgICAgUHJlZGljdGVkID0gcHJlZGljdGVkX2FyaW1hICU+JSByb3VuZCgwKSwgDQogICAgICAgICAgICAgICAgICAgICAgIEVycm9yID0gUHJlZGljdGVkIC0gQWN0dWFsLA0KICAgICAgICAgICAgICAgICAgICAgICBFcnJvcl9QZXJjZW50ID0gcm91bmQoRXJyb3IgLyBBY3R1YWwsIDIpKQ0KDQojIE1vZGVsIFBlcmZvcm1hbmNlOiANCmRmX2FyaW1hICU+JSBrbml0cjo6a2FibGUoY2FwdGlvbiA9ICJUYWJsZSAxOiBEZWF0aHMgYnkgQ29yb25hdmlydXMgUHJlZGljdGVkIHZzLiBBY3R1YWxzIikNCg0KDQoNCnRyYWluX2Z1bGwgPC0gaXRhbHlfZGVhdGhzJERlYXRocyU+JSB0cygpDQpteV9hcmltYV9mdWxsIDwtIGF1dG8uYXJpbWEodHJhaW5fZnVsbCkNCg0KIyBVc2UgdGhlIG1vZGVsIGZvciBmb3JlY2FzdGluZyBmb3IgbmV4dCAxMCBkYXlzIGZyb20gMjAyMC0wMy0xNzogDQpwcmVkaWN0ZWRfYXJpbWEgPC0gZm9yZWNhc3QobXlfYXJpbWFfZnVsbCwgaCA9IDEwKSRtZWFuICU+JSBhcy52ZWN0b3IoKQ0KDQoNCmRmX2FyaW1hX25leHQxMCA8LSBkYXRhX2ZyYW1lKERhdGUgPSBpdGFseV9kZWF0aHMkRGF0ZSAlPiUgbWF4KCkgKyAxOjEwLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIFByZWRpY3RlZCA9IHByZWRpY3RlZF9hcmltYSAlPiUgcm91bmQoMCkpDQoNCg0KZGZfYXJpbWFfbmV4dDEwICU+JSBrbml0cjo6a2FibGUoY2FwdGlvbiA9ICJUYWJsZSAyOiBEZWF0aHMgYnkgQ29yb25hdmlydXMgUHJlZGljdGVkIGZvciBOZXh0IDEwIERheXMiKQ0KDQoNCmRmX3ByZWRpY3RlZCA8LSBkZl9hcmltYV9uZXh0MTAgJT4lIA0KICByZW5hbWUoRGVhdGhzID0gUHJlZGljdGVkKSAlPiUgDQogIG11dGF0ZShUeXBlID0gIlByZWRpY3RlZCIpDQoNCg0KaXRhbHlfZGVhdGhzICU+JSANCiAgbXV0YXRlKFR5cGUgPSAiQWN0dWFsIikgLT4gZGZfYWN0dWFsDQoNCg0KbGlicmFyeShocmJydGhlbWVzKQ0KbXlfZm9udCA8LSAiVWJ1bnR1IENvbmRlbnNlZCINCg0KZ2dwbG90KCkgKyANCiAgZ2VvbV9wb2ludChkYXRhID0gZGZfYWN0dWFsLCBhZXMoRGF0ZSwgRGVhdGhzLCBjb2xvciA9IFR5cGUpLCBzaXplID0gMikgKyANCiAgZ2VvbV9wb2ludChkYXRhID0gZGZfcHJlZGljdGVkLCBhZXMoRGF0ZSwgRGVhdGhzLCBjb2xvciA9IFR5cGUpLCBzaXplID0gMikgKyANCiAgdGhlbWVfbW9kZXJuX3JjKGJhc2VfZmFtaWx5ID0gbXlfZm9udCkgKyANCiAgc2NhbGVfeV9jb250aW51b3VzKGxpbWl0cyA9IGMoMCwgNzAwMCksIGJyZWFrcyA9IHNlcSgwLCA3MDAwLCAxMDAwKSkgKyANCiAgc2NhbGVfY29sb3JfbWFudWFsKHZhbHVlcyA9IGMoImN5YW4iLCAicmVkIikpICsgDQogIGxhYnModGl0bGUgPSAiRGVhdGhzIGJ5IENvcm9uYXZpcnVzIFByZWRpY3RlZCBmb3IgbmV4dCAxMCBkYXlzIGluIEl0YWx5IikgKyANCiAgdGhlbWUobGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKSArIA0KICB0aGVtZShwbG90Lm1hcmdpbiA9IHVuaXQoYygxLjIsIDEuMiwgMS4yLCAxLjIpLCAiY20iKSkgKyANCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gInRvcCIpICsgDQogIHRoZW1lKGF4aXMudGV4dC55ID0gZWxlbWVudF90ZXh0KGNvbG9yID0gIndoaXRlIiwgc2l6ZSA9IDEyKSkgKyANCiAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoY29sb3IgPSAid2hpdGUiLCBzaXplID0gMTIpKSArIA0KICB0aGVtZShheGlzLnRpdGxlLnggPSBlbGVtZW50X2JsYW5rKCkpICsgDQogIHRoZW1lKGF4aXMudGl0bGUueSA9IGVsZW1lbnRfYmxhbmsoKSkgKyANCiAgdGhlbWUobGVnZW5kLnRleHQgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDExLCBmYWNlID0gImJvbGQiLCBjb2xvciA9ICJ3aGl0ZSIsIGZhbWlseSA9IG15X2ZvbnQpKSAtPiBmMSANCg0KDQpwbG90bHk6OmdncGxvdGx5KGYxKSANCmBgYA0KDQojIFJlZmVyZW5jZXMNCjEuIEtvdXJlbnR6ZXMsIE4uLCBCYXJyb3csIEQuIEsuLCAmIENyb25lLCBTLiBGLiAoMjAxNCkuIE5ldXJhbCBuZXR3b3JrIGVuc2VtYmxlIG9wZXJhdG9ycyBmb3IgdGltZSBzZXJpZXMgZm9yZWNhc3RpbmcuIEV4cGVydCBTeXN0ZW1zIHdpdGggQXBwbGljYXRpb25zLCA0MSg5KSwgNDIzNS00MjQ0Lg0KMi4gSHluZG1hbiwgUi4gSi4sIGFuZCBZLiBLaGFuZGFrYXIgKDIwMDgpOuKAnEF1dG9tYXRpYyB0aW1lIHNlcmllcyBmb3JlY2FzdGluZzogVGhlIGZvcmVjYXN0IHBhY2thZ2UgZm9yIFIs4oCdLiBKb3VybmFsIG9mIFN0YXRpc3RpY2FsIFNvZnR3YXJlLCAyNigzKS4NCg==

Deaths by Coronavirus Predicted for next 10 days in Italy

R Machine Learning Series

Author: Nguyen Chi Dung

Deaths Predicted for next 10 days in Italy

References