# Code - Dr Michael Tam
# m.tam@unsw.edu.au
# Projections using NLR to Gompertz equation of Australian COVID-19 outbreak in July 2020

# 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("remotes")
# remotes::install_github("OnofriAndreaPG/aomisc")


# Load libraries
library(drc)
library(aomisc)
library(ggplot2)
library(readxl)

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

View(au_covid)

today = 60
time = 2130
future_week = today + 7
all_days = 1:93

# Gompertz models (4 parameters)
# See: https://www.statforbiology.com/nonlinearregression/usefulequations
# 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/

model.au.cases <- au_covid$tot[1:today]
model.au.day <- au_covid$day[1:today]
model.au.gompertz <- drm(model.au.cases ~ model.au.day, fct = G.4())
model.au.summary <- summary(model.au.gompertz)

prediction_days <- data.frame(model.au.day = all_days)
model.au.pred <- data.frame(predict(model.au.gompertz, newdata = prediction_days, interval = 'prediction'))

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


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

plot_au <- ggplot(data=au_covid, aes(x=day, y=tot)) + 
	xlab("Days from June 1") + 
	ylab("Confirmed COVID-19 cases (Australia) starting June 1") + 
	xlim(1,today + 16) +
	ylim(-200*sf,1800*sf) +

	geom_vline(xintercept=today, color = "grey") +
	annotate("label", x=(today), y=50*sf, label=au_covid$date[today], size=3) +
	geom_vline(xintercept=(today + 7), color = "grey") +
	annotate("label", x=(today + 7), y=50*sf, label=au_covid$date[today+7], size=3) +
	geom_vline(xintercept=(today + 14), color = "grey") +
	annotate("label", x=(today + 14), y=50*sf, label=au_covid$date[today+14], size=3) +
    
    geom_hline(yintercept=-100*sf) + 
	annotate("text", x=1, y=-150*sf, label="Jun 1") +
	annotate("text", x=15, y=-150*sf, label="Jun 15") +
	annotate("text", x=30+1, y=-150*sf, label="Jul 1") +
	annotate("text", x=30+15, y=-150*sf, label="Jul 15") +   
	annotate("text", x=30+31+1, y=-150*sf, label="Aug 1") +

	annotate("label", x=((today + 16)/2), y=1800*sf, label="Projection of COVID-19 cases", size=5.5) +
	annotate("text", x=((today + 16)/2), y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10")) +
	annotate("text", x=((today + 16)/2), y=1625*sf, label="Assumption: fitted to Gompertz model") +

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

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

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


### Plots - new cases

model.au.pred.05 <- data.frame(predict(model.au.gompertz, 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.gompertz, 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.gompertz, 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.gompertz, 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.gompertz, 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.gompertz, 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.gompertz, 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.gompertz, 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.gompertz, 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.gompertz, 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.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.99[1:today+16]) / 1800
alphaPI = 0.08

plot_au_new <- ggplot(data=au_covid, aes(x=day, y=new)) + 
	xlab("Days from Jun 1") + 
	ylab("Confirmed COVID-19 cases (Australia)") +
	# xlim(1,today + 16) +
	# ylim(-250*sf,1800*sf) +
	geom_vline(xintercept=today, color = "grey") +
    geom_vline(xintercept=(today + 7), color = "grey") +
    geom_vline(xintercept=(today + 14), color = "grey") +

    geom_hline(yintercept=-150*sf) + 
	annotate("text", x=1, y=-200*sf, label="Jun 1") +
	annotate("text", x=15, y=-200*sf, label="Jun 15") +
	annotate("text", x=30+1, y=-200*sf, label="Jul 1") +
	annotate("text", x=30+15, y=-200*sf, label="Jul 15") +   
	annotate("text", x=30+31+1, y=-200*sf, label="Aug 1") +
	annotate("text", x=30+31+15, y=-200*sf, label="Aug 15") +
	annotate("text", x=30+31+31+1, y=-200*sf, label="Sep 1") +

	annotate("label", x=5, y=1800*sf, label="Projection of new daily cases of COVID-19", size=5.5, hjust=0) +
	annotate("text", x=5, y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) +
	annotate("text", x=5, y=1625*sf, label="Assumption: fitted to Gompertz model", hjust=0) +

	geom_ribbon(aes(ymin=model.au.lwr.deriv, ymax=model.au.upr.deriv), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.05, ymax=model.au.upr.deriv.05), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.15, ymax=model.au.upr.deriv.15), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.25, ymax=model.au.upr.deriv.25), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.35, ymax=model.au.upr.deriv.35), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.45, ymax=model.au.upr.deriv.45), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.55, ymax=model.au.upr.deriv.55), fill="blue", colour=NA, alpha=alphaPI) +
    geom_ribbon(aes(ymin=model.au.lwr.deriv.65, ymax=model.au.upr.deriv.65), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.75, ymax=model.au.upr.deriv.75), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.85, ymax=model.au.upr.deriv.85), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.99, ymax=model.au.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.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=(today + 7), y=model.au.fit.deriv[today + 7], label=round(model.au.fit.deriv[today + 7]), size=3) +    
	annotate("text", x=(today + 7), y=model.au.upr.deriv[today + 7], label=round(model.au.upr.deriv[today + 7]), size=3) +  
	annotate("text", x=(today + 7), y=model.au.lwr.deriv[today + 7], label=round(model.au.lwr.deriv[today + 7]), size=3) +  

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

    annotate("label", x=(today), y=-85*sf, label=au_covid$date[today], size=3) +
	annotate("label", x=(today + 7), y=-85*sf, label=au_covid$date[today+7], size=3) +
	annotate("label", x=(today + 14), y=-85*sf, label=au_covid$date[today+14], size=3) +
	coord_cartesian(ylim = c(-250*sf, 1800*sf), xlim = c(1,today+16)) 


# scaling factor for fixed elements on chart
sf = 1000 / 1800

# Plot for animation
plot_au_video <- ggplot(data=au_covid, aes(x=day, y=new)) + 
	xlab("Days from Jun 1") + 
	ylab("Confirmed COVID-19 cases (Australia)") + 
	# xlim(1,today+14) +
	# ylim(-250*sf,1800*sf) +

	geom_vline(xintercept=39, color = "red", size=1) +
	annotate("text", x=38.8, y=750, label="Lockdown starts (9/7)", hjust=1) +
	geom_vline(xintercept=today, color = "grey") +
    	geom_vline(xintercept=43, color = "grey") +
	geom_vline(xintercept=50, color = "grey") +
    	geom_vline(xintercept=57, color = "grey") +
	geom_vline(xintercept=64, color = "grey") +
	geom_vline(xintercept=71, color = "grey") +
	geom_vline(xintercept=78, color = "grey") +
	geom_vline(xintercept=85, color = "grey") +
	geom_vline(xintercept=92, color = "grey") +

    geom_hline(yintercept=-150*sf) + 
	annotate("text", x=1, y=-200*sf, label="Jun 1") +
	annotate("text", x=15, y=-200*sf, label="Jun 15") +
	annotate("text", x=30+1, y=-200*sf, label="Jul 1") +
	annotate("text", x=30+15, y=-200*sf, label="Jul 15") +   
	annotate("text", x=30+31+1, y=-200*sf, label="Aug 1") +
	annotate("text", x=30+31+15, y=-200*sf, label="Aug 15") +
	annotate("text", x=30+31+31+1, y=-200*sf, label="Sep 1") +

	annotate("label", x=5, y=1800*sf, label="Projection of new daily cases of COVID-19", size=5.5, hjust=0) +
	annotate("text", x=5, y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) +
	annotate("text", x=5, y=1625*sf, label="Assumption: fitted to Gompertz model", hjust=0) +

	geom_ribbon(aes(ymin=model.au.lwr.deriv, ymax=model.au.upr.deriv), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.05, ymax=model.au.upr.deriv.05), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.15, ymax=model.au.upr.deriv.15), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.25, ymax=model.au.upr.deriv.25), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.35, ymax=model.au.upr.deriv.35), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.45, ymax=model.au.upr.deriv.45), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.55, ymax=model.au.upr.deriv.55), fill="blue", colour=NA, alpha=alphaPI) +
    geom_ribbon(aes(ymin=model.au.lwr.deriv.65, ymax=model.au.upr.deriv.65), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.75, ymax=model.au.upr.deriv.75), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.85, ymax=model.au.upr.deriv.85), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model.au.lwr.deriv.99, ymax=model.au.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.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=43, y=model.au.fit.deriv[43], label=round(model.au.fit.deriv[43]), size=3) +    
	# annotate("text", x=43, y=model.au.upr.deriv[43], label=round(model.au.upr.deriv[43]), size=3) +  
	# annotate("text", x=43, y=model.au.lwr.deriv[43], label=round(model.au.lwr.deriv[43]), size=3) +  

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

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

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

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

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

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

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

    annotate("label", x=(today), y=-85*sf, label=au_covid$date[today], size=3) +
	annotate("label", x=43, y=-150*sf, label=au_covid$date[43], size=2.5) +
	annotate("label", x=50, y=-150*sf, label=au_covid$date[50], size=2.5) +
	annotate("label", x=57, y=-150*sf, label=au_covid$date[57], size=2.5) +
	annotate("label", x=64, y=-150*sf, label=au_covid$date[64], size=2.5) +
	annotate("label", x=71, y=-150*sf, label=au_covid$date[71], size=2.5) +
	annotate("label", x=78, y=-150*sf, label=au_covid$date[78], size=2.5) +
	annotate("label", x=85, y=-150*sf, label=au_covid$date[85], size=2.5) +
	annotate("label", x=92, y=-150*sf, label=au_covid$date[92], size=2.5) +
	coord_cartesian(ylim = c(-250*sf, 1800*sf), xlim = c(1,today+14))


# Model outputs

plot_au
plot_au_new
plot_au_video

plot_au_video + 
	annotate("label", x=5, y=1000*sf, label="INTERPRET WITH EXTREME CAUTION", colour="red", hjust=0, size=5) +
	coord_cartesian(ylim = c(-250*sf, 1800*sf), xlim = c(1,93)) 

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