##  1 Experiments and Results

Since June 2021, we have started our experiments and begun to submit results
of our HKUST-DNN model (trained on real data only, without DA) to the US Hub
[hub], and, since November 2021, results of our HKUST-DNN_DA model (trained on
real data and validated by augmented data) to the German Hub [hubde]. Since
then, based on the data collected up to every Saturday, our models are trained
to predict $cum\\_death$ in the US and $inc\\_case$ in Germany for the four
subsequent Saturdays, i.e., four different horizons spanning 1 to 4 weeks
ahead. Instead of predicting all the targets of the two countries, our scope
of experiments has focused on predicting each country with a specific target,
i.e., $cum\\_death$ for the US and $inc\\_case$ for Germany.

Our analysis below will focus on a 9-month window between June 2021 and
February 2022, in which our models were designed and results were closely
monitored. This window covers the period when the dominating variant of
COVID-19 was transiting from Delta to earlier versions of Omicron. For various
issues and the performance after this period of time, we will discuss them
later in Section LABEL:sec:discussion.

###  1.1 Evaluation Metrics

In this paper, the forecast performance of different models will mainly be
evaluated using the mean absolute error (MAE) of their point estimates (if
absent, the estimated medians) and the weighted interval score (WIS)
[bracher2021evaluating_wis] of their distributional forecasts in a 23-quantile
format.

##### Ground Truth

We adopt the death and case counts reported by Johns Hopkins University (JHU)
[jhu] and Robert Koch Institut (RKI) as ground truth values of the US and
Germany, respectively. The RKI’s data is retrieved through the German Hub
[hubde]. Please note that the JHU dataset is updated retrospectively, in which
the latest values may be different from earlier versions. As a matter of fact,
there are frequently corrections made to the historical record of JHU’s
$cum\\_death$. To isolate errors induced by such corrections from the
prediction errors, we opt to compare the predicted increments (since the day
of the last observation) to the ground truth increments instead of comparing
the actual predicted values to the corresponding ground truths directly.

##### Missing Data

Every week, the US Hub collects forecasts from dozens of models on various
prediction targets. Since there could be some missing forecasts (on certain
dates or in some locations), we will only compare models which have
contributed at least 90% of the forecasts falling in our scope of interest,
i.e., on $cum\\_death$, 51 states (including Washington DC), 23 quantiles, 4
horizons (1- to 4-week), and within our 9-month window. Furthermore, to
maintain a fair comparison between the selected models which might still have
slightly different missing values, our calculation of MAE and WIS will only
include events in which all models have participated. If there is any missing
forecast from any model, that forecast will not be included into our
calculation of MAE and WIS. This same rule applies to the German forecasts.

##### Mean Absolute Error (MAE)

For point estimates (or predicted medians), their performance is simply
measured by their MAE to the ground truth. However, since the reported state-
level $cum\\_death$ in the US can possibly be updated after our forecast, we
opt to compare the predicted increments to the ground truth increments,
instead of the absolute values, as explained above in Section 1.1.

##### Weighted Interval Score (WIS)

For probabilistic forecast evaluation (in 23-quantile format), we adopt the
WIS scoring method proposed by [bracher2021evaluating_wis] to compare forecast
submissions to the German Hub. Unlike traditional methods such as the
continuous ranked probability score (CRPS) for measuring the errors between
continuous distributions, WIS only requires certain quantiles of such
distributions to be evaluated, such as the 23-quantile format which is being
used in this work. It is also adopted by other researchers
[cramer2021evaluation_ensemble] to compare the results in the US Hub.

###  1.2 Results in the US (without DA)

Using our transformer-based deep learning model, we have participated in
forecasting 1 to 4 weeks ahead (4 horizons) the $cum\\_death$ in 51 states
(including Washington DC) of the US since June 2021. Our submitted results are
named under HKUST-DNN in the US Hub [hub].

Figure 1 shows the trends of MAE of the top 15 models in the US Hub based on
the criteria mentioned in Section 1.1. Its horizontal axis represents the
dates being forecasted, from July 2021 to February 2022. The MAE of every
forecasted date is an aggregated result of 4 different horizons (1 to 4 weeks
ahead).

In terms of the $cum\\_death$ forecast, our model HKUST-DNN has achieved a
relatively low MAE compared to other models submitted to the Hub. It also has
comparable performance to the Hub ensembles named under ”COVIDhub…ensemble”.
Note that they are the combined and balanced results of all the models from
various methods submitted to the Hub.

The WIS of forecasts from different models have very similar trends to their
corresponding MAE. It is also included here for completion.

Figure 1: Trends of forecast performance (MAE and WIS) in the US

Table LABEL:tab:mae_rank_wis_us reveals the overall MAE over this period of
time. They are the aggregated results from all 51 states and 4 horizons (1 to
4 weeks ahead). For each date, horizon and state, each model has a rank among
others, depending on who has a lower MAE. The average of these ranks is used
to sort the entries in the table in ascending order. In term of this average
rank, our model HKUST-DNN has ranked relatively high among others. Note that
the top-5 models are all Hub ensembles which are aggregated results of all
submitted models.

table/tab_us_mae-rank-wis

###  1.3 Results in Germany (with DA)

We have also tested our model with DA on forecasting the sub-national
$inc\\_case$ of the 16 states in Germany. Adopting the naming convention of
the Germany Hub [hubde], these 16 German states are also aliased as GM01 to
GM16.111https://github.com/cfong32/covid19-forecast-hub-
de/blob/master/template/state_codes_germany.csv

Compared to the US, Germany has a much smaller dataset, attributing to its
fewer states (also with a smaller population and territory). From our
observation, the trends of $case$ and $death$ in Germany exhibit less
variation, especially prior to the recent waves of Omicron variants.

Therefore, in order to avoid overfitting to the small training dataset, we
introduce our DA method in training the models for Germany. Empirically, we
found that the models trained with our DA method perform better than those
without. Based on the following divergence analysis, we will also show that
our augmented data has a distribution closer to the test data (to be
forecasted), rather than the real training data (historical data).

####  1.3.1 Augmented Data

##### Fitted parameters

Figure 2 shows the results of the fitted parameters of our data model. On
every day (or timestep) $t$, the parameters are obtained by fitting our SIRD
model to a recent history of $W\in\\{15,16,\dots,28\\}$ days in each of the 16
states. The parameters $\beta$, $\gamma$, $\omega$ and $\delta$ depict the
rates of infection, recovery, reinfection and death, respectively, $\eta$ is
the initial portion of infectious population among the infected ones, and the
ratio of $\eta:1-\eta$ represents the proposed initial ratio of $I$ to $R$.
Since $\eta$ is initialized to 0.5 before fitting, unless the L-BFGS-B
optimization needs a much lower or higher value to explain the trends, it
tends to converge to local optima around 0.5.

Figure 2: Fitted parameters of the 16 states in Germany.

##### Generating new data

Sampling from the observed or fitted distributions of $S$, $D$, $\beta$,
$\gamma$, $\omega$, $\delta$ and $\eta$, we obtained new parameters to
generate synthetic trends for DA. This process was described in Section
LABEL:sec:data-generation.

Figure 3 shows samples of the best and the worst fitted 28-day trends
($W=28$). They respectively represent the most easy and difficult cases that
our SIRD model has comprehended so far. These examples are shown primarily to
demonstrate the quality of the fitting processes.

A couple of points are worth noting. First of all, for visualization purpose,
the death counts are multiplied by a factor of 20, so that their trends can be
visually compared more easily. Green and red dashed lines are the inferred $I$
and $R$. Secondly, unlikely the real observations which are more noisy, these
fitted trends are smoother and lack of noise. Yet we hope that they can still
capture much information about the situations that are plausible to happen.

Figure 3: Best- and worst-fitted state-level trends in Germany

##### Divergence analysis

To verify the suitability of our synthetic data, we conducted a divergence
analysis between the original training data and our synthetic data. In
particular, we want to see if the synthetic data is close to the test data, so
as to enhance the predictability over the situations in the near future.

We first visualize the data using t-SNE [van2008visualizing_tsne] plots.
Training examples $[X_{i},Y_{i}]$ are first serialized into vectors
$Z_{i}\in\mathbb{R}^{L{\times}d_{\mathit{in}}+T{\times}d_{\mathit{out}}}$.
Figure 4 shows two t-SNE plots which visualize these vectors in a 2-D space.
Every plot corresponds to a particular date of forecast. In order to show the
test and recent observations more clearly, their point sizes are enlarged.

Two particular dates are shown. The first plot in Figure 4 visualizes various
available datasets when forecast was carried out on 2021-Jul-17, with $L=14$
and $T=28$. Its test data refers to a four-week horizon till 2021-Aug-14. The
recent observations ($L+T+7$ days by 2021-Jul-17) were supposed to form the
validation set in the usual practice without DA. All historical data
(approximately two years) was available for training. All synthetic datasets
(with $F\in{0,14,28}$) were generated based on observations as of Jul-17.
Based on our divergence analysis, as shown in Figure 5, this time reached the
lowest divergence so far between the test data and our generated data. For
comparison, the latter plot in Figure 4 depicts the circumstance on
2022-Jan-29. In contrast, It has the highest divergence so far in our
experiments.

Yet, both plots of different dates visually demonstrate that, compared to the
available training data (historical as well as recent observations), our
generated data could produce a smaller divergence from the tests. As the
pandemic has evolved so much since its emergence in 2019, its historical
course of observations (real training data) may become limited in capturing
the distribution of trends which we are trying to predict (test cases).

Figure 4: t-SNE plots of datasets on 2021-Jul-17 and 2022-Jan-29

Similar to the idea of [wang2009divergence_kNN], we estimate the KL divergence
using a $k$-nearest-neighbor (kNN) density estimator. In Equation 2, $z_{i}$
denotes every sample drawn from distribution $P$, $v_{i}$ is the spherical
volume of its $k$ nearest neighbors occupied in $Q$, and $n$ and $m$ are the
numbers of samples obtained from $P$ and $Q$, respectively.

$\displaystyle D_{KL}(P||Q)$
$\displaystyle=\mathbb{E}\Big{(}P(Z)log\frac{P(Z)}{Q(Z)}\Big{)}$ (1)
$\displaystyle\approx\frac{1}{n}\sum_{i=1}^{n}log\frac{P(z_{i})}{Q(z_{i})}$
(2) $\displaystyle\approx\frac{1}{n}\sum_{i=1}^{n}log\frac{1/n}{k/mv_{i}}$ (3)

In our analysis, $P$ corresponds to the distribution of the test data, which
is the target our models need to capture. Since every Saturday corresponds to
a new forecast with an updated model trained on the latest data, the test data
for each model normally contains just a few samples (i.e., one sample from
each state $i$), which are hardly formulated by a density function. We sought
to approximate it by a uniform weighting across its samples; therefore,
$P(z_{i})$ is simplified to $1/n$ in Equation (3), whereas $Q(z_{i})=k/mv_{i}$
is still estimated by its $k$-nearest-neighbors density.

Figure 5: Trend of the estimated divergence

Figure 5 shows a comparison in terms of the estimated divergence, from the
test data to the 1) historical data (training set), 2) recent data (usual
validation set), and 3) our augmented data. The densities of $Q$ are
approximated by the third nearest neighbor of the data points, i.e., $k=3$.
From the comparison, we can observe that the original training data has a
relatively larger divergence from the test set. We can see that our synthetic
data will likely be helpful for prediction because it has a smaller divergence
to the target data.

Please note that the negative values in the estimated divergence suggests that
our kNN approximation may produce a density $Q(z_{i})>1/n$. It means that the
third nearest neighbor may be so close to the query point (test data) that it
over-estimates the density around. Yet it is still a fair comparison, because
all datasets are estimated using the same set of $P(z_{i})$.

As a result, it shows that, over the whole period of time, the generated data
are closer to the test data, rather than the training or the usual validation
(recent) data. Consequently, it is expected that the generated data can be a
good supplement as an additional supply of training data.

####  1.3.2 Forecast Performance in Germany

##### Models included

Our results submitted to the German Hub are named HKUST-DNN_DA, which is
trained on the real data and validated using synthetic data. It will also be
aliased as HKUST-DNN_DA-val in our following discussion. We have also included
HKUST-DNN and HKUST-DNN_DA-train in this paper for comparison. They represent
our models trained without DA and using DA as training data, respectively.

Figure 6 shows a comparison of three of our trained models to other models
available in the German Hub. HKUST-DNN corresponds to our model without DA
(trained only on historical observations) and HKUST-DNN_DA-val is our model
trained using DA for validation. It is equivalent to our model HKUST-DNN_DA
submitted to the German Hub. HKUST-DNN_DA-train is our model trained using DA
as part of the training set, keeping the recent observations for validation.
Please note that the results of HKUST-DNN and HKUST-DNN_DA-train are only
available in our GitHub fork222https://github.com/cfong32/covid19-forecast-
hub-de of the German Hub.

Here, augmented data are generated using a 14-day forward ($F=14$) setting. To
supply data which is closer to the situation to be forecasted, only new
synthetic data which generated in previous seven days were included in each
training or validation.

Figure 6 shows the overall performance of our models compared to other models
submitted to the German Hub. They are measured in MAE and WIS, same as for the
US.

The models being compared are submissions from other research groups, such as
FIAS_FZJ-Epi1Ger from [barbarossa2020impact_hubde_FIAS] and USC-SIkJalpha from
[srivastava2020fast_alpha]. For details of these models, please refer to their
published meta data on the Hub [hubde]. To the best of our knowledge, ours are
the only deep learning models based on neural networks.

One of the Hub ensembles (KITCOVIDhub-median_ensemble) is also
compared.333Since the German Hub ceased publishing official KITCOVIDhub-
median_ensemble forecasts after 2021-Dec-06, we have continued calculating it
for subsequent dates. It is a simple quantile-wise median of all the submitted
models, including our HKUST-DNN_DA.

Comparing our own results, both our HKUST-DNN_DA-val (red) and HKUST-DNN_DA-
train (green) models could result in lower errors than the one without DA
(black). It shows that our augmented data could help our deep learning models
to capture and predict the trends. Especially on Jul-03, HKUST-DNN resulted in
a relatively large error due to overfitting to training data. With the help of
DA, in either validation or training, our models resulted in a much lower
error.

Compared with some other models in the Hub, especially during the plateaus in
2021 Sep and Dec, our model HKUST-DNN_DA-val was able to deliver more stable
performance.

Figure 6: Trends of forecast performance (MAE and WIS) in Germany

Table LABEL:tab:mae_rank_wis_de shows the overall MAE, WIS and rank of all the
models. In terms of MAE, HKUST-DNN_DA-val has achieved the second best among
all submissions. The average rank is calculated similarly to the US results.

table/tab_de_mae-rank-wis

###  1.4 Our DA on Other Deep Learning Models

####  1.4.1 Other Models

As a proof of effectiveness, we have tried to apply our generated data on two
other deep learning models. One of them is based on a simple LSTM encoder-
decoder architecture. Each of the encoder and decoder is simply a 2-layer
64-unit LSTM. Another model is the Adaptive Graph Convolutional Recurrent
Network (AGCRN) [bai2020adaptive_AGCRN]. It is a multivariate time series
prediction architecture proposed for traffic forecasting in road networks. By
applying the same quantile loss as described in Section LABEL:sec:loss-
function, on these two architectures, they are also able to produce
23-quantile outputs for each of the 16 states in Germany.

####  1.4.2 Results

Table 1 shows a comparison of our HKUST-DNN model with the two baselines with
and without DA, during the period from 2021-Jun-01 to 2022-Jan-01. With our DA
used for validation (DA-val) and saving the recent data for training, these
models are all improved with an overall reduction in MAE.

As the pandemic progresses, the number of cases sometimes can vary by orders
of magnitude. The overall MAE or WIS could be dominant by times when the
number of cases is high. Therefore, on top of the total reduction in
percentage, here we also calculate the per-date reduction which measures the
average of percentage reduction on different target dates. In terms of this
per-date MAE and WIS metrics, the performance of our model HKUST-DNN and AGCRN
can aslo be boosted by the augmented data.

However, the WIS evaluation seems to suggest that the improvement in
probabilistic forecasts (in 23-quantile format) is not as prominent as in the
point (median) forecasts. The smooth and simple generated data was not able to
promote performance for simple models like the standard LSTM architectures.

Table 1: MAE and WIS reduction using DA for validation MAE
---
model | no-DA | DA-val | Total reduction | Per-date reduction
HKUST-DNN | 3726.1 | 3349.6 | 10.1% | 11.3%
AGCRN | 3447.0 | 3298.7 | 4.3% | 3.5%
LSTM | 3412.2 | 3337.4 | 2.2% | -2.6%
WIS
model | no-DA | DA-val | Total reduction | Per-date reduction
HKUST-DNN | 2836.1 | 2901.6 | -2.3% | 6.1%
AGCRN | 2770.7 | 2633.3 | 5.0% | 4.5%
LSTM | 2669.5 | 3138.8 | -17.6% | -14.4%

###  1.5 Ablation Study on Dual-residual

To show the benefit of using our proposed dual-residual mechanism, we have
carried out the following ablation study.

Table 2 shows a comparison of our submitted models HKUST-DNN and HKUST-DNN_DA
with three other small ensemble models under different residual settings
($\alpha\in\\{0,0.5,1\\}$). The suffixes _R1, _R2 and _dualR indicate the
dual-residual settings of $\bm{R}^{1}$-only ($\alpha=0$), $\bm{R}^{2}$-only
($\alpha=1$) and balanced dual-residual ($\alpha=0.5$), respectively. These
three residual-specific models are in fact smaller ensembles combining five
composites of different random seeds, while the submitted ones are bigger
ensembles combining the best 10 of all these 15 composites. Please refer to
Section LABEL:sec:model-selection-and-ensemble for more details of the
selection process for the best.

To help focus on this comparison, the reported average ranks and standard
deviation of ranks are calculated among these four models only.

As HKUST-DNN and HKUST-DNN_DA are ensembles of the other models, no wonder
their average and variance of ranks are the lowest. Nonetheless, this
comparison clearly shows that balanced dual-residual has performed more
steadily (with low rank SD) and is on average better (with lower average rank)
than either one of the $\bm{R}^{1}$-only or $\bm{R}^{2}$-only setting.

Table 2: Comparison of difference residuals $cum\\_death$ forecast in the US
---
Model | Avg. rank↓ | Rank SD | MAE | WIS
HKUST-DNN | 2.43 | 0.83 | 147.3 | 109.6
HKUST-DNN_dualR | 2.47 | 1.06 | 142.6 | 103.5
HKUST-DNN_R1 | 2.51 | 1.22 | 169.7 | 124.8
HKUST-DNN_R2 | 2.59 | 1.30 | 139.3 | 107.3
$inc\\_case$ forecast in Germany
Model | Avg. rank↓ | Rank SD | MAE | WIS
HKUST-DNN_DA | 2.30 | 0.63 | 8555.8 | 7974.0
HKUST-DNN_DA_dualR | 2.44 | 0.87 | 7909.3 | 7357.6
HKUST-DNN_DA_R2 | 2.51 | 1.37 | 9000.9 | 8480.4
HKUST-DNN_DA_R1 | 2.74 | 1.37 | 8529.5 | 7597.1