# Code - Dr Michael Tam # m.tam@unsw.edu.au # Projections using NLR to Gompertz and Richards' growth curve of Australian (VIC) COVID-19 # outbreak in August 2021 (the "sixth" Melbourne lock down). # The data excludes cases known to have been acquired from overseas (for instance, detected # in hotel quarantine) # # I've released this code into the public domain. If you make use of this code, attribution # is welcome but not necessary. # # This version of the code was written on 13 October 2021. The original analysis undertook # the modelling starting from 4 August 2021, which was the day of the first case of the # outbreak. This new model now only starts from the day of a number of large anti-lock down # protests in Melbourne in September 2021. The reason for this is that after this time, there has # been a clear escalation in cases and a change in the trajectory. The case numbers and trajectory # have lesser direct mechanistic impact on the current and future numbers, and if anything, provide # a bias given the assumptions inherent in modelling using regression to the epidemiologic growth # curves. # Establish the date and time the code is run. # This sets up some of the date dependent variables in the analysis, and also provides the # date and time labels in the charts. Please update this daily so that the analyses run # correctly and so that the charts are date and timestamped. todaydate = as.Date("2021-10-22") time = 1000 # Install libraries # These are the libraries used in this code that are not part of base R. You'll need to # install this the first time you run the code. # To do this, remove the hash symbol to take it off the comments, and copy-paste into the R # console to run. Once these packages are installed, you won't need to do it again. # install.packages("drc") # install.packages("ggplot2") # install.packages("tidyquant") # install.packages("gridExtra") # install.packages("readxl") # install.packages("scales") # install.packages("remotes") # install.packages("DescTools") # remotes::install_github("OnofriAndreaPG/aomisc") # Load libraries library(drc) library(aomisc) library(ggplot2) library(tidyquant) library(gridExtra) library(readxl) library(scales) library(DescTools) # Loading the data # I'm assuming that the data file for today is an Excel spreadsheet named # "au_covid_vic.xlsx". # If it isn't, rename it before uploading it into RStudio. The spreadsheet basically # contains only five columns - Date, Day number (1 = 4 August 2021), Number of new local # cases of COVID-19 announced by Vic Health, and the Cumulative number of # local cases of COVID-19 in the Vic outbreak starting from 4 August 2021. The final # column contains potential comments and is not used in the analysis. au_covid <- read_excel("au_covid_vic.xlsx", col_types = c("date", "numeric", "numeric", "numeric", "text")) # View(au_covid) # Setting up some variables. # The day number starts at "1", and this is 4 August 2021. This was the date prior to the # reported local case in the Melbourne August 2021 outbreak. # all_days is simply a list of integers from 1 to 150. # Note: Day 150 is 31 December 2021. today = as.numeric(todaydate - as.Date("2021-8-4")) + 1 all_days = 1:150 # On 13 October 2021, I changed the models to only use the cumulative case data starting # on 20 August 2021 (Day 48). To not have to rewrite large sections of the code, I # created a new vector "tot2" that contains the cumulative totals from day 1 to "today", # and modified "tot" to only contain the cumulative totals from day 48 to "today". # Same with "new" and "new2". au_covid$tot2 <- au_covid$tot au_covid$tot[1:47] <- NA au_covid$new2 <- au_covid$new au_covid$new[1:47] <- NA # This analysis critically uses the following packages: # "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 4 October 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:150] <- NA model.history.g.7fit[1:150] <- NA model.history.g.7lwr[1:150] <- NA model.history.g.7upr[1:150] <- NA model.history.g.14fit[1:150] <- NA model.history.g.14lwr[1:150] <- NA model.history.g.14upr[1:150] <- NA model.history.g.7fit.error[1:150] <- NA model.history.g.7lwr.error[1:150] <- NA model.history.g.7upr.error[1:150] <- NA model.history.g.14fit.error[1:150] <- NA model.history.g.14lwr.error[1:150] <- NA model.history.g.14upr.error[1:150] <- 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:150] <- NA model.history.r.7fit[1:150] <- NA model.history.r.7lwr[1:150] <- NA model.history.r.7upr[1:150] <- NA model.history.r.14fit[1:150] <- NA model.history.r.14lwr[1:150] <- NA model.history.r.14upr[1:150] <- NA model.history.r.7fit.error[1:150] <- NA model.history.r.7lwr.error[1:150] <- NA model.history.r.7upr.error[1:150] <- NA model.history.r.14fit.error[1:150] <- NA model.history.r.14lwr.error[1:150] <- NA model.history.r.14upr.error[1:150] <- NA real.7 <- vector() real.14 <- vector() real.7[1:150] <- NA real.14[1:150] <- NA for(i in 62:today) { model.cases <- au_covid$tot[48:i] model.day <- au_covid$day[48: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) # Projections of CUMULATIVE case counts based on Gompertz/Richards' model # The code here is a little ugly, but there is a reason for this. For some reason, when # using the predict() function with a gompertz model, asking it to include the standard # error of the estimate in the output, se=TRUE, results in zero values being returned for # the upper and lower bounds. I have no idea why. So basically, I use the predict() # function twice, once to get the estimate and confidence interval, and once to get the # standard error. I then create a new data frame to clean this up and get rid of the # duplicated estimate column, and the useless zeroed columns for the CI bounds. model.cases <- au_covid$tot[1:150] model.day <- au_covid$day[1:150] model.g <- drm(model.cases ~ model.day, fct = G.4()) model.r <- drm(model.cases ~ model.day, fct = L.5()) model.au.gompertz.pred <- data.frame(predict(model.g, newdata = prediction_days, interval = 'confidence'), predict(model.g, newdata = prediction_days, interval = 'confidence', se=TRUE)) 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, se = model.au.gompertz.pred$SE) # Remove the projections backwards in time model.gompertz$fit[1:47] <- NA model.gompertz$lwr[1:47] <- NA model.gompertz$upr[1:47] <- NA model.gompertz$se[1:47] <- NA model.au.pred <- data.frame(predict(model.r, newdata = prediction_days, interval = 'confidence'), predict(model.r, newdata = prediction_days, interval = 'confidence', se=TRUE)) model.au <- data.frame( day = all_days, fit = model.au.pred$Prediction, lwr = model.au.pred$Lower, upr = model.au.pred$Upper, se = model.au.pred$SE) model.au$fit[1:47] <- NA model.au$lwr[1:47] <- NA model.au$upr[1:47] <- NA model.au$se[1:47] <- NA ### Plots ## Plot cumulative cases - Richards' growth curve model # Scaling factor for fixed elements on chart # Note: the value of 1800 as equal to 100% of the height of the chart is an arbitrary # choice made originally when I was hardcoding a chart in 2020. There is no special # significance to this value other than it became too tedious to recode all the # charts. sf = model.au$upr[today+10] * 1.05 / 1800 plot_au <- ggplot(data=au_covid, aes(x=as.Date(date), y=tot2)) + xlab("") + ylab("Cumulative Local COVID-19 cases (Vic) starting August 4") + 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=2.5) + geom_vline(xintercept=todaydate + 7, color = "grey") + annotate("label", x=todaydate+7, y=50*sf, label=au_covid$date[today+7], size=2.5) + annotate("label", x=as.Date("2021-8-6"), y=1800*sf, label="Projection of total COVID-19 cases", size=5.5, hjust=0) + annotate("text", x=as.Date("2021-8-6"), y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) + annotate("text", x=as.Date("2021-8-6"), y=1625*sf, label="Assumption: fitted to Richards' growth curve", hjust=0) + geom_point(colour="grey50", size = 1) + geom_point(aes(y=tot), 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] * 1.04, label=round(model.au$upr[today + 7]), size=2) + annotate("text", x=todaydate + 7, y=model.au$lwr[today + 7] * 0.96, label=round(model.au$lwr[today + 7]), size=2) + coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-8-4"), todaydate + 10) ) ## Plot cumulative cases - Gompertz # scaling factor for fixed elements on chart sf = model.gompertz$upr[today+10] * 1.05 / 1800 plot_au_gompertz <- ggplot(data=au_covid, aes(x=as.Date(date), y=tot2)) + xlab("") + ylab("Cumulative Local COVID-19 cases (Vic) starting August 4") + 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=2.5) + geom_vline(xintercept=todaydate + 7, color = "grey") + annotate("label", x=todaydate+7, y=50*sf, label=au_covid$date[today+7], size=2.5) + annotate("label", x=as.Date("2021-8-6"), y=1800*sf, label="Projection of total COVID-19 cases", size=5.5, hjust=0) + annotate("text", x=as.Date("2021-8-6"), y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) + annotate("text", x=as.Date("2021-8-6"), y=1625*sf, label="Assumption: fitted to Gompertz growth curve", hjust=0) + geom_point(colour="grey50", size = 1) + geom_point(aes(y=tot), 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] * 1.04, label=round(model.gompertz$upr[today + 7]), size=2) + annotate("text", x=todaydate + 7, y=model.gompertz$lwr[today + 7] * 0.94, label=round(model.gompertz$lwr[today + 7]), size=2) + coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-8-4"), todaydate + 10) ) ### Plots - new cases # Projections of NEW daily case counts based on Gompertz/Richards' model # Note: This is effectively the first derivative of the model. # # I've changed the way I compute the 95% confidence intervals on several occasions. # Initially, I couldn't think of how to easily compute this, so I resorted to a hack. # I basically used the confidence interval of the estimate from the model of the # cumulative cases as the estimate for the CI of the derivative. This would obviously # be an overestimate. # # On 7 August 2021, I changed the method to be more statistically defensible. # For each day's TOTAL case number projection in the fitted Gompertz/Richards' model, I use # the standard error of that fitted value to create a randomised normal distribution of # values. I compute a list of 1,000,000 values for each day. I then subtract a list of # these values from the list of the prior day, effectively creating a list of differences # between the days. This distribution of values will be overly conservative as in reality # it is not possible for the total number of cases to ever be less than a previous day's. # # On 9 September 2021, I changed the method once again, which I think is finally # acceptable. Both the Gompertz and Richards' models give standard errors for each of the # parameters. I use this to create a randomised normal distribution of these parameters, # and then use these parameters to simulate new curves in a bootstrapping process (a run of # 10,000 times). From each of these simulated growth curves, I extract the daily growth in # new cases. For each day, the median, or middle value of this list of differences is the # NEW case estimate as directly computed from the model. Specific percentiles, which are # the proportion of values of the distribution can be interpreted as the confidence # interval. For instance, the 95% confidence interval is the interval between the 2.5% and # 97.5% centiles. ## # Initial set up of the vectors model.gompertz.fit.deriv <- vector() model.gompertz.lwr.deriv <- vector() model.gompertz.upr.deriv <- vector() model.gompertz.lwr.deriv.05 <- vector() model.gompertz.upr.deriv.05 <- vector() model.gompertz.lwr.deriv.15 <- vector() model.gompertz.upr.deriv.15 <- vector() model.gompertz.lwr.deriv.25 <- vector() model.gompertz.upr.deriv.25 <- vector() model.gompertz.lwr.deriv.35 <- vector() model.gompertz.upr.deriv.35 <- vector() model.gompertz.lwr.deriv.45 <- vector() model.gompertz.upr.deriv.45 <- vector() model.gompertz.lwr.deriv.55 <- vector() model.gompertz.upr.deriv.55 <- vector() model.gompertz.lwr.deriv.65 <- vector() model.gompertz.upr.deriv.65 <- vector() model.gompertz.lwr.deriv.75 <- vector() model.gompertz.upr.deriv.75 <- vector() model.gompertz.lwr.deriv.85 <- vector() model.gompertz.upr.deriv.85 <- vector() model.gompertz.lwr.deriv.99 <- vector() model.gompertz.upr.deriv.99 <- vector() model.gompertz.fit.deriv[1:150] <- NA model.gompertz.lwr.deriv[1:150] <- NA model.gompertz.upr.deriv[1:150] <- NA model.gompertz.lwr.deriv.05[1:150] <- NA model.gompertz.upr.deriv.05[1:150] <- NA model.gompertz.lwr.deriv.15[1:150] <- NA model.gompertz.upr.deriv.15[1:150] <- NA model.gompertz.lwr.deriv.25[1:150] <- NA model.gompertz.upr.deriv.25[1:150] <- NA model.gompertz.lwr.deriv.35[1:150] <- NA model.gompertz.upr.deriv.35[1:150] <- NA model.gompertz.lwr.deriv.45[1:150] <- NA model.gompertz.upr.deriv.45[1:150] <- NA model.gompertz.lwr.deriv.55[1:150] <- NA model.gompertz.upr.deriv.55[1:150] <- NA model.gompertz.lwr.deriv.65[1:150] <- NA model.gompertz.upr.deriv.65[1:150] <- NA model.gompertz.lwr.deriv.75[1:150] <- NA model.gompertz.upr.deriv.75[1:150] <- NA model.gompertz.lwr.deriv.85[1:150] <- NA model.gompertz.upr.deriv.85[1:150] <- NA model.gompertz.lwr.deriv.99[1:150] <- NA model.gompertz.upr.deriv.99[1:150] <- NA model.au.fit.deriv <- vector() model.au.lwr.deriv <- vector() model.au.upr.deriv <- vector() model.au.lwr.deriv.05 <- vector() model.au.upr.deriv.05 <- vector() model.au.lwr.deriv.15 <- vector() model.au.upr.deriv.15 <- vector() model.au.lwr.deriv.25 <- vector() model.au.upr.deriv.25 <- vector() model.au.lwr.deriv.35 <- vector() model.au.upr.deriv.35 <- vector() model.au.lwr.deriv.45 <- vector() model.au.upr.deriv.45 <- vector() model.au.lwr.deriv.55 <- vector() model.au.upr.deriv.55 <- vector() model.au.lwr.deriv.65 <- vector() model.au.upr.deriv.65 <- vector() model.au.lwr.deriv.75 <- vector() model.au.upr.deriv.75 <- vector() model.au.lwr.deriv.85 <- vector() model.au.upr.deriv.85 <- vector() model.au.lwr.deriv.99 <- vector() model.au.upr.deriv.99 <- vector() model.au.fit.deriv[1:150] <- NA model.au.lwr.deriv[1:150] <- NA model.au.upr.deriv[1:150] <- NA model.au.lwr.deriv.05[1:150] <- NA model.au.upr.deriv.05[1:150] <- NA model.au.lwr.deriv.15[1:150] <- NA model.au.upr.deriv.15[1:150] <- NA model.au.lwr.deriv.25[1:150] <- NA model.au.upr.deriv.25[1:150] <- NA model.au.lwr.deriv.35[1:150] <- NA model.au.upr.deriv.35[1:150] <- NA model.au.lwr.deriv.45[1:150] <- NA model.au.upr.deriv.45[1:150] <- NA model.au.lwr.deriv.55[1:150] <- NA model.au.upr.deriv.55[1:150] <- NA model.au.lwr.deriv.65[1:150] <- NA model.au.upr.deriv.65[1:150] <- NA model.au.lwr.deriv.75[1:150] <- NA model.au.upr.deriv.75[1:150] <- NA model.au.lwr.deriv.85[1:150] <- NA model.au.upr.deriv.85[1:150] <- NA model.au.lwr.deriv.99[1:150] <- NA model.au.upr.deriv.99[1:150] <- NA # For the Gompertz model # We calculate the estimate and the 95% CI # Further 5%, 15%, 25%, 35%, 45%, 55%, 65%, 75%, 85%, 99% CIs are calculated. # These are used to create the coloured confidence bands in the charts. # This will be the data frame where the derivative data from each simulation will be # placed. # Create the initial empty data frame. model.g.deriv <- data.frame(day = double(), deriv = double()) # Create the matrix of the variations of the parameters from the Gompertz model that # we'll use for bootstrapping. Effectively, we use the standard error for each of the # coefficients in the model to create a normal distribution of each parameter. model.g.bootstrapcoef <- data.frame( b = rnorm(10000, mean=model.g$coefficients[1], sd=summary(model.g)$coefficients["b:(Intercept)","Std. Error"]), c = rnorm(10000, mean=model.g$coefficients[2], sd=summary(model.g)$coefficients["c:(Intercept)","Std. Error"]), d = rnorm(10000, mean=model.g$coefficients[3], sd=summary(model.g)$coefficients["d:(Intercept)","Std. Error"]), e = rnorm(10000, mean=model.g$coefficients[4], sd=summary(model.g)$coefficients["e:(Intercept)","Std. Error"])) # Start of the bootstrapping. # 10,000 simulations are run. for(i in 1:10000) { # Set up the model parameters for this instance of the simulation b = model.g.bootstrapcoef$b[i] c = model.g.bootstrapcoef$c[i] d = model.g.bootstrapcoef$d[i] e = model.g.bootstrapcoef$e[i] model.g.bootstrap.deriv <- vector() model.g.bootstrap.deriv[1:150] <- NA for(j in 49:150) { # The Gompertz equation is in the following form: # f(x) = c + (d-c) * ( exp( -exp( b * (x - e)))) # # Where x = day number # b, c, d, e are the parameters of the model # The calculation below uses the model to calculate the cumulative cases for the # current day, and subtracts it from that from the previous day. model.g.bootstrap.deriv[j] <- ( c + (d-c) * ( exp( -exp( b * (j - e)))) ) - ( c + (d-c) * ( exp( -exp( b * ((j-1) - e)))) ) } # Append the derivative data from this run of the simulation to those of all the previous # runs. bootstrap.deriv.temp <- data.frame(day = all_days, deriv = model.g.bootstrap.deriv) model.g.deriv <- data.frame( days = c(model.g.deriv$days, bootstrap.deriv.temp$day), deriv = c(model.g.deriv$deriv, bootstrap.deriv.temp$deriv)) } # Compute the confidence bands from the derivative data. for(i in 49:150) { day_diff <- model.g.deriv$deriv[model.g.deriv$days == i] model.gompertz.fit.deriv[i] <- median(day_diff) model.gompertz.lwr.deriv[i] <- quantile(day_diff, 0.025) model.gompertz.upr.deriv[i] <- quantile(day_diff, 0.975) model.gompertz.lwr.deriv.05[i] <- quantile(day_diff, 0.475) model.gompertz.upr.deriv.05[i] <- quantile(day_diff, 0.525) model.gompertz.lwr.deriv.15[i] <- quantile(day_diff, 0.425) model.gompertz.upr.deriv.15[i] <- quantile(day_diff, 0.575) model.gompertz.lwr.deriv.25[i] <- quantile(day_diff, 0.375) model.gompertz.upr.deriv.25[i] <- quantile(day_diff, 0.625) model.gompertz.lwr.deriv.35[i] <- quantile(day_diff, 0.325) model.gompertz.upr.deriv.35[i] <- quantile(day_diff, 0.675) model.gompertz.lwr.deriv.45[i] <- quantile(day_diff, 0.275) model.gompertz.upr.deriv.45[i] <- quantile(day_diff, 0.725) model.gompertz.lwr.deriv.55[i] <- quantile(day_diff, 0.225) model.gompertz.upr.deriv.55[i] <- quantile(day_diff, 0.775) model.gompertz.lwr.deriv.65[i] <- quantile(day_diff, 0.175) model.gompertz.upr.deriv.65[i] <- quantile(day_diff, 0.825) model.gompertz.lwr.deriv.75[i] <- quantile(day_diff, 0.125) model.gompertz.upr.deriv.75[i] <- quantile(day_diff, 0.875) model.gompertz.lwr.deriv.85[i] <- quantile(day_diff, 0.075) model.gompertz.upr.deriv.85[i] <- quantile(day_diff, 0.925) model.gompertz.lwr.deriv.99[i] <- quantile(day_diff, 0.005) model.gompertz.upr.deriv.99[i] <- quantile(day_diff, 0.995) } # For the Richards' growth model # This will be the data frame where the derivative data from each simulation will be # placed. # Create the initial empty data frame. model.r.deriv <- data.frame(day = double(), deriv = double()) # Create the matrix of the variations of the parameters from the Richards' model that # we'll use for bootstrapping. Effectively, we use the standard error for each of the # coefficients in the model to create a normal distribution of each parameter. model.r.bootstrapcoef <- data.frame( b = rnorm(10000, mean=model.r$coefficients[1], sd=summary(model.r)$coefficients["b:(Intercept)","Std. Error"]), c = rnorm(10000, mean=model.r$coefficients[2], sd=summary(model.r)$coefficients["c:(Intercept)","Std. Error"]), d = rnorm(10000, mean=model.r$coefficients[3], sd=summary(model.r)$coefficients["d:(Intercept)","Std. Error"]), e = rnorm(10000, mean=model.r$coefficients[4], sd=summary(model.r)$coefficients["e:(Intercept)","Std. Error"]), f = rnorm(10000, mean=model.r$coefficients[5], sd=summary(model.r)$coefficients["f:(Intercept)","Std. Error"])) # Start of the bootstrapping. # 10,000 simulations are run. for(i in 1:10000) { # Set up the model parameters for this instance of the simulation b = model.r.bootstrapcoef$b[i] c = model.r.bootstrapcoef$c[i] d = model.r.bootstrapcoef$d[i] e = model.r.bootstrapcoef$e[i] f = model.r.bootstrapcoef$f[i] model.r.bootstrap.deriv <- vector() model.r.bootstrap.deriv[1:150] <- NA for(j in 49:150) { # The Richards' growth curve is in the following form: # f(x) = c + ( (d-c) / (1 + exp( b * (x - e)))^f) # # Where x = day number # b, c, d, e, f are the parameters of the model # The calculation below uses the model to calculate the cumulative cases for the # current day, and subtracts it from that from the previous day. model.r.bootstrap.deriv[j] <- ( c + ( (d-c) / (1 + exp( b * (j - e)))^f) ) - (c + ( (d-c) / (1 + exp( b * ((j-1) - e)))^f) ) } # Append the derivative data from this run of the simulation to those of all the previous # runs. bootstrap.deriv.temp <- data.frame(day = all_days, deriv = model.r.bootstrap.deriv) model.r.deriv <- data.frame( days = c(model.r.deriv$days, bootstrap.deriv.temp$day), deriv = c(model.r.deriv$deriv, bootstrap.deriv.temp$deriv)) } # Compute the confidence bands from the derivative data. for(i in 49:150) { day_diff <- model.r.deriv$deriv[model.r.deriv$days == i] # exclude any values less than zero day_diff <- day_diff[day_diff > 0] model.au.fit.deriv[i] <- median(day_diff) model.au.lwr.deriv[i] <- quantile(day_diff, 0.025) model.au.upr.deriv[i] <- quantile(day_diff, 0.975) model.au.lwr.deriv.05[i] <- quantile(day_diff, 0.475) model.au.upr.deriv.05[i] <- quantile(day_diff, 0.525) model.au.lwr.deriv.15[i] <- quantile(day_diff, 0.425) model.au.upr.deriv.15[i] <- quantile(day_diff, 0.575) model.au.lwr.deriv.25[i] <- quantile(day_diff, 0.375) model.au.upr.deriv.25[i] <- quantile(day_diff, 0.625) model.au.lwr.deriv.35[i] <- quantile(day_diff, 0.325) model.au.upr.deriv.35[i] <- quantile(day_diff, 0.675) model.au.lwr.deriv.45[i] <- quantile(day_diff, 0.275) model.au.upr.deriv.45[i] <- quantile(day_diff, 0.725) model.au.lwr.deriv.55[i] <- quantile(day_diff, 0.225) model.au.upr.deriv.55[i] <- quantile(day_diff, 0.775) model.au.lwr.deriv.65[i] <- quantile(day_diff, 0.175) model.au.upr.deriv.65[i] <- quantile(day_diff, 0.825) model.au.lwr.deriv.75[i] <- quantile(day_diff, 0.125) model.au.upr.deriv.75[i] <- quantile(day_diff, 0.875) model.au.lwr.deriv.85[i] <- quantile(day_diff, 0.075) model.au.upr.deriv.85[i] <- quantile(day_diff, 0.925) model.au.lwr.deriv.99[i] <- quantile(day_diff, 0.005) model.au.upr.deriv.99[i] <- quantile(day_diff, 0.995) } # scaling factor for fixed elements on chart sf = max(model.au.upr.deriv[1:today+7], na.rm=TRUE) / 1800 alphaPI = 0.08 plot_au_new <- ggplot(data=au_covid, aes(x=as.Date(date), y=new2)) + xlab("") + ylab("Daily Local COVID-19 cases (Vic)") + 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=as.Date("2021-9-20"), color = "red", size=0.5) + annotate("text", x=as.Date("2021-9-20"), y=1400*sf, label="Large public protests (20/9) ", hjust=1) + annotate("label", x=as.Date("2021-8-6"), y=1800*sf, label="Projection of new daily cases of COVID-19", size=5.5, hjust=0) + annotate("text", x=as.Date("2021-8-6"), y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) + annotate("text", x=as.Date("2021-8-6"), 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="grey50", size = 2.8) + geom_point(aes(y=new), 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), y=85*sf, label=au_covid$date[today], size=3) + annotate("label", x=(todaydate + 7), y=85*sf, label=au_covid$date[today+7], size=2.5) + coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-8-4"),todaydate+10) ) # scaling factor for fixed elements on chart # sf = max(model.gompertz.upr.deriv[1:today+7], na.rm=TRUE) / 1800 # A single data point was unexpectedly high on 14 Oct 2021 sf = 2500 / 1800 alphaPI = 0.08 plot_au_new_gompertz <- ggplot(data=au_covid, aes(x=as.Date(date), y=new2)) + xlab("") + ylab("Daily Local COVID-19 cases (Vic)") + 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=as.Date("2021-9-20"), color = "red", size=0.5) + annotate("text", x=as.Date("2021-9-20"), y=1400*sf, label="Large public protests (20/9) ", hjust=1) + annotate("label", x=as.Date("2021-8-6"), y=1800*sf, label="Projection of new daily cases of COVID-19", size=5.5, hjust=0) + annotate("text", x=as.Date("2021-8-6"), y=1700*sf, label= paste(au_covid$date[today], "@", time, "GMT+10"), hjust=0) + annotate("text", x=as.Date("2021-8-6"), 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="grey50", size = 2.8) + geom_point(aes(y=new), 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), y=85*sf, label=au_covid$date[today], size=3) + annotate("label", x=(todaydate + 7), y=85*sf, label=au_covid$date[today+7], size=2.5) + coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-8-4"),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 = 2500 / 1800 p1 <- ggplot(data=au_covid, aes(x=as.Date(date), y=new2)) + xlab("") + ylab("Confirmed local COVID-19 cases (Vic)") + geom_vline(xintercept=todaydate, color = "grey") + annotate("label", x=as.Date("2021-8-6"), y=1800*sf, label="Moving Avg (7 days)", size=5.5, hjust=0) + annotate("text", x=as.Date("2021-8-6"), 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-8-4"),todaydate+10) ) p2 <- ggplot(data=au_covid, aes(x=as.Date(date), y=new2)) + xlab("") + ylab("") + geom_vline(xintercept=todaydate, color = "grey") + annotate("label", x=as.Date("2021-8-6"), 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-8-4"),todaydate+10) ) p3 <- ggplot(data=au_covid, aes(x=as.Date(date), y=new2)) + xlab("") + ylab("Confirmed local COVID-19 cases (Vic)") + geom_vline(xintercept=todaydate, color = "grey") + annotate("label", x=as.Date("2021-8-6"), y=1800*sf, label="Gompertz model projections", size=5.5, hjust=0) + geom_point(colour="white", size = 3.2) + geom_point(colour="grey50", size = 2.8) + geom_point(aes(y=new), 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-8-4"),todaydate+10) ) p4 <- ggplot(data=au_covid, aes(x=as.Date(date), y=new2)) + xlab("") + ylab("") + geom_vline(xintercept=todaydate, color = "grey") + annotate("label", x=as.Date("2021-8-6"), y=1800*sf, label="Richards' curve projections", size=5.5, hjust=0) + geom_point(colour="white", size = 3.2) + geom_point(colour="grey50", size = 2.8) + geom_point(aes(y=new), 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-8-4"),todaydate+10) ) # 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-10-2"), y=1800*sf, label="Gompertz - 7 day projection", size=5.5, hjust=0) + annotate("text", x=as.Date("2021-10-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=900*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-10-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-10-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=900*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-10-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-10-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=900*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-10-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-10-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=900*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-10-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("Date") + ylab("Total number of COVID-19 cases") + scale_x_date(date_breaks = "1 weeks", date_labels = "%b %d") + scale_y_continuous(labels = comma) + geom_vline(xintercept=todaydate, color = "grey") + annotate("label", x=as.Date("2021-10-2"), y=1800*sf, label="Projection of FINAL SIZE of the outbreak", size=5.5, hjust=0) + annotate("text", x=as.Date("2021-10-2"), y=1700*sf, label="Assumption: fitted to Gompertz model", hjust=0) + annotate("text", x=as.Date("2021-10-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-10-1"),todaydate+4) ) f2 <- ggplot(data=model.history, aes(x=as.Date(date), y=r.d)) + xlab("Date") + ylab("") + scale_x_date(date_breaks = "1 weeks", date_labels = "%b %d") + scale_y_continuous(labels = comma) + geom_vline(xintercept=todaydate, color = "grey") + annotate("label", x=as.Date("2021-10-2"), y=1800*sf, label="Projection of FINAL SIZE of the outbreak", size=5.5, hjust=0) + annotate("text", x=as.Date("2021-10-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-10-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) # Write model outputs into a file sink('vic_model_output.txt') cat("Date\n") au_covid$date[today] cat("\nDay number:\n") today cat("\n\nRichards' growth curve model - model.r\n") summary(model.r) cat("\nmodel.r projections\n") model.au cat("\nmodel.r bootstrap simulation derivatives\n") model.au.fit.deriv cat("\n\nGompertz model - model.g\n") summary(model.g) cat("\nmodel.g projections\n") model.gompertz cat("\nmodel.g bootstrap simulation derivatives\n") model.gompertz.fit.deriv cat("\n\nmodel.g 7-day projection fit error\n") model.history.g.7fit.error cat("\nmodel.g 14-day projection fit error\n") model.history.g.14fit.error cat("\n\nmodel.r 7-day projection fit error\n") model.history.r.7fit.error cat("\nmodel.r 14-day projection fit error\n") model.history.r.14fit.error sink()