sf = 2600 / 1800 alphaPI = 0.08 r0 <- 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 = "dark green") + geom_vline(xintercept=as.Date("2021-6-25"), color = "red", size=0.5) + annotate("text", x=as.Date("2021-6-25"), y=1200*sf, label="A ", hjust=1) + geom_vline(xintercept=as.Date("2021-7-9"), color = "red", size=0.5) + annotate("text", x=as.Date("2021-7-9"), y=1200*sf, label="B ", hjust=1) + geom_vline(xintercept=as.Date("2021-8-16"), color = "red", size=0.5) + annotate("text", x=as.Date("2021-8-16"), y=1200*sf, label="C ", hjust=1) + geom_vline(xintercept=as.Date("2021-8-23"), color = "red", size=0.5) + annotate("text", x=as.Date("2021-8-23"), y=1200*sf, label="D ", hjust=1) + annotate("text", x=as.Date("2021-6-20"), y=900*sf, label="A (25/6): Initial lockdown", hjust=0) + annotate("text", x=as.Date("2021-6-20"), y=775*sf, label="B (9/7): Intensifies in hot LGAs", hjust=0) + annotate("text", x=as.Date("2021-6-20"), y=650*sf, label="C (16/8): NSW-wide lockdown", hjust=0) + annotate("text", x=as.Date("2021-6-20"), y=525*sf, label="D (23/8): Curfews & permits in hot LGAs", hjust=0) + annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Projection of new daily cases of COVID-19 (Today)", 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) + annotate("text", x=as.Date("2021-6-20"), y=1475*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), y=85*sf, label=au_covid$date[today], size=3) + coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),todaydate+10) ) r0 # Code - r7 today = today - 7 todaydate = todaydate - 7 r7data <- vector() r7data.exclude <- vector() r7data[1:198] <- NA r7data.exclude[1:198] <- NA r7data[1:today] <- au_covid$new[1:today] r7data.exclude[(today+1):(today+7)] <- au_covid$new[(today+1):(today+7)] model.cases <- au_covid$tot[1:today] model.day <- au_covid$day[1:today] 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.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:198] <- NA model.au.lwr.deriv[1:198] <- NA model.au.upr.deriv[1:198] <- NA model.au.lwr.deriv.05[1:198] <- NA model.au.upr.deriv.05[1:198] <- NA model.au.lwr.deriv.15[1:198] <- NA model.au.upr.deriv.15[1:198] <- NA model.au.lwr.deriv.25[1:198] <- NA model.au.upr.deriv.25[1:198] <- NA model.au.lwr.deriv.35[1:198] <- NA model.au.upr.deriv.35[1:198] <- NA model.au.lwr.deriv.45[1:198] <- NA model.au.upr.deriv.45[1:198] <- NA model.au.lwr.deriv.55[1:198] <- NA model.au.upr.deriv.55[1:198] <- NA model.au.lwr.deriv.65[1:198] <- NA model.au.upr.deriv.65[1:198] <- NA model.au.lwr.deriv.75[1:198] <- NA model.au.upr.deriv.75[1:198] <- NA model.au.lwr.deriv.85[1:198] <- NA model.au.upr.deriv.85[1:198] <- NA model.au.lwr.deriv.99[1:198] <- NA model.au.upr.deriv.99[1:198] <- NA # For the Richards' growth model 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.boostrapcoef <- 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.boostrapcoef$b[i] c = model.r.boostrapcoef$c[i] d = model.r.boostrapcoef$d[i] e = model.r.boostrapcoef$e[i] f = model.r.boostrapcoef$f[i] model.r.bootstrap.deriv <- vector() model.r.bootstrap.deriv[1:198] <- NA for(j in 2:198) { # 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) ) } 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 2:198) { day_diff <- model.r.deriv$deriv[model.r.deriv$days == i] 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 = 2600 / 1800 alphaPI = 0.08 r7 <- 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 = "dark green") + geom_vline(xintercept=todaydate + 7, color = "grey") + geom_vline(xintercept=as.Date("2021-6-25"), color = "red", size=0.5) + geom_vline(xintercept=as.Date("2021-7-9"), color = "red", size=0.5) + geom_vline(xintercept=as.Date("2021-8-16"), color = "red", size=0.5) + geom_vline(xintercept=as.Date("2021-8-23"), color = "red", size=0.5) + annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Projection from 7 days ago", size=5.5, 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(aes(y=r7data), colour="grey20", size = 2.8) + geom_point(aes(y=r7data.exclude), colour="grey50", 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=2) + annotate("text", x=(todaydate + 7), y=model.au.lwr.deriv[today + 7], label=round(model.au.lwr.deriv[today + 7]), size=2) + annotate("label", x=(todaydate), y=1800*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-6-17"), todaydate+17) ) r7 # Code - r14 today = today - 7 todaydate = todaydate - 7 r14data <- vector() r14data.exclude <- vector() r14data[1:198] <- NA r14data.exclude[1:198] <- NA r14data[1:today] <- au_covid$new[1:today] r14data.exclude[(today+1):(today+14)] <- au_covid$new[(today+1):(today+14)] model.cases <- au_covid$tot[1:today] model.day <- au_covid$day[1:today] 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.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:198] <- NA model.au.lwr.deriv[1:198] <- NA model.au.upr.deriv[1:198] <- NA model.au.lwr.deriv.05[1:198] <- NA model.au.upr.deriv.05[1:198] <- NA model.au.lwr.deriv.15[1:198] <- NA model.au.upr.deriv.15[1:198] <- NA model.au.lwr.deriv.25[1:198] <- NA model.au.upr.deriv.25[1:198] <- NA model.au.lwr.deriv.35[1:198] <- NA model.au.upr.deriv.35[1:198] <- NA model.au.lwr.deriv.45[1:198] <- NA model.au.upr.deriv.45[1:198] <- NA model.au.lwr.deriv.55[1:198] <- NA model.au.upr.deriv.55[1:198] <- NA model.au.lwr.deriv.65[1:198] <- NA model.au.upr.deriv.65[1:198] <- NA model.au.lwr.deriv.75[1:198] <- NA model.au.upr.deriv.75[1:198] <- NA model.au.lwr.deriv.85[1:198] <- NA model.au.upr.deriv.85[1:198] <- NA model.au.lwr.deriv.99[1:198] <- NA model.au.upr.deriv.99[1:198] <- NA # For the Richards' growth model 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.boostrapcoef <- 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.boostrapcoef$b[i] c = model.r.boostrapcoef$c[i] d = model.r.boostrapcoef$d[i] e = model.r.boostrapcoef$e[i] f = model.r.boostrapcoef$f[i] model.r.bootstrap.deriv <- vector() model.r.bootstrap.deriv[1:198] <- NA for(j in 2:198) { # 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) ) } 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 2:198) { day_diff <- model.r.deriv$deriv[model.r.deriv$days == i] 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 = 2600 / 1800 alphaPI = 0.08 r14 <- 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 = "dark green") + geom_vline(xintercept=todaydate + 14, color = "grey") + geom_vline(xintercept=as.Date("2021-6-25"), color = "red", size=0.5) + geom_vline(xintercept=as.Date("2021-7-9"), color = "red", size=0.5) + geom_vline(xintercept=as.Date("2021-8-16"), color = "red", size=0.5) + geom_vline(xintercept=as.Date("2021-8-23"), color = "red", size=0.5) + annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Projection from 14 days ago", size=5.5, 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(aes(y=r14data), colour="grey20", size = 2.8) + geom_point(aes(y=r14data.exclude), colour="grey50", 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 + 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=2) + annotate("text", x=(todaydate + 14), y=model.au.lwr.deriv[today + 14], label=round(model.au.lwr.deriv[today + 14]), size=2) + annotate("label", x=(todaydate), y=1800*sf, label=au_covid$date[today], size=3) + 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+24) ) r14 # Code - r21 today = today - 7 todaydate = todaydate - 7 r21data <- vector() r21data.exclude <- vector() r21data[1:198] <- NA r21data.exclude[1:198] <- NA r21data[1:today] <- au_covid$new[1:today] r21data.exclude[(today+1):(today+21)] <- au_covid$new[(today+1):(today+21)] model.cases <- au_covid$tot[1:today] model.day <- au_covid$day[1:today] 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.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:198] <- NA model.au.lwr.deriv[1:198] <- NA model.au.upr.deriv[1:198] <- NA model.au.lwr.deriv.05[1:198] <- NA model.au.upr.deriv.05[1:198] <- NA model.au.lwr.deriv.15[1:198] <- NA model.au.upr.deriv.15[1:198] <- NA model.au.lwr.deriv.25[1:198] <- NA model.au.upr.deriv.25[1:198] <- NA model.au.lwr.deriv.35[1:198] <- NA model.au.upr.deriv.35[1:198] <- NA model.au.lwr.deriv.45[1:198] <- NA model.au.upr.deriv.45[1:198] <- NA model.au.lwr.deriv.55[1:198] <- NA model.au.upr.deriv.55[1:198] <- NA model.au.lwr.deriv.65[1:198] <- NA model.au.upr.deriv.65[1:198] <- NA model.au.lwr.deriv.75[1:198] <- NA model.au.upr.deriv.75[1:198] <- NA model.au.lwr.deriv.85[1:198] <- NA model.au.upr.deriv.85[1:198] <- NA model.au.lwr.deriv.99[1:198] <- NA model.au.upr.deriv.99[1:198] <- NA # For the Richards' growth model 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.boostrapcoef <- 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.boostrapcoef$b[i] c = model.r.boostrapcoef$c[i] d = model.r.boostrapcoef$d[i] e = model.r.boostrapcoef$e[i] f = model.r.boostrapcoef$f[i] model.r.bootstrap.deriv <- vector() model.r.bootstrap.deriv[1:198] <- NA for(j in 2:198) { # 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) ) } 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 2:198) { day_diff <- model.r.deriv$deriv[model.r.deriv$days == i] 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 = 2600 / 1800 alphaPI = 0.08 r21 <- 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 = "dark green") + geom_vline(xintercept=todaydate + 21, color = "grey") + geom_vline(xintercept=as.Date("2021-6-25"), color = "red", size=0.5) + geom_vline(xintercept=as.Date("2021-7-9"), color = "red", size=0.5) + geom_vline(xintercept=as.Date("2021-8-16"), color = "red", size=0.5) + geom_vline(xintercept=as.Date("2021-8-23"), color = "red", size=0.5) + annotate("label", x=as.Date("2021-6-20"), y=1800*sf, label="Projection from 21 days ago", size=5.5, 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(aes(y=r21data), colour="grey20", size = 2.8) + geom_point(aes(y=r21data.exclude), colour="grey50", 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 + 21), y=model.au.fit.deriv[today + 21], label=round(model.au.fit.deriv[today + 21]), size=3) + annotate("text", x=(todaydate + 21), y=model.au.upr.deriv[today + 21], label=round(model.au.upr.deriv[today + 21]), size=2) + annotate("text", x=(todaydate + 21), y=model.au.lwr.deriv[today + 21], label=round(model.au.lwr.deriv[today + 21]), size=2) + annotate("label", x=(todaydate), y=1800*sf, label=au_covid$date[today], size=3) + annotate("label", x=(todaydate + 21), y=85*sf, label=au_covid$date[today+21], size=2.5) + coord_cartesian(ylim = c(0, 1800*sf), xlim = c(as.Date("2021-6-17"),todaydate+31) ) r21