# Code - Dr Michael Tam
# m.tam@unsw.edu.au
# Projections using NLR to Richards growth curve of Australian (NSW) COVID-19 outbreak in June 2021
# The data excludes cases known to have been acquired from overseas (for instance, detected in hotel quarantine)

# Install libraries - might need to do this the very first time you run this code
# Delete the hash symbol to take it off the comments
# install.packages("drc")
# install.packages("ggplot2")
# install.packages("tidyquant")
# install.packages("gridExtra")
# install.packages("readxl")
# install.packages("scales")
# install.packages("remotes")
# remotes::install_github("OnofriAndreaPG/aomisc")


# Load libraries
library(drc)
library(aomisc)
library(ggplot2)
library(tidyquant)
library(gridExtra)
library(readxl)
library(scales)


# I'm assuming that the data file for today is an Excel spreadsheet named "au_covid.xlsx"
# If it isn't, rename it before uploading it into RStudio.

au_covid <- read_excel("au_covid.xlsx", col_types = c("date", 
    "numeric", "numeric", "numeric", "text"))

View(au_covid)


# Setting the "day". 17/6/2021: today=1  The rest should be obvious from context.
# all_days is just an integer list from 1 (17/6/2021) to 198 (31/12/2021).
# 17/6/2021 was the date of the first reported local cases in the outbreak
todaydate = as.Date("2021-7-15")
today = as.numeric(todaydate - as.Date("2021-6-17")) + 1
time = 1200
all_days = 1:198


# Richards' growth curve (generalised logistic function) (5 parameters)
# This uses:
# "drc": Ritz, C., Baty, F., Streibig, J. C., Gerhard, D. (2015) Dose-Response
# Analysis Using R PLOS ONE, 10(12), e0146021
# "aomisc": https://rdrr.io/github/OnofriAndreaPG/aomisc/


# This first section here isn't technically necessary. The code here is for the sake of interest.
# One of the things I want to do is to see how the model changes over time with each new daily
# data point.  The following code here generates the original Gompertz model for each day from
# 1 July 2021 to "today", and outputs the 4 parameters of the model into a data frame.  The interesting
# thing here for me is how the asymptote (co-efficient "d" in the model summary) changes over time.

model.history.b <- vector()
model.history.c <- vector()
model.history.d <- vector()
model.history.e <- vector()

model.history.b[1:198] <- NA
model.history.c[1:198] <- NA
model.history.d[1:198] <- NA
model.history.e[1:198] <- NA

for(i in 15:today) {
model.au.cases <- au_covid$tot[1:i]
model.au.day <- au_covid$day[1:i]
# model.au.richards <- drm(model.au.cases ~ model.au.day, fct = L.5())
model.au.gompertz <- drm(model.au.cases ~ model.au.day, fct = G.4())
model.au.summary <- summary(model.au.gompertz)
model.history.b[i] <- model.au.summary$coefficients[1]
model.history.c[i] <- model.au.summary$coefficients[2]
model.history.d[i] <- model.au.summary$coefficients[3]
model.history.e[i] <- model.au.summary$coefficients[4]
}

model.history <- data.frame(date = au_covid$date, day = all_days, b = model.history.b, c = model.history.c, d = model.history.d, e = model.history.e)


# The Richards' growth curve model
model.au.cases <- au_covid$tot[1:today]
model.au.day <- au_covid$day[1:today]
model.au.richards <- drm(model.au.cases ~ model.au.day, fct = L.5())
model.au.summary <- summary(model.au.richards)


# Projections based on Richards' growth model
prediction_days <- data.frame(model.au.day = all_days)
model.au.pred <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction'))
model.au.predCI <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'confidence'))

model.au <- data.frame(day = all_days, fit = model.au.pred$Prediction, lwr = model.au.pred$Lower, upr = model.au.pred$Upper)
model.auCI <- data.frame(day = all_days, fit = model.au.predCI$Prediction, lwr = model.au.predCI$Lower, upr = model.au.predCI$Upper)

# Replace any predicted values < 0 to be equal to zero.
# Obviously, there can't be less than 0 cases.
model.au$lwr[model.au$lwr < 0] <- 0
model.auCI$lwr[model.au$lwr < 0] <- 0


# Projections based on Gompertz model - deprecated
model.au.gompertz.pred <- data.frame(predict(model.au.gompertz, newdata = prediction_days, interval = 'confidence'))
model.gompertz <- data.frame(day = all_days, fit = model.au.gompertz.pred$Prediction, lwr = model.au.gompertz.pred$Lower, upr = model.au.gompertz.pred$Upper)
model.gompertz$lwr[model.gompertz$lwr < 0] <- 0


### Plots
## Plot cumulative cases
# scaling factor for fixed elements on chart
# sf = model.au$upr[today+14] * 1.05 / 1800
sf = model.au$upr[today+10] * 1.05 / 1800

plot_au <- ggplot(data=au_covid, aes(x=as.Date(date), y=tot)) + 
	xlab("") + 
	ylab("Cumulative Local COVID-19 cases (NSW) starting June 17") + 
    	scale_x_date(date_breaks = "1 weeks", date_labels = "%b %d") +

	geom_vline(xintercept=todaydate, color = "grey") +
	annotate("label", x=todaydate, y=50*sf, label=au_covid$date[today], size=3) +
	geom_vline(xintercept=todaydate + 7, color = "grey") +
	annotate("label", x=todaydate + 7, y=50*sf, label=au_covid$date[today+7], size=3) +
	geom_vline(xintercept=todaydate + 14, color = "grey") +
	annotate("label", x=todaydate + 14, y=50*sf, label=au_covid$date[today+14], size=3) +

	annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Projection of total COVID-19 cases", size=5.5, hjust=0) +
	annotate("text", x=as.Date("2021-6-20"), y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) +
	annotate("text", x=as.Date("2021-6-20"), y=1625*sf, label="Assumption: fitted to Richards' growth curve", hjust=0) +

	geom_point(colour="grey20", size = 1) + 
	geom_line(aes(y=model.au$fit), colour="grey20", alpha=0.5, size = 1, linetype=2) + 
	geom_ribbon(aes(ymin=model.au$lwr, ymax=model.au$upr), colour=NA, fill="dark green", alpha=0.1) +
    
	annotate("label", x=todaydate, y=au_covid$tot[today], label=au_covid$tot[today], size=3) +    

	annotate("label", x=todaydate + 7, y=model.au$fit[today + 7], label=round(model.au$fit[today + 7]), size=3) +    
	annotate("text", x=todaydate + 7, y=model.au$upr[today + 7] + 50, label=round(model.au$upr[today + 7]), size=2) +  
	annotate("text", x=todaydate + 7, y=model.au$lwr[today + 7] - 50, label=round(model.au$lwr[today + 7]), size=2) +  

	annotate("label", x=todaydate + 14, y=model.au$fit[today + 14], label=round(model.au$fit[today + 14]), size=3) +    
	annotate("text", x=todaydate + 14, y=model.au$upr[today + 14] + 50, label=round(model.au$upr[today + 14]), size=2) +  
	annotate("text", x=todaydate + 14, y=model.au$lwr[today + 14] - 50, label=round(model.au$lwr[today + 14]), size=2) +
    	# coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"), todaydate + 16))
    	coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"), todaydate + 10))


### Plots - new cases

model.au.pred.05 <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction', level=0.05))
model.au.05 <- data.frame(day = all_days, lwr = model.au.pred.05$Lower, upr = model.au.pred.05$Upper)
model.au.05$lwr[model.au.05$lwr < 0] <- 0

model.au.pred.15 <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction', level=0.15))
model.au.15 <- data.frame(day = all_days, lwr = model.au.pred.15$Lower, upr = model.au.pred.15$Upper)
model.au.15$lwr[model.au.15$lwr < 0] <- 0

model.au.pred.25 <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction', level=0.25))
model.au.25 <- data.frame(day = all_days, lwr = model.au.pred.25$Lower, upr = model.au.pred.25$Upper)
model.au.25$lwr[model.au.25$lwr < 0] <- 0

model.au.pred.35 <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction', level=0.35))
model.au.35 <- data.frame(day = all_days, lwr = model.au.pred.35$Lower, upr = model.au.pred.35$Upper)
model.au.35$lwr[model.au.35$lwr < 0] <- 0

model.au.pred.45 <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction', level=0.45))
model.au.45 <- data.frame(day = all_days, lwr = model.au.pred.45$Lower, upr = model.au.pred.45$Upper)
model.au.45$lwr[model.au.45$lwr < 0] <- 0

model.au.pred.55 <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction', level=0.55))
model.au.55 <- data.frame(day = all_days, lwr = model.au.pred.55$Lower, upr = model.au.pred.55$Upper)
model.au.55$lwr[model.au.55$lwr < 0] <- 0

model.au.pred.65 <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction', level=0.65))
model.au.65 <- data.frame(day = all_days, lwr = model.au.pred.65$Lower, upr = model.au.pred.65$Upper)
model.au.65$lwr[model.au.65$lwr < 0] <- 0

model.au.pred.75 <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction', level=0.75))
model.au.75 <- data.frame(day = all_days, lwr = model.au.pred.75$Lower, upr = model.au.pred.75$Upper)
model.au.75$lwr[model.au.75$lwr < 0] <- 0

model.au.pred.85 <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction', level=0.85))
model.au.85 <- data.frame(day = all_days, lwr = model.au.pred.85$Lower, upr = model.au.pred.85$Upper)
model.au.75$lwr[model.au.75$lwr < 0] <- 0

model.au.pred.99 <- data.frame(predict(model.au.richards, newdata = prediction_days, interval = 'prediction', level=0.99))
model.au.99 <- data.frame(day = all_days, lwr = model.au.pred.99$Lower, upr = model.au.pred.99$Upper)
model.au.99$lwr[model.au.99$lwr < 0] <- 0



model.au.fit.fn <- splinefun(all_days, model.au$fit)
model.au.fit.deriv <- model.au.fit.fn(all_days, deriv = 1)

model.au.upr.deriv <- model.au.fit.deriv + (model.au$upr - model.au$fit)
model.au.lwr.deriv <- model.au.fit.deriv - (model.au$fit - model.au$lwr)
model.au.lwr.deriv[model.au.lwr.deriv < 0] <- 0


model.auCI.fit.fn <- splinefun(all_days, model.auCI$fit)
model.auCI.fit.deriv <- model.auCI.fit.fn(all_days, deriv = 1)

model.auCI.upr.deriv <- model.auCI.fit.deriv + (model.auCI$upr - model.auCI$fit)
model.auCI.lwr.deriv <- model.auCI.fit.deriv - (model.auCI$fit - model.auCI$lwr)
model.auCI.lwr.deriv[model.auCI.lwr.deriv < 0] <- 0


model.gompertz.fit.fn <- splinefun(all_days, model.gompertz$fit)
model.gompertz.fit.deriv <- model.gompertz.fit.fn(all_days, deriv = 1)

model.gompertz.upr.deriv <- model.gompertz.fit.deriv + (model.gompertz$upr - model.gompertz$fit)
model.gompertz.lwr.deriv <- model.gompertz.fit.deriv - (model.gompertz$fit - model.gompertz$lwr)
model.gompertz.lwr.deriv[model.gompertz.lwr.deriv < 0] <- 0


model.au.upr.deriv.05 <- model.au.fit.deriv + (model.au.05$upr - model.au$fit)
model.au.lwr.deriv.05 <- model.au.fit.deriv - (model.au$fit - model.au.05$lwr)
model.au.lwr.deriv.05[model.au.lwr.deriv.05 < 0] <- 0

model.au.upr.deriv.15 <- model.au.fit.deriv + (model.au.15$upr - model.au$fit)
model.au.lwr.deriv.15 <- model.au.fit.deriv - (model.au$fit - model.au.15$lwr)
model.au.lwr.deriv.15[model.au.lwr.deriv.15 < 0] <- 0

model.au.upr.deriv.25 <- model.au.fit.deriv + (model.au.25$upr - model.au$fit)
model.au.lwr.deriv.25 <- model.au.fit.deriv - (model.au$fit - model.au.25$lwr)
model.au.lwr.deriv.25[model.au.lwr.deriv.25 < 0] <- 0

model.au.upr.deriv.35 <- model.au.fit.deriv + (model.au.35$upr - model.au$fit)
model.au.lwr.deriv.35 <- model.au.fit.deriv - (model.au$fit - model.au.35$lwr)
model.au.lwr.deriv.35[model.au.lwr.deriv.35 < 0] <- 0

model.au.upr.deriv.45 <- model.au.fit.deriv + (model.au.45$upr - model.au$fit)
model.au.lwr.deriv.45 <- model.au.fit.deriv - (model.au$fit - model.au.45$lwr)
model.au.lwr.deriv.45[model.au.lwr.deriv.45 < 0] <- 0

model.au.upr.deriv.55 <- model.au.fit.deriv + (model.au.55$upr - model.au$fit)
model.au.lwr.deriv.55 <- model.au.fit.deriv - (model.au$fit - model.au.55$lwr)
model.au.lwr.deriv.55[model.au.lwr.deriv.55 < 0] <- 0

model.au.upr.deriv.65 <- model.au.fit.deriv + (model.au.65$upr - model.au$fit)
model.au.lwr.deriv.65 <- model.au.fit.deriv - (model.au$fit - model.au.65$lwr)
model.au.lwr.deriv.65[model.au.lwr.deriv.65 < 0] <- 0

model.au.upr.deriv.75 <- model.au.fit.deriv + (model.au.75$upr - model.au$fit)
model.au.lwr.deriv.75 <- model.au.fit.deriv - (model.au$fit - model.au.75$lwr)
model.au.lwr.deriv.75[model.au.lwr.deriv.75 < 0] <- 0

model.au.upr.deriv.85 <- model.au.fit.deriv + (model.au.85$upr - model.au$fit)
model.au.lwr.deriv.85 <- model.au.fit.deriv - (model.au$fit - model.au.85$lwr)
model.au.lwr.deriv.85[model.au.lwr.deriv.85 < 0] <- 0

model.au.upr.deriv.99 <- model.au.fit.deriv + (model.au.99$upr - model.au$fit)
model.au.lwr.deriv.99 <- model.au.fit.deriv - (model.au$fit - model.au.99$lwr)
model.au.lwr.deriv.99[model.au.lwr.deriv.99 < 0] <- 0


# scaling factor for fixed elements on chart
# sf = max(model.au.upr.deriv[1:today+10]) / 1800
sf = max(model.au.upr.deriv[1:today+7]) / 1800
# sf = 800 / 1800
alphaPI = 0.08

plot_au_new <- ggplot(data=au_covid, aes(x=as.Date(date), y=new)) + 
	xlab("") + 
	ylab("Daily Local COVID-19 cases (NSW)") +
	scale_x_date(date_breaks = "1 weeks", date_labels = "%b %d") +
	geom_vline(xintercept=todaydate, color = "grey") +
    	geom_vline(xintercept=todaydate + 7, color = "grey") +
    	geom_vline(xintercept=todaydate + 14, color = "grey") +

	geom_vline(xintercept=as.Date("2021-6-25"), color = "red", size=1) +
	annotate("text", x=as.Date("2021-6-25"), y=1400*sf, label="Lockdown (25/6) ", hjust=1) +
	geom_vline(xintercept=as.Date("2021-7-9"), color = "red", size=1) +
	annotate("text", x=as.Date("2021-7-9"), y=1400*sf, label="Intensifies (9/7) ", hjust=1) +

	annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Projection of new daily cases of COVID-19", size=5.5, hjust=0) +
	annotate("text", x=as.Date("2021-6-20"), y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) +
	annotate("text", x=as.Date("2021-6-20"), y=1625*sf, label="Assumption: fitted to Richards' growth curve", hjust=0) +

	geom_ribbon(aes(ymin=model.au.lwr.deriv, ymax=model.au.upr.deriv), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.05, ymax=model.au.upr.deriv.05), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.15, ymax=model.au.upr.deriv.15), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.25, ymax=model.au.upr.deriv.25), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.35, ymax=model.au.upr.deriv.35), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.45, ymax=model.au.upr.deriv.45), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.55, ymax=model.au.upr.deriv.55), fill="dark green", colour=NA, alpha=alphaPI) +
    	geom_ribbon(aes(ymin=model.au.lwr.deriv.65, ymax=model.au.upr.deriv.65), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.75, ymax=model.au.upr.deriv.75), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.85, ymax=model.au.upr.deriv.85), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.99, ymax=model.au.upr.deriv.99), fill="dark green", colour=NA, alpha=0.03) +
    
	geom_point(colour="white", size = 3.2) + 
	geom_point(colour="grey20", size = 2.8) +     
	geom_line(aes(y=model.au.fit.deriv), colour="grey20", alpha=0.3, size = 1, linetype=2) + 
    	geom_line(aes(y=model.au.upr.deriv), colour="grey20", alpha=0.3, linetype=2) + 
    	geom_line(aes(y=model.au.lwr.deriv), colour="grey20", alpha=0.3, linetype=2) + 

	annotate("label", x=(todaydate + 7), y=model.au.fit.deriv[today + 7], label=round(model.au.fit.deriv[today + 7]), size=3) +    
	annotate("text", x=(todaydate + 7), y=model.au.upr.deriv[today + 7], label=round(model.au.upr.deriv[today + 7]), size=3) +  
	annotate("text", x=(todaydate + 7), y=model.au.lwr.deriv[today + 7], label=round(model.au.lwr.deriv[today + 7]), size=3) +  

	annotate("label", x=(todaydate + 14), y=model.au.fit.deriv[today + 14], label=round(model.au.fit.deriv[today + 14]), size=3) +    
	annotate("text", x=(todaydate + 14), y=model.au.upr.deriv[today + 14], label=round(model.au.upr.deriv[today + 14]), size=3) +  
	annotate("text", x=(todaydate + 14), y=model.au.lwr.deriv[today + 14], label=round(model.au.lwr.deriv[today + 14]), size=3) +

    	annotate("label", x=(todaydate), y=800, label=au_covid$date[today], size=3) +
	annotate("label", x=(todaydate + 7), y=85*sf, label=au_covid$date[today+7], size=2.5) +
	annotate("label", x=(todaydate + 14), y=85*sf, label=au_covid$date[today+14], size=2.5) +
	# coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),todaydate+16)) 
	coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),todaydate+10)) 


# scaling factor for fixed elements on chart
sf = max(model.au.upr.deriv[1:today+14]) / 1800
# sf = 300 / 1800

# Plot for animation
plot_au_video <- ggplot(data=au_covid, aes(x=as.Date(date), y=new)) + 
	xlab("") + 
	ylab("Daily Local COVID-19 cases (NSW)") + 
    	scale_x_date(date_breaks = "1 weeks", date_labels = "%b %d") +

	geom_vline(xintercept=as.Date("2021-6-25"), color = "red", size=1) +
	annotate("text", x=as.Date("2021-6-25"), y=1500*sf, label="Lockdown (25/6) ", hjust=1) +
	geom_vline(xintercept=as.Date("2021-7-9"), color = "red", size=1) +
	annotate("text", x=as.Date("2021-7-9"), y=1500*sf, label="Intensifies (9/7) ", hjust=1) +
	geom_vline(xintercept=todaydate, color = "grey") +

	annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Projection of new daily cases of COVID-19", size=5.5, hjust=0) +
	annotate("text", x=as.Date("2021-6-20"), y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) +
	annotate("text", x=as.Date("2021-6-20"), y=1625*sf, label="Assumption: fitted to Richards' growth curve", hjust=0) +

	geom_ribbon(aes(ymin=model.au.lwr.deriv, ymax=model.au.upr.deriv), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.05, ymax=model.au.upr.deriv.05), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.15, ymax=model.au.upr.deriv.15), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.25, ymax=model.au.upr.deriv.25), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.35, ymax=model.au.upr.deriv.35), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.45, ymax=model.au.upr.deriv.45), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.55, ymax=model.au.upr.deriv.55), fill="dark green", colour=NA, alpha=alphaPI) +
    	geom_ribbon(aes(ymin=model.au.lwr.deriv.65, ymax=model.au.upr.deriv.65), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.75, ymax=model.au.upr.deriv.75), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.85, ymax=model.au.upr.deriv.85), fill="dark green", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.99, ymax=model.au.upr.deriv.99), fill="dark green", colour=NA, alpha=0.03) +
    
	geom_point(colour="white", size = 3.2) + 
	geom_point(colour="grey20", size = 2.8) +     
	geom_line(aes(y=model.au.fit.deriv), colour="grey20", alpha=0.3, size = 1, linetype=2) + 
    	geom_line(aes(y=model.au.upr.deriv), colour="grey20", alpha=0.3, linetype=2) + 
    	geom_line(aes(y=model.au.lwr.deriv), colour="grey20", alpha=0.3, linetype=2) + 

	annotate("label", x=as.Date("2021-7-15"), y=model.au.fit.deriv[29], label=round(model.au.fit.deriv[29]), size=3) +    
	annotate("text", x=as.Date("2021-7-15"), y=model.au.upr.deriv[29], label=round(model.au.upr.deriv[29]), size=3) +  
	annotate("text", x=as.Date("2021-7-15"), y=model.au.lwr.deriv[29], label=round(model.au.lwr.deriv[29]), size=3) + 

	annotate("label", x=as.Date("2021-7-22"), y=model.au.fit.deriv[36], label=round(model.au.fit.deriv[36]), size=3) +    
	annotate("text", x=as.Date("2021-7-22"), y=model.au.upr.deriv[36], label=round(model.au.upr.deriv[36]), size=3) +  
	annotate("text", x=as.Date("2021-7-22"), y=model.au.lwr.deriv[36], label=round(model.au.lwr.deriv[36]), size=3) + 

    	annotate("label", x=todaydate, y=1000, label=au_covid$date[today], size=3) +
	coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),todaydate+16))


# Comparison of data trends (moving averages and generalised additive model), the Gompertz model, and the Richards'
# growth curve model
# scaling factor for fixed elements on chart
sf = 140 / 1800

p1 <- ggplot(data=au_covid, aes(x=as.Date(date), y=new)) + 
	xlab("") + 
	ylab("Confirmed local COVID-19 cases (NSW)") +

	geom_vline(xintercept=todaydate, color = "grey") +
	annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Moving Avg (7 days)", size=5.5, hjust=0) +
	annotate("text", x=as.Date("2021-6-20"), y=1600*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) +
  
	geom_point(colour="white", size = 3.2) + 
	geom_point(colour="grey20", size = 2.8) +     

    	geom_ma(n=7, colour="orange", linetype=1, size = 1) +

    	annotate("label", x=(todaydate), y=0*sf, label=au_covid$date[today], size=3) +
	coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),todaydate+4)) 

p2 <- ggplot(data=au_covid, aes(x=as.Date(date), y=new)) + 
	xlab("") + 
	ylab("") +

	geom_vline(xintercept=todaydate, color = "grey") +

	annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Generalised additive model", size=5.5, hjust=0) +

	geom_point(colour="white", size = 3.2) + 
	geom_point(colour="grey20", size = 2.8) +     

    	geom_smooth(colour="red3", fill="red", alpha=0.1, method="gam") + 

    	annotate("label", x=(todaydate), y=0*sf, label=au_covid$date[today], size=3) +
	coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),todaydate+4)) 

p3 <- ggplot(data=au_covid, aes(x=as.Date(date), y=new)) + 
	xlab("") + 
	ylab("Confirmed local COVID-19 cases (NSW)") +

	geom_vline(xintercept=todaydate, color = "grey") +

	annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Gompertz model projections", size=5.5, hjust=0) +

	geom_point(colour="white", size = 3.2) + 
	geom_point(colour="grey20", size = 2.8) +     

	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv, ymax=model.gompertz.upr.deriv), fill="blue", colour=NA, alpha=0.1) +
    	geom_line(aes(y=model.gompertz.fit.deriv), colour="blue", size = 1, linetype=1) + 

    	annotate("label", x=(todaydate), y=0*sf, label=au_covid$date[today], size=3) +
	coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),todaydate+4)) 

p4 <- ggplot(data=au_covid, aes(x=as.Date(date), y=new)) + 
	xlab("") + 
	ylab("") +

	geom_vline(xintercept=todaydate, color = "grey") +

	annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Richards' curve projections", size=5.5, hjust=0) +

	geom_point(colour="white", size = 3.2) + 
	geom_point(colour="grey20", size = 2.8) +     

	geom_ribbon(aes(ymin=model.auCI.lwr.deriv, ymax=model.auCI.upr.deriv), fill="dark green", colour=NA, alpha=0.1) +
    	geom_line(aes(y=model.au.fit.deriv), colour="dark green", size = 1, linetype=1) + 

    	annotate("label", x=(todaydate), y=0*sf, label=au_covid$date[today], size=3) +
	coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),todaydate+4)) 

# grid.arrange(p1, p2, p3, p4, nrow = 2)


# Gompertz projection asymptote (projected final size of outbreak)
# How the model changed with time. This is parameter "d" in the Gompertz model.

# Scaling factor for fixed chart elements
sf = max(model.history$d[15:today])*1.05 / 1800

plot.gompertz.asymptote <- ggplot(data=model.history, aes(x=as.Date(date), y=d)) + 
	xlab("") + 
	ylab("Total number of COVID-19 cases") +
    	scale_x_date(date_breaks = "1 weeks", date_labels = "%b %d") +

	geom_vline(xintercept=todaydate, color = "grey") +

	geom_vline(xintercept=as.Date("2021-6-25"), color = "red", size=1) +
	annotate("text", x=as.Date("2021-6-25"), y=1500*sf, label="Lockdown (25/6) ", hjust=1) +
	geom_vline(xintercept=as.Date("2021-7-9"), color = "red", size=1) +
	annotate("text", x=as.Date("2021-7-9"), y=1500*sf, label="Intensifies (9/7) ", hjust=1) +
	geom_vline(xintercept=todaydate, color = "grey") +

	annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Projection of FINAL SIZE of the outbreak (SPECULATIVE)", size=5.5, hjust=0) +
	annotate("text", x=as.Date("2021-6-20"), y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) +
	annotate("text", x=as.Date("2021-6-20"), y=1625*sf, label="Assumption: fitted to Gompertz model", hjust=0) +

	geom_point(colour="white", size = 3.2) + 
	geom_point(colour="grey20", size = 2.8) +     

	geom_smooth(method="gam", colour="blue", fill="blue", size=1, linetype=1, alpha=0.1) +

    	annotate("label", x=(todaydate), y=500, label=au_covid$date[today], size=3) +
	coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),todaydate)) 



# Model outputs

plot_au
plot_au_new
# plot_au_video

# plot_au_video + 
    annotate("label", x=as.Date("2021-6-20"), y=1000*sf, label="INTERPRET WITH EXTREME CAUTION", colour="red", hjust=0, size=5) + 
    scale_x_date(date_breaks = "2 weeks", date_labels = "%b %d") +
    coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),as.Date("2021-12-31")))

grid.arrange(p1, p2, p3, p4, nrow = 2)

plot.gompertz.asymptote


au_covid$date[today]
model.au.summary
model.au
model.au.fit.deriv