covid19.Rmd

---
title: "COVID-19 Epidemic: An Exploratory Data Analysis and Signal Extraction for Projecting Ultimate Counts"
author: "Daniel Suryakusuma"
date: "3/18/2019"
output: html_document
---


```{r setup, message= FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(httr)
library(rvest)
library(tidyverse)
library(lubridate)
library(pdftools)
library(pdfsearch)
library(kableExtra)
library(jsonlite)
library(plotly)



qkable <- function(x, height="360px") {
  x %>% kable(format = "html") %>% kable_styling(bootstrap_options = c("condensed", "responsive", "striped", "hover", "bordered"), font_size = 11, position = "center") %>% scroll_box(width="100%", height= height, fixed_thead =  list(enabled = TRUE, background = "lightgrey") ) 
}
```


In early stages of infection outbreaks, counts of confirmed cases suggest an underlying exponential growth trend. 

Let's take a closer look. 

<img align="center" width = "40%" src="/Users/dsury/covid-19/images/WHO-scrape.png">


## Web scraper to compile data from WHO (World Health Organization) daily reports. 

Noticing a pattern in the report URL structure, we might first try to exploit this to quickly get the files we need. 

```{r}
td <- Sys.Date()
firstreportdate <- as.Date("2020-01-21")
number.of.reports <- td - firstreportdate

# generate .pdf url 
# getReportUrl <- function(date = firstreportdate) {
#   return( paste0(urlstem, 
#                  gsub("-", x = date, replacement = ""), 
#                  "-sitrep-",
#                  td  - firstreportdate + 1,
#                  "-2019-ncov.pdf"
#                  ) )
# }

# proceeding would be pretty pointless, because it turns out the urls actually change structure mid-way

```


But of course, this is dependent upon the World Health Organization keeping up with this structure and not making errors or deviations. 

**It turns out that mid-way, the naming structure changes from 2/10 to 2/11, where `ncov` changes to `covid` in the URL. **

In general, we shouldn't rely on the web structure to have a neat coherent pattern behind its file naming system. 

Instead, let's look at the root page and pull all the urls from it. 


```{r}
# get the .pdf links from webpage
rooturl <- read_html("https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports") # all the .pdf links are here

situation.reports <- rooturl %>%
  html_nodes("a") %>% 
  html_attr("href") 
# %>%   filter(!grepl('.pdf'))

situation.reports <- situation.reports[grepl('.pdf',situation.reports)]
situation.reports <- situation.reports[grepl('coronaviruse/situation-reports/',situation.reports)] %>% unique()  # some links are duplicates
situation.reports <- paste0("https://www.who.int", situation.reports) # full link
```

We can check that the number of reports is correct and up to date. This difference should be 0, unless if there's a one-off error due to timezone / reporting lag on behalf of the WHO. 

```{r}
length(situation.reports) - number.of.reports
```


We can visualize our URLS and check that the dates match up correctly (especially the leap day February 29, thanks to `lubridate`). 

```{r}
sit.rep <- tibble(date = as.Date("2020-01-20") + 1:length(situation.reports), urls = rev(situation.reports)) 
sit.rep %>% qkable() 
```


Now we have the list of pdf URLS to parse through and clean up. Let's first download them offline so we can process them offline/faster. 


```{r}
# download offline for processing
if (! dir.exists("covid-19")) {
  dir.create("covid-19")
}
setwd("covid-19")

for (i in 1:nrow(sit.rep)) {
  filename <- paste0(getwd(),"/",sit.rep$date[i], ".pdf" )
  if (!file.exists(filename))
    download.file(url = sit.rep$urls[i], 
                  destfile = filename, 
                  mode = "wb")
}
```


```{r, warning = FALSE, message = FALSE}
setwd("covid-19")
pdf.paths <- paste0(sit.rep$date, ".pdf")
# getwd()

counts <- sit.rep %>% select(- "urls") %>% mutate(
              global = 0L,
              china = 0L,
              outside_china = 0L,
              # global_new = 0L,
              # china_new = 0L,
              global_deaths = 0L,
              china_deaths = 0L,
              outside_china_deaths = 0L
)


# REPORT TYPE A: Older, (actually the web scrape reveals this form was reused from the Zika virus)
# treat 1 ~ 6 MANUALLY; not worth the investment to automate this because this is an old format

editrow <- function(df, n = 0, global, china, outsidechina) {
  if ( (n <= nrow(df)) * (n >= 1) ) {
    df$global[n] <- global
    df$china[n] <- china
    df$outside_china[n] <- outsidechina
    }
  return(df)
}

# manually do first report's data
counts <- counts %>% editrow(n = 1, 
                             global = 282, 
                             china = 278, 
                             outsidechina = 3)
# "As of 20 January 2020, 282 confirmed cases of 2019-nCoV have been reported from four countries including China (278 cases), Thailand (2 cases), Japan (1 case) and the Republic of Korea (1 case);"


# report 2 : jan 22
counts <- counts %>% editrow(n = 2,
                             global = 314,
                             china = 309,
                             outsidechina = 4) # mismatch

# report 3: jan 23
counts <- counts %>% editrow(n = 3,
                             global = 581,
                             china = 571,
                             outsidechina = 10)

# report 4 : jan 24
counts <- counts %>% editrow(n = 4,
                             global = 846,
                             china = 830,
                             outsidechina = 11) # doesnt add up

# report 5 : jan 25
counts <- counts %>% editrow(n = 5,
                             global = 1320,
                             china = 1297,
                             outsidechina = 23)

# report 6 : jan 26
counts <- counts %>% editrow(n = 6,
                             global = 2014,
                             china = 1985,
                             outsidechina = 29)


# counts %>% qkable()


```



<!-- The data are in *serious* need of cleaning... From the above, we want to select out the "2798 confirmed", "2741 confirmed", "461 severe", "80 deaths", all while knowing what figures these represent exactly.  -->


## Data Cleaning

<!-- Let's do a custom tailored scrape as proof of concept.  -->

```{r, warning = FALSE, message = FALSE, eval = FALSE, include = FALSE}
tmp <- pdf_text(pdf = paste0("covid-19/", pdf.paths[7]))

result <- keyword_search(tmp[[1]], # first page only
                         keyword = c(" confirmed ", " suspected ", " severe ", " deaths "),
                         path = FALSE, 
                         surround_lines = 0
                         )
length(result)
head(result$line_text)
# result$line_text[1] %>% grep("^confirmed")
```

```{r, warning = FALSE, message = FALSE, eval = FALSE, include = FALSE}

bin <- c()
mess <- result$line_text[1]
# mess

# get global confirmed
bin <- c(bin, substr(
  x = mess,
  start = regexpr(pattern = " confirmed", text = mess) - 10,
  stop = regexpr(pattern = " confirmed", text = mess) + 9
) )

# get china confirmed
mess <- result$line_text[2]

bin <- c(bin, substr(
  x = mess,
  start = regexpr(pattern = " confirmed", text = mess) - 10,
  stop = regexpr(pattern = " confirmed", text = mess) + 9
) )

## get outsidechina confirmed
mess <- result$line_text[3]

bin <- c( bin, substr(
  x = mess,
  start = regexpr(pattern = " confirmed", text = mess) - 10,
  stop = regexpr(pattern = " confirmed", text = mess) + 9
) )


as.numeric(gsub(pattern =  ".*?([0-9]+).*", replacement =  "\\1", x = bin)) # alternatively just use tidyr's extract_numeric()
```


```{r, warning = FALSE, message = FALSE, eval = FALSE}
# head(result$line_text)
dict <- c(" confirmed", " deaths")
```

Let's write some functions to automate the data cleaning. 

```{r, warning = FALSE, message = FALSE}
getValues <- function(n, pattern) {
  # pattern : key from dict, such as " confirmed" or " deaths"
  # n : a particular situation report number (this function only support 7+; I manually filled out the first few)

  parsedsitrep <- pdf_text(pdf = paste0("covid-19/", pdf.paths[n]))

  result.n <- keyword_search(parsedsitrep[[1]], # care about first page only
                           keyword = pattern,
                           path = FALSE,
                           surround_lines = 0)

  txt <- result.n$line_text %>% str_split(pattern = "Globally", simplify = FALSE)
  txt <- txt %>% str_split(pattern = " China ", simplify = FALSE)
  txt <- txt %>% str_split(pattern = "•", simplify = FALSE, n = 2)


  bin <- c()
  
  for (i in 1:length(txt)) {
    for (j in 1:length(txt[[i]])) {
          mess <- txt[[i]][j]
          bin <- c(bin, substr(x = mess, 
                           start = regexpr(pattern = pattern, text = mess) - ifelse(pattern == " deaths ", 5, 8), # dont expect deaths to increase over 99,999 yet so minimize chance of picking up extraneous numbers
                           stop = regexpr(pattern = pattern, text = mess) + nchar(pattern) - 1))
    }
  }

  if ( pattern == " deaths " ) {
    count.vec <- bin %>%
      gsub(pattern = "\\.", replacement = "") %>%
      str_split(pattern = " ") %>%
      keep_numeric()
    count.vec <- count.vec[!is.na(count.vec)] %>% tail(3)
  } else if (pattern == " confirmed ") {
    count.vec <- bin %>%
      gsub(pattern = "\\.", replacement = "") %>%
      str_split(pattern = " ") %>%
      extract_numeric() # depreciated but better than parse_number()
    count.vec <- count.vec[!is.na(count.vec)] %>% tail(3)
  }
  
  return(count.vec)
}

keep_numeric <- function( obj ) {
  keep <- c()
  for (i in 1:length(obj)) {
    for (j in 1:length(obj[[i]])) {
      tmp <- try(parse_number(obj[[i]][j]))
      if (is.numeric(tmp)) {
        keep <- c(keep, tmp)
      }
    }
  }
  return(keep)
}


```

Here's a working example, where there are 1381 deaths in China by February 14 and 2 deaths outside of China.

```{r, warning = FALSE, message = FALSE}
getValues( n = 56, pattern = " confirmed " )
getValues( n = 56, pattern = " deaths " )
```


```{r, warning = FALSE, message = FALSE, eval = FALSE}

# troubleshooting & debugging for 3rd confirmed
n <- 28
pattern <- " confirmed "
bin <- c()

parsedsitrep <- pdf_text(pdf = paste0("covid-19/", pdf.paths[n]))

result.n <- keyword_search(parsedsitrep[[1]], # only first page
                           keyword = pattern,
                           path = FALSE,
                           surround_lines = 0)

txt <- result.n$line_text %>% str_split(pattern = "Globally", simplify = FALSE)
txt <- txt %>% str_split(pattern = "laboratory", simplify = FALSE)
txt <- txt %>% str_split(pattern = " China ", simplify = FALSE)


for (i in 1:length(txt)) {
  for (j in 1:length(txt[[i]])) {
        mess <- txt[[i]][j]
        bin <- c(bin, substr(x = mess,
                         start = regexpr(pattern = pattern, text = mess) - ifelse(pattern == " deaths ", 5, 10), # dont expect deaths to increase over 99,999 yet so minimize chance of picking up extraneous numbers
                         stop = regexpr(pattern = pattern, text = mess) + nchar(pattern)))
  }
}

if ( pattern == " deaths " ) {
  count.vec <- bin %>%
    gsub(pattern = "\\.", replacement = "") %>%
    str_split(pattern = " ") %>%
    keep_numeric()
} else if (pattern == " confirmed " ) {
  count.vec <- bin %>%
    gsub(pattern = "\\.", replacement = "") %>%
    str_split(pattern = " ") %>%
    extract_numeric() # depreciated but better than parse_number()
}

count.vec <- count.vec[!is.na(count.vec)]
count.vec

# txt
```


```{r, warning = FALSE, message = FALSE, output = FALSE}
# 

for ( i in 7:(nrow(counts)) ) {
  print(paste("Now parsing", counts$date[i]))
  
  current.confirmed <- getValues(i, " confirmed ")
  current.deaths <-  getValues(i, " deaths ")

  # CONFIRMED CASES - reports always have 3 instances of "confirmed" metrics
  instances.confirmed <- sum( !is.na( current.confirmed ))
  if ( instances.confirmed != 3) {
    warning('Cannot find 3 instances of \'confirmed\' metrics in report.')
  } else {
    counts$global[i] <- current.confirmed[1]
    counts$china[i] <- current.confirmed[2]
    counts$outside_china[i] <- current.confirmed[3]
  }
  
  # DEATHS - 
  instances.deaths <- sum( !is.na( current.deaths ))
  if ( instances.deaths == 1 ) { # only have china deaths
    counts$china_deaths[i] <- current.deaths[1]
  } else if ( instances.deaths == 2 ) {
    counts$china_deaths[i] <- current.deaths[1]
    counts$outside_china_deaths[i] <- current.deaths[2]
  } else if ( instances.deaths == 3 ) {
    counts$global_deaths[i] <- current.deaths[1]
    counts$china_deaths[i] <- current.deaths[2]
    counts$outside_china_deaths[i] <- current.deaths[3]
  }
  
}
```


```{r, warning = FALSE, message = FALSE, output = FALSE}
counts %>% qkable()

```

```{r, include = FALSE, eval = FALSE}
# case question regarding elasticity as given by Jie from Mercury Insurance

new_a <-  500 * 1000 * .40
new_b <- (500 * 6/5) * (1000 - 600) * .5


new_b - new_a


renewal_a <- 500 * 1000 * .6
renewal_b <- (500 * 4/5) * (1000 + 400) * .5

renewal_b - renewal_a

```


Most entries are correct, but we still need a considerable amount of cleaning. We find out this is because the WHO decide to use the term `"laboratory-confirmed"` instead of `"confirmed"` for a few, and generally, scraping from PDF sidebars like these can be difficult. 

Then just starting March 17, the reports stopped publishing metrics with `Global`, `China` and `Outside of China` metrics. Instead, the numbers are broken down into:

- Globally
- Western Pacific Region
- European Region
- South-East Asia
- Eastern Mediterranean Region
- Regions of the Americas
- African Region


At this point it really is faster to work in Excel. 


```{r, warning = FALSE, message = FALSE, eval = FALSE}
setwd("covid-19")

extract_MMDD <- function(input_date = Sys.Date() ) {
  # input_date : like Sys.Date() "2020-03-25" format
  return(substring(input_date, first = nchar(input_date) + 1, last = nchar(input_date) + 5))
}

# csvfilename <- paste0(getwd(), "/", extract_MMDD(), ".csv")

update_csv <- function(path_to_file = paste0(getwd(), "/", extract_MMDD(), ".csv")) {
  
  if (!file.exists(path_to_file)) {
    write_excel_csv(x = counts, 
                    path = path_to_file)
  } else {
    print( paste( ".csv for ",  ) )
  }
}

if (!file.exists(csvfilename)) {
  write_excel_csv(x = counts,
                  path = csvfilename)
} else {
  print( paste( ".csv for ", td, "already exists.") )
}
```


```{r, warning = FALSE, message = FALSE, message = FALSE}
covid <- as_tibble(read_csv(file = paste0(getwd(), "/covid-19/scrape-cleaned.csv")))
# covid <- head(scrape_cleaned, nrow(scrape_cleaned) - 2) # handle this with NA's instead 
covid$date <- mdy("1-20-2020") + 1:nrow(covid) # fix dates from csv

covid %>% qkable()
```


## Data Prep


```{r}
covid.confirmed <- covid %>% 
  select(c('date', 'global', 'china', 'outside_china')) %>% 
  pivot_longer( cols = - c('date'), names_to = 'region', values_to = 'count' )

covid.deaths <- covid %>%
  select(c('date', 'global_deaths', 'china_deaths', 'outside_china_deaths')) %>%
  pivot_longer( cols = - c('date'), names_to = 'region', values_to = 'count')

# set factor levels for plotting order
covid.confirmed$region <- factor(covid.confirmed$region, levels = c("global", "china", "outside_china"))
covid.deaths$region <- factor(covid.deaths$region, levels = c("global_deaths", "china_deaths", "outside_china_deaths"))


# change 0's in counts to NA to not be plotted (starting 3-17-2020, we only use global data)
covid.confirmed$count <- covid.confirmed$count %>% na_if(y = 0)
covid.deaths$count <- covid.deaths$count %>% na_if(y = 0)


covid.confirmed %>% qkable()
covid.deaths %>% qkable()
```

## Quick Data Visualization


```{r, warning = FALSE}
covid.confirmed %>% ggplot() + geom_jitter(aes(x = `date`, y = `count`, color = `region`)) + ggtitle('COVID-19 Confirmed Counts, Reported by WHO')

covid.confirmed %>% ggplot() + geom_point(aes(x = `date`, y = `count`, color = `region`)) + ggtitle('COVID-19 Confirmed Counts, Reported by WHO') + facet_grid( ~ `region`)


```

True to the plots, there was a huge discontinuous jump in the reported numbers of confirmed counts from February 16 to February 17. 

From Situation Report 28 on February 17, we see an explanation:

> "From today, WHO will be reporting all confirmed cases, including both
laboratory-confirmed as previously reported, and those reported as clinically
diagnosed (currently only applicable to Hubei province, China). From 13
February through 16 February, we reported only laboratory confirmed cases
for Hubei province as mentioned in the situation report published on 13
February. The change in reporting is now shown in the figures. This accounts
for the apparent large increase in cases compared to prior situation reports. 

It's important to understand this change in reporting before conducting estimations and projections. 

```{r}
covid.deaths %>% ggplot() + geom_point(aes(x = `date`, y = `count`, color = `region`)) + ggtitle('Coronavirus Death Counts, Reported by WHO')
```

There certainly are signs of stabilization in death counts within China, but it is of course still too early to say this conclusively. Notice that the cumulatie number of deaths outside of china has surpassed that from within China, and as such, the World Health Organization has changed its reporting structure to offer data at a much finer granularity. 


## Using Global Cross-verified Data

Most developers and researchers are using data gathered and provided by either the team at [Johns Hopkins CSSE](https://github.com/CSSEGISandData/COVID-19) or cleaned and normalized data from [DataHub](https://datahub.io/core/covid-19). 

Let's first take a look at DataHub's timeseries data. 


```{r, eval = FALSE}
#from DataHub
json_file <- 'https://datahub.io/core/covid-19/datapackage.json'
json_data <- fromJSON(paste(readLines(json_file), collapse=""))

# get list of all resources:
print(json_data$resources$name)

# print all tabular data(if exists any)
for (i in 1:length(json_data$resources$datahub$type)) {
  if (json_data$resources$datahub$type[i] == 'derived/csv') {
    path_to_file = json_data$resources$path[i]
    data <- read.csv(url(path_to_file))
    print(data)
  }
}

```

```{r, message = FALSE, warning = FALSE}
setwd("covid-19")
td <- Sys.Date()


covid_combined <- as_tibble(read_csv(file =  "https://datahub.io/core/covid-19/r/time-series-19-covid-combined.csv"))
covid_keycountries_counts <- as_tibble( read_csv("https://datahub.io/core/covid-19/r/key-countries-pivoted.csv") )



covid_combined <- as_tibble(read_csv(file = "https://datahub.io/core/covid-19/r/time-series-19-covid-combined.csv"))
covid_keycountries_counts <- as_tibble(read_csv("https://datahub.io/core/covid-19/r/key-countries-pivoted.csv"))


# normalize some data (US, UK)
covid_combined$`Country/Region`[which(covid_combined$`Country/Region` == "US")] <- "United States"
covid_keycountries_counts <- covid_keycountries_counts %>% rename(`United States` = US)
covid_keycountries_counts <- covid_keycountries_counts %>% rename(`United Kingdom` = `United_Kingdom`)


covid_keycountries_counts_pivoted <- covid_keycountries_counts %>%
  pivot_longer(cols = - c("Date"),
               names_to = "Regions",
               values_to = "Confirmed Counts")
                
# get death counts of key countries
keycountries <- covid_keycountries_counts_pivoted$Regions %>% unique()
    

covid_keycountries_deaths_pivoted <- covid_combined %>% 
  filter(`Country/Region` %in% keycountries) %>%
  select(c(Date, `Country/Region`, Deaths)) 

  covid_keycountries_deaths <- covid_keycountries_deaths_pivoted %>%
    pivot_wider(names_from = `Country/Region`,
                values_from = Deaths,
                values_fn = list(Deaths = max))

covid_keycountries_recovered_pivoted <- covid_combined %>% 
  filter(`Country/Region` %in% keycountries) %>%
  select(c(Date, `Country/Region`, Recovered)) 

  covid_keycountries_recovered <- covid_keycountries_recovered_pivoted %>%
    pivot_wider(names_from = `Country/Region`,
                values_from = Recovered,
                values_fn = list(Recovered = max))


mostrecent <- max(covid_keycountries_counts$Date)
ranking_confirmed <- covid_keycountries_counts_pivoted %>% 
  filter(Date == mostrecent) %>% 
  select(c(Regions, `Confirmed Counts`)) %>%
  arrange( desc(`Confirmed Counts`) ) %>%
  pull('Regions')

# reorder columns in decreasing order:
covid_keycountries_counts <- covid_keycountries_counts %>% select(c(Date, ranking_confirmed))





plotCases <- covid_keycountries_counts_pivoted %>% 
  plot_ly(
    x = ~ Date,
    y = ~ `Confirmed Counts`,
    color = ~ Regions,
    mode = 'lines'
  ) %>% layout(title = 'Confirmed Cases of COVID-19 by Region',
    legend = list(x = 0, y = 1, font = list(size = 12), bgcolor = "#F2F2F2"))

plotDeaths <- covid_keycountries_deaths_pivoted %>%
  plot_ly(
    x = ~ Date,
    y = ~ Deaths,
    color = ~ `Country/Region`,
    mode = 'lines'
  ) %>% layout(title = 'Deaths to COVID-19 by Region',
    legend = list(x = 0, y = 1, font = list(size = 12), bgcolor = "#F2F2F2"))

plotCases

plotDeaths

```



## Cautions against using premature predictions. 

There has been a lot of speculation regarding the future of the novel coronavirus, and for good reason, as the financial markets certainly adjust to news of outbreaks and as many peoples' lives and international plans are modified due to COVID-19 concerns. 

```{r}

```


According to this [The Atlantic](https://www.theatlantic.com/technology/archive/2020/03/how-many-americans-really-have-coronavirus/607348/) article, "every country's numbers are the result of a specific set of testing and accounting regimes... even though these inconsistencies are public and plain, people continue to rely on charts showing different numbers... it encourages dangerous behavior such as cutting back testing to bring a country's numbers down or slow-walking testing to keep a country's numbers low. 

Special attention needs to be paid towards policy shifts such as the recent lockdowns, curfews, and transitions. 


## Near-Future Direction


Because the total number of cases is only determined by actual tests, we should expect there to be a great deal of noise obscuring the true "signal". Any model we run should add a potential "noise" factor in the estimation, perhaps by adding a random number of cases at each data point that is relative to the population in consideration (i.e. add a number of possible cases in percentage to the population of the country). 

The following is a minor consideration, but paradoxically, if there are more cases than reported, then perhaps some areas are saturated with infected people. We can add a bias factor in the prediction that "underestimates" the value in order to add in this consideration.

We must of course add an "improvement factor" accounting for the discontinuity in reported numbers as (presumably) 

## Model Selection

Because we are only in the early phase of the outbreak on a global scale, we want to consider a few models:

**Generalized Logistic Model:**
$$\frac{d C(t)}{dt} = r C^p (t) \left( 1 - \frac{C(t)}{K}\right)$$

**Logistic Growth Model:**
$$ \frac{d C(t)}{dt} = r C(t) \left( 1 - \frac{C(t)}{K}\right)$$


**Generalized Growth Model**
$$\frac{d C(t)}{dt} = r C^p (t) $$



#### Assumptions and Remarks:

The generalized growth model can let us tune to sub-exponential growth of the COVID-19 outbreak in the early phase; however, it fails to describe the decay of the incedence rate. We can use this as a rough upper bound for our estimation, assuming that the outbreak continues to grow at a rate representated by historical data. 

Because the Logistic Growth Model and Generalized Logistic Model assume a logistic decay of the growth rate along the time horizon, we should expect these Logistic models to provide a lower bound estimate for future counts. The Generalized Logistic Model allows for an early sub-exponential growth and is better suited to explain asymmetries in the growth and decay (and may better capture government policy efforts with isolation, quarantine and lockdown).