# Understanding exponential growth: The corona pandemic
# R code from https://philippmasur.de/blog/2020-03-13-exponential-growth/

library(tidyverse)
library(papaja)

# Exponential function
exp_growth <- function(start = 1, 
                  steps = 10, 
                  factor = 2) {
  x <- vector()
  x[1] <- start
  for(i in 2:steps){
    x[i] = x[i-1]*factor
  }
  return(x)
}

# Where are we after only 10 steps?
exp_growth()

# Where are we after 20 steps?
exp_growth(steps = 20)

(plot1 <- data.frame(n = 1:30,
                     value = exp_growth(steps = 30)) %>%
  ggplot(aes(x = n, y = value)) +
  geom_point() +
  geom_line() +
  theme_bw())

plot1 +
  ylim(0,750) +
  scale_x_continuous(breaks = c(1:10), limits = c(1,10))

plot1 +
  scale_y_log10()

perc_growth <- function(start = 1, steps, perc) {
  x <- vector()
  x[1] <- start
  for(i in 2:steps){
    x[i] = x[i-1]+ perc*x[i-1]
  }
  return(x)
}

# Recreating the exponential growth from above
perc_growth(start = 1, steps = 10, perc = 1)

(plot2 <- data.frame (n = 1:30,
            "5" = perc_growth(start = 1, steps = 30, perc = 0.5),
            "10" = perc_growth(start = 1, steps = 30, perc = 0.10),
            "20" = perc_growth(start = 1, steps = 30, perc = 0.20)) %>%
  gather(key, value, -n) %>%
  ggplot(aes(x = n, y = value, color = key)) +
  geom_point() +
  geom_line() +
  theme_bw())

plot2 +
  xlim(0, 10) +
  ylim(0, 50)

# JHU CSSE data as of 12 March 2020 (pinned to that day's commit)
url <- paste0("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/",
              "82e1397911272be823cb93cbc19b62551608b8dc/",
              "csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv")

d <- read.csv(url) %>%
  select(-Lat, -Long)

d <- d %>%
  rename(country = Country.Region) %>%
  filter(country == "Germany") %>%
  gather(date, value, -country) %>%
  separate(date, c("x", "date"), sep = 1) %>%
  mutate(cases = as.numeric(value),
         date = lubridate::mdy(date)) %>%
  as_tibble

# plot
(plot3 <- d %>%
  ggplot(aes(x = date, y = cases)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  labs(x = "date", y = "cases"))

d %>%
  ggplot(aes(x = date, y = cases + 1)) +
  geom_point() +
  geom_smooth(se = F, color = "black") +
  theme_bw() +
  scale_y_log10() +
  labs(x = "date", y = "cases")

d %>%
  mutate(cases_day = cases - lag(cases),
         percent = cases_day/cases) %>%
  summarize(growthrate = mean(percent, na.rm = T),
            se = psych::describe(percent)$se,
            lower = growthrate - 1.96*se,
            upper = growthrate + 1.96*se)

data.frame(n = 1:30,
           cases = perc_growth(start = 3156, steps = 30, perc = 0.154),
           lower = perc_growth(start = 3156, steps = 30, perc = 0.0918),
           upper = perc_growth(start = 3156, steps = 30, perc = 0.216)) %>%
  gather(key, cases, -n) %>%
  ggplot(aes(x = n, y = cases, color = key)) +
  geom_line() +
  theme_bw() +
  ylim(0, 200000) +
  scale_color_manual(values=c("black", "grey", "grey")) +
  scale_x_continuous(breaks = seq(0,30, by = 5)) +
  labs(x = "days", title = "Simulated growth in Germany")

(plot4 <- d %>%
  filter(cases >= 100) %>%
  mutate("cases of corona" = cases,
         "33% growth curve" = perc_growth(start = 130, steps = 11, perc = .33)) %>%
  select(date, "cases of corona", "33% growth curve") %>%
  gather(key, value, -date) %>%
  ggplot(aes(x = date, y = value, color = key)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  scale_color_manual(values = c("grey", "black")) +
  labs(x = "date", y = "cases", color = ""))