Statistical Model of Covid-19 Vaccine Rollout in Australia, Feb to July 2021

vaxAus.R

# Data downloaded from CovidBaseAU and Our World in Data,
# https://github.com/owid/covid-19-data/
vax_aus <- read.csv("Australia.txt")
vax_aus$Date <- as.Date(vax_aus$date)
ix1 <- which.max(vax_aus$Date > "2021-05-02") - 1
idx <- which.max(vax_aus$Date > "2021-07-02") - 1
ix2 <- which.max(vax_aus$Date > "2021-09-02") - 1
ix3 <- which.max(vax_aus$Date > "2021-11-02") - 1
ix4 <- which.max(vax_aus$Date > "2021-10-02") - 1
vax_aus$Date[idx]
daily_vax <- diff(c(0,vax_aus$total_vaccinations[1:idx]))
plot.ts(daily_vax, xaxt='n', ylab="Doses per day", xlab="")
abline(h=0,lty=2)
vax_aus$date[c(3,40,70,98,126,157,188)]
axis(side=1, at=c(3,40,70,98,126,157,188),
     labels=c("Feb 23","Apr 1","May 1","June 1","July 1","Aug 1","Sept 1"))
x <- 1:ix1
y <- daily_vax[x]
fit1 <- lm(y ~ x - 1)
summary(fit1)
AIC(fit1)
BIC(fit1)
lines(c(0,ix2), c(0, fit1$coefficients * ix2), col=4, lty=2, lwd=2)
pred1 <- predict(fit1, newdata=list(x=ix1:ix2), interval="prediction")
#lines(ix1:ix2, pred1[,"lwr"], col=4, lty=3)
#lines(ix1:ix2, pred1[,"upr"], col=4, lty=3)

day <- 1:ix3
fit2 <- lm(log(y) ~ log(x))
summary(fit2)
AIC(fit2)
BIC(fit2)
yhat <- exp(fit2$coeff[1]) * day^fit2$coeff[2]
lines(day, yhat, col=2, lty=3, lwd=3)

ix7 <- which.max(vax_aus$Date > "2021-03-21") - 1
lines(c(ix7,ix2), c(318000/7, 318000/7), col=3, lwd=3, lty=4)

legend("topleft", legend=c("doses","flat","linear","exponential"),
       lty=c(1,4,2,3), col=c(1,3,4,2), lwd=c(1,3,2,3))
title(main="Vaccine Rollout in Australia, 2021")

## Weekly Totals ----
vax_weeklyMx <- matrix(nrow=19, ncol=6)
colnames(vax_weeklyMx) <- c("Index","Date","Observed","Flat","Linear","Exponential")
vax_weekly <- data.frame(vax_weeklyMx)
vax_weekly$Date <- as.Date(c("2021-02-27","2021-03-06","2021-03-13","2021-03-20","2021-03-27",
                     "2021-04-03","2021-04-10","2021-04-17","2021-04-24","2021-05-01",
                     "2021-05-08","2021-05-15","2021-05-22","2021-05-29","2021-06-05",
                     "2021-06-12","2021-06-19","2021-06-26","2021-07-03"))
vax_weekly$Flat <- 318000
iWk <- which.max(vax_aus$Date > vax_weekly$Date[1]) - 1
vax_weekly$Index[1] <- iWk
vax_weekly$Observed[1] <- vax_aus$total_vaccinations[iWk]
xPred <- 1:iWk
yPred <- predict(fit1, newdata=list(x=xPred), interval="prediction")
vax_weekly$Linear[1] <- sum(yPred[,"fit"])
yPred <- exp(fit2$coeff[1]) * xPred^fit2$coeff[2]
vax_weekly$Exponential[1] <- sum(yPred)
for (i in 2:nrow(vax_weekly)) {
  jWk <- which.max(vax_aus$Date > vax_weekly$Date[i]) - 1
  vax_weekly$Index[i] <- jWk
  vax_weekly$Observed[i] <- vax_aus$total_vaccinations[jWk] - vax_aus$total_vaccinations[iWk]
  xPred <- iWk + 1:7
  yPred <- predict(fit1, newdata=list(x=xPred), interval="prediction")
  vax_weekly$Linear[i] <- sum(yPred[,"fit"])
  yPred <- exp(fit2$coeff[1]) * xPred^fit2$coeff[2]
  vax_weekly$Exponential[i] <- sum(yPred)
  iWk <- jWk
}
sum(vax_weekly$Observed)
vax_aus$total_vaccinations[jWk]
sum(vax_weekly$Flat)
vax_aus$total_vaccinations[jWk] - sum(vax_weekly$Flat)
sum(vax_weekly$Linear)
vax_aus$total_vaccinations[jWk] - sum(vax_weekly$Linear)
sum(vax_weekly$Exponential)
vax_aus$total_vaccinations[jWk] - sum(vax_weekly$Exponential)

plot.ts(vax_weekly$Observed, xlab="", xaxt='n', ylab="Doses per week")
abline(h=318000, col=3, lwd=3, lty=4)
axis(side=1, at=c(1,2,6,10,15,19), labels=c("Feb 27","","Apr 3","May 1","June 5","July 3"))
lines(1:19, vax_weekly$Linear, col=4, lty=2, lwd=2)
lines(1:19, vax_weekly$Exponential, col=2, lty=3, lwd=3)
title(main="Weekly Totals")
legend("topleft", legend=c("doses","flat","linear","exponential"),
       lty=c(1,4,2,3), col=c(1,3,4,2), lwd=c(1,3,2,3))

## RMSPE
sqrt(sum((vax_weekly$Observed[11:19] - vax_weekly$Flat[11:19])^2)/9)
sqrt(sum((vax_weekly$Observed[11:19] - vax_weekly$Linear[11:19])^2)/9)
sqrt(sum((vax_weekly$Observed[11:19] - vax_weekly$Exponential[11:19])^2)/9)

iWon <- which.max(vax_aus$total_vaccinations > 4e7)
vax_aus$Date[iWon]
vax_aus$total_vaccinations[iWon]
iW2 <- which.max(vax_aus$people_fully_vaccinated > 2e7)
vax_aus$Date[iW2]
vax_aus$people_fully_vaccinated[iW2]

## now for two more months...
daily_vax <- diff(c(0,vax_aus$total_vaccinations[1:ix3]))
plot.ts(daily_vax, xaxt='n', ylab="Doses per day", xlab="")
axis(side=1, at=c(3,40,70,98,126,157,188,218,249),
     labels=c("Feb 23","Apr 1","May 1","Jun 1","July 1","Aug 1","Sept 1","Oct 1","Nov 1"))
lines(c(0,ix3), c(0, fit1$coefficients * ix3), col=4, lty=2, lwd=2)
lines(day, yhat, col=2, lty=3, lwd=3)
abline(h=0,lty=2)

x2 <- 1:ix4
y2 <- daily_vax[x2]
fit3 <- lm(log(y2) ~ log(x2))
summary(fit3)
yhat3 <- exp(fit3$coeff[1]) * day^fit3$coeff[2]
lines(day, yhat3, col=3, lty=4, lwd=3)
legend("topleft", legend=c("doses","linear","exp 1", "exp 2"),
       lty=1:4, col=c(1,4,2, 3), lwd=c(1:3,3))

plot(x2, fit3$residuals, ylab=expression(log(e)), xaxt='n', xlab="", type='l',col=4)
abline(h=0, lty=2)
axis(side=1, at=c(3,40,70,98,126,157,188,218,249),
     labels=c("Feb 23","Apr 1","May 1","Jun 1","July 1","Aug 1","Sept 1","Oct 1","Nov 1"))

## seasonality ----
ehat <- fit3$residuals[idx:ix4]
plot(idx:ix4, ehat - min(ehat), xaxt='n', xlab="", type='l',col=4)
abline(h=0, lty=2)

points(136 - 7, 0)
points(136, 0)
points(136 + 7, 0)

ix5 <- seq(3,ix4,by=7)
points(ix5, fit3$residuals[ix5] - min(ehat))

plot.ts(daily_vax, xaxt='n', ylab="Doses per day", xlab="")
points(ix5, daily_vax[ix5], col=4)
points(108, daily_vax[108], col=2, pch=2, cex=1.5)
points(109, daily_vax[109], col=2, pch=2, cex=1.5)
ix7 <- seq(1,100,by=7)
points(ix7, daily_vax[ix7])
ix5[16] <- 109
ix5[1:15] <- ix7
ix5[5] <- 28
ix5[9] <- 56
ix5[14] <- 96
ix5[15] <- 103
plot.ts(daily_vax, xaxt='n', ylab="Doses per day", xlab="")
points(ix5, daily_vax[ix5], col=4)
diff(ix5)
vax_aus$Date[85]
vax_aus$Date[96]

plot.ts(fit3$residuals, xaxt='n', ylab="Residuals", xlab="")
points(ix5, fit3$residuals[ix5], col=4)

ix6 <- (1:ix4-1) %% 7
ix6[86:95] <- 1:10 *7/11
ix6[96:102] <- 0:6
ix6[103:108] <- 0:5 * 7/6
ix6[109:114] <- 0:5 * 7/6
ix6[115:ix4] <- (115:ix4 -3) %% 7
qfit <- -4/49 * ix6^2 + 4/7 * ix6
lines(1:ix4, qfit - 3, lty=3, col=2, lwd=2)
points(ix5,qfit[ix5] - 3)

fit6 <- lm(log(y2) ~ log(x2) + qfit)
summary(fit6)
plot(x2, fit6$residuals, ylab=expression(log(e)), xaxt='n', xlab="", type='l',col=4)
points(ix5, fit6$residuals[ix5])
abline(h=0,lty=2)

plot.ts(daily_vax, xaxt='n', ylab="Doses per day", xlab="", ylim=c(0,4e5), lwd=2)
axis(side=1, at=c(3,40,70,98,126,157,188,218,249),
     labels=c("Feb 23","Apr 1","May 1","Jun 1","July 1","Aug 1","Sept 1","Oct 1","Nov 1"))
yhat6 <- exp(fit6$coef[1]) * x2^fit6$coef[2] * exp(qfit*fit6$coef[3])
lines(x2, yhat6, col=4, lty=3, lwd=3)
yhat7 <- exp(fit6$coef[1]) * x2^fit6$coef[2]
lines(x2, yhat7, col=2, lty=2, lwd=2)
lines(x2, yhat7 * exp(fit6$coef[3]), col=2, lty=2, lwd=2)
abline(h=0,lty=2)
#lines(x2, yhat7 * exp(0.5*fit6$coef[3]), col=2, lty=2, lwd=2)

legend("topleft",lty=c(1,3,2), lwd=c(2,3,2), col=c(1,4,2),
       legend=c("observed","fitted","bounds"))
title(main="Exponential Model with Seasonality")

yhat3 <- exp(fit3$coeff[1]) * x2^fit3$coeff[2]
plot(x2, y2 - yhat3, xaxt='n', xlab="", ylab="Residuals", type='l', col=4)
abline(h=0, lty=2)

plot(x2, abs(y2 - yhat3), xaxt='n', xlab="", ylab=expression(abs(e)), type='l', col=4)
abline(h=0, lty=2)
title(main="Absolute Residuals")

y3 <- abs(y2 - yhat3)
fit4 <- lm(y3 ~ x2 - 1)
lines(c(0,ix3),c(0,fit4$coefficients * ix3), lty=2, lwd=2)

wt <- 1/(x2 * fit4$coefficients)
wt <- wt/sum(wt)
fit5 <- lm(log(y2) ~ log(x2), weights = wt)
BIC(fit5)
BIC(fit3)

plot(x2, fit3$residuals, ylab=expression(log(e)), xaxt='n', xlab="", type='l',col=4)
abline(h=0, lty=2)

## whole dataset ----
daily_vax <- diff(c(0,vax_aus$total_vaccinations))
plot.ts(daily_vax, xaxt='n', ylab="Doses per day", xlab="Month", ylim=c(0,4e5))
axis(side=1, at=c(9,40,70,98,126,157,188,218,249,279),
     labels=c("M","A","M","J","J","A","S","O","N","D"))
ixB <- which.max(vax_aus$Date > "2021-12-01") - 1
ixB0 <- which.max(vax_aus$Date > "2021-03-01") - 1
ixB1 <- which.max(vax_aus$Date > "2022-01-01") - 1
ixB2 <- which.max(vax_aus$Date > "2022-02-01") - 1
ixB3 <- which.max(vax_aus$Date > "2022-03-01") - 1
ixB4 <- which.max(vax_aus$Date > "2022-04-01") - 1
axis(side=1, at=c(ixB1,ixB2,ixB3,ixB4), labels=c("2022","F","M","A"))
abline(h=0, lty=2)

lines(x2, yhat6, col=5, lty=3, lwd=3)
lines(x2, yhat7, col=2, lty=2, lwd=2)
lines(x2, yhat7 * exp(fit6$coef[3]), col=2, lty=2, lwd=2)
title(main="Vaccine Rollout in Australia")

plot.ts(diff(vax_aus$total_vaccinations), xaxt='n', ylab="Doses per day (logarithmic)", xlab="", log='y')
axis(side=1, at=c(9,40,70,98,126,157,188,218,249,279),
     labels=c("M","A","M","J","J","A","S","O","N","D"))
axis(side=1, at=c(ixB1,ixB2,ixB3,ixB4), labels=c("2022","F","M","A"))
lines(x2, yhat7, col=2, lty=2, lwd=2)
lines(x2, yhat7 * exp(fit6$coef[3]), col=2, lty=2, lwd=2)
title(main="logarithmic plot")

plot.ts(diff(vax_aus$total_vaccinations), xaxt='n', ylab="Doses per day (logarithmic)", xlab="", log='xy')
axis(side=1, at=c(9,40,70,98,126,157,188,218,249,279),
     labels=c("March","April","M","J","J","","","","",""))
axis(side=1, at=c(ixB1,ixB2,ixB3,ixB4), labels=c("2022","F","M","A"))
lines(x2, yhat7, col=2, lty=2, lwd=2)
lines(x2, yhat7 * exp(fit6$coef[3]), col=2, lty=2, lwd=2)
title(main="doubly-logarithmic plot")