# Code for Model 1, Model 2, Model 3 and the Linear plot
# Michael Tam m.tam@unsw.edu.au
# I'm a newb at R. For the ease of adding new data, I do it manually into an Excel spreadsheet.
# I use that spreadsheet as the data source for the below code.
# The data can be found at facebook.com/vitualis or just email me.

# Step 0 - load libraries

library(drc)
library(aomisc)
library(ggplot2)
library(splines)



# Step 1 - define the "today" and current "time" variable.  1 = March 1, 2 = March 2, 32 = April 1, etc.

today = 33
time = 1230
chart_day = today - 9
index_day = chart_day + 45
all_days = -44:83


# Model 3 - Gompertz model (4 parameters) from 25 January to 2 April 2020
# 
# 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/

model3.cases <- COVID_AU_world$cases[1:index_day]
model3.day <- COVID_AU_world$day[1:index_day]
model3.gompertz <- drm(model3.cases ~ model3.day, fct = G.4())
model3.summary <- summary(model3.gompertz)

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

model3 <- data.frame(day = all_days, fit = model3.pred$Prediction, lwr = model3.pred$Lower, upr = model3.pred$Upper)
model3$lwr[model3$lwr < 5] <- 5



# Step 3 - Plots

# Australian COVID-19 cases projection vs comparators (Model 3)
# Projections are reported for up to 14 days.

#   scaling for fixed chart elements of y-axis in the log-plots
    sf = model3$upr[index_day + 16] * 100 / 50000

ggplot(data=COVID_AU_world, aes(x=day, y=cases)) + 
    xlab("Days from 100 confirmed cases") + 
    ylab("Confirmed COVID-19 cases [log scale]") + 
	geom_vline(xintercept=22.5, color="white", size=2) +
	geom_hline(yintercept=100, color="white", size=5) +
	geom_hline(yintercept=1000, color="white", size=5) +
	geom_hline(yintercept=10000, color="white", size=5) +
	geom_hline(yintercept=100000, color="white", size=5) +
	annotate("text", x=(chart_day + 19), y=100, label="100") +
	annotate("text", x=(chart_day + 19), y=1000, label="1,000") +
	annotate("text", x=(chart_day + 19), y=10000, label="10,000") +
	annotate("text", x=(chart_day + 19), y=100000, label="100,000") +
	geom_hline(yintercept=10, color="white", size=5) +
	geom_hline(yintercept=10) +
	geom_segment(aes(x=-8, y=9, xend = -8, yend = 11)) +
	geom_segment(aes(x=-1, y=9, xend = -1, yend = 11)) +
	geom_segment(aes(x=6, y=9, xend = 6, yend = 11)) +
	geom_segment(aes(x=13, y=9, xend = 13, yend = 11)) +
	geom_segment(aes(x=20, y=9, xend = 20, yend = 11)) +
	geom_segment(aes(x=27, y=9, xend = 27, yend = 11)) +
	geom_segment(aes(x=34, y=9, xend = 34, yend = 11)) +
	geom_segment(aes(x=41, y=9, xend = 41, yend = 11)) +
	annotate("text", x=-8, y=6, label="Mar 1") +
	annotate("text", x=-1, y=6, label="Mar 8") +
	annotate("text", x=6, y=6, label="Mar 15") +
	annotate("text", x=13, y=6, label="Mar 22") +
	annotate("text", x=20, y=6, label="Mar 29") +
	annotate("text", x=27, y=6, label="Apr 5") +
	annotate("text", x=34, y=6, label="Apr 12") +
	annotate("text", x=41, y=6, label="Apr 19") +
	annotate("text", x=48, y=6, label="Apr 26") +
	annotate("text", x=55, y=6, label="May 3") +
	annotate("text", x=62, y=6, label="May 10") +
	geom_segment(aes(x=48, y=9, xend = 48, yend = 11)) +
	geom_segment(aes(x=55, y=9, xend = 55, yend = 11)) +
	geom_segment(aes(x=62, y=9, xend = 62, yend = 11)) +

	annotate("label", x=-4, y=10000*sf, label="Australia", colour="blue", size=4) +
	annotate("label", x=-4, y=5000*sf, label="South Korea", colour="aquamarine4") +
	annotate("label", x=-4, y=2500*sf, label="Italy", colour="red") +
	annotate("label", x=-4, y=1250*sf, label="US", colour="purple") +
	annotate("label", x=-4, y=650*sf, label="UK", colour="coral") +
	annotate("label", x=-4, y=375*sf, label="Norway", colour="orange") +
	annotate("label", x=-4, y=187*sf, label="Singapore", colour="green4") +
	annotate("label", x=-4, y=94*sf, label="New Zealand", colour="grey20") +     
	scale_y_log10(limits = c(5,50000*sf), labels = scales::comma) +
	annotation_logticks(sides = "lr") +
    
    geom_point(aes(y=sing_cases), colour="green4", alpha=0.5) + 
    geom_line(aes(y=sing_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="green1") + 
    geom_point(aes(y=ita_cases), colour="red", alpha=0.5) + 
	geom_line(aes(y=ita_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="red1") + 
    geom_point(aes(y=nor_cases), colour="orange4", alpha=0.5) + 
	geom_line(aes(y=nor_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="orange1") + 
    geom_point(aes(y=SK_cases), colour="aquamarine4", alpha=0.5) + 
	geom_line(aes(y=SK_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="aquamarine1") + 
    geom_point(aes(y=UK_cases), colour="coral4", alpha=0.5) + 
	geom_line(aes(y=UK_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="coral1") + 
    geom_point(aes(y=US_cases), colour="purple4", alpha=0.5) + 
	geom_line(aes(y=US_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="purple1") + 
    geom_point(aes(y=nz_cases), colour="grey20", alpha=0.5) + 
	geom_line(aes(y=nz_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="grey20") + 

	geom_line(aes(y=model3$fit), colour="blue", size=1, alpha=0.5, linetype=2) + 
    geom_ribbon(aes(ymin=model3$lwr, ymax=model3$upr), colour=NA, fill="blue", alpha=0.1) +
    
    geom_point(size=3, colour="blue", alpha=0.6) + 
	geom_segment(aes(x=1, y=10, xend = 1, yend = 100), color="grey") +
	annotate("text", x=1, y=15, label="D1: 10/3/20") +
	annotate("text", x=1, y=20, label="100", colour="blue") +

	annotate("label", x=-8, y=40000*sf, label="Projection of Australian COVID-19 cases (Model 3)", size=5.5, hjust=0) +
	annotate("text", x=-8, y=25000*sf, label= paste(COVID_AU_world$date[index_day], "@", time, "GMT+10"), hjust=0) +
	annotate("text", x=-8, y=18000*sf, label="Assumption: trend based on Gompertz model", hjust=0) +
	geom_segment(aes(x=chart_day, y=10, xend = chart_day, yend = COVID_AU_world$cases[index_day]), color="grey") +
	annotate("text", x=chart_day, y=15, label= paste("D", chart_day, ": ", COVID_AU_world$date[index_day], sep="", size=3)) +
	annotate("text", x=chart_day, y=20, label=COVID_AU_world$cases[index_day], colour="blue") +
	geom_segment(aes(x=(chart_day + 14), y=10, xend = (chart_day + 14), yend = model3$fit[index_day + 14]), color="grey") +
	annotate("text", x=(chart_day + 14), y=15, label= paste("D", chart_day + 14, ": ", COVID_AU_world$date[index_day+14], " 2 weeks", sep=""), size=3) +
	annotate("label", x=(chart_day + 14), y=60, label= paste("~",round(model3$fit[index_day + 14]), " cases", sep=""), size=5.5, colour="blue") +
	annotate("text", x=(chart_day + 14), y=40, label="95% prediction interval:", colour="blue") +
	annotate("text", x=(chart_day + 14), y=30, label=paste(round(model3$lwr[index_day + 14]), "to", round(model3$upr[index_day + 14])), colour="blue") +
	xlim(-8,chart_day + 20) -> plot_model3


# Australian COVID-19 cases projection - linear plot of Model 3
#   scaling for fixed chart elements of y-axis of model 3
    sf = model3$upr[index_day + 14] * 1.5 / 30000
    
ggplot(data=COVID_AU_world, aes(x=day, y=cases)) + 
    xlab("Days from 100 confirmed cases") + 
    ylab("Confirmed COVID-19 cases [linear scale]") + 
	ylim(-4000*sf,30000*sf) +
	xlim(-8,chart_day + 16) +
	geom_hline(yintercept=0) +
	geom_hline(yintercept=-2000*sf, color="white", size=5) +
	geom_hline(yintercept=-2000*sf) +
	geom_segment(aes(x=-8, y=-2200*sf, xend = -8, yend = -1800*sf)) +
	geom_segment(aes(x=-1, y=-2200*sf, xend = -1, yend = -1800*sf)) +
	geom_segment(aes(x=6, y=-2200*sf, xend = 6, yend = -1800*sf)) +
	geom_segment(aes(x=13, y=-2200*sf, xend = 13, yend = -1800*sf)) +
	geom_segment(aes(x=20, y=-2200*sf, xend = 20, yend = -1800*sf)) +
	geom_segment(aes(x=27, y=-2200*sf, xend = 27, yend = -1800*sf)) +
	geom_segment(aes(x=34, y=-2200*sf, xend = 34, yend = -1800*sf)) +
	geom_segment(aes(x=41, y=-2200*sf, xend = 41, yend = -1800*sf)) +
	geom_segment(aes(x=48, y=-2200*sf, xend = 48, yend = -1800*sf)) +
	geom_segment(aes(x=55, y=-2200*sf, xend = 55, yend = -1800*sf)) +
	geom_segment(aes(x=62, y=-2200*sf, xend = 62, yend = -1800*sf)) +

	annotate("text", x=-8, y=-3000*sf, label="Mar 1") +
	annotate("text", x=-1, y=-3000*sf, label="Mar 8") +
	annotate("text", x=6, y=-3000*sf, label="Mar 15") +
	annotate("text", x=13, y=-3000*sf, label="Mar 22") +
	annotate("text", x=20, y=-3000*sf, label="Mar 29") +
	annotate("text", x=27, y=-3000*sf, label="Apr 5") +
	annotate("text", x=34, y=-3000*sf, label="Apr 12") +
	annotate("text", x=41, y=-3000*sf, label="Apr 19") +
	annotate("text", x=48, y=-3000*sf, label="Apr 26") +
	annotate("text", x=55, y=-3000*sf, label="May 3") +
	annotate("text", x=62, y=-3000*sf, label="May 10") +


	annotate("label", x=-3, y=18000*sf, label="Australia", colour="blue", size=4) +
	annotate("label", x=-3, y=16000*sf, label="South Korea", colour="aquamarine4") +
	annotate("label", x=-3, y=14000*sf, label="Italy", colour="red") +
	annotate("label", x=-3, y=12000*sf, label="US", colour="purple") +
	annotate("label", x=-3, y=10000*sf, label="UK", colour="coral") +
	annotate("label", x=-3, y=8000*sf, label="Norway", colour="orange") +
	annotate("label", x=-3, y=6000*sf, label="Singapore", colour="green4") +
	annotate("label", x=-3, y=3000*sf, label="New Zealand", colour="grey20") +   
 
    geom_point(aes(y=sing_cases), colour="green4", alpha=0.5) + 
    geom_line(aes(y=sing_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="green1") + 
    geom_point(aes(y=ita_cases), colour="red", alpha=0.5) + 
	geom_line(aes(y=ita_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="red1") + 
    geom_point(aes(y=nor_cases), colour="orange4", alpha=0.5) + 
	geom_line(aes(y=nor_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="orange1") + 
    geom_point(aes(y=SK_cases), colour="aquamarine4", alpha=0.5) + 
	geom_line(aes(y=SK_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="aquamarine1") + 
    geom_point(aes(y=UK_cases), colour="coral4", alpha=0.5) + 
	geom_line(aes(y=UK_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="coral1") + 
    geom_point(aes(y=US_cases), colour="purple4", alpha=0.5) + 
	geom_line(aes(y=US_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="purple1") + 
    geom_point(aes(y=nz_cases), colour="grey20", alpha=0.5) + 
	geom_line(aes(y=nz_cases), stat="smooth", alpha=0.7, arrow = arrow(type="closed", angle=15), colour="grey20") + 

	geom_line(aes(y=model3$fit), colour="blue", size=1, alpha=0.5, linetype=2) + 
    geom_ribbon(aes(ymin=model3$lwr, ymax=model3$upr), colour=NA, fill="blue", alpha=0.1) +
    
    geom_point(size=3, colour="blue", alpha=0.6) + 
	geom_segment(aes(x=1, y=-2000*sf, xend = 1, yend = 100), color="grey") +
	annotate("text", x=1, y=-1100*sf, label="D1: 10/3") +
	annotate("label", x=-8, y=30000*sf, label="Projection of Australian COVID-19 cases (Model 3)", size=5.5, hjust=0) +
	annotate("text", x=-8, y=28000*sf, label="Assumption: trend based on Gompertz model", hjust=0) +
	annotate("text", x=-8, y=26000*sf, label= paste(COVID_AU_world$date[index_day], "@", time, " GMT+10"), hjust=0) +
	geom_segment(aes(x=chart_day, y=-2000*sf, xend = chart_day, yend = COVID_AU_world$cases[index_day]), color="grey") +
	annotate("label", x=chart_day, y=-1100*sf, label=paste(COVID_AU_world$date[index_day],": ", COVID_AU_world$cases[index_day], sep=""), colour="blue", size=3) +

	geom_segment(aes(x=(chart_day + 14), y=-2000*sf, xend = (chart_day + 14), yend = model3$fit[index_day + 14]), color="grey") +
	annotate("label", x=(chart_day + 14), y=4000*sf, label= paste("~",round(model3$fit[index_day + 14]), " cases", sep=""), size=5.5, colour="blue") +
	annotate("text", x=(chart_day + 14), y=2600*sf, label="95% prediction interval:", colour="blue") +
	annotate("text", x=(chart_day + 14), y=1600*sf, label=paste(round(model3$lwr[index_day + 14]), "to", round(model3$upr[index_day + 14])), colour="blue") +

	annotate("label", x=(chart_day + 14), y=-1100*sf, label= paste(COVID_AU_world$date[index_day+14], ": 2 weeks", sep=""), colour="blue", size=3) -> plot_model3_linear
    


### Plot - Projection of new cases from Model 3

prediction_days <- data.frame(model3.day = all_days)
model3.pred <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction'))
model3 <- data.frame(day = all_days, fit = model3.pred$Prediction, lwr = model3.pred$Lower, upr = model3.pred$Upper)
model3$lwr[model3$lwr < 0] <- 0

model3.pred.05 <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction', level=0.05))
model3.05 <- data.frame(lwr = model3.pred.05$Lower, upr = model3.pred.05$Upper)
model3.05$lwr[model3.05$lwr < 0] <- 0

model3.pred.15 <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction', level=0.15))
model3.15 <- data.frame(lwr = model3.pred.15$Lower, upr = model3.pred.15$Upper)
model3.15$lwr[model3.15$lwr < 0] <- 0

model3.pred.25 <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction', level=0.25))
model3.25 <- data.frame(lwr = model3.pred.25$Lower, upr = model3.pred.25$Upper)
model3.25$lwr[model3.25$lwr < 0] <- 0

model3.pred.35 <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction', level=0.35))
model3.35 <- data.frame(lwr = model3.pred.35$Lower, upr = model3.pred.35$Upper)
model3.35$lwr[model3.35$lwr < 0] <- 0

model3.pred.45 <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction', level=0.45))
model3.45 <- data.frame(lwr = model3.pred.45$Lower, upr = model3.pred.45$Upper)
model3.45$lwr[model3.45$lwr < 0] <- 0

model3.pred.55 <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction', level=0.55))
model3.55 <- data.frame(lwr = model3.pred.55$Lower, upr = model3.pred.55$Upper)
model3.55$lwr[model3.55$lwr < 0] <- 0

model3.pred.65 <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction', level=0.65))
model3.65 <- data.frame(lwr = model3.pred.65$Lower, upr = model3.pred.65$Upper)
model3.65$lwr[model3.65$lwr < 0] <- 0

model3.pred.75 <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction', level=0.75))
model3.75 <- data.frame(lwr = model3.pred.75$Lower, upr = model3.pred.75$Upper)
model3.75$lwr[model3.75$lwr < 0] <- 0

model3.pred.85 <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction', level=0.85))
model3.85 <- data.frame(lwr = model3.pred.85$Lower, upr = model3.pred.85$Upper)
model3.85$lwr[model3.85$lwr < 0] <- 0

model3.pred.99 <- data.frame(predict(model3.gompertz, newdata = prediction_days, interval = 'prediction', level=0.99))
model3.99 <- data.frame(lwr = model3.pred.99$Lower, upr = model3.pred.99$Upper)
model3.99$lwr[model3.99$lwr < 0] <- 0


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


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

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

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

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

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

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

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

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

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

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

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



# scaling factor for fixed elements on chart
sf = 700 / 1800
alphaPI = 0.08

plot_model3_new <- ggplot(data=COVID_AU_world, aes(x=day, y=new)) + 
	xlab("Days from 100 confirmed cases") + 
	ylab("New daily COVID-19 cases") + 
	# xlim(-8,chart_day + 16) +
	# ylim(-250*sf,1800*sf) +

	geom_vline(xintercept=chart_day, color = "blue", size=1) +
	annotate("text", x=(chart_day - 0.2), y=550, label="+10 days (2/4)", hjust=1) +

	geom_vline(xintercept=14, color = "red", size=1) +
	annotate("text", x=13.8, y=550, label="Lockdown starts (23/3)", hjust=1) +

    geom_vline(xintercept=(chart_day + 7), color = "grey") +
    geom_vline(xintercept=(chart_day + 14), color = "grey") +

    geom_hline(yintercept=-150*sf) + 
	annotate("text", x=-37, y=-200*sf, label="Feb 1") +
	annotate("text", x=-23, y=-200*sf, label="Feb 15") +   
	annotate("text", x=-8, y=-200*sf, label="Mar 1") +
	annotate("text", x=6, y=-200*sf, label="Mar 15") +
	annotate("text", x=23, y=-200*sf, label="Apr 1") +
	annotate("text", x=37, y=-200*sf, label="Apr 15") +
	annotate("text", x=53, y=-200*sf, label="May 1") +
	annotate("text", x=67, y=-200*sf, label="May 15") +
    
	geom_ribbon(aes(ymin=model3.lwr.deriv, ymax=model3.upr.deriv), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model3.lwr.deriv.05, ymax=model3.upr.deriv.05), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model3.lwr.deriv.15, ymax=model3.upr.deriv.15), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model3.lwr.deriv.25, ymax=model3.upr.deriv.25), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model3.lwr.deriv.35, ymax=model3.upr.deriv.35), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model3.lwr.deriv.45, ymax=model3.upr.deriv.45), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model3.lwr.deriv.55, ymax=model3.upr.deriv.55), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model3.lwr.deriv.65, ymax=model3.upr.deriv.65), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model3.lwr.deriv.75, ymax=model3.upr.deriv.75), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model3.lwr.deriv.85, ymax=model3.upr.deriv.85), fill="blue", colour=NA, alpha=alphaPI) +
	geom_ribbon(aes(ymin=model3.lwr.deriv.99, ymax=model3.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=model3.fit.deriv), colour="grey20", alpha=0.3, size = 1, linetype=2) + 
    geom_line(aes(y=model3.upr.deriv), colour="grey20", alpha=0.3, linetype=2) + 
    geom_line(aes(y=model3.lwr.deriv), colour="grey20", alpha=0.3, linetype=2) + 

	annotate("label", x=(chart_day + 7), y=model3.fit.deriv[index_day + 7], label=round(model3.fit.deriv[index_day + 7]), size=3) +    
	annotate("text", x=(chart_day + 7), y=model3.upr.deriv[index_day + 7], label=round(model3.upr.deriv[index_day + 7]), size=3) +  
	# annotate("text", x=(chart_day + 7), y=model3.lwr.deriv[index_day + 7], label=round(model3.lwr.deriv[index_day + 7]), size=3) +  

	annotate("label", x=(chart_day + 14), y=model3.fit.deriv[index_day + 14], label=round(model3.fit.deriv[index_day + 14]), size=3) +    
	annotate("text", x=(chart_day + 14), y=model3.upr.deriv[index_day + 14], label=round(model3.upr.deriv[index_day + 14]), size=3) +  
	# annotate("text", x=(chart_day + 14), y=model3.lwr.deriv[index_day + 14], label=round(model3.lwr.deriv[index_day + 14]), size=3) +

    annotate("label", x=(chart_day), y=-85*sf, label=COVID_AU_world$date[index_day], size=3) +
	annotate("label", x=(chart_day + 7), y=-85*sf, label=COVID_AU_world$date[index_day+7], size=3) +
	annotate("label", x=(chart_day + 14), y=-85*sf, label=COVID_AU_world$date[index_day+14], size=3) +

	annotate("label", x=-8, y=1800*sf, label="Australian COVID-19 projection based on data from 2 April 2020", size=5.5, hjust=0) +
	annotate("text", x=-8, y=1700*sf, label= "2020-07-19 @ 2130 GMT+10", hjust=0) +
	annotate("text", x=-8, y=1625*sf, label="Fitted to Gompertz model from data up to 2/4/2020", hjust=0) +
	coord_cartesian(ylim = c(-250*sf, 1800*sf), xlim = c(-8,chart_day + 16)) 



# Step 4 - Output summaries and stats into the console      
#

    plot_model3
    plot_model3_linear
    plot_model3_new

    model3.summary
    model3