# Code - Dr Michael Tam
# m.tam@unsw.edu.au
# Projections using NLR to Gompertz and 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-8-1")
today = as.numeric(todaydate - as.Date("2021-6-17")) + 1
time = 1445
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 is for model evaluation as the model changes over time with each new daily
# data point.  It generates the model for each day from 1 July 2021 to "today".  The parameters of the
# Gompertz and Richards' models are saved into a data frame.  The 7 and 14 days estimates and 95% CIs
# for each of the model are also saved, which allows comparison to subsequent reality.

model.history.g.d <- vector()
model.history.g.7fit <- vector()
model.history.g.7lwr <- vector()
model.history.g.7upr <- vector()
model.history.g.14fit <- vector()
model.history.g.14lwr <- vector()
model.history.g.14upr <- vector()
model.history.g.7fit.error <- vector()
model.history.g.7lwr.error <- vector()
model.history.g.7upr.error <- vector()
model.history.g.14fit.error <- vector()
model.history.g.14lwr.error <- vector()
model.history.g.14upr.error <- vector()

model.history.g.d[1:198] <- NA
model.history.g.7fit[1:198] <- NA
model.history.g.7lwr[1:198] <- NA
model.history.g.7upr[1:198] <- NA
model.history.g.14fit[1:198] <- NA
model.history.g.14lwr[1:198] <- NA
model.history.g.14upr[1:198] <- NA
model.history.g.7fit.error[1:198] <- NA
model.history.g.7lwr.error[1:198] <- NA
model.history.g.7upr.error[1:198] <- NA
model.history.g.14fit.error[1:198] <- NA
model.history.g.14lwr.error[1:198] <- NA
model.history.g.14upr.error[1:198] <- NA

model.history.r.d <- vector()
model.history.r.7fit <- vector()
model.history.r.7lwr <- vector()
model.history.r.7upr <- vector()
model.history.r.14fit <- vector()
model.history.r.14lwr <- vector()
model.history.r.14upr <- vector()
model.history.r.7fit.error <- vector()
model.history.r.7lwr.error <- vector()
model.history.r.7upr.error <- vector()
model.history.r.14fit.error <- vector()
model.history.r.14lwr.error <- vector()
model.history.r.14upr.error <- vector()

model.history.r.d[1:198] <- NA
model.history.r.7fit[1:198] <- NA
model.history.r.7lwr[1:198] <- NA
model.history.r.7upr[1:198] <- NA
model.history.r.14fit[1:198] <- NA
model.history.r.14lwr[1:198] <- NA
model.history.r.14upr[1:198] <- NA
model.history.r.7fit.error[1:198] <- NA
model.history.r.7lwr.error[1:198] <- NA
model.history.r.7upr.error[1:198] <- NA
model.history.r.14fit.error[1:198] <- NA
model.history.r.14lwr.error[1:198] <- NA
model.history.r.14upr.error[1:198] <- NA

real.7 <- vector()
real.14 <- vector()
real.7[1:198] <- NA
real.14[1:198] <- NA

for(i in 15:today) {
model.cases <- au_covid$tot[1:i]
model.day <- au_covid$day[1:i]
model.r <- drm(model.cases ~ model.day, fct = L.5())
model.g <- drm(model.cases ~ model.day, fct = G.4())
model.summary.g <- summary(model.g)
model.summary.r <- summary(model.r)

model.history.g.d[i] <- model.summary.g$coefficients[3]
model.history.r.d[i] <- model.summary.r$coefficients[3]

prediction_days <- data.frame(day = all_days)
model.g.pred <- data.frame(predict(model.g, newdata = prediction_days, interval = 'confidence'))
model.history.g.7fit[i] <- model.g.pred$Prediction[i+7]
model.history.g.7lwr[i] <- model.g.pred$Lower[i+7]
model.history.g.7upr[i] <- model.g.pred$Upper[i+7]
model.history.g.14fit[i] <- model.g.pred$Prediction[i+14]
model.history.g.14lwr[i] <- model.g.pred$Lower[i+14]
model.history.g.14upr[i] <- model.g.pred$Upper[i+14]

model.r.pred <- data.frame(predict(model.r, newdata = prediction_days, interval = 'confidence'))
model.history.r.7fit[i] <- model.r.pred$Prediction[i+7]
model.history.r.7lwr[i] <- model.r.pred$Lower[i+7]
model.history.r.7upr[i] <- model.r.pred$Upper[i+7]
model.history.r.14fit[i] <- model.r.pred$Prediction[i+14]
model.history.r.14lwr[i] <- model.r.pred$Lower[i+14]
model.history.r.14upr[i] <- model.r.pred$Upper[i+14]

real.7[i] <- au_covid$tot[i+7]
real.14[i] <- au_covid$tot[i+14]

model.history.g.7fit.error[i] <- (model.g.pred$Prediction[i+7] - real.7[i]) / real.7[i]
model.history.g.7lwr.error[i] <- (model.g.pred$Lower[i+7] - real.7[i]) / real.7[i]
model.history.g.7upr.error[i] <- (model.g.pred$Upper[i+7] - real.7[i]) / real.7[i]
model.history.g.14fit.error[i] <- (model.g.pred$Prediction[i+14] - real.14[i]) / real.14[i]
model.history.g.14lwr.error[i] <- (model.g.pred$Lower[i+14] - real.14[i]) / real.14[i]
model.history.g.14upr.error[i] <- (model.g.pred$Upper[i+14] - real.14[i]) / real.14[i]

model.history.r.7fit.error[i] <- (model.r.pred$Prediction[i+7] - real.7[i]) / real.7[i]
model.history.r.7lwr.error[i] <- (model.r.pred$Lower[i+7] - real.7[i]) / real.7[i]
model.history.r.7upr.error[i] <- (model.r.pred$Upper[i+7] - real.7[i]) / real.7[i]
model.history.r.14fit.error[i] <- (model.r.pred$Prediction[i+14] - real.14[i]) / real.14[i]
model.history.r.14lwr.error[i] <- (model.r.pred$Lower[i+14] - real.14[i]) / real.14[i]
model.history.r.14upr.error[i] <- (model.r.pred$Upper[i+14] - real.14[i]) / real.14[i]
}

model.history <- data.frame(date = au_covid$date, day = all_days, g.d = model.history.g.d, g.7fit = model.history.g.7fit, g.7lwr = model.history.g.7lwr, g.7upr = model.history.g.7upr, g.14fit = model.history.g.14fit, g.14lwr = model.history.g.14lwr, g.14upr = model.history.g.14upr, r.d = model.history.r.d, r.7fit = model.history.r.7fit, r.7lwr = model.history.r.7lwr, r.7upr = model.history.r.7upr, r.14fit = model.history.r.14fit, r.14lwr = model.history.r.14lwr, r.14upr = model.history.r.14upr)

model.history.g.error <- data.frame(date = au_covid$date, day = all_days, g.7fit = model.history.g.7fit.error, g.7lwr = model.history.g.7lwr.error, g.7upr = model.history.g.7upr.error, g.14fit = model.history.g.14fit.error, g.14lwr = model.history.g.14lwr.error, g.14upr = model.history.g.14upr.error)

model.history.r.error <- data.frame(date = au_covid$date, day = all_days, r.7fit = model.history.r.7fit.error, r.7lwr = model.history.r.7lwr.error, r.7upr = model.history.r.7upr.error, r.14fit = model.history.r.14fit.error, r.14lwr = model.history.r.14lwr.error, r.14upr = model.history.r.14upr.error)


# 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 = 'confidence'))
model.au <- data.frame(day = all_days, fit = model.au.pred$Prediction, lwr = model.au.pred$Lower, upr = model.au.pred$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


# Projections based on Gompertz model
model.au.gompertz.pred <- data.frame(predict(model.g, 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 - Richards'
# 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))


## Plot cumulative cases - Gompertz
# 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_gompertz <- 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 Gompertz growth curve", hjust=0) +

	geom_point(colour="grey20", size = 1) + 
	geom_line(aes(y=model.gompertz$fit), colour="grey20", alpha=0.5, size = 1, linetype=2) + 
	geom_ribbon(aes(ymin=model.gompertz$lwr, ymax=model.gompertz$upr), colour=NA, fill="blue", 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.gompertz$fit[today + 7], label=round(model.gompertz$fit[today + 7]), size=3) +    
	annotate("text", x=todaydate + 7, y=model.gompertz$upr[today + 7] + 50, label=round(model.gompertz$upr[today + 7]), size=2) +  
	annotate("text", x=todaydate + 7, y=model.gompertz$lwr[today + 7] - 50, label=round(model.gompertz$lwr[today + 7]), size=2) +  

	annotate("label", x=todaydate + 14, y=model.gompertz$fit[today + 14], label=round(model.gompertz$fit[today + 14]), size=3) +    
	annotate("text", x=todaydate + 14, y=model.gompertz$upr[today + 14] + 50, label=round(model.gompertz$upr[today + 14]), size=2) +  
	annotate("text", x=todaydate + 14, y=model.gompertz$lwr[today + 14] - 50, label=round(model.gompertz$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 = 'confidence', 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 = 'confidence', 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 = 'confidence', 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 = 'confidence', 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 = 'confidence', 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 = 'confidence', 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 = 'confidence', 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 = 'confidence', 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 = 'confidence', 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 = 'confidence', 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.gompertz.pred.05 <- data.frame(predict(model.g, newdata = prediction_days, interval = 'confidence', level=0.05))
model.gompertz.05 <- data.frame(day = all_days, lwr = model.gompertz.pred.05$Lower, upr = model.gompertz.pred.05$Upper)
model.gompertz.05$lwr[model.gompertz.05$lwr < 0] <- 0

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

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

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

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

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

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

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

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

model.gompertz.pred.99 <- data.frame(predict(model.g, newdata = prediction_days, interval = 'confidence', level=0.99))
model.gompertz.99 <- data.frame(day = all_days, lwr = model.gompertz.pred.99$Lower, upr = model.gompertz.pred.99$Upper)
model.gompertz.99$lwr[model.gompertz.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.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

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

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

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

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

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

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

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

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

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

model.gompertz.upr.deriv.99 <- model.gompertz.fit.deriv + (model.gompertz.99$upr - model.gompertz$fit)
model.gompertz.lwr.deriv.99 <- model.gompertz.fit.deriv - (model.gompertz$fit - model.gompertz.99$lwr)
model.gompertz.lwr.deriv.99[model.gompertz.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.gompertz.upr.deriv[1:today+10]) / 1800
sf = max(model.gompertz.upr.deriv[1:today+7]) / 1800
# sf = 800 / 1800
alphaPI = 0.08

plot_au_new_gompertz <- 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 Gompertz curve", hjust=0) +

	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv, ymax=model.gompertz.upr.deriv), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv.05, ymax=model.gompertz.upr.deriv.05), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv.15, ymax=model.gompertz.upr.deriv.15), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv.25, ymax=model.gompertz.upr.deriv.25), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv.35, ymax=model.gompertz.upr.deriv.35), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv.45, ymax=model.gompertz.upr.deriv.45), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv.55, ymax=model.gompertz.upr.deriv.55), fill="blue", colour=NA, alpha=alphaPI) +
    	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv.65, ymax=model.gompertz.upr.deriv.65), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv.75, ymax=model.gompertz.upr.deriv.75), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv.85, ymax=model.gompertz.upr.deriv.85), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.gompertz.lwr.deriv.99, ymax=model.gompertz.upr.deriv.99), fill="blue", colour=NA, alpha=0.03) +
    
	geom_point(colour="white", size = 3.2) + 
	geom_point(colour="grey20", size = 2.8) +     

	geom_line(aes(y=model.gompertz.fit.deriv), colour="grey20", alpha=0.3, size = 1, linetype=2) + 
    	geom_line(aes(y=model.gompertz.upr.deriv), colour="grey20", alpha=0.3, linetype=2) + 
    	geom_line(aes(y=model.gompertz.lwr.deriv), colour="grey20", alpha=0.3, linetype=2) + 

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

	annotate("label", x=(todaydate + 14), y=model.gompertz.fit.deriv[today + 14], label=round(model.gompertz.fit.deriv[today + 14]), size=3) +    
	annotate("text", x=(todaydate + 14), y=model.gompertz.upr.deriv[today + 14], label=round(model.gompertz.upr.deriv[today + 14]), size=3) +  
	annotate("text", x=(todaydate + 14), y=model.gompertz.lwr.deriv[today + 14], label=round(model.gompertz.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)) 



# 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 = 250 / 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.au.lwr.deriv, ymax=model.au.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)



# Evaluation of the model estimates with reality.
sf = 0.8 / 1800
e1 <- ggplot(data=au_covid, aes(x=as.Date(date), y=new)) + 
	xlab("") + 
	ylab("Error (%) from reality") +

	geom_vline(xintercept=todaydate, color = "grey") +
	geom_hline(yintercept=0, color = "grey") +
	annotate("label", x=as.Date("2021-7-2"), y=1800*sf, label="Gompertz - 7 day projection", size=5.5, hjust=0) +
	annotate("text", x=as.Date("2021-7-2"), y=1500*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) +
  
	geom_ribbon(aes(ymin=model.history.g.error$g.7lwr, ymax=model.history.g.error$g.7upr), fill="blue", colour=NA, alpha=0.1) +
    	geom_line(aes(y=model.history.g.error$g.7fit), colour="blue", size = 1, linetype=1) + 

    	annotate("label", x=(todaydate), y=0*sf, label=au_covid$date[today], size=3) +
 	scale_y_continuous(labels = scales::percent) +
	coord_cartesian(ylim = c(-0.8,0.8), xlim = c(as.Date("2021-7-1"),todaydate+4)) 

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

	geom_vline(xintercept=todaydate, color = "grey") +
	geom_hline(yintercept=0, color = "grey") +
	annotate("label", x=as.Date("2021-7-2"), y=1800*sf, label="Gompertz - 14 day projection", size=5.5, hjust=0) +
  
	geom_ribbon(aes(ymin=model.history.g.error$g.14lwr, ymax=model.history.g.error$g.14upr), fill="blue", colour=NA, alpha=0.1) +
    	geom_line(aes(y=model.history.g.error$g.14fit), colour="blue", size = 1, linetype=1) + 

    	annotate("label", x=(todaydate), y=0*sf, label=au_covid$date[today], size=3) +
 	scale_y_continuous(labels = scales::percent) +
	coord_cartesian(ylim = c(-0.8,0.8), xlim = c(as.Date("2021-7-1"),todaydate+4)) 

e3 <- ggplot(data=au_covid, aes(x=as.Date(date), y=new)) + 
	xlab("Date") + 
	ylab("Error (%) from reality") +

	geom_vline(xintercept=todaydate, color = "grey") +
	geom_hline(yintercept=0, color = "grey") +
	annotate("label", x=as.Date("2021-7-2"), y=1800*sf, label="Richards' - 7 day projection", size=5.5, hjust=0) +
  
	geom_ribbon(aes(ymin=model.history.r.error$r.7lwr, ymax=model.history.r.error$r.7upr), fill="dark green", colour=NA, alpha=0.1) +
    	geom_line(aes(y=model.history.r.error$r.7fit), colour="dark green", size = 1, linetype=1) + 

    	annotate("label", x=(todaydate), y=0*sf, label=au_covid$date[today], size=3) +
 	scale_y_continuous(labels = scales::percent) +
	coord_cartesian(ylim = c(-0.8,0.8), xlim = c(as.Date("2021-7-1"),todaydate+4)) 

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

	geom_vline(xintercept=todaydate, color = "grey") +
	geom_hline(yintercept=0, color = "grey") +
	annotate("label", x=as.Date("2021-7-2"), y=1800*sf, label="Richards' - 14 day projection", size=5.5, hjust=0) +
  
	geom_ribbon(aes(ymin=model.history.r.error$r.14lwr, ymax=model.history.r.error$r.14upr), fill="dark green", colour=NA, alpha=0.1) +
    	geom_line(aes(y=model.history.r.error$r.14fit), colour="dark green", size = 1, linetype=1) + 

    	annotate("label", x=(todaydate), y=0*sf, label=au_covid$date[today], size=3) +
 	scale_y_continuous(labels = scales::percent) +
	coord_cartesian(ylim = c(-0.8,0.8), xlim = c(as.Date("2021-7-1"),todaydate+4)) 

# grid.arrange(e1, e2, e3, e4, nrow = 2)


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

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

f1 <- ggplot(data=model.history, aes(x=as.Date(date), y=g.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") +

	annotate("label", x=as.Date("2021-7-2"), y=1800*sf, label="Projection of FINAL SIZE of the outbreak", size=5.5, hjust=0) +
	annotate("text", x=as.Date("2021-7-2"), y=1700*sf, label="Assumption: fitted to Gompertz model", hjust=0) +
	annotate("text", x=as.Date("2021-7-2"), y=1625*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_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-7-1"),todaydate+4)) 


f2 <- ggplot(data=model.history, aes(x=as.Date(date), y=r.d)) + 
	xlab("Date") + 
	ylab("Date") +
    	scale_x_date(date_breaks = "1 weeks", date_labels = "%b %d") +

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

	annotate("label", x=as.Date("2021-7-2"), y=1800*sf, label="Projection of FINAL SIZE of the outbreak", size=5.5, hjust=0) +
	annotate("text", x=as.Date("2021-7-2"), y=1700*sf, label="Assumption: fitted to Richards' growth model", hjust=0) +

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

	geom_smooth(method="gam", colour="dark green", fill="dark green", 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-7-1"),todaydate+4)) 

# grid.arrange(f1, f2, nrow = 1)


# Model outputs

plot_au
plot_au_gompertz
plot_au_new
plot_au_new_gompertz


grid.arrange(p1, p2, p3, p4, nrow = 2)
grid.arrange(e1, e2, e3, e4, nrow = 2)
grid.arrange(f1, f2, nrow = 1)


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

summary(model.g)
model.g
model.gompertz.fit.deriv