# The Variations of SIkJalpha Model for COVID-19 Forecasting and Scenario
Projections

Ajitesh Srivastava  University of Southern California
<EMAIL_ADDRESS>
###### Abstract

We proposed the SIkJalpha model at the beginning of the COVID-19 pandemic
(early 2020). Since then, as the pandemic evolved, more complexities were
added to capture crucial factors and variables that can assist with projecting
desired future scenarios. Throughout the pandemic, multi-model collaborative
efforts have been organized to predict short-term outcomes (cases, deaths, and
hospitalizations) of COVID-19 and long-term scenario projections. We have been
participating in five such efforts. This paper presents the evolution of the
SIkJalpha model and its many versions that have been used to submit to these
collaborative efforts since the beginning of the pandemic. Specifically, we
show that the SIkJalpha model is an approximation of a class of
epidemiological models. We demonstrate how the model can be used to
incorporate various complexities, including under-reporting, multiple
variants, waning of immunity, and contact rates, and to generate probabilistic
outputs.

##  1 Introduction

Since the beginning of the COVID-19 pandemic, many models have been proposed
to predict the trajectory of the outcomes including cases, deaths, and
hospitalizations. To standardize the prediction tasks and communicate the
results to stakeholders, several collaborative efforts (hubs) were initiated.
These include US/CDC [16], Europe/ECDC [17] and Germany/Poland [20] forecast
hubs, and the US and the European scenario modeling hubs [4, 8]. The forecast
hubs ask for short-term predictions (up to 4 weeks ahead), while scenario
modeling hubs ask for long-term projections (months to years ahead) under
various scenarios. A collection of teams participate in these efforts, each
one independently producing predictions and submitting regularly – weekly for
forecast hubs, and once in 1-2 months for scenario modeling hubs.

During the early phase of the pandemic, we proposed the SIkJalpha model [42,
44] for cases and deaths forecasts (see Section 2.1). The approach has evolved
with time, depending on dataset availability and new factors that we believed
to have a significant impact. Further, for scenario modeling hubs, scenarios
are defined on a certain decision or a variable of uncertainty. To account for
the given scenario, that variable also needed to be included in the model. The
modeling philosophy has been attempting to design the simplest model that can
account for all such variables. We also avoid estimations where all parameters
are estimated simultaneously. In prior work, we have shown that simultaneously
estimating parameters of a simple non-linear version of our model can still
lead to overfitting [41]. Instead, we split the estimation problem into
multiple sub-problems, most of which are solved by linear regressions. Some
parameters are borrowed from the literature and some are provided by the
specifying scenario.

A key design feature of our approach is that all regions that a version of the
model is applied to receive the same treatment. There is no region-specific
manual tuning of the code. As an example, reporting rates differ significantly
between the states in the US – cases in Florida may be reported only once in
two weeks, while some other states may report on all weekdays or certain days
of the week. Yet, the same data pre-processing code applies to both. This is
done to improve the generalizability of the approach, i.e., the same code can
be applied to other regions in the world, as long as the data inputs the
version of the model needs are available for those regions. A second key
feature of the approach used for forecast hub submissions is that the
forecasts are generated without any human interventions. Occasionally,
hyperparameter-tuning is performed, when debugging or when new variables are
to be included creating a new version of the model. The relevant hyper-
parameters are presented in Table 1. The script is automated that generates
all forecasts for submission. As a result, sometimes, data anomalies can cause
unusual forecasts. However, we take the automated approach for scalability.
The philosophy is that the approach should scale to a large number of regions
if compute power and memory are available. Reliance on human assessment would
limit the scalability. Further, the method should stand on its own without
relying on the developer, and thus the open-source code can be used by anyone,
including a non-expert. While these are desirable design features, it is
crucial for the model forecasts to be accurate to be meaningful. Several works
have evaluated our model demonstrating that it has performed well (see Section
5).

###  1.1 Contributions

Our key contributions are the following. (i) We prove that the basic SIkJalpha
model can approximate a class of epidemiological models (Theorem 1). (ii) We
present various versions of the SIkJalpha model that have been used over time
over different forecast hubs (Section 2).(iii) We present the implementation
details including data sources, parameter estimation methods, runtimes, and
uncertainty quantification (Section 3 and 4). Note that, over time, many
modeling decisions were changed and new variables were included. The details
of the impact of each modeling decision are beyond the scope of this paper.
Instead, we focus on the breadth of modeling decisions. Evaluations of our
models are publicly available on our dashboard as well as on others’
dashboards, articles, and reports (see Section 5).

###  1.2 Related Work

During the COVID-19 pandemic, many models have been proposed for short-term
and long-term forecasts. For the models that participated in the Forecast and
Scenario Modeling Hubs, one way to categorize the model is the following.

Mechanistic models simulate the spread of COVID-19 through populations by
explicitly representing the transmission dynamics of the virus. These models
typically take into account factors such as the rate of contact between
individuals, the probability of infection given contact, and the duration of
the infectious period [46, 32]. These also include “meta-population” models
that incorporate various demographics, spatial resolutions, and interactions
between them [27, 31]. Agent-based models may be considered a special case
where each individual is modeled in the network of contacts [24].

Statistical models use historical data to predict future COVID-19 outcomes.
These include non-parametric spline-based models [47], and weighted regression
models that estimate growth rates [23]. A statistical model can be more
sophisticated than mechanistic models but they are also more dependent on the
quality and completeness of the data on which they are trained.

Machine Learning models are also data-driven, but they can be distinguished
from statistical models based on the fact they do not necessarily require an
explicit model of the relationship between the variables. Most such models
used some form of deep learning [48, 38] to ingest data from multiple sources
to predict future cases and mortality.

Ensemble models have also been used to generate forecasts and projections
based on multiple model outputs. Besides the “Hubs” that generate ensembles of
submitted models, some submitting models also generate an ensemble of their
own models. For instance, one submission combined auto-regressive models, Long
short-term memory (LSTM) models, ensemble Kalman filter, and an SEIR model
[18].

Scenario projections are more naturally performed by mechanistic models as the
desired real-world variable can be directly incorporated. However, there have
been Machine Learning-based models for scenario projections as well [25]. Our
approach is a combination of Mechanistic and Statistical category, as we model
susceptibility as a mechanistic process. Given the susceptibility, the
estimation of transmission parameters is done using a linear regression over
new infections.

##  2 Methodology

###  2.1 The Basic SIkJalpha Model

The key idea behind the SIkJalpha model [44] is to approximate the disease
dynamics as a discrete heterogeneous rate model. For instance, suppose the new
infections are created by the following true dynamics:

$\Delta I(t)=\frac{S(t)}{N}\mathcal{R}_{0}\sum_{\tau}\alpha(\tau)\Delta
I(t-\tau)=\frac{S(t)}{N}\sum_{\tau}\alpha^{\prime}(\tau)\Delta I(t-\tau)$ (1)

Here $S(t)$ is the number of susceptible individuals, $N$ is the population,
and $\mathcal{R}_{0}$ is the reproduction number [28] that represents the
expected number of infections created by each infected individual in a fully
susceptible population. The probability distribution of the serial interval,
i.e., the delay between an infected individual infecting another, is given by
$\alpha(\tau)$. As a result, someone who was infected at time $t-\tau$ on
average will infect $\alpha^{\prime}(\tau)=\mathcal{R}_{0}\alpha(\tau)$
individuals at time $t$.

Note that Equation 1 represents the new infections rather than active
infections which is common in variations of SIR models [19]. This choice has
been made as new infections (daily/weekly) are more easily observed (adjusting
for reporting rates and noise) compared to active infections. It can be shown
that by proper choice of $\alpha(\tau)$ we can mimic the SIR model.

######  Theorem 1.

The SIR model is a special case of SIkJalpha model.

###### Proof.

Suppose $\mathscr{I}(t)$ represents the number of active infections and
$\mathscr{R}(t)$ represents the number of recovered at time $t$ in the SIR
model. Then in the general model, total cumulative infections are the sum of
active and recovered infections, i.e., $I(t)=\mathscr{I}(t)+\mathscr{R}(t)$.
Therefore, new infections are the sum of new active infections and new
recoveries.

$\Delta
I(t)=\Delta\mathscr{I}(t)+\Delta\mathscr{R}(t)=\left(\beta\frac{S(t)}{N}\mathscr{I}-\Delta\mathscr{R}(t)\right)+\Delta\mathscr{R}(t)=\beta\frac{S(t)}{N}\mathscr{I}(t).$
(2)

Suppose $\chi$ is the recovery rate in the SIR model such that
$\Delta\mathscr{R}(t)=\chi\mathscr{I}(t)$. Also, the number of active
infections at time $t$ is the sum of all new infections in the past that have
not recovered by the time $t$. Therefore, under the SIR model:

$\displaystyle\mathscr{I}(t)$ $\displaystyle=\sum_{t^{\prime}=0}^{t}\Delta
I(t^{\prime})P(\mbox{infected at $t^{\prime}$ not recovered until $t$})$
$\displaystyle=\sum_{t^{\prime}=0}^{t}\Delta
I(t^{\prime})(1-\chi)^{t-t^{\prime}}=\sum_{\tau=0}^{t}\Delta
I(t-\tau)(1-\chi)^{\tau}\,.$ (3)

Putting this value in Equation 2

$\Delta
I(t)=\beta\frac{S(t)}{N}\mathscr{I}(t)=\beta\frac{S(t)}{N}\sum_{\tau=0}^{t}\Delta
I(t-\tau)(1-\chi)^{\tau}=\frac{S(t)}{N}\sum_{\tau=0}^{t}\left(\beta(1-\chi)^{\tau}\right)\Delta
I(t-\tau)\,.$ (4)

Setting $\alpha^{\prime}(\tau)=\beta(1-\chi)^{\tau}$, this is equivalent to
Equation 1. ∎

Figure 1: Demonstration of how the SIkJalpha approximates the distribution of
the serial interval with $k=4$ and $J=5$.

We approximate the dependence on past infections by $k$ temporal bins of size
$J$ as shown in Figure 1. Mathematically,

$\Delta I(t)=\frac{S(t)}{N}\sum_{\tau}\beta^{\prime}(\tau)\Delta
I(t-\tau)\approx\frac{S(t)}{N}\sum_{i=1}^{k}\beta^{\prime}_{i}\left(I\left(t-(i-1)J\right)-I(t-iJ)\right)\,.$
(5)

This enables the representation of arbitrary dependence with few learnable
parameters. Further, in real-world setting, obtaining reliable daily data can
be infeasible and may induce noise due to reporting patterns. Having $J>1$ has
a smoothing effect.

We note that only reported cases are observed, that can be represented by
$\Delta R(t)=\sum_{\tau}\gamma(\tau)\Delta I(t-\tau)$, which when combined
with Equation 5 results in a further approximation:

$\Delta
R(t)=\frac{S_{t}}{N}\sum_{i=1}^{k}\beta_{i}\left(R\left(t-(i-1)J\right)-R(t-iJ)\right)\,.$
(6)

Deaths (and hospitalizations) are considered to be outcomes of reported cases:

$\Delta D_{t}=\sum_{i=1}^{k_{D}}\theta_{i}(R_{t-(i-1)J_{D}}-R_{t-iJ_{D}}).$
(7)

The assumption of independence of deaths and hospitalization enables
independent learning of the parameters. Theoretically, the transition from
cases to hospitalizations, and that from hospitalizations to deaths are
conditionally independent, so they can be seen as a single transition from
cases to deaths. However, if the data on cases is of poor quality, it will
affect the prediction of deaths. In that scenario, it is better to estimate
deaths from hospitalizations rather than cases in Equation 7. The lag between
the reported hospitalizations and reported deaths gives a portion of a clearer
signal that is a better predictor of deaths.

The “alpha” in SIkJalpha represents the learning technique. We use a weighted
linear regression to emphasize recently seen data to adapt for dynamically
changing behavior. Suppose we have observed the reported cases for T time
steps. Then the regression is performed by minimizing the weighted least
square error

$LSE=\sum_{t=1}^{T}\alpha^{T-t}\left(\Delta R_{obs}(t)-\Delta
R_{fit}(t)\right)^{2}\,.$ (8)

Next, we present how various features were included and evolved over time. For
a summary of changes over time across the Hubs, please see Table 3.

###  2.2 True Infections

Since not all infections are reported, we used several strategies over time to
estimate the true number of infections. Our objective is to determine true new
infections at a given time.

Estimated.  We assume that there exists a reporting probability $\gamma$ for
each region, such that cumulative reported cases at time $t$ is given by
$R(t)=\gamma I(t-t_{0})$ for some constant $t_{0}$. In our prior work [41], we
showed that it is often not possible to reliably estimate under-reporting even
in simple models. Whenever it is possible (certain periods and certain
regions), we use the methods described in [41]. Then we assume that similar
reporting rates are observed in all other regions. We also assume that this
reporting rate is constant in the past and will remain constant in the future.
For short-term forecasts, in the early phases of the pandemic, the choice of
this parameter did not impact the performance [45].

Seroprevalence.  Wherever available (for US states), we use seroprevalence
data [12] that provides an estimated total number of individuals who had
infection-induced immunity at the frequency of once per 2-3 months. Each of
these points provides an estimate of reporting probability $\gamma(t)$. We use
linear interpolation to fill in the missing days within each of these 2-3
month periods. Finally, we take a weighted average to find the average
reporting probability:

$\gamma=\frac{\sum_{t}\gamma(t)\Delta R(t)}{R(T)},$ (9)

where $T$ is the length of the time series of reported cases. The
seroprevalence dataset also provides an interval for the percentage population
with infection-induced immunity using which we get an interval
$(\gamma_{low},\gamma_{high})$ of reporting probability. Note that this
assumes a constant reporting rate in the past, and we use the same rate to
project in the future. While a constant reporting rate in the past may not be
an accurate reflection of the truth, we consider it a reasonable modeling
choice for forecasting and scenario projection as it still provides a good
estimate of total infections. Also, note that the number of individuals with
infection-induced immunity is less than the number of total infections as one
individual may have multiple infections. However, before the rise of the
omicron variant, the number of reinfections was orders of magnitude smaller
[13] than uncertainty in true infections introduced due to the range
$(\gamma_{low},\gamma_{high})$. For the regions where seroprevalence data was
not available we used similar range for $(\gamma_{low},\gamma,\gamma_{high})$.

Wastewater.  Due to the low frequency of updates in seroprevalence data and
increase in reinfections, we switched to a wastewater dataset in early 2022
for US states [1]. Instead of a constant reporting rate assumption, we use the
following model

$\displaystyle\Delta I(t)$ $\displaystyle\propto C(t)=c_{0}C(t)\,$ (10)
$\displaystyle\gamma(t)$ $\displaystyle=I(t)/R(t)$ (11)
$\displaystyle\gamma(t)$ $\displaystyle=\gamma(T)\forall t>T\,,$ (12)

where $T$ is the last observed time and $C(t)$ is the average effective
concentration at time $t$. We also apply a cubic spline smoothing [37] on
$\gamma(t)$. To estimate $c_{0}$, we assume that seroprevalence data between
August 2020 ($t_{1}$) to July 2021 ($t_{2}$) provides a good estimate of
incident infections. Then,

$c_{0}=\frac{\sum_{t=t_{1}}^{t_{2}}\Delta
I(t)}{\sum_{t=t_{1}}^{t_{2}}C(t)}\,.$ (13)

This enables estimation of $\Delta I(t)$ over the entire period. While this
results in variable reporting probability in the past (Equation 11), for
forecasting and projections, we still use a fixed reporting rate based on most
recent rate (Equation 12).

###  2.3 Variants

During the emergence of the Alpha variant, we incorporated it into the model.
Our general approach is to split the cases into multiple time series
representing the cases caused by each variant. While genomic sampling data is
available that can provide proportions over time, a major challenge is that
the sampling data is updated with a lag of 1-2 weeks [2]. Instead of using the
variants proportions directly, we fit a model to identify the proportions as a
function of time. This enables us to fill in for the days with missing data
and removes the noise due to sampling. Specifically, we use the last 2-3
months of data to fit a logistic model, so that:

$\displaystyle\log\frac{\Delta I^{(i)}(t)}{\Delta
I^{(1)}(t)}=s_{i}t+c_{i}\,,\sum_{i}\Delta I^{(i)}(t)=\Delta I(t)\,,$ (14)

where $\Delta I^{(i)}$ represents the infections created by the variant $i$.
We may interpret $s_{i}$ as the advantage in the growth rate of variant $i$
over the base variant $1$. The model assumes a constant growth rate advantage
within the short period under consideration. In the long-term, the variants
dynamics can significantly deviate from the logistic model. However, in the
short term, it can be shown that the deviation is small. A detailed analysis
of the errors due to this violation is beyond the scope of this paper. Various
implementations were used over time:

2-strain.  In this model, we only focus on two strains - two most dominant
strains, or the current dominant strain, and one strain known to be rapidly
increasing in prevalence. E.g., wild-type vs Alpha, and non-Omicron vs
Omicron. A disadvantage of this model is that if there are more than two
variants with significantly different transmission rates, then forcing them
into two groups may lead to variable growth advantage deviating from the
logistic model, even in the short term. Another drawback is the human
intervention requirement – two strains need to be manually selected and
updated when a new highly transmissible variant appears.

All PANGO.  In this model, we use all the variants that appear as different
PANGO lineages [14]. This approach does not require manual selection of
variants and automatically includes new variants that appear as the dataset
[9, 11] is updated. However, this approach also has some drawbacks as the
number of different lineages active in a 3-month window can be of the order of
a few 100s. As a result, the computation time increases. Also, distributing
the samples across so many lineages results in higher uncertainty due to the
small number of samples per lineage.

Selected.  In this model, we take the list of all PANGO lineages, and filter
it to select those of interest. For example, during the Omicron wave, we
extracted all lineages with names starting with “BA”. All other lineages that
appeared before omicron are grouped into one variant called “pre-omicron”.

###  2.4 Vaccines and Waning Immunity

Initially, a simple all-or-nothing model model was used, where the acquired
immunity (through vaccines or prior infection) against variant $i$ provides
full protection to $(1-\epsilon_{A}(i,j))$ fraction of the total population
$N_{j}$ in age group $j$ at time $t$, and no protection to the rest. As a
result, the susceptible population is given by

$S_{j}(t)=N_{j}-I_{j}(t)-\epsilon_{A}(i,j)V(t)\left(1-I_{j}(t)/N_{j}\right)\,.$
(15)

This assumes that infection provides full immunity and those who are
vaccinated but not infected are susceptible with probability
$\epsilon_{A}(i,j)$. The value of the parameter $\epsilon_{A}(i,j)$ was
borrowed from the vaccine trials [36, 4]. For simplicity, we assumed
$\epsilon_{A}(i,j)=\epsilon_{A}(j)\forall i$. However, since July 2021, a
time-waning model was used in accordance with the discussions with the
Scenario Modeling Hub. In this model, immunity is lost with time based on some
function $w_{j}^{(i)}(\tau)$, where $\tau$ is the time since
vaccination/infection. The expected susceptible population is given by

$S_{j}(t)\approx
N_{j}-I_{j}(t)-V(t)\left(1-\frac{I_{j}(t)}{N_{j}}\right)+\left(\sum_{\tau}w_{j}(\tau)\left(\Delta
V_{j}(t-\tau)+\Delta I_{j}(t-\tau)\right)\right)\,.$ (16)

From Aug 2021 - Oct 2021, the function $w_{j}^{(i)}$ was assumed to be a
function representing the expected susceptibility based on the following
model: (i) An individual transitions to a “partially susceptible state at a
time given by exponential distribution with mean $\lambda$. (ii) In this
state, the individual has a residual protection $p_{inf}$ against infection
and $p_{death}$ and $p_{hosp}$ against death and hospitalization,
respectively. This means that they will only be susceptible with probability
$1-p_{inf}$. The parameters were determined by the scenarios in Scenario
Modeling Hub. Since then, we have switched to using a gamma distribution to
represent the transition into the partially susceptible state. The gamma
distribution requires two parameters to be specified. These are obtained from
two equations: 1) The mean/median of transition probability should match the
specification of the scenarios. 2) The expected immunity after two months
should match the vaccine trial results obtained at two months.

Note that both Equations 15 and 16, are approximations in that they only
support one breakthrough of acquired immunity and that total infections is
roughly equal to the total individuals infected. With increasing reinfections,
this approximation does not hold, and a more sophisticated model is used later
(see Section 2.8).

###  2.5 Vaccine Adoption

During the early phase of vaccinations, from Jan 2021 to July 2021, a linear
approximation was used for short-term forecasts (up to 4-weeks ahead). We use
linear extrapolation determined by the last two weeks of new vaccinations to
get the future new vaccinations. For longer-term projections, a sigmoid model
was used:

$V(t)=\frac{a}{1+\exp(-b(t-c))}\,$ (17)

where $b$ and $c$ are learnable parameters based on early observed data. The
parameter $a$ is fixed based on vaccine adoption surveys.

When a few month of data was available, we modeled vaccine adoption as a
contagious behavior:

$\displaystyle\Delta V(t)$
$\displaystyle=\frac{E_{V}(t)-V_{E}(t)}{N}\sum_{i=1}^{k}\beta_{i}\left(V(t-(i-1)J)-V(t-iJ)\right)$
(18)

where $V(t)$ is the cumulative number of people vaccinated by time $t$,
$E_{V}(t)$ is the number of individuals eligible for the vaccine at $t$, and
$V_{E}(t)\leq V(t)$ is the number of individuals who are eligible at time $t$
but are already vaccinated. The above model may be applied to different age
groups independently. For the first dose, everyone may be eligible in a given
age group. We first predict the future-first doses. This determines the
eligibility for the second dose, which in turn determines the eligibility for
the booster. In the implementations so far, $V_{E}(t)=V(t)$ as the eligibility
is cumulative - someone who is eligible at $t$ will remain eligible at
$t^{\prime}>t$.

The hyperparameters are set to $k=2,J=7$. The parameters $\beta_{i}$ are
estimated using a weighted regression similar to the one described in Section
2.1 with a weight of 0.9.

To simulate target adoption scenarios, such as targeting a fraction $u$ of
eligible population, a scaling function $\alpha(t)$ is introduced, after
estimating $\beta_{i}$ in Equation 18.

$\displaystyle\Delta V(t)$
$\displaystyle=\frac{uE_{V}(t)-V_{E}(t)}{N}\sum_{i=1}^{k}\alpha(t)\beta_{i}\left(V(t-(i-1)J)-V(t-iJ)\right).$
(19)

We pick $\alpha(t)=\exp(\frac{H}{H-t})$, where $H$ is the time before which
the target adoption is to be reached. The impact of this function is that it
smoothly and rapidly increases the rate with time, thus making the term
$\frac{uE_{V}(t)-V_{E}(t)}{N}$ decisive in the saturation of adoption.

###  2.6 Age-groups

When age-specific targets (cases, deaths, and hospitalizations) are available,
we incorporate interactions between age-groups in a “contact matrix”. We
assume different age groups may have different susceptibility and there is a
transmission rate that is scaled by the contact matrix $C=[c_{g,g^{\prime}}]$,
where $c_{g,g^{\prime}}$ represents a relative measure of contact rates
between age groups $g$ and $g^{\prime}$. In other words, $c_{g,g^{\prime}}$
represents the impact of age-group $g^{\prime}$ on $g$. For $G$ age groups,
this may introduce $G(G-1)/2$ parameters if we assume symmetric impacts. To
avoid overfitting, we reduce this to $G$ parameters, assuming that for some
parameters $c_{1},c_{2},\dots,c_{G}$,
$c_{g,g^{\prime}}=c_{g}{c_{g}^{\prime}}$. One way to interpret this is to
assume that age-group $g$ independently participates in a randomly selected
contact with probability $\propto c_{g}$. Therefore, probability of an
individual in age group $g$ being in contact with one in $g^{\prime}$ is
proportional to $c_{g}c_{g^{\prime}}$. This is similar to the formulation of
the contact matrix in Ma et. al. [34].

Since contact is a physical process, it is assumed to be the same for all
variants. We rewrite Equation 6 to incorporate age-groups as:

$\Delta
R^{(i)}(t,g)=\frac{S^{(i)}(t,g)}{N(t,g)}\sum_{l=1}^{k}\sum_{g^{\prime}=1}^{G}\beta_{l}^{(i)}c_{g}c_{g^{\prime}}\left(R^{(i)}(t-(l-1)J,g^{\prime})-R^{(i)}(t-lJ,g^{\prime})\right)\,.$
(20)

To fit the model to the data, first, we fix the variant $i$ by choosing the
one which has the most cumulative cases in the recent weeks (in our
implementation, we used the last four weeks). Then, we set $c_{G}=1$, for the
last age group, and learn $\beta_{l}^{(i)}$ and $c_{1},\dots,c_{G-1}$. Then,
we fix the learned values of $c_{g}$ and learn $\beta_{l}^{(i)}\,\forall l,i$.
The purpose of setting $c_{G}=1$ is to avoid multiple equivalent solutions
obtained by arbitrarily scaling up all $c_{g}$ and scaling down
$\beta_{l}^{(i)}$.

As an example, the following matrices are obtained from this method based on
data as of June 6, 2021 for California and Washington, respectively.

$\begin{bmatrix}0.88&0.88&1.15&2.13&0.94\\\ 0.88&0.89&1.16&2.15&0.94\\\ 1.15&1.16&1.51&2.80&1.23\\\ 2.13&2.15&2.80&5.19&2.28\\\ 0.94&0.94&1.23&2.28&1.00\end{bmatrix}$ |  | $\begin{bmatrix}1.02&1.01&1.32&2.40&1.01\\\ 1.01&1.01&1.32&2.40&1.01\\\ 1.32&1.32&1.73&3.13&1.31\\\ 2.40&2.40&3.13&5.66&2.38\\\ 1.01&1.01&1.31&2.38&1.00\end{bmatrix}$
---|---|---

The age groups considered here are 0-4, 5-11, 12-17, 18-65, 65+. The entry in
the $i^{th}$ row and $j^{th}$ column represent the relative contact rate
between age groups $i$ and $j$. The matrix is normalized so that the contact
rate within the oldest age group is 1. Figure 2 shows the fitted time-series
by age-group using these matrices. Since the contact matrix and the infection
rates vary over time, these fits are generated by providing a higher weight to
the more recently seen data. Therefore, the method shows a better fit towards
the end of the calibration period.

Figure 2: Result of fitting age-specific infection trends by learning contact
matrix and infection rates. The markers indicate the pre-processed age-wise
data and the solid lines indicate the fitted values. More weight is given to
the end of the calibration period to identify recent rates.

###  2.7 Immune Escape Properties

To incorporate immune escape variants like Omicron, a “cross-protection”
matrix $C_{P}$ was introduced, such that $C_{P}(i,i^{\prime})$ represents the
protection against variant $i$ due to an infection by variant $i^{\prime}$.
This is the protection before waning of immunity takes place. As a result,
each variant has its own view of the population. In the absence of vaccines
and waning immunity, susceptibility to variant $i$ will be given by:

$S^{(i)}(t)=\sum_{i^{\prime}}\sum_{t^{\prime}}(1-C_{P}(i,i^{\prime}))I^{(i^{\prime})}(t)$
(21)

The matrix $C_{P}$ is populated based on the scenario requirements. For
instance, if Omicron is required to have an immune escape of 80%, then
$C_{i,i^{\prime}}=1-0.8=2,$ where $i$ is a lineage of Omicron and $i^{\prime}$
is a lineage of past variants. Otherwise, $C_{i,i^{\prime}}=1$. Additionally,
a vector $\mathbf{r}_{VE(j)}$ is introduced such that the $i^{th}$ element
represents the initial protection against variant $i$ from the vaccine type
$j$. In absence of any waning immunity and infection,

$S^{(i)}(t)=\sum_{j}(1-\mathbf{r}_{VE(j)}(i))V_{j}(t)$ (22)

where $j$ is a vaccine type (first dose, 2nd dose, booster, …) and $V_{j}(t)$
is the total individuals with vaccine type $j$. Exact computation of
susceptibility across variants in presence of waning immunity requires keeping
track of various states and is discussed in Section 2.8.

###  2.8 All-state Model

Figure 3: The “all-state” version of the SIkJalpha model. The transmission is
determined by those who were reported infected $1,\dots,k$ time-steps ago in
groups $1,\dots,g$ (denoted by $R_{i,j}$) to create the new reported cases
$R$. The values of $R_{i,j}$ also determine the new deaths $D$ and any other
target ($X$) such as hospitalizations. The states on the right indicate
various susceptibility states representing the combination of prior infection
and vaccines, time of these events, variant that infected them, and their age
group.

To incorporate multiple reinfections, and various susceptibility groups that
arise due to multiple rounds of vaccine, we define the following types of
state: first-time infections ($I_{1}$), two or more infections ($I_{2}+$),
first dose ($F$), infection with a prior first dose ($I_{f}$), first dose
after a prior infection ($F_{i}$), second dose before boosting ($VnB$),
infected with a prior second dose ($I_{v}$), second dose after a prior
infection ($VnB_{i}$), booster dose ($B$), infected with a prior booster dose
($I_{b}$), and booster dose after a prior infection ($F_{i}$). Figure 3 shows
the allowed transitions between states. Let state $0$ denote those who are not
in any of the above states (naive state). This is the initial state of
everyone in the population before they get their first infection or vaccine.

The state of an individual is determined by the $4$-tuple $(j,i,g,t)$
determined by the state-type $j$, age group $g$, time $t$ and variant $i$. Let
$N(j,t,i,g)$ represent the number of individuals in state $(j,i,g,t)$. The
association with the variant index represents the last variant that affected
that individual. For instance, $N_{I2+}(g,i,t)$ would represent that the last
transition they made on or before time $t$ was to a reinfection state due to
variant $i$. Those in the vaccination-only state, i.e., $F$, $VnB$, and $B$
exist for all variants, i.e., $N_{j}(g,i,t)=N_{j}(g,i^{\prime},t)\forall
i,i^{\prime}$ when $j\in\\{F,VnB,B\\}$. There are at most two ways of
transitioning from a state either due to an infection or due to vaccination.

Infection-induced transition.  At $t=0$, $N_{j}(g,i,t)=0,\forall j,g,i,t$.
Given how immunity wanes over time, we can compute the expected number of
susceptible individuals in each state $(j,g,i^{\prime},t^{\prime})$ as
perceived by each variant $i$ at time $t$.

$\displaystyle S_{i,t}(j,g,i^{\prime},t^{\prime})$
$\displaystyle=N(j,g,i^{\prime},t^{\prime}),j\in\\{0\\}$ (23) $\displaystyle
S_{i,t}(j,g,i^{\prime},t^{\prime})$
$\displaystyle=N(j,g,i^{\prime},t^{\prime})\bigg{(}1-\mathbf{r}_{VE(j)}(i)\left(1-w_{j}(i,t-t^{\prime})\right)\bigg{)}\,,j\in\\{F,VnB,B\\}$
(24) $\displaystyle S_{i,t}(j,g,i^{\prime},t^{\prime})$
$\displaystyle=N(j,g,i^{\prime},t^{\prime})\bigg{(}1-C_{p}(i,i^{\prime})\left(1-w_{j}(i,t-t^{\prime})\right)\bigg{)}\,,j\in\\{I1,I2+\\}$
(25) $\displaystyle S_{i,t}(j,g,i^{\prime},t^{\prime})$
$\displaystyle=N(j,g,i^{\prime},t^{\prime})\bigg{(}1-\max\\{C_{p}(i,i^{\prime}),\mathbf{r}_{VE(j)}(i)\\}\left(1-w_{j}(i,t-t^{\prime})\right)\bigg{)}\,,$
$\displaystyle j\in\\{I_{f},I_{v},I_{b},F_{i},VnB_{i},B_{i}\\}$ (26)

In the naive state, everyone is susceptible (Equations 23). In a state-type
with vaccination but no infection, probability of immunity is given by the
initial vaccine efficacy and lack of waning immunity (Equations 24). In a
state-type with prior infection but no vaccine probability of immunity is
given by cross-protection and lack of waning immunity (Equations 25). In a
state-type with prior infection and vaccination the probability of immunity is
given by the maximum of cross-protection and the initial vaccine efficacy, if
the protection has not waned (Equations 26).

Let $\Delta I^{(i)}(t,g)$ be the new infections at time $t$ in age group $g$
due to variant $i$. Then the new infections due to this variant in this age-
group that will transition from state $(j,g,i,t)$ can be estimated as

$\Delta_{I}N_{t}(j,g,i^{\prime},t^{\prime})=\Delta
I^{(i)}(t,g)\frac{S_{i,t}(j,g,i^{\prime},t^{\prime})}{\sum_{j^{\prime},t^{\prime}}S_{i,t}(j^{\prime},g,i^{\prime},t^{\prime})}$
(27)

where $\Delta_{I}$ represents the change due to infections at time $t$.

Vaccination-induced transition.  In infection induced transitions we
calculated susceptibility in each state to identify where the transitions are
coming from. Similarly, in vaccine induced transitions we need to compute
eligibility to identify where the transitions are coming from. The eligible
population to go through vaccine-induced transition from the state
$(j^{\prime},g,i,t)$ due to vaccine type determined by $j$ (i.e., $F,VnB,B$),
at time $t$ is given by:

$\displaystyle E_{j,t}(j^{\prime},g,i,t^{\prime})$
$\displaystyle=N_{t}(j^{\prime},g,i,t^{\prime}),\mbox{ s.t.,
}j^{\prime}\rightarrow j,\,t^{\prime}<t-\delta(j^{\prime}\rightarrow
j),j\in\\{F,VnB,B\\}$ (28)

where, $j^{\prime}\rightarrow j$ represents a valid transition, and
$\delta(j^{\prime}\rightarrow j)$ denotes the minimum delay allowed between
the transitions. For instance, $\delta(j^{\prime}\rightarrow j)$ may be set to
six months for transitioning from 2nd dose only state ($VnB$) to the “boosted
without prior infections” state $B$. Let $V_{j}(t,g)$ denote the number of
vaccinations of type $j$ (first dose, second dose, or booster) given at time
$t$ to the age-group $g$, then the population in state
$(j^{\prime},g,i,t^{\prime})$ transitioning out induced by vaccinating
$V_{j}(t,g)$ individuals at time $t$ is

$\Delta_{V}N_{t}(j^{\prime},g,i,t^{\prime})=V_{j}(t,g)\frac{E_{j,t}(j^{\prime},g,i,t^{\prime})}{\sum_{j^{\prime},t^{\prime}}E_{j,t}(j^{\prime},g,i,t^{\prime})}$
(29)

Finally, the change in number of individuals in each state is given by

Initially, $\displaystyle N_{0}(0,g,i,t^{\prime})=N(0,g)\forall
g,i,t^{\prime}$ (30) $\displaystyle N_{0}(j,g,i,t^{\prime})=0\,\forall
j>0,g,i,t^{\prime}$ (31) Transitions to state $j$ at time $t$ $\displaystyle
N_{t}(j,g,i,t)=\sum_{i^{\prime},t^{\prime},j\rightarrow
j^{\prime}}\Delta_{I}N_{t}(j^{\prime},g,i^{\prime},t^{\prime})+\sum_{i^{\prime},t^{\prime},j\rightarrow
j^{\prime}}\Delta_{V}N_{t}(j^{\prime},g,i^{\prime},t^{\prime})$ (32)
Transitions coming from state $j$ at time $t$ $\displaystyle
N_{t}(j,g,i,t^{\prime})=N_{t-1}(j,g,i,t^{\prime})-\sum_{i^{\prime},t^{\prime},,j\rightarrow
j^{\prime}}\Delta_{I}N_{t}(j^{\prime},g,i^{\prime},t^{\prime})-\sum_{i^{\prime},t^{\prime},,j\rightarrow
j^{\prime}}\Delta_{V}N_{t}(j^{\prime},g,i^{\prime},t^{\prime})$ (33)

$\displaystyle S_{j}^{(i)}(t,g)$
$\displaystyle=\sum_{i^{\prime},t^{\prime}}S_{i,t}(j,g,i^{\prime},t^{\prime})$
(34)

As before, we assume that all severe outcomes at time $t$ in age-group $g$
($X_{j}^{(i)}(t,g)$), including hospitalizations and deaths are simply linear
functions of the temporal sequence of infections.

$\Delta
X_{j}^{(i)}(t,g)=\sum_{l=1}^{k_{X}}\chi_{l}^{(i)}(g)(1-P_{j}^{(i)}(t,g))\left(I_{j}^{(i)}(t-(l-1)J,g)-I_{j}^{(i)}(t-lJ_{X},g)\right).$
(35)

By this decoupling of the dependency, instead of sequential dependency between
cases to hospitalizations to deaths, we are able to independently produce
estimate parameters. Here $\chi_{l}^{(i)}(g)$ are rates of severe outcome for
a “naive” individual - one who has no protection. We assume that the ratio of
severity of any two variants is a constant. Therefore, we set
$\chi_{l}^{(i)}(g)=\rho^{(i)}\chi_{l}(g)$, where $\rho^{(i)}$ represents
relative severity of variant $i$. This parameter makes a significant
difference only if there exist variants with significantly different severity
in the last few days used for training ( 100 days). For example, in the US,
this parameter was used from December 2021 to March 2022, when the last 100
days were expected to have significant proportions of both omicron and delta
variants. We treat $\rho^{i}$ as a known parameter when provided in the
scenarios [4]. Otherwise, we perform a grid search to identify the $\rho^{i}$
that best fits a held-out validation data (new deaths/hospitalizations) of the
last week. When all the circulating variants are expected to have the same
severity, we set $\rho^{(i)}=1\forall i$.

The parameter $P_{j}^{(i)}(t,g)$ is the protection against severe outcome of
an individual in state-type $j$ in age group $g$ at time $t$ from variant $i$.
We treat $P_{j}^{(i)}(t,g)$ as known parameters borrowed from the provided
scenarios [4] for the scenario modeling hub. In absence of recommended values
of protection (for forecast hubs), we use the various values suggested by the
scenario modeling hub one at a time as “sampled” hyperparameters to generate
“sub-scenarios” (see Section 3). It should be noted that by 2022, the naive
population (no past vaccinations or infections) had become a small fraction of
the population as suggested by the seroprevalence data [12]. Therefore, this
parameter has a negligible impact on the projections. Since we do not observe
severe outcomes per variant, we perform regression on:

$\Delta
X(t,g)=\sum_{l=1}^{k_{X}}\chi_{l}(g)\left(\sum_{i,j}\rho^{(i)}(1-P_{j}^{(i)}(t,g))\left(I_{j}^{(i)}(t-(l-1)J,g)-I_{j}^{(i)}(t-lJ_{X},g)\right)\right).$
(36)

The quantity in the parenthesis is computed before the regression resulting in
a linear regression only in $\chi_{l}(g)$, thus producing extremely fast
results free from overfitting common in non-linear settings.

###  2.9 Other Capabilities

In addition to the above-described factors that are accounted for, the model
also supports other capabilities as described below.

Anomaly detection and Smoothing.  To address the reporting noise and delayed
dumps of backlogged cases and deaths, various strategies were used. Until July
2020, only moving average smoothing was used. However, delayed dumps caused
large spikes, thus misleading the model. Therefore, an anomaly detection
module was added which considered a daily value (case or death) to be an
outlier if it exceeded 4 times the median value of the last few ($\sim$70)
days. The outlier is then replaced with the linearly interpolated value from
the neighboring non-outliers. Then we apply a moving average smoothing. With
changes to reporting schedule (once a week, reporting only on certain
weekdays, etc.) the anomaly detection module had to be changed. We first
accumulate the reported targets into a weekly time series instead of daily.
The anomaly detection is then applied to this weekly time series. Then we
convert this time series back into daily time series by distributing the
reports in a week equally among its days. Finally, we apply to smooth on this
time series.

Future Contact Behavior.  Future changes in Non-Pharmaceutical Interventions
(NPIs) for scenario modeling are implemented as linear scaling of contact
matrix 2.6. In absence of age stratification, this results in a linear scaling
of the transmission rates in the future. In earlier rounds, when NPIs were
asked to be lifted linearly, the scaling was changed linearly over time. In
rounds where future NPI modeling was left up to the teams, the scaling was
based on contact rates of the previous year obtained from Cuebiq [5] to
account for seasonal contact changes.

Arbitrary variant.  Arbitrary variants in the future can be incorporated into
the model along with the existing ones. The model requires as input the day of
introduction (or distribution over days), the number of introductions, and
properties of the new variant (transmissibility and immune escape advantages).

More vaccine doses.  While we present the model with three doses (first,
second, and booster), the model supports an arbitrary number of doses. This is
done by repeating Equations 24, and 28 to more vaccine doses beyond
$\\{F,VnB,B\\}$.

Learnable Parameters (estimated by regression) |  Transmission rates vector per variant $i$ $\beta_{l}^{(i)}$, contact matrix $C$, severe outcome rates $\chi_{l}(g)$ for deaths and hospitalizations
---|---
Selected Hyperparameters (selected by grid search) |  For cases $k=2,J=7,\alpha\in\\{0.9\\}$, for deaths $k_{D}\in\\{3,4,\dots,7\\},J_{D}=7,\alpha_{D}\in\\{0.95,0.98\\}$, for hospitalizations $k_{H}=\\{2,3\\},J_{H}=\\{2,3,7\\},\alpha_{H}\in\\{0.95,0.98\\}$, relative variant severity $\rho^{(i)}$
Provided/Sampled Hyperparameters (sampled from a reasonable set) |  under-reporting/true cases interval, interval estimate of variant proportions (affecting $\Delta I^{(i)}(t)$), waning immunity parameters $\epsilon_{A}(i,j),w_{j}(i,\tau)$, protection against severe disease $P^{(i)}_{j}(t,g)$, Relative vaccine efficacy per variant per group $\mathbf{r}_{VE(j)}(i)$, cross protection matrix $C_{P}$, future vaccine coverage $u$
Table 1: Categorization of parameters involved in model training.

##  3 Model Uncertainty

Two methods have been used to model uncertainty: sub-scenarios and Random
Forest. The latter has only been used for the Germany/Poland Forecast Hub.

###  3.1 Sub-scenario Based

Here, our goal is to generate multiple future trajectories from which we can
sample desired quantiles. We divide the parameters involved in the model into
three categories (Table 1): (i) Learnable parameters include those that are
estimated using regression. (ii) Selected hyperparameters are those that are
identified via grid search based on a held-out validation set (combination of
hyperparameters that best predicts the unseen data), and (iii)
Provided/Sampled hyperparameters are those that are either provided by the
scenarios or are considered to take one of several different values. Among the
provided/sampled hyperparameters that can take several values, each setting
corresponds to a “sub-scenario”. Once we fix a sub-scenario, learnable
parameters can be estimated and selected hyperparameters can be identified by
grid-search, leading to one trajectory. We form all possible combinations of
sub-scenarios thus typically getting 50-300 trajectories. The desired
quantiles are then sampled from these trajectories.

###  3.2 Random Forest

In this approach, we first generate the “mean” prediction. This is the time-
series outcome of regression which we convert to a vector of four elements
$(y_{1},y_{2},y_{3},y_{4})$. The $i^{th}$ element represents the $i$-week
ahead forecast. We assume that the true observation $i$-week ahead is
$\bar{y}_{i}=y_{i}+e_{i}$. For each $i$, to estimate $e_{i}$, we train a
Random Forest [35] with 100 trees. The inputs to the Random Forest are current
and last week’s incident data (cases or deaths), the population of the region,
and the prediction from SIkJalpha, $y_{i}$. The output is the deviation from
the true observed data $\bar{y}_{i}-y_{i}$. The training is performed by
considering recently seen (last 35 days) data, and the difference between
prediction and observed data over all the regions. The desired quantiles are
then sampled from the Random Forest.

##  4 Real-time Prediction

###  4.1 Participation in Real-time Prediction Efforts

Various versions of SIkJalpha model were submitted to US/CDC Forecast Hub (US
FH), Europe/ECDC Forecast Hub (Euro FH), Germany/Poland Forecast Hub (G-P FH),
US/CDC Scenario Modeling Hub (US SMH), and Europe/ECDC Scenario Modeling Hub
(EU SMH). Submissions to the forecast hubs were made every week on Sunday.
Scenario Modeling Hubs have submissions once in 1-2 months. The details of
targets and modeling choices are presented in Table 3.

###  4.2 Real-time Publicly Available Forecasts

We have maintained an online dashboard since April 2020 [40] to visualize the
generated forecasts for US states and countries around the world. These
forecasts are generated using the same modeling choices as mentioned in Table
3 for US FH and Euro FH, respectively. Since December 2020, we have also added
forecasts for all the locations around the world [39] for which Google
provides case and death data [3]. This includes 214 Admin-0 locations, around
1000 Admin-1 locations, and around 16,000 Admin-2 locations. This
implementation includes modeling choices identical to G-P FH with the
exception of the inclusion of 2-dose vaccines which are assumed to provide
all-or-nothing protection. When vaccination data is not available from Google
for Admin-1 or Admin-2 locations, we assume that the vaccination uptake per
unit population is the same as that of the higher Admin-level for which data
is available. We do not generate uncertainty bounds for these forecasts.

###  4.3 Implementation

All code was written in MATLAB and is publicly available on Github [40]. Daily
forecasts have been generated through an automated script at 5 a.m. Pacific
Time every day since May 2020 on an Intel 2-core desktop. The scenario
projections were developed and generated on the same machine as per the
submission schedule until September, after which we moved to a Dell PowerEdge
R540 server. It has two Intel Xeon Gold 5218 processors each running @2.3G,
16C/32T, 10.4GT/s, 22M Cache (total 32 core/64 threads, 44M Cache). It has
64GB RDIMM, 3200MT/s, Dual Rank RAM, and 1.92TB SSD SATA Read Intensive 6Gbps
drive. MATLAB implicitly converts many matrix operations and functions into
multi-threading. Starting in September, we also introduced parallelism across
sub-scenarios. The runtimes of the latest implementation are presented in
Table 2. The most time-consuming aspect is the estimation of the susceptible
population which requires evaluating all states from the beginning of the
epidemic, which in our model is set to Jan 23, 2020. To accelerate the
computation, we update the states at a weekly granularity rather than daily.

Targets | “all-state”? | Runtime/subscenario
---|---|---
US states / US FH | Yes | 25s on server
$\sim$190 countries | Yes |  $\sim$200s on desktop
Euro FH/SMH | Yes | 15s on server
GP FH | No |  $<$ 10s on desktop (incl. quantiles)
$>$17,000 locations (Google) | No |  $\sim$400s on desktop
Table 2: Approximate runtimes of various implementations.

##  5 Evaluations

(a) US Forecast Hub - cases

(b) US Forecast Hub - deaths

(c) ECDC Forecast Hub - cases

(d) ECDC Forecast Hub - deaths

Figure 4: Comparison of our submitted 4-week against Hub ensemble in terms of
MAE

(a) Cases in the US states

(b) Cases in the European countries

Figure 5: Reported cases after applying the same smoothing to both datasets:
States in the US have smoother trends than the European countries. Figure 6:
The model often over-predicts just before the surge: Here the purple line is
the model prediction for California before the Omicron surge. The grey line is
the ground truth. Due to the inclusion of variants-specific data, the model
can predict an impending surge, however, it over-predicts the intensity of the
surge.

Our models have been evaluated by us [43] as well as several other studies
[20, 26, 22, 29, 30], reports [15, 7, 10] and dashboards [6]. Despite the lack
of human involvement in generating forecasts (except for debugging and adding
new features), our methods have been in the top 25% during most of the
pandemic based on mean absolute error and weighted interval score [21]. For
case forecasts in the US, among all models that were submitted consistently in
2020-2021, our submission ranks highest in terms of Weighted Interval Score
[33]. For country-level death forecasts over multiple countries, our forecasts
performed the best in terms of median absolute percentage error [29] among 7
models that fit the inclusion criteria among 386 models analyzed. This
analysis included long-term forecasts up to 12 weeks ahead generated in 2020
to mid 2021.

However, the US hospitalization forecasts have not been good. This is
particularly because it was assumed that there is at least one one-week lag
between the reported cases and hospitalizations. However, this was not
consistently true and in the majority of the weeks, the reported cases were
well-aligned (without a lag) with reported cases due to delay in reporting. It
was fixed in February of 2022, by removing the enforced one-week lag. This is
a drawback of the “automated” approach and could have been fixed earlier with
more human involvement – by visualizing the results every week. However, due
to the involvement of the lone author in multiple collaborative efforts, this
was not feasible and would have required multiple personnel.

Figure 4 shows the comparison of our submissions and the Hub ensemble in terms
of mean absolute error for US and ECDC 4-week ahead case and death
forecasts.We observe that until the end of 2021, the performance was close to
the Hub ensembles. The case forecasts had poor performance just before and
during the rise of the Omicron variant. The errors are higher in the US
compared to those in the ECDC. This could be due to sparser state-level data
on variants leading to poor estimates of variant proportions. Death forecasts
are impacted less in the US, with the exception of one outlier. On the other
hand, ECDC death forecasts are impacted more during this period. This may be
due to the poorer quality of case reporting during this period in Europe
compared to the US states (see Figure 5). Since, during this time, our death
forecasts were based on a function of reported cases, death forecasts were of
relatively poor quality as well.

Our models often correctly predict the existence of an upcoming surge, but
they tend to over-predict its intensity (Figure 6). While qualitatively,
predicting the existence of a surge is useful, with traditional metrics such
as mean absolute error, the models are penalized for overpredicting. How to
improve the forecast based on these metrics during a time of variant-drive
surge is not clear. Better quality of data on genomic surveillance and case
reporting can help with better estimation of variant proportions, resulting in
better forecasts, at least in the short term.

##  6 Conclusions

We presented various versions of the SIkJalpha model that have been used over
time over different forecast hubs and scenario modeling hubs. The models have
evolved into one capable of accounting for various age groups, an arbitrary
number of variants, imperfect protection from natural and vaccine-induced
immunity, and immunity states formed by multiple vaccine doses and infections.
We have presented the implementation details including data sources, parameter
estimation methods, runtimes, and uncertainty quantification. The goal of this
paper was to provide an overview of the modeling choices over time. The
details of the impact of each modeling decision are beyond the scope of this
paper. Instead, we focused on the breadth of modeling decisions. Our forecasts
for the forecast hubs were generated through automated scripts. The advantage
of our approach is scalability due to automation and fast implementation. That
is, we do not perform any manual tuning or analyze individual forecasts for
each region. The same advantage has the drawback that the models may produce
poor forecasts if there is an anomaly or high noise in data that have not been
already accounted for in the scripts. A particular period where the models
perform poorly is just before a surge due to a new variant. The model is able
to predict the existence of a surge, however, it often over-predicts the
intensity.

In future work, we will perform retrospective analysis to identify which
modeling choices lead to better performance. Particularly, how to predict the
existence of a surge without over-predicting its intensity. This will include
identification of the transmissibility advantage, the immune escape advantage,
and better estimates of population immunity. Estimating population immunity
will require exact identification of the dynamics of waning immunity, which so
far has been borrowed from the suggestions of the Scenario Modeling Hub.

## Acknowledgments

This work was supported by the Centers for Disease Control and Prevention and
the National Science Foundation under the awards no. 2027007, 2135784, and
2223933. Any opinions, findings, and conclusions or recommendations expressed
in this material are those of the author and do not necessarily reflect the
views of the National Science Foundation and the Center for Disease Control
and Prevention. The author would like to thank Tianjuan Xu and Majd Al Aawar
for assistance with the submission of forecasts to the forecast hubs. The
author would also like to thank the US Scenario Modeling Hub for useful
discussions that have helped with the evolution of the models.

Table 3: Modeling choices for various Hub submissions.

|  | US FH
---
(start: 06/20)
| G+P FH
---
(start: 11/20)
| Euro FH
---
(start: 02/21)
| US SMH
---
(start: 01/21)
| Euro SMH
---
(start: 03/22)
Cases | Yes | Yes | Yes | Yes | Yes
Deaths | Yes | Yes | Yes | Yes | Yes
Hospitalizations | 11/20-2/23 | - | - | Yes | -
True Infections |  | Estimated: start-12/20
---
Seroprevalence: 1/21-3/22
Wastewater: 3/22-2/23
Extrapolated | Extrapolated |  | Seroprevalence: start-3/22
---
Wastewater: 3/22-2/23
Extrapolated
Vaccine Adoption |  | Linear: 1/21 – 7/21
---
Contagious: 7/21 – 2/23
- |  | Linear: start – 7/21
---
Contagious: 7/21 – 2/23
| Sigmoid: start – 7/21
---
Contagious: 7/21 – 2/23
| Contagious:
---
start – 2/23
Vaccine Dosage |  2-dose: 1/21 – 7/21 All-doses: 07/21 – 2/23 | - |  2-dose: start – 7/21 2-dose+ boost: 07/21 – 2/23 | 2-dose: start – 7/21 |  2-dose+booster: start - 2/23
2-dose+boost: 07/21 – 2/23
Vaccine Model |  | AoN: 1/21 – 7/21
---
Waning: 7/21 – 2/23
- |  | AoN: start – 7/21
---
Waning: 7/21 – 2/23
| AoN: start – 7/21
---
Waning: 7/21 – 2/23
-
| Natural Immunity
---
Waning
7/21 - 2/23 | - | 7/21 - 2/23 | 7/21 – 2/23 | start - 2/23
Variants |  | 2-strain: 03/21 - 07/21
---
All PANGO: 07/21 - 12/21
Selected: 12/21 - 2/23
- |  | All PANGO: 07/21 - 12/21
---
Selected: 12/21 - 2/23
| 2-strain: 03/21 - 07/21
---
All PANGO: 07/21 - 12/21
Selected: 12/21 - 2/23
start - 2/23
Immune Escape | Since 12/21 | - | Since 12/21 | Since 12/21 | Yes
Age groups | - | - | - | Yes | -
| All Susceptibility
---
States
Since 12/21 | - | Since 12/21 | Since 12/21 | Yes
Quantiles |  | Sub-scenarios:
---
11/20-2/23
Random Forest |  | Sub-scenarios:
---
11/20-2/23
| Sub-scenarios:
---
1/21-2/23
| Sub-scenarios:
---
start-2/23

## References

  * [1] Biobot analytics: Population health analytics powered by sewage.  https://github.com/biobotanalytics/covid19-wastewater-data.
  * [2] CoVariants.  https://covariants.org/.
  * [3] COVID-19 Open-Data.  https://github.com/GoogleCloudPlatform/covid-19-open-data.
  * [4] COVID-19 Scenario Modeling Hub.  https://github.com/midas-network/covid19-scenario-modeling-hub.
  * [5] Cuebiq: Contact Index Dashboard.  https://www.cuebiq.com/visitation-insights-contact-index/.
  * [6] DELPHI Group: Forecast Evaluation Dashboard.  https://delphi.cmu.edu/forecast-eval/.
  * [7] European COVID-19 Forecast Hub Evaluation Report.  https://covid19forecasthub.eu/reports.html.
  * [8] European COVID-19 Scenario Hub.  https://github.com/covid19-forecast-hub-europe/covid19-scenario-hub-europe.
  * [9] Extension of SARS-CoV-2 data processing to incorporate Open Data through GenBank.  https://nextstrain.org/blog/2021-07-08-ncov-open-announcement.
  * [10] Forecast evaluations - german and polish covid-19 forecast hub.  https://jobrac.shinyapps.io/app_evaluation/.
  * [11] GISAID: Tracking of Variants.  https://www.gisaid.org/hcov19-variants/.
  * [12] Nationwide Commercial Laboratory Seroprevalence Survey.  https://covid.cdc.gov/covid-data-tracker/#national-lab.
  * [13] New York State Statewide COVID-19 Reinfection Data.  https://health.data.ny.gov/Health/New-York-State-Statewide-COVID-19-Reinfection-Data/7aaj-cdtu.
  * [14] Pango lineages: Latest epidemiological lineages of sars-cov-2.  https://cov-lineages.org/.
  * [15] US COVID-19 Forecast Hub: Forecast Evaluations.  https://covid19forecasthub.org/eval-reports/.
  * [16] The COVID-19 Forecast Hub.  https://covid19forecasthub.org/, 2020.
  * [17] European COVID-19 Forecast Hub.  https://covid19forecasthub.eu/index.html, 2020.
  * [18] Aniruddha Adiga, Lijing Wang, Benjamin Hurt, Akhil Peddireddy, Przemyslaw Porebski, Srinivasan Venkatramanan, Bryan Leroy Lewis, and Madhav Marathe.  All models are useful: Bayesian ensembling for robust high resolution covid-19 forecasting.  In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pages 2505–2513, 2021.
  * [19] Linda JS Allen.  Some discrete-time si, sir, and sis epidemic models.  Mathematical biosciences, 124(1):83–105, 1994.
  * [20] Johannes Bracher, Jannik Deuschel, Tilmann Gneiting, Konstantin Görgen, and Melanie Schienle.  Assembling forecasts of COVID19 cases and deaths in Germany and Poland.  https://jobrac.shinyapps.io/app_forecasts_de/.
  * [21] Johannes Bracher, Evan L Ray, Tilmann Gneiting, and Nicholas G Reich.  Evaluating epidemic forecasts in an interval format.  PLoS computational biology, 17(2):e1008618, 2021.
  * [22] Johannes Bracher, Daniel Wolffram, Jannik Deuschel, Konstantin Görgen, Jakob L Ketterer, Alexander Ullrich, Sam Abbott, Maria Vittoria Barbarossa, Dimitris Bertsimas, Sangeeta Bhatia, et al.  A pre-registered short-term forecasting study of covid-19 in germany and poland during the second wave.  Nature communications, 12(1):1–16, 2021.
  * [23] Lauren Castro, Geoffrey Fairchild, Isaac Michaud, and Dave Osthus.  Coffee: Covid-19 forecasts using fast evaluations and estimation.  arXiv preprint arXiv:2110.01546, 2021.
  * [24] Jiangzhuo Chen, Anil Vullikanti, Joost Santos, Srinivasan Venkatramanan, Stefan Hoops, Henning Mortveit, Bryan Lewis, Wen You, Stephen Eubank, Madhav Marathe, et al.  Epidemiological and economic impact of covid-19 in the us.  Scientific reports, 11(1):20451, 2021.
  * [25] Shi Chen, Rajib Paul, and Jean-Claude Janies, Danieland Thill.  UNCC-hierbin (Github).  https://github.com/midas-network/covid19-scenario-modeling-hub/blob/master/data-processed/UNCC-hierbin/metadata-UNCC-hierbin.txt.
  * [26] Estee Y Cramer, Evan L Ray, Velma K Lopez, Johannes Bracher, Andrea Brennen, Alvaro J Castro Rivadeneira, Aaron Gerding, Tilmann Gneiting, Katie H House, Yuxin Huang, et al.  Evaluation of individual and ensemble probabilistic forecasts of covid-19 mortality in the united states.  Proceedings of the National Academy of Sciences, 119(15):e2113561119, 2022.
  * [27] Jessica T Davis, Matteo Chinazzi, Nicola Perra, Kunpeng Mu, Ana Pastore y Piontti, Marco Ajelli, Natalie E Dean, Corrado Gioannini, Maria Litvinova, Stefano Merler, et al.  Cryptic transmission of sars-cov-2 and the first covid-19 wave.  Nature, 600(7887):127–132, 2021.
  * [28] Klaus Dietz.  The estimation of the basic reproduction number for infectious diseases.  Statistical methods in medical research, 2(1):23–41, 1993.
  * [29] Joseph Friedman, Patrick Liu, Christopher E Troeger, Austin Carter, Robert C Reiner Jr, Ryan M Barber, James Collins, Stephen S Lim, David M Pigott, Theo Vos, et al.  Predictive performance of international covid-19 mortality forecasting models.  Nature communications, 12(1):2609, 2021.
  * [30] Emily Howerton, Lucie Contamin, Luke C Mullany, Michelle M Qin, Nicholas G Reich, Samantha J Bents, Rebecca K Borchering, Sung-mok Jung, Sara L Loo, Claire P Smith, et al.  Informing pandemic response in the face of uncertainty. an evaluation of the us covid-19 scenario modeling hub.  medRxiv, pages 2023–06, 2023.
  * [31] Joseph C Lemaitre, Kyra H Grantz, Joshua Kaminsky, Hannah R Meredith, Shaun A Truelove, Stephen A Lauer, Lindsay T Keegan, Sam Shah, Josh Wills, Kathryn Kaminsky, et al.  A scenario modeling pipeline for covid-19 emergency planning.  Scientific reports, 11(1):7534, 2021.
  * [32] Benjamin P Linas, Jade Xiao, Ozden O Dalgic, Peter P Mueller, Madeline Adee, Alec Aaron, Turgay Ayer, and Jagpreet Chhatwal.  Projecting covid-19 mortality as states relax nonpharmacologic interventions.  In JAMA Health Forum, volume 3, pages e220760–e220760. American Medical Association, 2022.
  * [33] Velma Lopez, Estee Y Cramer, Robert Pagano, John M Drake, Eamon B O’Dea, Benjamin P Linas, Turgay Ayer, Jade Xiao, Madeline Adee, Jagpreet Chhatwal, et al.  Challenges of covid-19 case forecasting in the us, 2020-2021.  medRxiv, pages 2023–05, 2023.
  * [34] Kevin C Ma, Tigist F Menkir, Stephen Kissler, Yonatan H Grad, and Marc Lipsitch.  Modeling the impact of racial and ethnic disparities on covid-19 epidemic dynamics.  Elife, 10:e66601, 2021.
  * [35] Nicolai Meinshausen.  Quantile regression forests.  Journal of Machine Learning Research, 7(Jun):983–999, 2006.
  * [36] Sara E Oliver, Julia W Gargano, Mona Marin, Megan Wallace, Kathryn G Curran, Mary Chamberland, Nancy McClung, Doug Campos-Outcalt, Rebecca L Morgan, Sarah Mbaeyi, et al.  The advisory committee on immunization practices’ interim recommendation for use of pfizer-biontech covid-19 vaccine—united states, december 2020.  Morbidity and mortality weekly report, 69(50):1922, 2020.
  * [37] Christian H Reinsch.  Smoothing by spline functions.  Numerische mathematik, 10(3):177–183, 1967.
  * [38] Alexander Rodriguez, Anika Tabassum, Jiaming Cui, Jiajia Xie, Javen Ho, Pulak Agarwal, Bijaya Adhikari, and B Aditya Prakash.  Deepcovid: An operational deep learning-driven framework for explainable real-time covid-19 forecasting.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 15393–15400, 2021.
  * [39] Ajitesh Srivastava.  Forecasts for “Almost” Everywhere.  https://scc-usc.github.io/ReCOVER-COVID-19/#/row.
  * [40] Ajitesh Srivastava.  ReCOVER: Accurate Predictions and Scenario Projections for COVID-19 Epidemic Response.  https://github.com/scc-usc/ReCOVER-COVID-19.
  * [41] Ajitesh Srivastava and Viktor K Prasanna.  Data-driven Identification of Number of Unreported Cases for COVID-19: Bounds and Limitations.  ACM SIGKDD international conference on Knowledge discovery and data mining (Health Day), arXiv preprint arXiv:2006.02127, 2020.
  * [42] Ajitesh Srivastava and Viktor K Prasanna.  Learning to Forecast and Forecasting to Learn from the COVID-19 Pandemic.  arXiv preprint arXiv:2004.11372, 2020.
  * [43] Ajitesh Srivastava and Tianjian Xu.  COVID19 Forecast Bench.  https://github.com/scc-usc/covid19-forecast-bench, 2020.
  * [44] Ajitesh Srivastava, Tianjian Xu, and Viktor K Prasanna.  Fast and Accurate Forecasting of COVID-19 Deaths Using the SIkJalpha Model.  arXiv preprint arXiv:2007.05180, 2020.
  * [45] Ajitesh Srivastava, Tianjian Xu, and Viktor K Prasanna.  The epibench platform to propel ai/ml-based epidemic forecasting: A prototype demonstration reaching human expert-level performance.  5th International Workshop on Health Intelligence, AAAI, 2021.
  * [46] Brad Suchoski, Steve Stage, Heidi Gurung, and Prasith Baccam.  Gpu accelerated parallel processing for large-scale monte carlo analysis: Covid-19 parameter estimation and new case forecasting.  Frontiers in Applied Mathematics and Statistics, 8:818016, 2022\.
  * [47] Yueying Wang, Myungjin Kim, Shan Yu, Xinyi Li, Guannan Wang, and Li Wang.  Nonparametric estimation and inference for spatiotemporal epidemic models.  Journal of Nonparametric Statistics, 34(3):683–705, 2022.
  * [48] Shun Zheng, Zhifeng Gao, Wei Cao, Jiang Bian, and Tie-Yan Liu.  Hierst: A unified hierarchical spatial-temporal framework for covid-19 trend forecasting.  In Proceedings of the 30th ACM International Conference on Information & Knowledge Management, pages 4383–4392, 2021.