Comparing the Nowcast Estimates of the Poisson and negative binomial Models
Now let us plot and compare the estimates obtained from each of these nowcasts:
plot_estimates <- function(nowcast1, nowcast2, cases_per_date, now) {
# Ensure input data is valid
if (is.null(nowcast1$estimates) || is.null(nowcast2$estimates)) {
stop("Nowcast estimates are missing.")
}
# Extract estimates and credible intervals for nowcast without DoW effect
onset_dates1 <- nowcast1$estimates$onset_date
estimates1 <- nowcast1$estimates$estimate
lower1 <- nowcast1$estimates$q_0.025 # 2.5% quantile (lower bound of 95% CI)
upper1 <- nowcast1$estimates$q_0.975 # 97.5% quantile (upper bound of 95% CI)
# Extract estimates and credible intervals for nowcast with DoW effect
onset_dates2 <- nowcast2$estimates$onset_date
estimates2 <- nowcast2$estimates$estimate
lower2 <- nowcast2$estimates$q_0.025
upper2 <- nowcast2$estimates$q_0.975
# Extract eventual case counts
case_dates <- cases_per_date$onset_week
case_counts <- cases_per_date$count
# Calculate plot range
min_val <- min(c(lower1, lower2, case_counts), na.rm = TRUE)
max_val <- max(c(upper1, upper2, case_counts), na.rm = TRUE) * 1.35
# Create the plot
plot(
onset_dates1, estimates1, col = 'blue', type = 'l',
xlab = 'Onset Date', ylab = 'Cases',
ylim = c(min_val, max_val), lwd = 2,
main = paste0('Incidence Estimates for ', now)
)
lines(onset_dates2, estimates2, col = 'red', lwd = 2)
# Add 95% PI shaded regions for both nowcasts
polygon(c(onset_dates1, rev(onset_dates1)), c(lower1, rev(upper1)),
col = rgb(0, 0, 1, 0.2), border = NA)
polygon(c(onset_dates2, rev(onset_dates2)), c(lower2, rev(upper2)),
col = rgb(1, 0, 0, 0.2), border = NA)
# Add true case counts as points
points(case_dates, case_counts, col = 'black', pch = 20)
# Add a legend
legend(
'topleft',
legend = c('Estimates Poisson model',
'Estimates negative binomial model',
'95% PI (Poisson model)',
'95% PI (negative binomial model)',
'Eventual cases'),
col = c('blue', 'red', rgb(0, 0, 1, 0.2), rgb(1, 0, 0, 0.2), 'black'),
lty = c(1, 1, NA, NA, NA), lwd = c(2, 2, NA, NA, NA),
pch = c(NA, NA, 15, 15, 20),
pt.cex = c(NA, NA, 1.5, 1.5, 1),
cex = 0.9
)
}
cases_per_week <- denguedat %>% group_by(onset_week) %>% summarize(count=n())
for (t in seq_along(test_dates)) {
now <- test_dates[t]
nowcast_pois <- nowcasts_pois[[t]]
nowcast_nb <- nowcasts_nb[[t]]
cases_per_week1 <- cases_per_week[cases_per_week$onset_week <= now, ]
plot_estimates(nowcast_pois, nowcast_nb, cases_per_week1, now)
}
As seen in these plots, the mean estimates of the two models are very similar, but the negative binomial model has wider prediction intervals, which appear to better capture the actual number of eventual cases.
Now, let us calculate the WIS for these nowcasting estimates. To do
this, we first need to extract the quantile estimates from the nowcast
results and use them to construct a data frame compatible with the
scoringutils package.
