# Regularised B-splines projected Gaussian Process priors to estimate time-
trends of age-specific COVID-19 deaths related to vaccine roll-out

Mélodie<EMAIL_ADDRESS>[ Alexandra Blenkinsop
Andrea Brizzi  Yu Chen  Carlos Cardoso Correia Perello  Vidoushee Jogarah
Yuanrong Wang  Seth Flaxman  Samir Bhatt  Oliver
<EMAIL_ADDRESS>[ Corresponding authors e2

(2021)

# Supplementary Materials

Mélodie<EMAIL_ADDRESS>[ Alexandra Blenkinsop
Andrea Brizzi  Yu Chen  Carlos Cardoso Correia Perello  Vidoushee Jogarah
Yuanrong Wang  Seth Flaxman  Samir Bhatt  Oliver
<EMAIL_ADDRESS>[ Corresponding authors e2

(2021)

###### Abstract

The COVID-19 pandemic has caused severe public health consequences in the
United States. In this study, we use a hierarchical Bayesian model to estimate
the age-specific COVID-19 attributable deaths over time in the United States.
The model is specified by a novel non-parametric spatial approach, a low-rank
Gaussian Process (GP) projected by regularised B-splines. We show that this
projection defines a new GP with attractive smoothness and computational
efficiency properties, derive its kernel function, and discuss the penalty
terms induced by the projected GP. Simulation analyses and benchmark results
show that the spatial approach performs better than standard B-splines and
Bayesian P-splines and equivalently well as a standard GP, for considerably
lower runtimes. The B-splines projected GP priors that we develop are likely
an appealing addition to the arsenal of Bayesian regularising priors. We apply
the model to weekly, age-stratified COVID-19 attributable deaths reported by
the US Centers for Disease Control, which are subject to censoring and
reporting biases. Using the B-splines projected GP, we can estimate
longitudinal trends in COVID-19 associated deaths across the US by 1-year age
bands. These estimates are instrumental to calculate age-specific mortality
rates, describe variation in age-specific deaths across the US, and for
fitting epidemic models. Here, we couple the model with age-specific
vaccination rates to show that lower vaccination rates in younger adults aged
18-64 are associated with significantly stronger resurgences in COVID-19
deaths, especially in Florida and Texas. These results underscore the critical
importance of medically able individuals of all ages to be vaccinated against
COVID-19 in order to limit fatal outcomes.

62F15,

60G15,

65D07,

62P10,

Bayesian Inference,

Gaussian Processes,

Splines,

Non-parametric estimation,

COVID-19,

###### keywords:

[class=MSC]

###### keywords:

††volume: TBA††issue: TBA

, , , , , , , , and

t1First supporter of the project t2Second supporter of the project
a1Department of Mathematics, Imperial College London a2MRC Centre for Global
Infectious Disease Analysis, Jameel Institute, School of Public Health,
Imperial College London a3Section of Epidemiology, Department of Public
Health, University of Copenhagen

##  1 Introduction

A new pathogen, Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2),
emerged in the Wuhan region of China in December 2019 and continues to spread
worldwide. The resulting disease, COVID-19, is severe, with overall infection
fatality rates (IFRs) between 0.1% and 1% (Meyerowitz-Katz and Merone, 2020;
Brazeau et al., 2020), which increase exponentially with age (Levin et al.,
2020). Developed vaccines are highly effective to prevent deaths (Haas et al.,
2021; Baden et al., 2021), and initially have been prioritised to older age
groups. In the United States (US), vaccines are now offered to all adults
since May 2021, but uptake has been variable across states (Centers for
Disease Control and Prevention, 2021a). Despite increasing vaccine coverage,
several US states, most notably Florida and Texas, are reporting resurgent
COVID-19 death waves that are on par or exceed in magnitude the death waves
seen in the prevaccination period. A key question is if these death waves
since summer 2021 are linked to low vaccination coverage in younger adults,
which are known to drive transmission (Monod et al., 2021).

Here, we provide methods for estimating the age composition of COVID-19
attributable deaths over time, and to characterise the potential shifts in the
age composition of deaths. Doing so requires a statistical model because data
are partially censored, reported with delays, and reported in age bands that
can be inconsistent across locations or do not match those of other data
streams such as reported vaccinations. Numerous methods have been developed to
interpret the time evolution of overall COVID-19 attributable deaths (e.g.,
Lavezzo et al. (2020); Institute for Health Metrics and Evaluation (2020);
Chen et al. (2020); Blangiardo et al. (2020); Zheng et al. (2021); Flaxman et
al. (2020)), but few are suitable to analyze longitudinal data stratified by
age brackets or other discrete strata (Monod et al., 2021), let alone at high
resolution such as 1-year age bands. One reason for this paucity of methods is
that in SIR-type models, calculating age-specific next generation operators at
a time resolution of days becomes computationally slow over observation
periods that span years, especially when calculations need to be repeated
millions of times within a Bayesian framework (Wikle et al., 2020; Monod et
al., 2021). These considerations are prompting us and others to consider non-
mechanistic, flexible estimation approaches (Shah et al., 2020; Pokharel and
Deardon, 2021). As our objective is to reconstruct spatiotemporal trends in
age-specific COVID-19 attributable deaths, it does not require invoking many
of the assumptions or complexities that underlie mechanistic or semi-
mechanistic COVID-19 transmission dynamics models and Bayesian non-parametric
models are sufficient. We present a fully Bayesian, computationally scalable
approach to estimate a two-dimensional (2D) surface over ages and weeks that
describes the time evolution of COVID-19 deaths by $1$-year age bands at US
state level. To impute missing entries on the surface and estimate global
trends over ages and weeks, we borrow information across neighboring entries
by using a non-parametric smoothing method.

A natural starting point for modelling a surface is a 2D Gaussian process (GP)
(Rasmussen and Williams, 2005). However, their computational complexity makes
the use of 2D GPs in a fully Bayesian framework difficult when the surface
dimension becomes large, even after using a Kronecker decomposition of the
kernel function (Saatçi, 2011; Wilson et al., 2014). Here we adopt a low-rank
approximation via a tensor product of B-splines for which the parameters
follow a 2D GP. The resulting approach is equivalent to a 2D GP defined by a
low-rank covariance matrix projected by B-splines. B-splines are a popular
choice for non-parametric modelling, due to their continuity properties,
ensuring smoothness of the fitted surface, and their easy implementation. But
choosing the optimal number and position of knots—the defining grid segments
where the surface is expected to change its behavior—on the space of ages and
weeks is a complex task. Some approaches have focused on adding a penalty to
restrict the flexibility of the fitted surface in a frequentist framework
(O’Sullivan, 1986, 1988; Eilers and Marx, 1996; Eilers et al., 2006).
Following this idea, we regularise the fitted surface by using a kernel
function with a free complexity parameter to define the covariance matrix of
the low-rank 2D GP. We qualitatively compare the penalty induced by this
choice to that of related regularisation methods. We benchmark the proposed
regularised B-splines projected GP against several other popular smoothing
methods, and we demonstrate that our approach results in substantive
computational gains over a standard 2D GP for similar estimation accuracy.

This paper focus on estimating age-specific COVID-19 attributable deaths in
the context of vaccine roll-out on publicly available, age-specific COVID-19
death and vaccination data from the $4$ most populated US states, California,
Florida, New York (state) and Texas, that are reported by the Centers for
Disease Control and Prevention (CDC) (2021a; 2021c). It is structured as
follows. Section 2 introduces the data and their limitations. Section 3
describes the proposed methodology to model trends in the age-composition of
COVID-19 deaths, including a theoretical characterization of the penalties
introduced by B-splines projected GPs. Section 4 presents a comparison of the
proposed method to related smoothing approaches on simulated data, and Section
5 on real world data used for benchmarking. In Section 6 we estimate the time
evolution of age-specific COVID-19 deaths for the four most populated US
states. We document the marked variation in the summer 2021 resurgence in age-
specific COVID-19 deaths across states, and show strong resurgences in deaths
are associated low vaccination coverage in younger adults. Lastly, Section 7
closes with a discussion.

The proposed modelling framework can be further applied to estimate COVID-19
mortality rates in arbitrary age groups, investigate differences in the age
composition of deaths across locations or identify time shifts in the age
composition of deaths. The code to use our approach and to reproduce the
results is available at
https://github.com/ImperialCollegeLondon/B-SplinesProjectedGPs. Median and 95%
credible interval of the state- and age-specific COVID-19 attributable deaths
predicted by our approach are available in the GitHub repository. Posterior
samples are also available upon request. Our results suggest regularised
B-splines projected GP priors are likely useful for other surface estimation
problems, and we provide templated Stan model files in Supplement S12.

##  2 COVID-19 deaths and vaccination data

###  2.1 COVID-19 deaths data

The Center for Disease Control (CDC) and National Center for Health Statistics
(NCHS) report each week on the total number of deaths involving COVID-19 that
have been received and tabulated through the National Vital Statistics
database (Centers for Disease Control and Prevention, 2021c) for every US
state across the age groups

$\displaystyle\mathcal{B}$
$\displaystyle=\big{\\{}0,1-4,5-14,15-24,25-34,35-44,45-54,$ (2.1)
$\displaystyle\quad\quad 55-64,65-74,75-84,85+\big{\\}}.$

We refer to the latter as cumulative reported COVID-19 attributable deaths,
$D_{m,b,w}$, for state $m$, in age group $b$ and on week $w$. Historical
records are made available by Rearc (2021). To simplify the longitudinal
dependence and create a time series from the cumulative counts (King et al.,
2015), we obtain the weekly COVID-19 attributable deaths by differencing,

$d_{m,b,w}=D_{m,b,w+1}-D_{m,b,w},$ (2.2)

in each location and age group for all but the last available week (daily
deaths, $d$, for state $m$, in age group $b$, on week $w$). We index the
weekly deaths by $w=1,\dotsc,W$.

Reflecting the reporting nature of the age-stratified CDC data, the weekly
COVID-19 attributable $d_{m,b,w}$ are subject to reporting delays and do not
necessarily correspond to the number of individuals who died of COVID-19 in
state $m$ and week $w$. The CDC does not report the cumulative deaths count if
the count is between $1$ and $9$. Thus, it is not possible to retrieve all
weekly deaths, and we refer to the set of weeks that are retrievable in state
$m$ for age group $b$ through first order differencing by
$\mathcal{W}_{m,b}^{\text{WR}}$, and to the set of weeks that were non-
retrievable because of censoring by $\mathcal{W}_{m,b}^{\text{WNR}}$. The
censored cumulative deaths are bounded, such that the sum of non-retrievable
weekly deaths is also bounded. The exact computations are described in
Appendix S10.1. In addition, there was no update on July 4, 2020 resulting in
missing cumulative deaths in that week. The weekly deaths of that week and the
preceding week are declared missing. Note that the missing weekly deaths are
not equivalent to the non-retrievable weekly deaths because we cannot bound
their sums. In the main text, we focus on data reported from May 2, 2020 to
September 25, 2021. We thus have $W=74$. Figure 1 show the COVID-19
attributable weekly deaths data for the four most populated US states,
California, Florida, New York and Texas.

###  2.2 COVID-19 vaccination data

The CDC reports weekly time series of the proportions of individuals aged
$18$-$64$ and $65+$ who are fully vaccinated (Centers for Disease Control and
Prevention, 2021b), referred to as vaccination rate in the following. Fully
vaccinated individuals are defined as having received the second dose of a
two-dose vaccine or one dose of a single-dose vaccine. We use records from
March 02, 2021 to September 25, 2021, and denote by $\text{v}_{m,c,w}$ the
vaccination rate for state $m$ in age group $c\in\mathcal{C}$ at the start of
week $w$,

$\displaystyle\mathcal{C}=\\{18-64,65+\\}.$ (2.3)

Figure 2 presents the vaccination rates among individuals aged $18$-$64$ and
$65+$ reported in the four most populated US states, California, Florida, New
York and Texas.

Figure 1: Retrievable weekly COVID-19 attributable deaths in four US states.
The retrievable weekly deaths are computed using the first order difference of
reported cumulative deaths by the CDC (Centers for Disease Control and
Prevention, 2021c). Grey cells indicate non-retrievable weekly death counts
due to the censoring of cumulative deaths. Orange bars show the missing weekly
deaths due to non-reported cumulative deaths. Figure 2: Age-specific
proportion of fully vaccinated individuals in four US states. The proportion
of fully vaccinated individuals over time were report by the CDC (Centers for
Disease Control and Prevention, 2021b). Fully vaccinated individuals are
defined as having receive the second dose of a two-dose vaccine or one dose of
a single-dose vaccine. Column facets show the proportion of fully vaccinated
individuals aged $18$-$64$ and $65+$.

###  2.3 Other data sets

We use two additional data sets. Firstly, we retrieved the weekly COVID-19
attributable deaths regardless of age from John Hopkins University (JHU) for
all US states from May 2, 2020 to September 25, 2021 (John Hopkins University
of Medicine, 2020), which we denote by $d^{\text{JHU}}_{m,w}$ for state $m$
and week $w$. In these data, reporting-delayed deaths are back-distributed
where possible, and for this reason we use the overall deaths as a reference
to mitigate the reporting delays in the age-stratified CDC data (2.2).
Secondly, age-specific weekly COVID-19 attributable deaths are also directly
reported by state Departments of Health (DoH) on their websites, or data
repositories. For California, Florida and Texas, data records are available up
to April 1, 2021 at https://github.com/ImperialCollegeLondon/US-
covid19-agespecific-mortality-data. The DoH do not censor low death counts,
and the all-ages sum of the DoH data correspond well to JHU data. This prompts
us to use the DoH data as an independent data set to assess the accuracy of
our model estimates that are derived from the CDC data.

##  3 Methods

###  3.1 Modelling weekly COVID-19 attributable deaths by age

For simplicity we suppress the state index in what follows, with all equations
being analogous. Our aim is to estimate the weekly deaths by $1$-year age
bands, $a\in\mathcal{A}=\\{0,1,\dots,104,105\\}$ in week
$w\in\mathcal{W}=\\{1,\dotsc,W\\}$, and we denote their expectation by
$\mu_{a,w}$. We first decompose $\mu_{a,w}$ as the product of the weekly
deaths for all ages $\lambda_{w}$ with the relative contribution $\pi_{a,w}$
of age $a$ to weekly deaths, where $\sum_{a}\pi_{a,w}=1$ for all $w$. Second,
we model the longitudinal age composition of deaths $\pi_{a,w}$ as a bivariate
random function $(a,w)\to f(a,w)\in\mathbb{R}$, that is exponentiated,
normalised, and evaluated on the 2D grid $\mathcal{A}\times\mathcal{W}$. Our
basic latent model structure is thus

$\displaystyle\begin{split}\mu_{a,w}&=\lambda_{w}\pi_{a,w}\\\
\pi_{a,w}&=\text{softmax}\Big{(}[f(a,w)]_{a\in\mathcal{A}}\Big{)}\\\
&=\bigg{(}\frac{\exp f(a,w)}{\sum_{\tilde{a}\in\mathcal{A}}\exp
f(\tilde{a},w)}\bigg{)}.\end{split}$ (3.1)

To link the expected weekly deaths by $1$-year age band, $\mu_{a,w}$, to the
data (2.2), we aggregate them over the age groups specified by the CDC,
$\mu_{b,w}=\sum_{a\in b}\mu_{a,w}$, for all $b\in\mathcal{B}$. Then, we model
the observed, weekly deaths in age $b$ and week $w$ for all ages and weeks
through Negative Binomial distributions in the shape-scale parameterisation,

$\displaystyle d_{b,w}\>|\>\alpha_{b,w},\theta$
$\displaystyle\sim\text{NegBin}(\alpha_{b,w},\theta)$ (3.2a)
$\displaystyle\mu_{b,w}$ $\displaystyle=\alpha_{b,w}\frac{\theta}{1-\theta}$
(3.2b) $\displaystyle\sigma_{b,w}^{2}$
$\displaystyle=\alpha_{b,w}\frac{\theta}{(1-\theta)^{2}}=\mu_{b,w}(1-\theta)^{-1}=\mu_{b,w}\>(1+\nu),$
(3.2c)

with mean $\mu_{b,w}$ and variance $\sigma_{b,w}^{2}$, and where $\nu>0$ is
interpreted as an overdispersion parameter. The purpose of the shape-scale
parameterisation with identical $\theta$ (3.2a) is that then, the weekly
deaths conditional on their total follow a Dirichlet-Multinomial distribution
with parameters $\alpha_{b,w}$, resulting in the succinct identity
$\alpha_{b,w}=\sum_{a\in b}\alpha_{a,w}$ (Townes, 2020). For implementation,
notice that the shape and scale parameters can be rewritten as
$\alpha_{b,w}=\mu_{b,w}/\nu$ and $\theta=\nu/(1+\nu)$.

Our Bayesian model comprises as log likelihood the sum of the log Negative
Binomial densities (3.2a) over retrievable weekly deaths parameterised by
$\mu_{b,w}$ and $\nu$. Considering data censoring, similar terms involving log
cumulative density functions of the Negative Binomial that bound sums of the
non-retrievable weekly deaths are added in the log likelihood. Taken together,
for the collection of weekly deaths
$\mathbf{d}=\big{(}d_{b,w},b\in\mathcal{B},w\in\\{\mathcal{W}_{b}^{\text{WR}},\mathcal{W}_{b}^{\text{WNR}}\\}\big{)}$,
our log posterior density is

$\displaystyle\log
p(\boldsymbol{\lambda},\boldsymbol{f},\boldsymbol{\phi},\nu\>|\>\mathbf{d})$
$\displaystyle\propto\sum_{b\in\mathcal{B}}\sum_{w\in\mathcal{W}_{b}^{\text{WR}}}\log\text{NegBin}(d_{b,w}\>|\>\mu_{b,w}/\nu,\frac{\nu}{1+\nu})$
(3.3a) $\displaystyle\quad\quad\quad\quad+\>\log
p\big{(}\\{d_{b,w}\\}_{w\in\mathcal{W}_{b}^{\text{WNR}}}\>|\>\dotsc\big{)}$
(3.3b)
$\displaystyle\quad+\>\sum_{w\in\mathcal{W}}\log\text{Gamma}(\lambda_{w}\>|\>T_{w},T_{w}/(2\eta))$
(3.3c) $\displaystyle\quad+\>\log\mathcal{N}^{+}(\nu^{-1/2}\>|\>0,1)\>+\>\log
p(\boldsymbol{f}\>|\>\boldsymbol{\phi})\>+\>\log p(\boldsymbol{\phi}),$ (3.3d)

where for brevity the term
$p\big{(}\\{d_{b,w}\\}_{w\in\mathcal{W}_{b}^{\text{WNR}}}|\dotsc\big{)}$ in
(3.3b) is detailed in Appendix S10.2. The priors are specified as follows.
First, the prior on the total number of weekly deaths $\lambda_{w}$, in
(3.3c), is specified in the mean-standard deviation parametrisation. The prior
expectation on the total deaths is found by summing the retrievable deaths,
such that
$T_{w}=\sum_{b\in\mathcal{B},w\in\mathcal{W}_{b}^{\text{WR}}}d_{b,w}$. Next,
we assume that the standard deviation is equal to twice the first order
difference in the total deaths. We find the empirical ratio of the total death
relative to their first order difference $\eta$, using ordinary least square
method without intercept, and specify the standard deviation accordingly.
Second, the inverse of the square root of the overdispersion parameter is
given a standard normal distribution truncated on the positive support, the
recommended prior in this context. Lastly, the model is completed with a prior
on the random function $\boldsymbol{f}$ and its hyperparameters
$\boldsymbol{\phi}$, which we investigate in detail in the next section.

Considering reporting delays in the data, we explicitly allow for a rescaling
of the sum of deaths across age groups in the model according to other curated
data sets, which effectively re-distributes reporting delayed deaths in the
CDC data set to earlier dates. Specifically, we require that our posterior
predictions of the COVID-19 attributable deaths in age group $a$ and week $w$,
$d^{\star}_{a,w}$, sum to $\sum_{a}d^{\star}_{a,w}=d_{w}^{\text{JHU}}$. We
achieve this by exploiting the probabilistic relationship between Negative
Binomial distributions in the shape-scale parameterisation and the Dirichlet-
Multinomial, through

$\begin{split}&p(d^{\star}_{\cdot,w}\>|\>\text{{d}},d_{w}^{\text{JHU}})=\\\
&\quad\quad\int\text{DirMult}\big{(}d^{\star}_{\cdot,w}\>|\>d^{\text{JHU}}_{w},(\alpha_{a,w})_{a\in\mathcal{A}}\big{)}\>p(\boldsymbol{\lambda},\nu,\boldsymbol{f},\boldsymbol{\phi}\>|\>\text{{d}})\>d(\boldsymbol{\lambda},\nu,\boldsymbol{f},\boldsymbol{\phi}),\end{split}$
(3.4)

where as before $\alpha_{a,w}=\mu_{a,w}/\nu$. We denote by $\mu^{\star}_{a,w}$
the expected predictive COVID-19 attributable deaths in age group $a$ at week
$w$,

$\displaystyle\mu^{\star}_{a,w}\>|\>\text{{d}},d_{w}^{\text{JHU}}=\mathbb{E}\big{[}d^{\star}_{a,w}\>|\>\text{{d}},d_{w}^{\text{JHU}}\big{]}.$
(3.5)

###  3.2 Modelling age-specific contributions to COVID-19 weekly deaths

In (3.1), we introduced a 2D function $f(a,w)$, for which we shall find a
prior in this section. Let the number of points of the age axis be
$n=|\mathcal{A}|$ and on the week axis be $m=|\mathcal{W}|$. The total number
of points on the grid is $N=n\times m$. The ensemble of pairs of points is
$\boldsymbol{X}=\big{(}\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{N}\big{)}=\mathcal{A}\times\mathcal{W}$.
We now investigate different modelling approaches for the function $f(a,w)$.

####  3.2.1 Two-dimensional Gaussian Process

Given observations
$(\boldsymbol{X},\boldsymbol{f})=\big{\\{}\big{(}\boldsymbol{x}_{1},f(\boldsymbol{x}_{1})\big{)},\dots,\big{(}\boldsymbol{x}_{N},f(\boldsymbol{x}_{N})\big{)}\big{\\}}$,
we start by considering as model of $\boldsymbol{f}$ a zero-mean 2D GP,

$\displaystyle\boldsymbol{f}\>|\>\boldsymbol{\phi}$
$\displaystyle\sim\mathcal{GP}\big{(}\boldsymbol{0},\boldsymbol{K}\big{)}.$
(3.6)

The covariance matrix $\boldsymbol{K}$ is evaluated at all pairs of points in
$\boldsymbol{X}$ and has entries
$\boldsymbol{K}_{\boldsymbol{x},\boldsymbol{x}^{\prime}}=\text{Cov}\big{(}f(\boldsymbol{x}),f(\boldsymbol{x}^{\prime})\big{)}=k(\boldsymbol{x},\boldsymbol{x}^{\prime})$
with $\boldsymbol{x},\boldsymbol{x}^{\prime}\in\boldsymbol{X}$, where $k(.,.)$
is a kernel function. For computational efficiency and because our output is
on a multidimensional grid, we decompose the kernel function,

$\displaystyle
k\Big{(}(a,w),(a^{\prime},w^{\prime})\Big{)}=k^{1}(a,a^{\prime})\>k^{2}(w,w^{\prime})$
(3.7)

where $k^{1}(.,.)$ and $k^{2}(.,.)$ are kernel functions over ages and weeks,
respectively. The corresponding covariance matrix $\boldsymbol{K}$ is
calculated with the Kronecker product
$\boldsymbol{K}=\boldsymbol{K}^{2}\otimes\boldsymbol{K}^{1}$, and the number
of operations to evaluate the covariance matrix reduces from
$\mathcal{O}(N^{3})$ to $\mathcal{O}(2N^{3/2})$ (Gönen and Alpaydin, 2011;
Saatçi, 2011). For implementation, we note that the order of the product
follows from the fact that the matrix’s entries are stacked columnwise. We
rely on the efficient Kronecker product implementation via matrix-vector
products proposed in the Supplementary Material, Section 4 of Wilson et al.
(2014). In our applications, we set $k^{1}$ and $k^{2}$ to be squared
exponential kernel functions with variance scale $\zeta^{2}$ and specific
lengthscales for the rows and the columns, $\gamma_{1}$ and $\gamma_{2}$. Our
priors on the length scales are independent $\gamma_{i}\sim\text{Inv-
Gamma}(5,5)$ for $i=1,2$, and $\zeta\sim\text{Cauchy}^{+}(0,1)$. For brevity,
we denote $\boldsymbol{\phi}=\\{\zeta,\gamma_{1},\gamma_{2}\\}$.

####  3.2.2 B-splines surface

B-splines basis functions - or, more simply, B-splines - are constructed from
polynomial pieces that are joined at certain values over the input space,
called knots, and defined by a polynomial degree, $d$, and a non-decreasing
sequence of knots, $(t_{1}\leq t_{2}\leq\dots\leq t_{K})$, where $K$ is the
number of knots. Given those, the total number of B-splines is $I=K+d-1$. We
show how a B-spline is constructed in Appendix S9.1. In this paper, we use
cubic B-splines, with $d=3$. Moreover, we use equally spaced knots on both
dimensions, such that the only tuning parameter is the number of knots. A
fundamental property of B-splines, that constitutes most of their
attractiveness, is that they are smooth. More rigorously, a cubic B-spline
defined on strictly increasing knots is piecewise infinitely differentiable
between the knots, and of continuity $C^{2}$ on the knots (Goldman (2002),
chap. 7). $\boldsymbol{f}$ can be modelled with a tensor product of B-splines
given by

$\displaystyle
f(a,w)=\sum_{i=1}^{I}\sum_{j=1}^{J}\beta_{i,j}\>B^{1}_{i}(a)\>B^{2}_{j}(w),$
(3.8)

where $B^{1}_{i}(.)$ is the $i$th B-spline defined on the knot vector
$(t^{1}_{1}\leq t^{1}_{2}\leq\dots\leq t^{1}_{K_{1}})$ over the space
$\mathcal{A}$ with $K_{1}$ being the number of knots,
$i\in\mathcal{I}=\\{1,\dots,I\\}$ and $I=K_{1}+2$. Similarly, $B^{2}_{j}(.)$
is the $j$th B-spline defined on the knot vector $(t^{2}_{1}\leq
t^{2}_{2}\leq\dots\leq t^{2}_{K_{2}})$ over the space $\mathcal{W}$ with
$K_{2}$ being the number of knots, $j\in\mathcal{J}=\\{1,\dots,J\\}$ and
$J=K_{2}+2$. The ensemble of B-splines form a matrix basis denoted by
$\boldsymbol{B}^{1}$ and $\boldsymbol{B}^{2}$ of size $I\times n$ and $J\times
m$, respectively. The pairs of B-splines indices is denoted by
$\boldsymbol{U}=\big{(}\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{M}\big{)}=\mathcal{I}\times\mathcal{J}$
with $M=I\times J$. $[\beta_{\boldsymbol{u}}]_{u\in\boldsymbol{U}}$ is a
rectangular set of coefficients. In our applications, we obtain standard
B-splines surfaces by placing independent standard normal priors on all
B-splines coefficients $\beta_{\boldsymbol{u}}$. To keep notation streamlined,
notice that there are no hyperparameters $\boldsymbol{\phi}$.

####  3.2.3 Regularised B-splines projected Gaussian Process

Given the B-splines indices and corresponding coefficients
$(\boldsymbol{U},\boldsymbol{\beta})=\big{\\{}(\boldsymbol{u}_{1},{\beta}_{\boldsymbol{u}_{1}}),\dots,(\boldsymbol{u}_{M},{\beta}_{\boldsymbol{u}_{M}})\big{\\}}$,
we place a 2D GP on the coefficients,

$\displaystyle\boldsymbol{\beta}\>|\>\boldsymbol{\phi}\sim\mathcal{GP}(0,\boldsymbol{K}_{\beta}).$
(3.9)

The covariance matrix has entries
$(\boldsymbol{K}_{\beta})_{\boldsymbol{u},\boldsymbol{u}^{\prime}}=\text{Cov}\big{(}\beta_{\boldsymbol{u}},\beta_{\boldsymbol{u}^{\prime}}\big{)}=k_{\beta}\big{(}\boldsymbol{u},\boldsymbol{u}^{\prime}\big{)}$
where $k_{\beta}(.,.)$ is a kernel function depending on unknown
hyperparameters $\boldsymbol{\phi}$. This approach is equivalent to directly
placing a GP on $\boldsymbol{f}$ with a covariance matrix projected by
B-splines. To show this, we rewrite (3.8) as a matrix calculation,
$\boldsymbol{f}=(\boldsymbol{B}^{1})^{T}\boldsymbol{\beta}\boldsymbol{B}^{2}$,
and find its vectorized form

$\displaystyle\text{vec}(\boldsymbol{f})$
$\displaystyle=\text{vec}\Big{(}(\boldsymbol{B}^{1})^{T}\boldsymbol{\beta}\boldsymbol{B}^{2}\Big{)}=(\boldsymbol{B}^{2}\otimes\boldsymbol{B}^{1})^{T}\text{vec}(\boldsymbol{\beta}).$
(3.10)

The linear operator produces functions from $\mathbb{R}^{IJ}$ to
$\mathbb{R}^{nm}$. GPs are closed under linear operations (Papoulis and Pillai
(2002), chap.10), and so $\boldsymbol{f}$ specified by (3.8-3.9) is the GP

$\displaystyle\boldsymbol{f}\>|\>\boldsymbol{\phi}\sim\mathcal{GP}\Big{(}0,\big{(}\boldsymbol{B}^{2}\otimes\boldsymbol{B}^{1}\big{)}^{T}\boldsymbol{K}_{\beta}\big{(}\boldsymbol{B}^{2}\otimes\boldsymbol{B}^{1}\big{)}\Big{)}.$
(3.11)

In our applications, we decompose the kernel function $k_{\beta}$ into two
kernel functions over the rows and columns of the B-spline coefficients matrix
as in (3.7). We use squared exponential kernel functions with variance scale
$\zeta^{2}$ and specific lengthscales for the rows and the columns,
$\gamma_{1}$ and $\gamma_{2}$. Our priors on the length scale priors are
independent $\gamma_{i}\sim\text{Inv-Gamma}(5,5)$ for $i=1,2$, and
$\zeta\sim\text{Cauchy}^{+}(0,1)$. Under this prior,
$\boldsymbol{\phi}=\\{\zeta,\gamma_{1},\gamma_{2}\\}$.

The kernel function obtained by projecting a base kernel function with cubic
B-splines is $C^{2}$ as it can be represented as a linear combination of
$C^{2}$ functions (proof in Appendix S9.2). Therefore the surface obtained
when modelling a surface with the B-spline tensor product (3.11), unlike a
standard 2D GP (3.6), has inherent smoothness properties not matter the base
kernel functions $k_{\beta}$ used. We can also interpret (3.11) as a
regularised splines method. Indeed, existing regularised spline methods such
as smoothing splines (O’Sullivan, 1986, 1988) and P-splines (Eilers and Marx,
1996; Eilers et al., 2006) aim to minimise in a frequentist setting the loss
function

$\displaystyle\underbrace{\sum_{a\in\mathcal{A}}\sum_{w\in\mathcal{W}}\Big{(}f(a,w)-\sum_{i=1}^{I}\sum_{j=1}^{J}\hat{\beta}_{i,j}B^{1}_{i}(a)B^{2}_{j}(w)\Big{)}^{2}}_{\text{data
fit}}+\underbrace{P}_{\text{penalty}},$ (3.12)

where $\hat{\boldsymbol{\beta}}$ are the estimated coefficients and $P$ is a
penalty applied on the second derivative of the fitted curve for smoothing
splines and on finite differences of adjacent B-splines coefficients for
P-splines. Now adopting a Bayesian approach, our aim is to maximise the
probability of the posterior parameter conditional upon the data, which is
proportional to the likelihood multiplied by the prior. By placing a GP on the
B-splines coefficients, the prior of $\boldsymbol{f}$ modelled with the
B-spline tensor product (3.11) can be decomposed in the same way as (3.12),

$\displaystyle\begin{split}\text{log
p}(\text{vec}(\boldsymbol{f})\>|\>\boldsymbol{\phi})&\propto-\frac{1}{2}\Big{(}\text{vec}(\boldsymbol{f})^{T}\big{(}\boldsymbol{B}^{2}\otimes\boldsymbol{B}^{1}\big{)}^{T}\boldsymbol{K}_{\beta}\big{(}\boldsymbol{B}^{2}\otimes\boldsymbol{B}^{1}\big{)}\text{vec}(\boldsymbol{f})\\\
&\quad\quad\quad+\text{log}\>\Big{|}\big{(}\boldsymbol{B}^{2}\otimes\boldsymbol{B}^{1}\big{)}^{T}\boldsymbol{K}_{\beta}\big{(}\boldsymbol{B}^{2}\otimes\boldsymbol{B}^{1}\big{)}\Big{|}\Big{)}\\\
&\propto-\frac{1}{2}\Big{(}\underbrace{\text{vec}(\boldsymbol{f})^{T}\big{(}\boldsymbol{B}^{2}\otimes\boldsymbol{B}^{1}\big{)}^{T}\boldsymbol{K}_{\beta}\big{(}\boldsymbol{B}^{2}\otimes\boldsymbol{B}^{1}\big{)}\text{vec}(\boldsymbol{f})}_{\text{data
fit}}+\underbrace{\text{log}\>|\boldsymbol{K}_{\beta}|}_{\text{penalty}}\Big{)},\end{split}$
(3.13)

where $|\boldsymbol{A}|$ denotes the determinant of matrix $\boldsymbol{A}$.
Thus, the log determinant of the covariance matrix in (3.9) can be interpreted
as a complexity penalty, or regulariser, which comes into play if the kernel
has a free parameter to control the model complexity (MacKay (2003), chap.
28).

####  3.2.4 Bayesian P-splines

Bayesian P-splines have been developed as an extension of the frequentist
P-splines, in Lang and Brezger (2004); Brezger and Lang (2006), and impose a
spatial dependency on the B-splines coefficients to avoid overfitting.
Bayesian P-splines are obtained by placing the Gaussian Markov Random Field
prior on the B-splines coefficients,

$\displaystyle\beta_{i,j}\>|\>\boldsymbol{\beta}_{-(i,j)},\boldsymbol{\phi}\sim\mathcal{N}\Big{(}\sum_{i^{\prime}=1}^{I}\sum_{j^{\prime}=1}^{J}\frac{\delta_{(i,j),(i^{\prime},j^{\prime})}\>\beta_{i^{\prime},j^{\prime}}}{\delta_{(i,j)}^{+}},\frac{\tau^{2}}{\delta_{(i,j)}^{+}}\Big{)},$
(3.14)

where $\boldsymbol{\beta}_{-(i,j)}$ is the matrix $\boldsymbol{\beta}$ without
entry $(i,j)$, $\delta_{(i,j),(i^{\prime},j^{\prime})}$ evaluates to one if
$(i,j)$ and $(i^{\prime},j^{\prime})$ are neighboring coordinates vertically
or horizontally, and otherwise to zero, and
$\delta_{(i,j)}^{+}=\sum_{i^{\prime}=1}^{I}\sum_{j^{\prime}=1}^{J}\delta_{(i,j),(i^{\prime},j^{\prime})}$.
Here, $\tau$ is a spatially varying precision parameter. This model implies
that each $\beta_{i,j}$ is normally distributed with a mean equal to the
average of its neighbors. We use the Bayesian P-splines prior in benchmark
comparisons, due to their close structural relationship with B-splines
projected GP. In our applications, we accelerated the evaluation of the
likelihood by simplifying the joint likelihood to the pairwise difference
formulation and we placed a sum-to-zero constrain on the parameters to ensure
the identifiability of the pairwise differences as proposed by (Morris et al.,
2019). We used the prior $\tau\sim\text{Cauchy}^{+}(0,1)$. To keep notation
streamlined, we write $\boldsymbol{\phi}=\\{\tau\\}$.

###  3.3 Numerical inference

All inferences using model (3.1-3.2) coupled with different priors on the
random surface $\boldsymbol{f}$ were fitted with RStan version 2.21.0, using
an adaptive Hamiltonian Monte Carlo (HMC) sampler (Stan Development Team,
2020). 8 HMC chains were run in parallel for 1,500 iterations, of which the
first 500 iterations were specified as warm-up. Inferences in simulation
analysis and benchmarking were performed similarly.

##  4 Simulation results

Prior to application to the CDC data, we compared the performance of the
spatial models defined in Section 3.2 at retrieving the mean surface of 2D
count data. In the simulations, we considered observations on a 2D grid
$\\{0,0.02\dots,0.98,1\\}\times\\{0,0.02\dots,0.98,1\\}$, such that there are
$50\times 50$ entries. The observations are generated through a Negative
Binomial model using as mean the exponential of a 2D GP with a squared
exponential kernel. The variance scale of the 2D GP is fixed to $1$ and the
length scale is varied ($\gamma\in\\{0.05,0.25,1\\}$) to generate weakly,
mildly and strongly correlated observations. The training set includes 40%
uniformly sampled observations (i.e., number of observations in the training
set $N=640$). Here, we compare the performance of four models with different
prior specifications on the mean surface, a standard 2D GP, a standard
B-splines surface, Bayesian P-splines and regularised B-splines projected 2D
GP. For the methods using B-splines, the same number of knots are placed along
the two axes.

Figure 3A shows the simulated mean surface in the mildly correlation scenario
obtained with $\gamma=0.25$, Figure 3B the simulated observations and Figure
3C the simulated observations included in the training set. The mean surface
estimated by the standard 2D GP model is shown in Supplementary Figure S8.
Figure 4 presents the estimated mean surface when a standard B-spline surface,
Bayesian P-splines, and regularised B-splines projected 2D GP priors are used
for inference with different numbers of knots.

Figure 3: Midly correlated 2D simulated observations. (A) The mean surface is
defined by the exponential of a 2D GP with a squared exponential kernel
implying mild correlations ($\gamma=0.25$). (B) Count data were simulated from
a negative Binomial model. (C) 40% of the simulated observations were included
in the training set.

We find that the standard 2D GP and the regularised B-splines projected 2D GP
recover well the simulated true mean surface (Figures 4 and Supplementary
Figure S8). Moreover, the regularised B-splines projected 2D GP obtains
equivalent or better predictive performance as the standard 2D GP for a
shorter running time, 73.95% faster on average for $30$ knots (Table 1 and
Supplementary Table S2). In contrast, the standard B-splines and the Bayesian
P-splines approaches overfit the data when the number of knots increase
(Figures 4 and Table 1). Similar results were obtained for the other
correlation scenarios (Figures S9-S12, Table 1 and Supplementary Table S2).

Figure 4: Estimated mean surfaces using non-regularized and regularized
B-splines spatial models. Posterior median of the mean surface estimated by
(A) standard B-splines, (B) Bayesian P-splines and (C) regularised B-splines
projected 2D GP with different number of knots given mildly correlated
simulated data ($\gamma=0.25$). The simulated data are illustrated in Figure
3.

Method | Simulation scenarios
---|---
Weakly correlated | Mildly correlated | Strongly correlated
|  | Mean squared error of the mean surface
Standard 2D GP |  |  |
|  | 0.12 |  0.04$\star$ |  0.05$\star$
Standard B-splines surface
| Number of knots |  |  |
| $10$ | 0.11 | 0.15 | 0.25
| $30$ | 0.46 | 0.61 | 0.91
| $40$ | 0.67 | 0.93 | 1.35
Bayesian P-splines surface
| Number of knots |  |  |
| $10$ | 0.11 | 0.09 | 0.10
| $30$ | 0.14 | 0.13 | 0.12
| $40$ | 0.15 | 0.13 | 0.12
Regularised B-splines projected 2D GP
| Number of knots |  |  |
| $10$ | 0.10 | 0.06 | 0.06
| $30$ | 0.09 | 0.05 | 0.06
| $40$ |  0.08$\star$ | 0.05 | 0.06
|  | runtime in minutes | runtime in minutes | runtime in minutes
|  | (-%longest runtime) | (-%longest runtime) | (-%longest runtime)
Standard 2D GP |  |  |
|  | 24 (-0.00%) | 47 (-0.00%) | 21 (-0.00%)
Standard B-splines surface
| Number of knots |  |  |
| $10$ | 1 (-97.23%) | 1 (-98.24%) | 1 (-96.96%)
| $30$ | 1 (-93.92%) | 1 (-96.98%) | 1 (-93.22%)
| $40$ | 2 (-90.97%) | 2 (-95.46%) | 2 (-88.18%)
Bayesian P-splines
| Number of knots |  |  |
| $10$ | 7 (-69.18%) | 6 (-86.44%) | 4 (-78.81%)
| $30$ | 3 (-85.95%) | 3 (-93.04%) | 4 (-81.05%)
| $40$ | 3 (-86.95%) | 4 (-90.49%) | 4 (-82.95%)
Regularised B-splines projected 2D GP
| Number of knots |  |  |
| $10$ | 5 (-80.58%) | 3 (-93.26%) | 4 (-79.78%)
| $30$ | 6 (-76.60%) | 8 (-83.36%) | 8 (-61.89%)
| $40$ | 6 (-76.61%) | 8 (-82.00%) | 8 (-60.02%)

  * $\star$

Best predictive performance.

Table 1: Performance of four spatial models on simulated data. (Top) Means and
standard deviations of the squared errors between the posterior median
estimates of the mean surface and the true values, evaluated on the test data
withheld from fitting. (Bottom) Proportional decreases in runtimes compared to
the slowest model (for which the proportional decrease is null).

##  5 Benchmark results

We next benchmarked the performance of the spatial models defined in Section
3.2 on data unrelated to COVID-19, which has previously been used to benchmark
spatial models. The data set consists of more than 83,000 locations across
East Africa with measured deviations in land surface temperatures (Ton et al.,
2018). We use the same training data as in Mishra et al. (2020), consisting of
6,000 uniformly sampled locations, and fit them using a Gaussian likelihood
with mean surface specified by the various spatial models, and a free
observation variance parameter. We compared standard 2D B-splines, Bayesian
P-splines, regularised B-splines projected 2D GPs, low rank Gaussian Markov
random fields, neural network models, and the $\pi$VAE model (Mishra et al.,
2020) as a prior for the mean surface. For the spline methods, $125$
equidistant knots were placed along both axes. The results are summarised in
Supplementary Table S3 and Figure S13. We find that the regularised B-splines
projected 2D GP model (testing MSE: 2.96) obtained better predictive
performance than the standard B-splines surface model (testing MSE: 4.45), a
low rank Gaussian Markov random field (testing MSE: 4.36), and a neural
network model (testing MSE: 14.94) and similar predictive performance than
Bayesian P-splines (testing MSE: 2.57) and a standard 2D GP (testing MSE:
2.47). It had a worse predictive performance than the $\pi$VAE model (testing
MSE: 0.38), but training $\pi$VAE is not trivial and requires considerable
computational runtimes.

##  6 Time trends of age-specific COVID-19 deaths in the US

Based on the simulation and benchmark results, we used the regularised
B-splines projected 2D GP model to reconstruct the time trends in age-specific
COVID-19 attributable deaths across the US. We placed $12$ equidistant knots
over the age axis and $10$ over the week axis. This choice was determined such
that the predictive performance did not significantly increase with more
knots. Overall, fitting the regularised B-splines projected 2D GP model took
48 hours on a high performance computing environment, with typical effective
sample sizes between 342 to 22127, and there were no reported divergences in
Stan’s Hamiltonian Monte Carlo algorithm. Here we focus on studying the
recovered trends in age-specific COVID-19 deaths during the summer 2021
resurgence in relation to vaccine coverage. Vaccinate coverage is reported by
age strata that do not match the age strata in which the weekly deaths are
reported by the CDC, which renders our Bayesian semi-parametric modelling
approach on 1-year age groups advantageous.

Figure 5: Predicted age-specific COVID-19 attributable deaths in the four most
populated US states. Shown are the posterior median (line) and 95% credible
intervals (ribbon) of the predicted weekly COVID-19 attributable deaths
obtained with (3.4). The CDC data are shown with dots. The start of the 2021
summer resurgence period is indicated as dashed vertical line.

###  6.1 Smooth estimates of age-specific COVID-19 attributable deaths over
time

Figure 5 shows the predicted weekly attributable COVID-19 deaths over time
that occurred in the four most populated US states. The predictions for each
$1$-year age are aggregated into the age bands $0$-$54$, $55$-$74$ and $75+$
for direct comparison to the CDC data, shown as black dots. The data and
estimates clearly reflect the multiple COVID-19 waves across the four states,
with high numbers of deaths in the elderly.

The CDC data contain many anomalies, including censored and delayed entries,
such as the unrealistically high number of deaths reported for the week
starting on January 23, 2021 (Figure 5). We used the curated all-ages weekly
COVID-19 deaths records reported by JHU to re-distribute the delayed deaths to
earlier weeks (6.1). The posterior median and 95% credible intervals in Figure
5 illustrate how the reporting delays apparent in the CDC data are adjusted by
using the JHU data as external calibration, and how the discretized CDC data
informs our estimates of the age composition of deaths. Using our model in
conjunction with the curated JHU data on all COVID-19 deaths, we can estimate
cumulative mortality for any age band. We find that tragically, as of
September 25, 2021, the cumulative mortality rates in individuals aged 85+ now
exceed 3% in Texas and New York (Supplementary Figure S14).

We validated our estimates of the age profile of COVID-19 attributable deaths
using an external data set, i.e. that was not used to inform the model. The
external data set includes age-specific COVID-19 deaths reported directly by
US states DoH (presented in Section 2). We computed the empirical age-specific
contribution to weekly COVID-19 attributable deaths from the DoH data and
checked if the latter lied inside the 95% posterior credible interval
estimated by our model. We find a 93.82% prediction using the regularised
B-splines projected 2D GP model across the four states, California, Florida,
New York and Texas. State stratification of the prediction performance are
presented in Supplementary Table S4.

For completeness, we also fitted the model using alternative smoothing priors.
Supplementary Section S11.1 shows the estimated age-specific composition of
the COVID-19 attributable deaths in three weeks in Florida, as obtained with a
standard 2D GP smoothing prior, standard 2D B-splines, Bayesian P-splines and
the regularised B-splines projected 2D GP. Moreover, Supplementary Table S4
presents the accuracy of their estimates compared to the external data set,
ranging from 94.31% to 95.29%. Under the model specified with a Bayesian
P-splines and standard 2D B-Splines surface, the estimated 2D surface over
ages and weeks was wiggly and had a large uncertainty. Moreover, better
predictive performances, quantified by the expected log pointwise predictive
density Verity et al. (2020), were obtained for the B-Splines projected 2D GP
compared to all other methods (Supplementary Table S5). Those results motivate
the use of our proposed prior, the B-Splines projected 2D GP.

###  6.2 Strong summer 2021 resurgence in age-specific COVID-19 deaths are
associated with limited vaccine coverage

Since July 2021, COVID-19 cases started to increase substantially (John
Hopkins University of Medicine, 2020) and with them COVID-19 deaths in all age
groups (Figure 5). We define the summer 2021 resurgence as the upward trend in
COVID-19 attributable deaths spanning from July 03, 2021 to September 25, 2021
across the United States, and denote by $w^{\text{start-r}}_{m}$ the week
index corresponding to the first week included in the summer 2021 resurgence
period for state $m$. The week index $w^{\text{start-r}}_{m}$ is defined as
the first week from July 01 2021 for which a 4-week central moving average on
the all-age weekly COVID-19 deaths was increasing. The state-specific start of
the resurgence period is indicated by a dashed vertical line in Figure 5. As
is evident from the empirical all-age data and confirmed in our reporting-
adjusted age-specific estimates, the increase in COVID-19 deaths was highly
heterogeneous across US states, with small increases seen in California and
New York, and large increases seen in Florida and Texas, relative to the
magnitude of COVID-19 deaths in the previous waves.

We hypothesised that the variation in the magnitudes of the summer 2021
resurgence in COVID-19 deaths could be associated with differences in how
comprehensively age-specific populations in each state took up the COVID-19
vaccine offer. To test this hypothesis, we attached as transformed parameter
to the model presented in Section 3.1, age-specific predicted COVID-19 deaths
relative to previous waves over the summer 2021. The relative COVID-19 deaths
for state $m$, in age group $c$, at week $w$ in summer 2021, denoted by
$r_{m,c,w}$, are defined as the ratio of the expected posterior predictive
weekly deaths in week $w$ to the maximum of the expected posterior predictive
weekly deaths attained before the summer 2021 resurgence period in the same
state and age group,

$\displaystyle\begin{split}r_{m,c,w}&=\mu^{\star}_{m,c,w}\>\big{/}\>\mu^{\star,\text{max-
pre-r}}_{m,c}\end{split}$ (6.1)

for $c\in\mathcal{C}$ and $w\geq w^{\text{start-r}}_{m}$, where
$\mu^{\star}_{m,c,w}$ are the expected posterior predictive weekly deaths for
state $m$, in age group $c$, at week $w$ defined in (3.5), and
$\mu^{\star,\text{max-
pre-r}}_{m,c}=\text{max}\big{(}\\{\mu^{\star}_{m,c,w}\\}_{w<w^{\text{start-r}}_{m}}\\}\big{)}$.
For simplicity we here suppress in the notation the conditional dependence on
the data. By dividing the mean posterior predictive weekly deaths by their
maximum value attained on previous waves, the resulting relative COVID-19
deaths are adjusted for state-specific factors that are known to influence the
magnitude of SARS-CoV-2 infections such as population density, stringency of
non-pharmaceutical interventions, and comprehensiveness of behaviour changes.

We then assessed if the relative resurgences in COVID-19 deaths across age
groups, weeks and states could be associated with differences in pre-
resurgence vaccination rates, within a random-effects meta-regression model
across states. The pre-resurgence vaccination rates are defined as the
proportion of fully vaccinated individuals 14 days before the start of the
summer 2021 resurgence period and denoted by $\text{v}_{m,c}^{\text{pre-r}}$
for state $m$ and age group $c$. We formulated the model as follows,

$\displaystyle\begin{split}r_{m,c,w}&\>|\>\xi_{m,c,w},\kappa_{m}\sim\text{Gamma}\big{(}\xi_{m,c,w},\kappa_{m}^{2}\big{)}\\\
\xi_{m,c,w}&=\exp\big{(}\chi_{m,c}+\psi_{m,c}\>(w-w^{\text{start-r}}_{m})\big{)}\\\
\begin{bmatrix}\chi_{m,18-64}\\\
\chi_{m,65+}\end{bmatrix}&=\begin{bmatrix}\chi_{18-64}^{\text{base}}+\chi_{m,18-64}^{\text{state}}\\\
\chi_{65+}^{\text{base}}+\chi_{m,65+}^{\text{state}}\\\
\end{bmatrix}+\begin{bmatrix}\chi^{\text{vacc}}&\chi^{\text{vacc-
cross}}_{18-64}\\\ \chi^{\text{vacc-cross}}_{65+}&\chi^{\text{vacc}}\\\
\end{bmatrix}\begin{bmatrix}\text{v}_{m,18-64}^{\text{pre-r}}\\\
\text{v}_{m,65+}^{\text{pre-r}}\\\ \end{bmatrix}\\\
\begin{bmatrix}\psi_{m,18-64}\\\
\psi_{m,65+}\end{bmatrix}&=\begin{bmatrix}\psi_{18-64}^{\text{base}}\\\
\psi_{65+}^{\text{base}}\\\
\end{bmatrix}+\begin{bmatrix}\psi^{\text{vacc}}&\psi^{\text{vacc-
cross}}_{18-64}\\\ \psi^{\text{vacc-cross}}_{65+}&\psi^{\text{vacc}}\\\
\end{bmatrix}\begin{bmatrix}\text{v}_{m,18-64}^{\text{pre-r}}\\\
\text{v}_{m,65+}^{\text{pre-r}}\\\ \end{bmatrix}\end{split}$ (6.2)

for $c\in\\{[18-64],[65+]\\}$. The $\chi_{m,c}$ correspond to state- and age-
specific baseline terms at the start of the summer resurgence, and the
$\psi_{m,c}$ capture state- and age-specific growth rates. Note that the model
also allows for indirect effects of vaccination rates in the other age group.
Random effects and fixed effects (in terms of vaccination rate) were specified
to the baseline terms, and fixed effects to the growth rates,

$\displaystyle\begin{split}\chi_{c}^{\text{base}},\psi_{c}^{\text{base}}&\sim\mathcal{N}(0,0.5),\\\
\chi_{m,c}^{\text{state}}&\sim\mathcal{N}(0,\sigma_{c,\chi}^{2}),\\\
\chi^{\text{vacc}},\psi^{\text{vacc}}&\sim\mathcal{N}(0,0.5),\\\
\chi^{\text{vacc-cross}}_{c},\psi^{\text{vacc-
cross}}_{c}&\sim\mathcal{N}(0,0.5),\\\
\sigma_{c,\chi},\kappa_{m}&\sim\text{Cauchy}^{+}(0,1).\end{split}$ (6.3)

for $c\in\\{[18-64],[65+]\\}$. The model structure was motivated by the fact
that by the start of the summer resurgence, the epidemics in each state had
widely different magnitudes, likely depending on several state-specific
unobserved factors, such as non-pharmaceutical intervention measures, mask
use, or mobility, which we aimed to capture with the
$\chi_{m,c}^{\text{state}}$ baseline random effects. Next, we hypothesized
that after accounting for baseline differences, deaths would grow at smaller
rates in populations with higher vaccine coverage at the start of the summer
resurgence. To summarize, in the meta-regression model, the growth in deaths
during the summer 2021 resurgence is expressed with initial values in terms of
baseline epidemic magnitude and initial vaccination rates, which in turn
characterise state-level exponential growth rates of this highly transmissible
virus over the resurgence period. We found that model (6.2) without state-
specific random effects on the growth rate already fitted the data very well
(Supplementary Figures S20-S21), suggesting a clear negative association of
growth rates in deaths and vaccination coverage (after accounting for baseline
difference at the start of resurgence, and after accounting for long-term
differences in epidemic magnitude across states through $r_{m,c,w}$).

Figure 6: Resurgence in COVID-19 deaths and pre-resurgence vaccination
coverage during the summer 2021 in four states. Posterior median estimates
(line) along with 95% credible intervals (ribbon) of the relative COVID-19
deaths (6.1) over the summer 2021 resurgence period are shown on the y-axis
for four states. The relative death estimates for each age group were
aggregated to the reporting strata of the vaccination data. Week indices of
the summer 2021 resurgence are shown on the x-axis, and correspond to
different calendar weeks for each state (see text). The pre-resurgence
vaccination rate are shown in color. Row facets show resurgences in relative
COVID-19 deaths in individuals aged 18-64 and 65+. Column facets show vaccine
coverage in individuals aged 18-64 and 65+.

Figure 6 shows on the x-axis the week indices of the summer 2021 resurgence
period and in color the state-specific pre-resurgence vaccination rate among
individuals aged $18$-$64$ (left) and $65+$ (right). The estimated relative
COVID-19 deaths (6.1) are shown on the y-axis, now aggregated by the reporting
strata of the vaccination data. We show results for four states for visual
clarity; please see Supplementary Figures S15-S19 for results across the 10
most populated US states California, Florida, Georgia, Illinois, Michigan, New
York, North Carolina, Ohio, Pennsylvania and Texas. Across US states, we
observe that the resurgences in the weekly, relative COVID-19 attributable
deaths were stronger in individuals aged $18$-$64$ compared to individuals
aged $65+$, in line with higher vaccination rates in the $65+$ age group
compared to the $18$-$64$ age group (note the scales on the y-axes).

We evaluated model (6.2) on the relative COVID-19 deaths of the ten most
populated US states, California, Florida, Georgia, Illinois, Michigan, New
York, North Carolina, Ohio, Pennsylvania and Texas. We find that the age-
specific vaccination rates in individuals aged $18$-$64$ and $65+$ were
statistically significantly associated with a protective effect on averting a
strong summer 2021 resurgence in COVID-19 deaths. These relationships are
visualized in Figure 6 and further in Supplementary Figure S22, which shows
effect sizes as a forest plot, indicating that higher proportions of fully
vaccinated individuals prior to resurgence were statistically significantly
associated with smaller relative COVID-19 deaths among the same age group.
Specifically, the baseline effect of vaccination in the same age group
($\chi^{\text{vacc}}$) was strongly negatively associated with resurgent
deaths, with a Bayesian p-value of being non-negative of 0.025%. Similarly,
the effect of vaccination in the same age group on the rate of exponential
increase in resurgent deaths ($\psi^{\text{vacc}}$) was strongly negatively
associated with resurgent deaths, with a Bayesian p-value of being non-
negative of 0%.

With regard to indirect effects, all the terms $\chi^{\text{vacc-
cross}}_{18-64}$, $\chi^{\text{vacc-cross}}_{65+}$, $\psi^{\text{vacc-
cross}}_{18-64}$, and $\psi^{\text{vacc-cross}}_{65+}$ were negative but not
statistically significantly so in terms of the 95% posterior credible
intervals. Yet, such negative associations are consistent with previous
observations that identified young adults as being the main spreaders of
COVID-19 across all age groups (Monod et al., 2021; Wikle et al., 2020). In
this context, infectious disease theory predicts that reducing infections in
the younger adults, has substantial indirect benefits in averting COVID-19
deaths both in younger adults and the elderly. This rationale has been widely
demonstrated for other highly transmissible pathogens such as seasonal
influenza (Baguelin et al., 2013), which is primarily transmitted through
children and their parents, and so achieving high vaccine coverage in these
populations is optimal to minimize deaths and hospitalisations, and maximize
other cost-benefit metrics. The effect sizes associated with higher vaccine
coverage are thus consistent with infectious disease theory, and underscore
the tremendous importance that younger adults get vaccinated against COVID-19
if they are medically fit to halt future COVID-19 resurgences among their age
group.

Figure 7: Projected avoidable deaths under higher vaccine coverage in
individuals aged 18-64 during the summer 2021 resurgence period in four US
states. In counterfactual scenarios, we predicted the number of COVID-19
deaths assuming that the vaccine coverage in individuals aged $18$-$64$ in all
states had been the same as the maximum vaccine coverage in individuals aged
$18$-$64$ across states ($46$%). Posterior median estimates of the predicted
weekly deaths (line) and 95% credible intervals (ribbon) are shown for the
observed COVID-19 deaths (black) and the counterfactual (purple). The start of
the 2021 summer resurgence is shown as dashed vertical line. In Florida, an
estimated 57.15% [40.96, 70.47] of deaths in the $18$-$64$ age group, and
20.09% [-29.54, 54.04] of deaths in the $65+$ age group could have been
avoided. In Texas, an estimated 39.80% [20.81, 54.85] of deaths in the
$18$-$64$ age group, and 13.41% [-41.91, 49.61] of deaths in the $65+$ age
group could have been avoided.

###  6.3 Projected impact of higher vaccination rates in younger adults

To illustrate the potential impact of the inferred associations between
vaccine coverage in younger adults and resurgent COVID-19 deaths, we
considered counterfactual analyses. Specifically, we project the weekly
COVID-19 deaths in all age groups assuming that the pre-resurgence vaccine
coverage among $18$-$64$ year olds in all states would have been the same as
the maximum pre-resurgence vaccine coverage among $18$-$64$ year olds across
the ten most populated states, and assuming that the estimated associations
are causal. The maximum pre-resurgence vaccine coverage among $18$-$64$ year
olds was reported in New York, and was equal to $46$%. We consider this
counterfactual scenario realistic, as it is a retrospective evaluation rather
than a forecast, and within the observed range of vaccination rates, i. e. not
an extrapolation. In the counterfactual, we did not consider higher
vaccination rates in individuals aged 65+ because the differences in
vaccination rates across states were much smaller in this age group, and we
here aim to highlight the potential consequences of the more substantial
variation in vaccination rates in younger individuals.

Figure 7 summarises the counterfactual projections for the four most populated
states and Supplementary Figures S23 for the remaining six of the ten most
populated US states. We estimate that significantly fewer deaths would have
occurred if vaccine coverage in $18$-$64$ year old in all states would have
been the same as the maximum vaccine coverage in $18$-$64$ year old across
states. Overall, from July 03, 2021 to September 25, 2021 and across the ten
most populous US states, we project that 41.96% [31.57, 50.89] of all deaths
in individuals aged 1 could have been avoided if pre-resurgence vaccine
coverage rates in $18$-$64$ year olds had been $46$%, and 13.41% [-41.91,
49.61] of deaths in individuals aged $65+$. Stratified results by US state are
shown in Table S6.

##  7 Discussion

We present a Bayesian non-parametric modelling approach to estimate and report
longitudinal trends in COVID-19 attributable deaths by 1-year age bands across
the United States. The model is informed by COVID-19 attributable weekly
deaths reported by the CDC, the official agency providing such data at the
finer spatial resolution in the US. In comparison to the crude CDC data, the
predictions are for all weeks, and account for missing and censoring present
in the crude data. The predictions of the model also adjust for reporting
delays as e. g. seen in November to December 2020 in Florida and for under-
reporting as e. g. seen in December 2020 to March 2021 in New York (Figure 5),
by allowing the use of an external dataset for calibration, for which we used
the curated JHU all-age COVID-19 attributable deaths. Importantly, the
Bayesian model predicts COVID-19 attributable deaths by $1$-year age band,
which for example allows reporting of cumulative mortality rates in any age
category of interest. Finally, the model implementation is portable, freely
available and computationally efficient through its implementation in Stan and
the use of B-splines projected Gaussian process priors. For these reasons, we
find that the model is a useful tool to improve the reporting of COVID-19
deaths and facilitate demographic and epidemiologic analyses at finer
granularity.

Our analyses of the US summer 2021 resurgence provide further support for the
hypothesis that lower vaccine coverages in several states, most notably
Florida and Texas, are associated with larger COVID-19 death counts (Sah et
al., 2021; Moghadas et al., 2021). However, there are a large number of
behavioural, governmental and state specific differences that modify this
effect size. We have here attempted to agnostically account for these via
random effect terms, but it should be noted that we do not explicitly model
these factors.

We developed B-splines projected Gaussian process priors to describe the time
trends in the age profile of COVID-19 deaths through a non-parametric 2D
random surface. We placed GP priors on the B-splines coefficients, showed that
the model is again a GP, and derived the B-splines projected kernel function
of the GP. We also showed that the model can be interpreted as regularised
B-splines, and find that the negative determinant of the covariance matrix in
the log-likelihood, often referred to as Occam’s factor (MacKay (2003), chap.
28), plays the same role as the penalty of smoothing splines and P-splines
(O’Sullivan, 1986; Eilers and Marx, 1996). Compared to a standard 2D GP,
improvements in runtimes lie both in the fact that fewer parameters on the 2D
surface are estimated (from $n\times m$ to $I\times J$), and that the
computational complexity in calculating covariance matrices is reduced (from
$\mathcal{O}(2(n\times m)^{3/2})$ to $\mathcal{O}(2(I\times J)^{3/2})$), where
$I$ and $J$ are the number of B-splines on the $x$ and $y$ axes, and $n$ and
$m$ are the number of points on the $x$ and $y$ axes and with $I<n$ and $J<m$.
On simulations and on real data, we find the model induces better smoothness
properties than standard B-splines (d. Boor, 1978) and Bayesian P-splines
(Lang and Brezger, 2004) where these are warranted, and so we find the
B-splines projected Gaussian process priors are likely an appealing addition
to the arsenal of Bayesian regularising priors.

Supplementary Figures and Tables

B-splines and properties of the regularised B-splines projected Gaussian
Process prior 2 parts: (1) Construction of a B-spline curve and a B-spline
surface and, (2) Proof that a baseline kernel function projected with cubic
B-splines is a $\mathcal{C}^{2}$ function.

Modelling COVID-19 weekly deaths 2 parts: (1) Finding weekly death count from
censored and missing cumulative death count from the CDC data and, (2)
Specifying the likelihood form.

Additional analyses (1) Comparison of random surface priors for one US state
and (2) Differences in the pre-vaccination age profile of COVID-19
attributable deaths across states.

Templated Stan model files Stan model to fit regularised B-splines projected
Gaussian Process on the 2D mean surface of count or continuous data.

State level summary Detailed fit summary for every US state

## References

  * Baden et al. (2021) Baden, L. R., Sahly, H. M. E., Essink, B., et al. (2021).  “Efficacy and Safety of the mRNA-1273 SARS-CoV-2 Vaccine.”  New England Journal of Medicine, 384(5): 403–416.
  * Baguelin et al. (2013) Baguelin, M., Flasche, S., Camacho, A., Demiris, N., Miller, E., and Edmunds, W. J. (2013).  “Assessing Optimal Target Populations for Influenza Vaccination Programmes: An Evidence Synthesis and Modelling Study.”  PLoS Medicine, 10(10): e1001527.
  * Blangiardo et al. (2020) Blangiardo, M., Cameletti, M., Pirani, M., Corsetti, G., Battaglini, M., and Baio, G. (2020).  “Estimating weekly excess mortality at sub-national level in Italy during the COVID-19 pandemic.”  PLOS ONE, 15(10): e0240286.
  * Brazeau et al. (2020) Brazeau, N. F., Verity, R., Jenks, S., Fu, H., Whittaker, C., Winskill, P., Dorigatti, I., Walker, P., Riley, S., Schnekenberg, R. P., Hoeltgebaum, H., Mellan, T. A., Mishra, S., Whittles, L., Bhatt, S., Ghani, A. C., Ferguson, N. M., and Okell, L. C. (2020).  “Report 34 - COVID-19 Infection Fatality Ratio Estimates from Seroprevalence.”  Imperial College London.
  * Brezger and Lang (2006) Brezger, A. and Lang, S. (2006).  “Generalized structured additive regression based on Bayesian P-splines.”  Computational Statistics & Data Analysis, 50(4): 967–991.
  * Centers for Disease Control and Prevention (2021a) Centers for Disease Control and Prevention (2021a).  “”COVID-19 Vaccinations in the United States”.”  https://covid.cdc.gov/covid-data-tracker.
  * Centers for Disease Control and Prevention (2021b) — (2021b).  “”COVID-19 Vaccinations in the United States, Jurisdiction”.”  https://data.cdc.gov/Vaccinations/COVID-19-Vaccinations-in-the-United-States-Jurisdi/unsk-b7fc.
  * Centers for Disease Control and Prevention (2021c) — (2021c).  “”Weekly Updates by Select Demographic and Geographic Characteristics”.”  https://www.cdc.gov/nchs/nvss/vsrr/covid_weekly/index.htm.
  * Chen et al. (2020) Chen, Y.-C., Lu, P.-E., Chang, C.-S., and Liu, T.-H. (2020).  “A Time-Dependent SIR Model for COVID-19 With Undetectable Infected Persons.”  IEEE Transactions on Network Science and Engineering, 7(4): 3279–3294.
  * d. Boor (1978) d. Boor, C. (1978).  A Practical Guide to Splines.  New York: Springer Verlag.
  * Eilers et al. (2006) Eilers, P. H., Currie, I. D., and Durbán, M. (2006).  “Fast and compact smoothing on large multidimensional grids.”  Computational Statistics & Data Analysis, 50(1): 61–76.
  * Eilers and Marx (1996) Eilers, P. H. C. and Marx, B. D. (1996).  “Flexible smoothing with B-splines and penalties.”  Statistical Science, 11(2).
  * Flaxman et al. (2020) Flaxman, S., Mishra, S., Gandy, A., Unwin, H. J. T., Mellan, T. A., Coupland, H., Whittaker, C., Zhu, H., Berah, T., Eaton, J. W., Monod, M., Ferguson, N. M., Okell, L. C., Bhatt, S., et al. (2020).  “Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe.”  Nature, 584(7820): 257–261.
  * Goldman (2002) Goldman, R. (2002).  Pyramid Algorithms.  Morgan Kaufmann.
  * Gönen and Alpaydin (2011) Gönen, M. and Alpaydin, E. (2011).  “Multiple Kernel Learning Algorithms.”  Journal of Machine Learning Research, 12(64): 2211–2268.
  * Haas et al. (2021) Haas, E. J., Angulo, F. J., McLaughlin, J. M., Anis, E., Singer, S. R., Khan, F., Brooks, N., Smaja, M., Mircus, G., Pan, K., Southern, J., Swerdlow, D. L., Jodar, L., Levy, Y., and Alroy-Preis, S. (2021).  “Impact and effectiveness of mRNA BNT162b2 vaccine against SARS-CoV-2 infections and COVID-19 cases, hospitalisations, and deaths following a nationwide vaccination campaign in Israel: an observational study using national surveillance data.”  The Lancet, 397(10287): 1819–1829.
  * Institute for Health Metrics and Evaluation (2020) Institute for Health Metrics and Evaluation (2020).  “Modeling COVID-19 scenarios for the United States.”  Nature Medicine, 27(1): 94–105.
  * John Hopkins University of Medicine (2020) John Hopkins University of Medicine (2020).  “Coronavirus Resource Center.”  Available at https://coronavirus.jhu.edu/data.
  * King et al. (2015) King, A. A., de Cellès, M. D., Magpantay, F. M. G., and Rohani, P. (2015).  “Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to Ebola.”  Proceedings of the Royal Society B: Biological Sciences, 282(1806): 20150347.
  * Lang and Brezger (2004) Lang, S. and Brezger, A. (2004).  “Bayesian P-Splines.”  Journal of Computational and Graphical Statistics, 13(1): 183–212.
  * Lavezzo et al. (2020) Lavezzo, E., Franchin, E., Ciavarella, C., Cuomo-Dannenburg, G., Barzon, L., Vecchio, C. D., Rossi, L., Manganelli, R., Ferguson, N. M., et al. (2020).  “Suppression of a SARS-CoV-2 outbreak in the Italian municipality of Vo’.”  Nature, 584(7821): 425–429.
  * Lax and Terrell (2017) Lax, P. D. and Terrell, M. S. (2017).  Multivariable Calculus with Applications.  Springer International Publishing.
  * Levin et al. (2020) Levin, A. T., Hanage, W. P., Owusu-Boaitey, N., Cochran, K. B., Walsh, S. P., and Meyerowitz-Katz, G. (2020).  “Assessing the age specificity of infection fatality rates for COVID-19: systematic review, meta-analysis, and public policy implications.”  European Journal of Epidemiology, 35(12): 1123–1138.
  * MacKay (2003) MacKay, D. J. C. (2003).  Information Theory, Inference and Learning Algorithms.  Cambridge University Press, 1st edition.
  * Meyerowitz-Katz and Merone (2020) Meyerowitz-Katz, G. and Merone, L. (2020).  “A systematic review and meta-analysis of published research data on COVID-19 infection fatality rates.”  International Journal of Infectious Diseases, 101: 138–148.
  * Mishra et al. (2020) Mishra, S., Flaxman, S., Zhu, H., and Bhatt, S. (2020).  “$\pi$VAE: Encoding stochastic process priors with variational autoencoders.”  arXiv.
  * Moghadas et al. (2021) Moghadas, S. M., Sah, P., Fitzpatrick, M. C., Shoukat, A., Pandey, A., Vilches, T. N., Singer, B. H., Schneider, E. C., and Galvani, A. P. (2021).  “COVID-19 deaths and hospitalizations averted by rapid vaccination rollout in the United States.”  MedRxiv.
  * Monod et al. (2021) Monod, M., Blenkinsop, A., Xi, X., Hebert, D., Bershan, S., Tietze, S., Baguelin, M., Bradley, V. C., Chen, Y., Coupland, H., Filippi, S., Ish-Horowicz, J., McManus, M., Mellan, T., Gandy, A., Hutchinson, M., Unwin, H. J. T., van Elsland, S. L., Vollmer, M. A. C., Weber, S., Zhu, H., Bezancon, A., Ferguson, N. M., Mishra, S., Flaxman, S., Bhatt, S., and Ratmann, O. (2021).  “Age groups that sustain resurging COVID-19 epidemics in the United States.”  Science, 371(6536): eabe8372.
  * Morris et al. (2019) Morris, M., Wheeler-Martin, K., Simpson, D., Mooney, S. J., Gelman, A., and DiMaggio, C. (2019).  “Bayesian hierarchical spatial models: Implementing the Besag York Mollié model in stan.”  Spatial and Spatio-temporal Epidemiology, 31: 100301.
  * O’Sullivan (1986) O’Sullivan, F. (1986).  “A Statistical Perspective on Ill-Posed Inverse Problems.”  Statistical Science, 1(4): 502–518.
  * O’Sullivan (1988) — (1988).  “Fast Computation of Fully Automated Log-Density and Log-Hazard Estimators.”  SIAM Journal on Scientific and Statistical Computing, 9(2): 363–379.
  * Papoulis and Pillai (2002) Papoulis, A. and Pillai, S. U. (2002).  Probability, Random Variables, and Stochastic Processes.  McGraw-Hill Higher Education.
  * Perperoglou et al. (2019) Perperoglou, A., Sauerbrei, W., Abrahamowicz, M., and Schmid, M. (2019).  “A review of spline function procedures in R.”  BMC Medical Research Methodology, 19(1).
URL https://doi.org/10.1186/s12874-019-0666-3

  * Pokharel and Deardon (2021) Pokharel, G. and Deardon, R. (2021).  “Emulation-based inference for spatial infectious disease transmission models incorporating event time uncertainty.”  Scandinavian Journal of Statistics.
  * Rasmussen and Williams (2005) Rasmussen, C. E. and Williams, C. K. I. (2005).  Gaussian Processes for Machine Learning.  The MIT Press.
  * Rearc (2021) Rearc (2021).  “”Provisional COVID-19 Death Counts by Sex, Age, and State — CDC / NCHS”.”  https://github.com/rearc-data/covid-19-death-counts-sex-age-state.
  * Saatçi (2011) Saatçi, Y. (2011).  “Scalable Inference for Structured Gaussian Process Models.”  PhD thesis, University of Cambridge.
  * Sah et al. (2021) Sah, P., Moghadas, S. M., Vilches, T. N., Shoukat, A., Singer, B. H., Hotez, P. J., Schneider, E. C., and Galvani, A. P. (2021).  “Implications of suboptimal COVID-19 vaccination coverage in Florida and Texas.”  The Lancet Infectious Diseases, 21(11): 1493–1494.
  * Shah et al. (2020) Shah, C., Dehmamy, N., Perra, N., Chinazzi, M., Barab’asi, A. L., Vespignani, A., and Yu, R. (2020).  “Finding Patient Zero: Learning Contagion Source with Graph Neural Networks.”  ArXiv.
  * Stan Development Team (2020) Stan Development Team (2020).  “RStan: the R interface to Stan.”  R package version 2.21.2.
URL http://mc-stan.org/

  * Ton et al. (2018) Ton, J.-F., Flaxman, S., Sejdinovic, D., and Bhatt, S. (2018).  “Spatial mapping with Gaussian processes and nonstationary Fourier features.”  Spatial Statistics, 28: 59–78.
  * Townes (2020) Townes, F. W. (2020).  “Review of Probability Distributions for Modeling Count Data.”  arXiv.
  * United States Census Bureau (2018) United States Census Bureau (2018).  “Age and Sex American Community Survey by U.S States, 2018 Census.”  Available at https://www.census.gov/data.html.
  * Vehtari et al. (2016) Vehtari, A., Gelman, A., and Gabry, J. (2016).  “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.”  Statistics and Computing, 27(5): 1413–1432.
  * Verity et al. (2020) Verity, R., Okell, L. C., Dorigatti, I., Winskill, P., Whittaker, C., Imai, N., Cuomo-Dannenburg, G., Thompson, H., Walker, P. G. T., Fu, H., Dighe, A., Griffin, J. T., Baguelin, M., Bhatia, S., Boonyasiri, A., Cori, A., Cucunubá, Z., FitzJohn, R., Gaythorpe, K., Green, W., Hamlet, A., Hinsley, W., Laydon, D., Nedjati-Gilani, G., Riley, S., van Elsland, S., Volz, E., Wang, H., Wang, Y., Xi, X., Donnelly, C. A., Ghani, A. C., and Ferguson, N. M. (2020).  “Estimates of the severity of coronavirus disease 2019: a model-based analysis.”  The Lancet Infectious Diseases, 20(6): 669–677.
URL https://doi.org/10.1016/s1473-3099(20)30243-7

  * Wikle et al. (2020) Wikle, N. B., Tran, T. N.-A., Gentilesco, B., Leighow, S. M., Albert, E., Strong, E. R., Brinda, K., Inam, H., Yang, F., Hossain, S., Chan, P., Hanage, W. P., Messick, M., Pritchard, J. R., Hanks, E., and Boni, M. F. (2020).  “SARS-CoV-2 epidemic after social and economic reopening in three US states reveals shifts in age structure and clinical characteristics.”  MedRxiv.
  * Wilson et al. (2014) Wilson, A. G., Gilboa, E., Nehorai, A., and Cunningham, J. P. (2014).  “Fast Kernel Learning for Multidimensional Pattern Extrapolation.”  Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2, 3626–3634.
  * Zheng et al. (2021) Zheng, J., Wu, X., Fang, F., Li, J., Wang, Z., Xiao, H., Zhu, J., Pain, C., Linden, P., and Xiang, B. (2021).  “Numerical study of COVID-19 spatial–temporal spreading in London.”  Physics of Fluids, 33(4): 046605.

[Acknowledgments] We thank Sam Abbott and two anonymous reviewers for their
helpful comments. This work was supported by the the Imperial College Research
Computing Service, DOI: 10.14469/hpc/2232, the Imperial College COVID-19
Response Fund and the EPSRC through the EPSRC Centre for Doctoral Training in
Modern Statistics and Statistical Machine Learning. S. Bhatt acknowledges The
UK Research and Innovation (MR/V038109/1), the Academy of Medical Sciences
Springboard Award (SBF004/1080), The MRC (MR/R015600/1), The BMGF
(OPP1197730), Imperial College Healthcare NHS Trust-BRC Funding (RDA02), The
Novo Nordisk Young Investigator Award (NNF20OC0059309) and The NIHR Health
Protection Research Unit in Modelling Methodology. O. Ratmann acknowledges the
Bill & Melinda Gates Foundation (OPP1175094).

Supplementary Materials to Regularised B-splines projected Gaussian Process
priors to estimate time-trends of age-specific COVID-19 deaths related to
vaccine roll-out by Monod et al.

###### Contents

  1. 1 Introduction
  2. 2 COVID-19 deaths and vaccination data
    1. 2.1 COVID-19 deaths data
    2. 2.2 COVID-19 vaccination data
    3. 2.3 Other data sets
  3. 3 Methods
    1. 3.1 Modelling weekly COVID-19 attributable deaths by age
    2. 3.2 Modelling age-specific contributions to COVID-19 weekly deaths
      1. 3.2.1 Two-dimensional Gaussian Process
      2. 3.2.2 B-splines surface
      3. 3.2.3 Regularised B-splines projected Gaussian Process
      4. 3.2.4 Bayesian P-splines
    3. 3.3 Numerical inference
  4. 4 Simulation results
  5. 5 Benchmark results
  6. 6 Time trends of age-specific COVID-19 deaths in the US
    1. 6.1 Smooth estimates of age-specific COVID-19 attributable deaths over time
    2. 6.2 Strong summer 2021 resurgence in age-specific COVID-19 deaths are associated with limited vaccine coverage
    3. 6.3 Projected impact of higher vaccination rates in younger adults
  7. 7 Discussion
  8. S8 Supplementary Figures and Tables
  9. S9 B-splines and properties of the regularised B-splines projected Gaussian Process prior
    1. S9.1 Construction of a B-spline curve and a B-spline surface
    2. S9.2 Proof that a base kernel function projected with cubic B-splines is a $\mathcal{C}^{2}$ function
  10. S10 Modelling COVID-19 weekly deaths
    1. S10.1 Finding weekly death count from censored and missing cumulative death count
    2. S10.2 Likelihood form
  11. S11 Additional analyses
    1. S11.1 Comparison of random surface priors for one US state
    2. S11.2 Differences in the pre-vaccination age profile of COVID-19 attributable deaths across states
  12. S12 Templated Stan model files
    1. S12.1 Stan model to fit regularised B-splines projected Gaussian Process on the mean surface of 2D count data
    2. S12.2 Stan model to fit regularised B-splines projected Gaussian Process on the mean surface of 2D continuous data
    3. S12.3 R code to obtain B-splines
  13. S13 State level summary

##  S8 Supplementary Figures and Tables

(a) $\gamma=0.05$

(b) $\gamma=0.25$

(c) $\gamma=1$

Figure S8: Estimated mean surfaces using a 2D GP model. Posterior median of
the mean surface estimated with a 2D GP model given weakly correlated
($\gamma=0.05$), mildly correlated ($\gamma=0.25$) and strongly correlated
($\gamma=1$) simulated data.

Figure S9: Weakly correlated simulated 2D observations. (A) The mean surface
is defined by the exponential of a 2D GP with a squared exponential kernel
implying weak correlations ($\gamma=0.05$). (B) Count data were simulated from
a negative Binomial model. (C) 40% of the simulated observations were included
in the training set.  Figure S10: Strongly correlated simulated 2D
observations. (A) The mean surface is defined by the exponential of a 2D GP
with a squared exponential kernel implying strong correlations ($\gamma=1$).
(B) Count data were simulated from a negative Binomial model. (C) 40% of the
simulated observations were included in the training set.  Figure S11:
Estimated mean surfaces using non-regularized and regularized B-splines
spatial models (weak correlation, $\gamma=0.05$). Posterior median of the mean
surface estimated by (A) standard B-splines, (B) Bayesian P-splines and (C)
regularised B-splines projected 2D GP with different number of knots given
weakly correlated simulated data ($\gamma=0.05$). The simulated data are
illustrated in Figure S9. Figure S12: Estimated mean surfaces using non-
regularized and regularized B-splines spatial models (strong correlation,
$\gamma=1$). Posterior median of the mean surface estimated by (A) standard
B-splines, (B) Bayesian P-splines and (C) regularised B-splines projected 2D
GP with different number of knots given strongly correlated simulated data
($\gamma=1$). The simulated data are illustrated in Figure S10. Figure S13:
Deviation in land surface temperature for East Africa trained on $6000$ random
uniformly chosen locations. (A) Complete data from Ton et al. (2018). (B)
Training data, similar as the one used in Mishra et al. (2020). The remaining
subfigures show the estimated posterior mean surface by (C) a standard
B-splines surface, (D) Bayesian P-Splines, and (E) regularised B-splines
projected 2D GP.

Figure S14: Predicted COVID-19 attributable mortality rates. Posterior median
estimates of the mortality rate as of September 25, 2021 (barplot) and 95%
credible intervals (error bars) are presented. Column facets indicate the age
group among which the COVID-19 mortality rates are predicted.

Figure S15: Retrievable weekly COVID-19 attributable deaths in ten US states.
The retrievable weekly deaths are computed using the first order difference of
reported cumulative deaths by the CDC (Centers for Disease Control and
Prevention, 2021c). Grey cells indicate non-retrievable weekly death counts
due to the censoring of cumulative deaths. Orange bars show the missing weekly
deaths due to non-reported cumulative deaths. Figure S16: Age-specific
proportion of fully vaccinated individuals in ten US states. The proportion of
fully vaccinated individuals over time were retrieved from Centers for Disease
Control and Prevention (2021b). Fully vaccinated individuals are defined as
having received the second dose of a two-dose vaccine or one dose of a single-
dose vaccine. Row facets show the proportion of fully vaccinated individuals
aged $18$-$64$ and $65+$. A vertical line indicates the 50% threshold. Figure
S17: Predicted age-specific COVID-19 attributable deaths in the six US states.
Shown are the posterior median (line) and 95% credible intervals (ribbon) of
the predicted weekly COVID-19 attributable deaths obtained with (3.4). The CDC
data are shown with dots. The start of the 2021 summer resurgence period is
indicated as dashed vertical line. Figure S18: Resurgence in COVID-19 deaths
and pre-resurgence vaccination coverage during the summer 2021 in the ten most
populated states (part 1). Posterior median estimates (line) along with 95%
credible intervals (ribbon) of the relative COVID-19 deaths (6.1) over the
summer 2021 resurgence period are shown on the y-axis for ten states. The
relative death estimates for each age group were aggregated to the reporting
strata of the vaccination data. Week indices of the summer 2021 resurgence are
shown on the x-axis, and correspond to different calendar weeks for each state
(see main text). The pre-resurgence vaccination rate among individuals aged
18-64 are shown in color. The states are ordered from left to right by the
pre-resurgence vaccination rate among individuals aged 18-64. Row facets show
resurgences in relative COVID-19 deaths in individuals aged 18-64 and $65+$.
Figure S19: Resurgence in COVID-19 deaths and pre-resurgence vaccination
coverage during the summer 2021 in the ten most populated states (part 2).
Posterior median estimates (line) along with 95% credible intervals (ribbon)
of the relative COVID-19 deaths (6.1) over the summer 2021 resurgence period
are shown on the y-axis for ten states. The relative death estimates for each
age group were aggregated to the reporting strata of the vaccination data.
Week indices of the summer 2021 resurgence are shown on the x-axis, and
correspond to different calendar weeks for each state (see main text). The
pre-resurgence vaccination rate among individuals aged $65+$ are shown in
color. The states are ordered from left to right by the pre-resurgence
vaccination rate among individuals aged $65+$. Row facets show resurgences in
relative COVID-19 deaths in individuals aged 18-64 and $65+$.  Figure S20:
Posterior predictive checks of model (6.2) (part 1). Posterior median
estimates (line) and 95% credible intervals (ribbon) of the relative COVID-19
deaths during the summer 2021 resurgence period estimated with pre-resurgence
vaccination coverage (6.2) are shown in blue. Posterior median estimates
(line) and 95% credible intervals (ribbon) of the relative COVID-19 deaths
during the summer 2021 resurgence period estimated with weekly deaths data
(3.4), on which (6.2) was fitted, are shown in red. Row facets show
resurgences in relative COVID-19 deaths in individuals aged 18-64 and $65+$.
Figure S21: Posterior predictive checks of model (6.2) (part 2). Posterior
median estimates (line) and 95% credible intervals (ribbon) of the relative
COVID-19 deaths during the summer 2021 resurgence period estimated with pre-
resurgence vaccination coverage (6.2) are shown in blue. Posterior median
estimates (line) and 95% credible intervals (ribbon) of the relative COVID-19
deaths during the summer 2021 resurgence period estimated with weekly deaths
data (3.4), on which (6.2) was fitted, are shown in red. Row facets show
resurgences in relative COVID-19 deaths in individuals aged 18-64 and $65+$.
Figure S22: Forest plot of the effects in model (6.2). Posterior median
estimates (dot) and 95% credible intervals (error bars) of model (6.2)’s
parameters are shown. The parameters predicting the relative COVID-19 deaths
among $18$-$64$ are shown on the left and among $65+$ on the right. Figure
S23: Projected avoidable deaths under higher vaccine coverage in individuals
aged 18-64 during the summer 2021 resurgence period in six US states. In
counterfactual scenarios, we predicted the number of COVID-19 deaths assuming
that the vaccine coverage in individuals aged $18$-$64$ in all states had been
the same as the maximum vaccine coverage in individuals aged $18$-$64$ across
states ($46$%). Posterior median estimates of the predicted weekly deaths
(line) and 95% credible intervals (ribbon) are shown for the observed COVID-19
deaths (black) and the counterfactual (purple). The start of the 2021 summer
resurgence is shown as dashed vertical line.

Method | Simulation scenarios
---|---
Weakly correlated | Mildly correlated | Strongly correlated
|  | Proportion of data inside the 95% credible interval
Standard 2D GP |  |  |
|  | 98.27% | 99.07% | 98.67%
Standard B-splines surface
| Number of knots |  |  |
| $10$ | 98.93% | 98.87% | 98.60%
| $30$ | 99.13% | 98.53% | 97.87%
| $40$ | 98.73% | 98.00% | 97.40%
Bayesian P-splines surface
| Number of knots |  |  |
| $10$ | 99.13% | 99.00% | 98.27%
| $30$ | 99.67% | 98.87% | 98.93%
| $40$ | 99.80% | 98.60% | 99.00%
Regularised B-splines projected 2D GP
| Number of knots |  |  |
| $10$ | 98.47% | 98.93% | 98.47%
| $30$ | 98.53% | 99.00% | 98.53%
| $40$ | 98.73% | 99.00% | 98.60%
|  |  |  |
Table S2: Predictive accuracy of the four spatial models on simulated data.
Proportion of simulated observations from the test data withheld from fitting,
lying inside the estimated 95% credible interval.

Method | Testing MSE
---|---
Results extracted from Mishra et al. (2020)
|  Full rank GP with Matérn $\frac{3}{2}$ kernel | 2.47
|  Low rank SPDE approximation with 1046 basis function and a Matérn $\frac{3}{2}$ kernel | 4.34
| Neural processes | 14.94
Results obtained by the authors
| Standard B-splines surface | 4.45
| Bayesian P-splines surface | 2.57
| Regularised B-splines projected 2D GP | 2.96

Table S3: Benchmark results obtained on 83,000 locations across East Africa
with measured deviations in land surface temperatures.

| Method
---|---
| Standard 2D GP | Standard B-splines | Bayesian P-splines | Regularised B-splines
| surface |  | projected 2D GP
State |  $\%$DoH observations |  $\%$DoH observations |  $\%$DoH observations |  $\%$DoH observations
inside 95% CI | inside 95% CI | inside 95% CI | inside 95% CI
California | 90.91% | 86.36% | 84.09% | 86.36%
Florida | 94.57% | 96.52% | 95.87% | 94.35%
Texas | 97.14% | 97.92% | 97.14% | 96.61%
Average | 94.90% | 95.29% | 94.31% | 93.82%

Table S4: Predictive accuracy of the four spatial prior on age-specific weekly
COVID-19 attributable deaths data. Proportion of the empirical age-specific
contribution to COVID-19 weekly deaths reported by DoH inside the 95% credible
interval estimated by the model presented in Section 3.1 under four prior
specifications. The four prior are placed on the random surface which
specifies the age composition of deaths over time, a standard 2D GP, a
standard B-splines surface and our regularised B-splines projected GP. Note
that age-specific COVID-19 attributable deaths were not reported by the DoH of
New York.

Method
---
Standard 2D GP | Standard B-splines | Bayesian P-splines | Regularised B-splines
surface | projected 2D GP
$\Delta$ELPD (SD) |  $\Delta$ELPD (SD) |  $\Delta$ELPD (SD) |  $\Delta$ELPD (SD)
-22.03 (11.99) | -295.11 (54.45) | -16.44 (13.05) | 0.00 (0.00)

Table S5: Predictive performance of the four spatial models on age-specific
weekly COVID-19 attributable deaths. The differences in the expected log
pointwise predictive density (ELPD) (Vehtari et al., 2016), and standard error
of the difference, of each model compared to the best performing model (for
which the difference is null) are presented. A higher expected log pointwise
predictive density suggests a better predictive capacity. Four prior placed on
the random surface which specifies the age composition of deaths over time are
compared, a standard 2D GP, a standard B-splines surface and our regularised
B-splines projected GP.

State | Avoidable COVID-19 | Percent of COVID-19
---|---|---
attributable deaths | attributable deaths avoidable
Among age group | Among age group
$18$-$64$ | $65+$ |  $18$-$64$ | $65+$
California | 185 [-47, 419] | 86 [-666, 862] | 9.93% [-2.58, 20.90] | 3.09% [-27.64, 26.80]
Florida | 2218 [1531, 2887] | 1495 [-1902, 4613] | 57.15% [40.96, 70.47] | 20.09% [-29.54, 54.04]
Georgia | 810 [597, 999] | 433 [-450, 1222] | 60.28% [45.68, 71.24] | 22.13% [-25.12, 55.42]
Illinois | 123 [51, 195] | 101 [-329, 540] | 29.90% [13.42, 43.82] | 9.33% [-35.78, 41.21]
Michigan | 116 [69, 165] | 101 [-150, 352] | 34.97% [22.15, 46.75] | 10.62% [-17.29, 33.86]
North Carolina | 387 [261, 503] | 336 [-356, 958] | 58.34% [40.99, 71.50] | 21.10% [-24.97, 53.94]
New York | - | - | - | -
Ohio | 234 [145, 319] | 224 [-570, 995] | 48.69% [31.90, 61.81] | 16.30% [-50.94, 56.41]
Pennsylvania | 134 [63, 207] | 122 [-474, 734] | 27.87% [13.96, 40.11] | 8.58% [-40.71, 42.66]
Texas | 1631 [790, 2452] | 765 [-2037, 3500] | 39.80% [20.81, 54.85] | 13.41% [-41.91, 49.61]
Total | 5850 [4341, 7251] | 3730 [-2958, 9585] | 41.96% [31.57, 50.89] | 14.39% [-11.86, 35.30]

Table S6: Projected avoidable deaths under higher vaccine coverage in
individuals aged 18-64 during the summer 2021 resurgence period. In
counterfactual scenarios, we predicted the number of COVID-19 deaths assuming
that the vaccine coverage in individuals aged $18$-$64$ in all states had been
the same as the maximum vaccine coverage in individuals aged $18$-$64$ across
states ($46$% reported in New York). These counterfactuals consider that the
inferred associations between resurgent deaths and vaccine coverage are
causal, which may be incorrect. Posterior median estimates and 95% credible
interval of the projected avoidable deaths (left) and proportion of deaths
that were avoidable (right) over the summer 2021 resurgence period in the 10
most populated US states are shown.

##  S9 B-splines and properties of the regularised B-splines projected
Gaussian Process prior

###  S9.1 Construction of a B-spline curve and a B-spline surface

We start this section by constructing a B-spline basis function, also called
B-spline segment, which we refer to as B-spline. A B-spline is entirely
defined by a polynomial degree, $d$, and a non-decreasing sequence of knots
$(t_{1}\leq\dots\leq t_{K})$ spanning in $\mathcal{X}\subseteq\mathbb{R}$,
where $K$ is the number of knots. For the B-spline to cover the entire span of
knots, the sequence of knots has to be extended on the left ($d$ times) and on
the right ($d$ times),
$(\xi_{1}\leq\dots\leq\xi_{d}\leq\xi_{d+1}\leq\dots\leq\xi_{d+K}\leq\xi_{d+K+1}\leq\dots\leq\xi_{2d+K+1})$,
where $\xi_{d+1}:=t_{1},\dots,\xi_{d+K}:=t_{K}$. The value of the additional
knots is generally set to the boundary knots (Perperoglou et al., 2019). The
order of a B-spline, denoted by $p$, is equal to $d-1$. A B-spline of order
$p$, written $B_{i,p}(x)$, is defined by the recursive formula,

$\displaystyle\begin{split}B_{i,0}(x)&=\begin{cases}1&\text{if }\xi_{i}\leq
x<\xi_{i+1},\\\ 0&\text{otherwise}\end{cases}\\\
B_{i,p}(x)&=\frac{x-\xi_{i}}{\xi_{i+p}-\xi_{i}}B_{i,p-1}(x)+\frac{\xi_{i+p+1}-x}{\xi_{i+p+1}-\xi_{i+1}}B_{i+1,p-1}(x)\end{split}$
(S9.1)

and $B_{i,0}(x)=0$ if $\xi_{i}=\xi_{i+1}$, for $x\in\mathcal{X}$ and
$i=1,\dots,I$ with $I=K+d-1$. The above is usually referred to as the Cox-de
Boor recursion formula (d. Boor, 1978).

Let us define $\mathcal{C}^{k}(A;B)$ as the set of $\mathcal{C}^{k}$ functions
from $A$ to $B$.

Proposition 1.  Let a B-spline of degree $d$, constructed in (S9.1) given a
vector of strictly increasing sequence knots (no knot occurs more than $1$
time), spanning over $\mathcal{X}$. Then, the B-spline is a
$C^{d-1}(\mathcal{X};\mathbb{R}^{+})$.

###### Proof.

Recall from Goldman (2002) that a B-spline of degree $d$, constructed in
(S9.1) given a vector of strictly increasing sequence knots, spanning over
$\mathcal{X}$ is of continuity $C^{\infty}(\mathcal{X};\mathbb{R}^{+})$
piecewise, and of continuity $C^{d-1}(\mathcal{X};\mathbb{R}^{+})$ on the
knots. Therefore, the same B-spline is $C^{d-1}(\mathcal{X};\mathbb{R}^{+})$
everywhere. ∎

Next, we construct a B-spline curve of degree $d$, a linear combination of
coefficients and B-splines of degree $d$, given by,

$\displaystyle\sum_{i=1}^{I}\beta_{i}\>B_{i,d-1}(x)$ (S9.2)

where $\boldsymbol{\beta}=[\beta_{i}]_{i=1,\dots,I}$ are the associated
B-spline coefficients, also referred to as control points or de Boor points.

The surface analogue to the B-spline curve is the B-spline surface. The
surface is defined on the space $\mathcal{X}_{1}\times\mathcal{X}_{2}$ with
$\mathcal{X}_{1}\subseteq\mathbb{R}$ and $\mathcal{X}_{2}\subseteq\mathbb{R}$.
We start by defining two non decreasing sequences of knots,
$\boldsymbol{t}^{1}=(t^{1}_{1}\leq\dots\leq t^{1}_{K_{1}})$ and
$(t^{2}_{1}\leq\dots\leq t^{2}_{K_{2}})$ placed respectively on the space
$\mathcal{X}_{1}$ and $\mathcal{X}_{2}$, where $K_{1}$ is the number of knots
along the first axis, and $K_{2}$ along the second axis. The corresponding
integral of the B-spline surface evaluated at $x_{1}\in\mathcal{X}_{1}$ and
$x_{2}\in\mathcal{X}_{2}$ is given by,

$\displaystyle\sum_{i=1}^{I}\sum_{j=1}^{J}\>\beta_{(i,j)}\>B^{1}_{i,d_{1}-1}(x_{1})\>B^{2}_{j,d_{2}-1}(x_{2}),$
(S9.3)

where $B^{1}_{i,d_{1}-1}(.)$ is the $i$th B-spline of degree $d_{1}$ defined
over the space $\mathcal{X}_{1}$ on the knot vector $\boldsymbol{t}^{1}$.
Similarly, $B^{2}_{j,d_{2}-1}(.)$ is the $j$th B-spline of degree $d_{2}$
defined over the space $\mathcal{X}_{2}$ on the knot vector
$\boldsymbol{t}^{2}$. The total number of B-splines defined over the space
$\mathcal{X}_{1}$ is $I\triangleq K_{1}+d_{1}-1$ and over the space
$\mathcal{X}_{2}$ it is $J\triangleq K_{2}+d_{2}-1$. We denote the ensemble of
pairs of B-spline’s indices by

$\displaystyle\boldsymbol{U}\triangleq\big{(}\boldsymbol{u}_{1},\dots,\boldsymbol{u}_{M}\big{)}=\\{1,\dots,I\\}\times\\{1,\dots,J\\},$
(S9.4)

with $M=I\times J$.

###  S9.2 Proof that a base kernel function projected with cubic B-splines is
a $\mathcal{C}^{2}$ function

We follow the notations introduced in Section S9.1 and denote the $i$th and
$j$th cubic B-spline by $B^{1}_{i}(x_{1})$ and $B^{2}_{j}(x_{2})$, with
$i=1,\dots I$ and $j=1,\dots,J$111Notice that we dropped the polynomial degree
index as it is fixed to $3$. and where
$x_{1}\in\mathcal{X}_{1}\subseteq\mathbb{R}$ and
$x_{2}\in\mathcal{X}_{2}\subseteq\mathbb{R}$. They are respectively defined on
two vectors of strictly increasing knots of lengths $K_{1}$ and $K_{2}$. For
ease of notation, we denote the ensemble of B-spline indices by
$\boldsymbol{U}$ defined in (S9.4).

Let us define a low-rank covariance matrix $\boldsymbol{K}_{\beta}$ of size
$IJ\times IJ$ obtained by evaluating at every pair of points in
$\boldsymbol{u},\boldsymbol{u}^{\prime}\in\boldsymbol{U}$ a base kernel
function $k_{\beta}$,

$\displaystyle\big{(}\boldsymbol{K}_{\beta}\big{)}_{\boldsymbol{u},\boldsymbol{u}^{\prime}}=k_{\beta}\big{(}\boldsymbol{u},\boldsymbol{u}^{\prime}\big{)}.$
(S9.5)

Let us define a high-rank covariance matrix projected by cubic B-splines by
$\boldsymbol{W}$, obtained by evaluating at every pair of points a kernel
function $w$,
$\boldsymbol{W}_{(x_{1},x_{2}),(x_{1}^{\prime},x_{2}^{\prime})}=w(x_{1},x_{2},x_{1}^{\prime},x_{2}^{\prime})$
with,

$w(x_{1},x_{2},x_{1}^{\prime},x_{2}^{\prime})=\\\
\sum_{(u_{1},u_{2})\in\boldsymbol{U}}\sum_{(u_{1}^{\prime},u_{2}^{\prime})\in\boldsymbol{U}}D(x_{1},x_{2},u_{1},u_{2})\>k_{\beta}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})\>D(x_{1}^{\prime},x_{2}^{\prime},u_{1}^{\prime},u_{2}^{\prime}).$
(S9.6)

where
$(x_{1},x_{2}),(x_{1}^{\prime},x_{2}^{\prime})\in\mathcal{X}_{1}\times\mathcal{X}_{2}$,
and

$D(x_{1},x_{2},u_{1},u_{2})=B_{u_{2}}^{2}(x_{2})\>B_{u_{1}}^{1}(x_{1}).$
(S9.7)

We repeat here the definition presented in S9.1: we define
$\mathcal{C}^{k}(A;B)$ as the set of $\mathcal{C}^{k}$ functions from $A$ to
$B$.

Proposition 2.  Let $w(\cdot)$ be defined as in (S9.6) and
$\mathcal{S}=\mathcal{X}_{1}\times\mathcal{X}_{2}\times\mathcal{X}_{1}\times\mathcal{X}_{2}\subseteq\mathbb{R}^{4}$.
Then, $w(\cdot)$ is a $\mathcal{C}^{2}(\mathcal{S};\mathbb{R})$ function.

###### Proof.

In order to show that
$w(x_{1},x_{2},x_{1}^{\prime},x_{2}^{\prime})\in\mathcal{C}^{2}(\mathcal{S};\mathbb{R})$
it is sufficient to show that all its second partial derivatives are
continuous functions on $\mathcal{S}$. We show this for the case of the
partial derivative $\frac{\partial^{2}w}{\partial x_{1}^{2}}$ and the mixed
derivative $\frac{\partial^{2}w}{\partial x_{1}\partial x_{2}}$, as the proof
for the remaining second partial derivatives is identical. The proof is
constructed as follows: (1) we compute $\frac{\partial^{2}w}{\partial
x_{1}^{2}}$ and $\frac{\partial^{2}w}{\partial x_{1}\partial x_{2}}$ and (2)
we show that $\frac{\partial^{2}w}{\partial x_{1}^{2}}$ and
$\frac{\partial^{2}w}{\partial x_{1}\partial x_{2}}$ are continuous.

Step (1): A direct computation shows,

$\frac{\partial^{2}w}{\partial x_{1}^{2}}=\\\
\sum_{(u_{1},u_{2})\in\boldsymbol{U}}\sum_{(u_{1}^{\prime},u_{2}^{\prime})\in\boldsymbol{U}}B^{2}_{u_{2}}(x_{2})\frac{\partial^{2}B^{1}_{u_{1}}}{\partial
x_{1}^{2}}(x_{1})k_{\beta}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})D(x_{1}^{\prime},x_{2}^{\prime},u_{1}^{\prime},u_{2}^{\prime}),$
(S9.8) $\frac{\partial w}{\partial x_{1}\partial x_{2}}=\\\
\sum_{(u_{1},u_{2})\in\boldsymbol{U}}\sum_{(u_{1}^{\prime},u_{2}^{\prime})\in\boldsymbol{U}}\frac{\partial
B^{2}_{u_{2}}}{\partial x_{2}}(x_{2})\frac{\partial B^{1}_{u_{1}}}{\partial
x_{1}}(x_{1})k_{\beta}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})D(x_{1}^{\prime},x_{2}^{\prime},u_{1}^{\prime},u_{2}^{\prime}),$
(S9.9)

where we have used the expression of $D$ as stated in (S9.7) and
$\boldsymbol{U}$ is defined in (S9.4).

Step (2): We show that (S9.8) is a continuous function on $\mathcal{S}$. As
stated in Proposition 1, the functions $B^{1}_{u_{1}}(\cdot)$ and
$B^{2}_{u_{2}}(\cdot)$ are $\mathcal{C}^{2}(\mathcal{X}_{1};\mathbb{R}^{+})$
and $\mathcal{C}^{2}(\mathcal{X}_{2};\mathbb{R}^{+})$ functions, respectively.
Therefore, $\frac{\partial^{2}B^{1}_{u_{1}}}{\partial x_{1}^{2}}(\cdot)$ and
$B^{2}_{u_{2}}(x_{2})$ are continuous function on $\mathcal{X}_{1}$ and
$\mathcal{X}_{2}$, respectively. Moreover, recall that $k_{\beta}$ is a kernel
function, so for any fixed
$(u_{1},u_{2}),(u_{1}^{\prime},u_{2}^{\prime})\in\boldsymbol{U}$,
$k_{\beta}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})\in\mathbb{R}$. By
standard results of Multivariate Calculus (Lax and Terrell (2017) Theorem
2.5), for any fixed $u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime}$ the term

$B^{2}_{u_{2}}(x_{2})\frac{\partial^{2}B^{1}_{u_{1}}}{\partial
x_{1}^{2}}(x_{1})k_{\beta}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})B^{2}_{u_{2}^{\prime}}(x_{2}^{\prime})B^{1}_{u_{1}^{\prime}}(x_{1}^{\prime})$

is a continuous function from $\mathcal{S}$ to $\mathbb{R}$. Moreover, the RHS
of (S9.8) is a linear combination of continuous function from $\mathcal{S}$ to
$\mathbb{R}$, therefore it is continuous from $\mathcal{S}$ to $\mathbb{R}$.
Hence, $\frac{\partial^{2}w}{\partial x_{1}^{2}}$ is a continuous function
from $\mathcal{S}$ to $\mathbb{R}$.

Now, we show that (S9.9) is a continuous function on $\mathcal{S}$. As stated
in Proposition 1, the functions $B^{1}_{u_{1}}(\cdot)$ and
$B^{2}_{u_{2}}(\cdot)$ are $\mathcal{C}^{2}(\mathcal{X}_{1};\mathbb{R}^{+})$
and $\mathcal{C}^{2}(\mathcal{X}_{2};\mathbb{R}^{+})$ functions, respectively.
Therefore, $\frac{\partial B^{1}_{u_{1}}}{\partial x_{1}}(\cdot)$ and
$\frac{\partial B^{2}_{u_{2}}}{\partial x_{2}}(\cdot)$ are continuous
functions on $\mathcal{X}_{1}$ and $\mathcal{X}_{2}$, respectively. Moreover,
recall that $k_{\beta}$ is a kernel function, so for any fixed
$(u_{1},u_{2}),(u_{1}^{\prime},u_{2}^{\prime})\in\boldsymbol{U}$,
$k_{\beta}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})\in\mathbb{R}$. By
standard results of Multivariate Calculus (Lax and Terrell (2017) Theorem
2.5), for any fixed $u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime}$ the term

$\frac{\partial B^{2}_{u_{2}}}{\partial x_{2}}(x_{2})\frac{\partial
B^{1}_{u_{1}}}{\partial
x_{1}}(x_{1})k_{\beta}(u_{1},u_{2},u_{1}^{\prime},u_{2}^{\prime})B^{2}_{u_{2}^{\prime}}(x_{2}^{\prime})B^{1}_{u_{1}^{\prime}}(x_{1}^{\prime})$

is a continuous function from $\mathcal{S}$ to $\mathbb{R}$. Moreover, the RHS
of (S9.9) is a linear combination of continuous functions on from
$\mathcal{S}$ to $\mathbb{R}$, therefore it is continuous from $\mathcal{S}$
to $\mathbb{R}$. Hence, $\frac{\partial w}{\partial x_{1}\partial x_{2}}$ is a
continuous function from $\mathcal{S}$ to $\mathbb{R}$.

The proof for the other terms follows from a similar argument, with the
necessary replacement of sets and functions. Therefore $w(\cdot)$ is a
$\mathcal{C}^{2}(\mathcal{S};\mathbb{R})$ function. ∎

##  S10 Modelling COVID-19 weekly deaths

###  S10.1 Finding weekly death count from censored and missing cumulative
death count

We denote by $\mathcal{W}^{\text{cum}}=\\{1,\dots,W^{\text{cum}}\\}$ the set
of week indices from the first to the last reported cumulative deaths by the
CDC. The reports span from May 2, 2020 to October 02, 2021 such that
$W^{\text{cum}}=75$. We denote by $D_{m,b,w}$ the reported cumulative deaths
by the CDC in state $m$ for age group $b$ in week $w$. The cumulative deaths
count $D_{m,b,w}$ is reported if the count count is 0 or strictly greater than
9, and otherwise it is censored. For simplicity we suppress the state index in
what follows, with all equations being analogous. We denote by
$\mathcal{W}_{b}^{\text{censored}}$ the set of weeks where $D_{b,w}$ has been
censored.

For a given age group $b$, if $D_{b,w}$ is reported on weeks $w$ and $w+1$, we
find the weekly deaths among age group $b$ in week $w$, $d_{b,w}$, with the
first order difference,

$\displaystyle d_{b,w}=D_{b,(w+1)}-D_{b,w}$ (S10.1)

The weekly deaths that are retrievable from the first order difference are
defined as “retrievable” and the set of week indices for which we can retrieve
the weekly deaths is denoted by $\mathcal{W}_{b}^{\text{WR}}$ (WR := Week
Retrievable). However, if $D_{b,w}$ is censored for $w$ or $w+1$, the first-
order difference is not obtainable and the weekly deaths is said to be “non-
retrievable”. We denote by $\mathcal{W}_{b}^{\text{WNR}}$ the set of week
indices for which we cannot retrieve the weekly deaths (WNR := Week non-
retrievable), with
$\mathcal{W}_{b}^{\text{WR}}\cap\mathcal{W}_{b}^{\text{WNR}}=\varnothing$.
Notice that for a fixed age group $b$, weekly deaths are non-retrievable for a
joint sequence of weeks over time, and importantly this period can happen only
once (because the cumulative deaths are strictly increasing). Therefore, we
consider the sequence of non-retrievable weekly deaths for a fixed age group
together. Because the boundaries of the censored cumulative death are known,
we can obtain the boundaries of the sum of the non-retrievable weekly deaths
sequence. Let us denote by $w_{b}^{\text{FC}}$ the first week when the
cumulative death count is censored (FC := First Censored). Such that,

$\displaystyle
w_{b}^{\text{FC}}=\begin{cases}\text{min}(\mathcal{W}_{b}^{\text{censored}})&\text{if
}\mathcal{W}_{b}^{\text{censored}}\neq\varnothing\\\ \text{does not
exist}&\text{otherwise}.\end{cases}$ (S10.2)

Similarly, we denote by $w_{b}^{\text{FNC}}$ the first week when the
cumulative death count is observed after being censored (FNC := First Non
Censored). Such that,

$\displaystyle
w_{b}^{\text{FNC}}=\begin{cases}\text{max}(\mathcal{W}_{b}^{\text{censored}})+1&\text{if
}\text{max}(\mathcal{W}_{b}^{\text{censored}})<W^{\text{cum}}\\\ \text{does
not exist}&\text{otherwise}.\end{cases}$ (S10.3)

Four different scenarios may apply.

  1. 1.

If $D_{b,w}$ is censored for some weeks $w\in\mathcal{W}^{\text{cum}}$ but is
observed on the first and the last weeks (i.e., $w=1$ and $w=W^{\text{cum}}$),
the non-retrievable weekly deaths must sum to the first positive cumulative
deaths, such that,

$\displaystyle\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}d_{b,w}=D_{b,w_{b}^{\text{FNC}}}.$
(S10.4)

where
$\mathcal{W}_{b}^{\text{WNR}}=\\{(w_{b}^{\text{FC}}-1),\mathcal{W}_{b}^{\text{censored}}\\}$.
In this case,
$\mathcal{W}_{b}^{\text{WR}}=\\{1,\dots,(w_{b}^{\text{FC}}-2),w_{b}^{\text{FNC}},\dots,(W^{\text{cum}}-1)\\}$.
We show an example of this scenario in Table S7.

  2. 2.

If $D_{b,w}$ is censored for some weeks $w\in\mathcal{W}^{\text{cum}}$,
including the first week (i.e., $w=1$), but is observed for the last week
(i.e., $w=W^{\text{cum}}$), the sum of non-retrievable weekly deaths must be
between,

$\displaystyle
D_{b,w_{b}^{\text{FNC}}}-9\leq\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}d_{b,w}\leq
D_{b,w_{b}^{\text{FNC}}}-1.$ (S10.5)

where $\mathcal{W}_{b}^{\text{WNR}}=\mathcal{W}_{b}^{\text{censored}}$. In
this case,
$\mathcal{W}_{b}^{\text{WR}}=\\{w_{b}^{\text{FNC}},\dots,(W^{\text{cum}}-1)\\}$.
We show an example of this scenario in Table S8.

  3. 3.

If $D_{b,w}$ is censored for some weeks $w\in\mathcal{W}^{\text{cum}}$,
including the last week (i.e., $w=W^{\text{cum}}$), but is observed for the
first week (i.e., $w=1$), the sum of the non-retrievable weekly deaths must be
between,

$\displaystyle 1\leq\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}d_{b,w}\leq 9$
(S10.6)

where
$\mathcal{W}_{b}^{\text{WNR}}=\\{(w_{b}^{\text{FC}}-1),\mathcal{W}_{b}^{\text{censored}}\setminus
W^{\text{cum}}\\}$. In this case,
$\mathcal{W}_{b}^{\text{WR}}=\\{1,\dots,(w_{b}^{\text{FC}}-2)\\}$. Notice that
$D_{b,w_{b}^{\text{FNC}}}$ does not exist. We show an example of this scenario
in Table S9.

  4. 4.

Lastly, if $D_{b,w}$ is censored for all weeks $w\in\mathcal{W}^{\text{cum}}$,
including the first week (i.e., $w=1$) and the last week (i.e.,
$w=W^{\text{cum}}$), the sum of the non-retrievable weekly deaths must be
between,

$\displaystyle 0\leq\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}d_{b,w}\leq 8$
(S10.7)

where
$\mathcal{W}_{b}^{\text{WNR}}=\\{\mathcal{W}_{b}^{\text{censored}}\setminus
W^{\text{cum}}\\}$. In this case, $\mathcal{W}_{b}^{\text{WR}}=\varnothing$.
Notice that $D_{b,w_{b}^{\text{FNC}}}$ does not exist. We show an example of
this scenario in Table S10.

The CDC did not publish a report on July 4 2020, such that the cumulative
deaths in the corresponding week are not available. We denote this week by
$w^{\text{missing}}$. The weekly deaths could not be obtained for the missing
week or the week before, such that,
$\mathcal{W}_{b}^{\text{WR}}\cup\mathcal{W}_{b}^{\text{WNR}}=\mathcal{W}^{\text{cum}}\setminus\\{w^{\text{missing}}-1,w^{\text{missing}},W^{\text{cum}}\\}$.
Weekly deaths that occurred on those $w^{\text{missing}}$ and
$w^{\text{missing}}-1$ are defined to be missing. It is important to note the
difference between non-retrievable weekly deaths, for which we have
information through the censoring of cumulative deaths, and missing weekly
deaths for which no information has been provided. Note that if the index of
the missing week $w^{\text{missing}}$ and the week before
$w^{\text{missing}}-1$, are in between the range of
$\mathcal{W}_{b}^{\text{censored}}$, or are before or after the set’s limits,
they are included in $\mathcal{W}_{b}^{\text{censored}}$, with all inequations
still holding.

$w$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$D_{w}$ | 0 | 0 | NA | NA | 11
$d_{w}$ | 0 | 11 | -
Table S7: Finding weekly death when $D_{w}$ is censored as in scenario 1. In this example $\mathcal{W}^{\text{censored}}=\\{3,4\\}$, $w^{\text{FC}}=3$, $w^{\text{FNC}}=5$, $\mathcal{W}^{\text{WNR}}=\\{2,3,4\\}$ and $\mathcal{W}^{\text{WR}}=\\{1\\}$. From (S10.4), we know that $\sum_{w\in\mathcal{W}^{\text{WNR}}}d_{w}=11$. $w$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$D_{w}$ | NA | NA | NA | 11 | 11
$d_{w}$ | [2-10] | 0 | -
Table S8: Finding weekly death when $D_{w}$ is censored as in scenario 2. In this example $\mathcal{W}^{\text{censored}}=\\{1,2,3\\}$, $w^{\text{FC}}=1$, $w^{\text{FNC}}=4$, $\mathcal{W}^{\text{WNR}}=\\{1,2,3\\}$ and $\mathcal{W}^{\text{WR}}=\\{4\\}$. From (S10.5), we know that $2\leq\sum_{w\in\mathcal{W}^{\text{WNR}}}d_{w}\leq 10$. $w$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$D_{w}$ | 0 | 0 | NA | NA | NA
$d_{w}$ | 0 | [1-9] | -
Table S9: Finding weekly death when $D_{w}$ is censored as in scenario 3. In this example $\mathcal{W}^{\text{censored}}=\\{3,4,5\\}$, $w^{\text{FC}}=3$, $w^{\text{FNC}}$ does not exist, $\mathcal{W}^{\text{WNR}}=\\{2,3,4\\}$ and $\mathcal{W}^{\text{WR}}=\\{1\\}$. From (S10.6), we know that $1\leq\sum_{w\in\mathcal{W}^{\text{WNR}}}d_{w}\leq 9$. $w$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$D_{w}$ | NA | NA | NA | NA | NA
$d_{w}$ | [0-8] | -
Table S10: Finding weekly death when $D_{w}$ is censored as in scenario 4. In
this example $\mathcal{W}^{\text{censored}}=\\{1,2,3,4,5\\}$,
$w^{\text{FC}}=1$, $w^{\text{FNC}}$ does not exist,
$\mathcal{W}^{\text{WNR}}=\\{1,2,3,4\\}$ and
$\mathcal{W}^{\text{WR}}=\varnothing$. From (S10.7), we know that
$0\leq\sum_{w\in\mathcal{W}^{\text{WNR}}}d_{w}\leq 8$.

###  S10.2 Likelihood form

In this section, we use the notation introduced in Section S10.1. The log
likelihood is,

$\displaystyle\log\mathcal{L}$
$\displaystyle=\sum_{b\in\mathcal{B}}\sum_{w\in\mathcal{W}_{b}^{\text{WR}}}\log
p(d_{b,w}\>;\>\alpha_{b,w},\theta)$ (S10.8a)
$\displaystyle+\sum_{b\in\mathcal{B}}\log
p\big{(}\\{d_{b,w}\\}_{w\in\mathcal{W}_{b}^{\text{WNR}}}\>;\>\\{\alpha_{b,w}\\}_{w\in\mathcal{W}_{b}^{\text{WNR}}},\theta\big{)},$
(S10.8b)

The log likelihood for the retrievable weekly deaths in (S10.8a) is

$\displaystyle\log
p(d_{b,w}\>;\>\alpha_{b,w},\theta)=\log\text{NegBin}(d_{b,w}\>;\>\alpha_{b,w},\theta).$
(S10.9)

The log likelihood for the non-retrievable weekly deaths in (S10.8b) is

$\displaystyle\begin{split}&p\big{(}\\{d_{b,w}\\}_{w\in\mathcal{W}_{b}^{\text{WNR}}}\>;\>\\{\alpha_{b,w}\\}_{w\in\mathcal{W}_{b}^{\text{WNR}}},\theta\big{)}=p\Big{(}\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}d_{b,w}\>;\>\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}\alpha_{b,w},\theta\Big{)}\\\
=&\begin{cases}\text{NegBin}\Big{(}D_{b,w_{b}^{\text{FNC}}}\>;\>\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}\alpha_{b,w},\theta\Big{)}&\text{if
}1,W^{\text{cum}}\notin\mathcal{W}_{b}^{\text{censored}}\\\\[17.07182pt]
\text{Cdf
NegBin}\Big{(}D_{bw_{b}^{\text{FNC}}}-1\>;\>\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}\alpha_{b,w},\theta\Big{)}&\text{if
}1\in\mathcal{W}_{b}^{\text{censored}}\\\\[5.69046pt] -\text{Cdf
NegBin}\Big{(}D_{bw_{b}^{\text{FNC}}}-9\>;\>\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}\alpha_{b,w},\theta\Big{)}&\text{and
}W^{\text{cum}}\notin\mathcal{W}_{b}^{\text{censored}}\\\\[17.07182pt]
\text{Cdf
NegBin}\Big{(}9\>;\>\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}\alpha_{b,w},\theta\Big{)}&\text{if
}1\notin\mathcal{W}_{b}^{\text{censored}}\\\\[5.69046pt] -\text{Cdf
NegBin}\Big{(}1\>;\>\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}\alpha_{b,w},\theta\Big{)}&\text{and
}W^{\text{cum}}\in\mathcal{W}_{b}^{\text{censored}}\\\\[17.07182pt] \text{Cdf
NegBin}\Big{(}8\>;\>\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}\alpha_{b,w},\theta\Big{)}&\text{if
}1,W^{\text{cum}}\in\mathcal{W}_{b}^{\text{censored}}\\\\[5.69046pt]
-\text{Cdf
NegBin}\Big{(}0\>;\>\sum_{w\in\mathcal{W}_{b}^{\text{WNR}}}\alpha_{b,w},\theta\Big{)}&\end{cases}\end{split}$
(S10.10)

where the boundaries for the different scenarios have been discussed in
Section S10.1. Note that to find (S10.10), we used the fact that the sum of
independent negative-binomially distributed random variables with shape
parameter $\alpha_{1}$ and $\alpha_{2}$, and with the same scale parameter
$\theta$ is negative-binomially distributed with a scale parameter $\theta$
and with a shape parameter $\alpha_{1}+\alpha_{2}$.

##  S11 Additional analyses

###  S11.1 Comparison of random surface priors for one US state

We used the non-parametric model (3.1-3.2) to fit age-specific weekly COVID-19
attributable deaths from the CDC data. We compare results, for one state,
Florida, across the four prior specifications on the random surface which
specifies the age composition of deaths over time, a standard 2D GP, a
standard B-splines surface and our regularised B-splines projected GP.

Figure S24 shows the estimated age-specific composition of weekly COVID-19
attributable deaths in three weeks for the four prior specifications. The
posterior distribution of the age profile of COVID-19 attributable deaths
estimated with the standard B-splines and Bayesian P-splines prior is rough
because the model’s complexity is not penalized, or not enough. Such a penalty
is introduced in our regularised B-splines projected GP prior, which results
in a less wiggly posterior distribution with less uncertainty. To quantify the
predictive performance of the model, we used the expected log pointwise
predictive density (ELPD) (Verity et al., 2020) (Supplementary Table S5). The
best performance was obtained by the regularised B-splines projected GP. In
particular, the expected log pointwise predictive density (ELPD) difference,
between a standard GP and a regularised B-splines projected GP was -22.03 (sd
11.99), suggesting significantly better predictive performance with the
regularised B-splines projected GP prior. In addition, a standard GP requires
a longer running time (2.34 times more). Similar results were obtained for all
other states, with the regularised B-splines projected GP prior consistently
outperforming the other standard methods.

Figure S24: Estimated age-specific contribution to weekly COVID-19
attributable deaths by four priors. Posterior median (black line) and 95%
credible interval (ribbon) of the age-specific contribution to COVID-19
attributable deaths in three weeks. Four prior placed on the random surface
which specifies the age composition of deaths over time are compared, a
standard 2D GP, a standard B-splines surface and our regularised B-splines
projected GP.

###  S11.2 Differences in the pre-vaccination age profile of COVID-19
attributable deaths across states

We investigate whether the age composition of COVID-19 attributable deaths
differs across states. To address this question, we focus on the baseline
period of the first $3$ months after the $10$th cumulated death in each state,
and aggregate the predictions of the share of $1$-year age bands from our non-
parametric model into the age strata $0$-$24$, $25$-$54$, $55$-$74$,
$75$-$84$, $85+$ years for each state during the baseline period. In these
calculations, we adjusted for age differences in state populations with the
age composition of the US population. Pre-pandemic age-specific population
counts were drawn from the U.S. Census Bureau in 2018 (United States Census
Bureau, 2018). Strikingly, the resulting posterior estimates of the share of
age groups among COVID-19 attributable deaths were statistically significantly
different across states even after standardising the age composition of state
populations (Figure S25). For example, the contributions of individuals aged
$85+$ to COVID-19 attributable deaths ranged from 26.17% (24.56%- 27.77%) in
Texas to 38.99% (37.49%- 40.53%) in Pennsylvania. High heterogeneity was also
found in the contributions from individuals age $55$-$74$, which varied from
28.85% (27.56%- 30.2%) in Pennsylvania to 37.71% (36.16%- 39.31%) in Texas.

These discrepancies across states are clearly supported in the empirical data,
albeit without an assessment of their statistical significance (Figure S26).
Notice that the contribution from individuals aged $55+$ estimated by the
model is always lower or equal to the empirical contribution while the one in
$25$-$54$ is always greater or equal. This is because the empirical
contributions were computed considering that unretrievable weekly deaths were
null while they could, in fact, be positive.

Figure S25: Estimated age-specific contribution to COVID-19 attributable
deaths adjusted for population size in the first 3 months of state epidemics.
Posterior median estimates of each age group to COVID-19 attributable deaths
are shown for each state (barplot), along with 95% credible intervals (error
bars). Estimates are for the $3$ months after the $10$th cumulative death in
each state, and are adjusted for differences in the age composition of state
populations (see text). Brighter colors indicate higher values of the baseline
contribution and darker colors indicate lower values.  Figure S26: Estimated
and empirical age-specific contribution to COVID-19 attributable deaths in the
first 3 months of state epidemics. Posterior median estimates of each age
group to COVID-19 attributable deaths are shown for each state (barplot),
along with 95% credible intervals (error bars). Estimates are for the $3$
months after the $10$th cumulative death in each state, and are not adjusted
for differences in the age composition of state populations. Brighter colors
indicate higher values of the baseline contribution and darker colors indicate
lower values. The empirical estimates (red points) of contribution are
calculated by setting the non-retrievable weekly deaths to zero.

##  S12 Templated Stan model files

###  S12.1 Stan model to fit regularised B-splines projected Gaussian Process
on the mean surface of 2D count data

⬇

1functions {

2 matrix kron_mvprod(matrix A, matrix B, matrix V)

3 {

4 return transpose(A*transpose(B*V));

5 }

6

7 matrix gp(int N_rows, int N_columns, real[] rows_idx, real[] columns_index,

8 real delta0,

9 real alpha_gp,

10 real rho_gp1, real rho_gp2,

11 matrix z1)

12 {

13

14 matrix[N_rows,N_columns] GP;

15

16 matrix[N_rows, N_rows] K1;

17 matrix[N_rows, N_rows] L_K1;

18

19 matrix[N_columns, N_columns] K2;

20 matrix[N_columns, N_columns] L_K2;

21

22 K1 = cov_exp_quad(rows_idx, alpha_gp, rho_gp1) +
diag_matrix(rep_vector(delta0, N_rows));

23 K2 = cov_exp_quad(columns_index, alpha_gp, rho_gp2) +
diag_matrix(rep_vector(delta0, N_columns));

24

25 L_K1 = cholesky_decompose(K1);

26 L_K2 = cholesky_decompose(K2);

27

28 GP = kron_mvprod(L_K2, L_K1, z1);

29

30 return(GP);

31 }

32}

33

34data {

35 int<lower=1> n; // number of rows

36 int<lower=1> m; // number of columns

37 int<lower=1> N; // number of entries observed

38 int coordinates[N,2]; // coordinate of entries observed

39 int y[N]; // data on entries observed

40

41 //splines

42 int num_basis_rows; // number of B-Splines basis functions rows

43 int num_basis_columns; // number of B-Splines basis functions columns

44 matrix[num_basis_rows, n] BASIS_ROWS; // B-splines basis functions on the
rows

45 matrix[num_basis_columns, m] BASIS_COLUMNS; // B-splines basis functions on
the columns

46

47 // GP

48 real IDX_BASIS_ROWS[num_basis_rows]; // index of the B-splines basis
functions rows

49 real IDX_BASIS_COLUMNS[num_basis_columns]; // index of the B-splines basis
functions columns

50}

51

52transformed data {

53 real delta = 1e-9;

54}

55

56parameters {

57 real<lower=0> rho_1; // length scale rows

58 real<lower=0> rho_2; // length scale columns

59 real<lower=0> alpha_gp; // output variance

60 matrix[num_basis_rows,num_basis_columns] eta; // GP variables

61 real<lower=0> nu_unscaled; // overdispersion parameter

62}

63

64transformed parameters {

65 real<lower=0> nu = (1/nu_unscaled)^2;

66 real<lower=0> theta = (1 / nu);

67 matrix[num_basis_rows,num_basis_columns] beta = gp(num_basis_rows,

68 num_basis_columns,

69 IDX_BASIS_ROWS,

70 IDX_BASIS_COLUMNS,

71 delta,

72 alpha_gp,

73 rho_1, rho_2,

74 eta);

75

76 matrix[n,m] f = exp( (BASIS_ROWS’) * beta * BASIS_COLUMNS );

77 matrix[n,m] alpha = f / nu;

78}

79

80model {

81 nu_unscaled ~ normal(0,1);

82

83 rho_1 ~ inv_gamma(5, 5);

84 rho_2 ~ inv_gamma(5, 5);

85 alpha_gp ~ cauchy(0,1);

86

87 for(i in 1:num_basis_rows){

88 for(j in 1:num_basis_columns){

89 eta[i,j] ~ std_normal();

90 }

91 }

92

93 for(k in 1:N){

94 y[k] ~ neg_binomial(alpha[coordinates[k,1],coordinates[k,2]], theta);

95 }

96

97}

Listing 1: Stan model to fit regularised B-splines projected Gaussian Process
on the mean surface of 2D count data

###  S12.2 Stan model to fit regularised B-splines projected Gaussian Process
on the mean surface of 2D continuous data

⬇

1functions {

2 matrix kron_mvprod(matrix A, matrix B, matrix V)

3 {

4 return transpose(A*transpose(B*V));

5 }

6

7 matrix gp(int N_rows, int N_columns, real[] rows_idx, real[] columns_index,

8 real delta0,

9 real alpha_gp,

10 real rho_gp1, real rho_gp2,

11 matrix z1)

12 {

13

14 matrix[N_rows,N_columns] GP;

15

16 matrix[N_rows, N_rows] K1;

17 matrix[N_rows, N_rows] L_K1;

18

19 matrix[N_columns, N_columns] K2;

20 matrix[N_columns, N_columns] L_K2;

21

22 K1 = cov_exp_quad(rows_idx, alpha_gp, rho_gp1) +
diag_matrix(rep_vector(delta0, N_rows));

23 K2 = cov_exp_quad(columns_index, alpha_gp, rho_gp2) +
diag_matrix(rep_vector(delta0, N_columns));

24

25 L_K1 = cholesky_decompose(K1);

26 L_K2 = cholesky_decompose(K2);

27

28 GP = kron_mvprod(L_K2, L_K1, z1);

29

30 return(GP);

31 }

32}

33

34data {

35 int<lower=1> n; // number of rows

36 int<lower=1> m; // number of columns

37 int<lower=1> N; // number of entries observed

38 int coordinates[N,2]; // coordinate of entries observed

39 vector[N] y; // data on entries observed

40

41 //splines

42 int num_basis_rows; // number of B-Splines basis functions rows

43 int num_basis_columns; // number of B-Splines basis functions columns

44 matrix[num_basis_rows, n] BASIS_ROWS; // B-splines basis functions on the
rows

45 matrix[num_basis_columns, m] BASIS_COLUMNS; // B-splines basis functions on
the columns

46

47 // GP

48 real IDX_BASIS_ROWS[num_basis_rows]; // index of the B-splines basis
functions rows

49 real IDX_BASIS_COLUMNS[num_basis_columns]; // index of the B-splines basis
functions columns

50}

51

52transformed data {

53 real delta = 1e-9;

54}

55

56parameters {

57 real<lower=0> rho_1; // length scale rows

58 real<lower=0> rho_2; // length scale columns

59 real<lower=0> alpha_gp; // output variance

60 matrix[num_basis_rows,num_basis_columns] eta; // GP variables

61 real<lower=0> sigma; // observational noise

62}

63

64transformed parameters {

65 matrix[num_basis_rows,num_basis_columns] beta = gp(num_basis_rows,

66 num_basis_columns,

67 IDX_BASIS_ROWS,

68 IDX_BASIS_COLUMNS,

69 delta,

70 alpha_gp,

71 rho_1, rho_2,

72 eta);

73

74 matrix[n,m] f = (BASIS_ROWS’) * beta * BASIS_COLUMNS;

75}

76

77model {

78 rho_1 ~ inv_gamma(5, 5);

79 rho_2 ~ inv_gamma(5, 5);

80 alpha_gp ~ cauchy(0,1);

81

82 sigma ~ cauchy(0,1);

83

84 for(i in 1:num_basis_rows){

85 for(j in 1:num_basis_columns){

86 eta[i,j] ~ std_normal();

87 }

88 }

89

90 for(k in 1:N){

91 y[k] ~ normal(f[coordinates[k,1],coordinates[k,2]], sigma);

92 }

93

94}

Listing 2: Stan model to fit regularised B-splines projected Gaussian Process
on the mean surface of 2D continuous data

###  S12.3 R code to obtain B-splines

⬇

1bsplines = function(data, knots, degree)

2 {

3 # data: support

4 # knots: value of the knots

5 # spline_degree: spline degree

6

7 K = length(knots)

8 num_basis = K + degree - 1

9

10 intervals = find_intervals(knots, degree)

11

12 m = matrix(nrow = num_basis, ncol = length(data), 0)

13

14 for(k in 1:num_basis)

15 {

16 m[k,] = bspline(data, k, degree + 1, intervals)

17 }

18

19 m[num_basis,length(data)] = 1

20

21 return(m)

22}

23

24bspline = function(x, k, order, intervals)

25 {

26

27 if(order == 1){

28 return(x >= intervals[k] & x < intervals[k+1])

29 }

30

31 w1 = 0; w2 = 0

32

33 if(intervals[k] != intervals[k+order-1])

34 w1 = (x - intervals[k]) / (intervals[k+order-1] - intervals[k])

35 if(intervals[k+1] != intervals[k+order])

36 w2 = 1 - (x - intervals[k+1]) / (intervals[k+order] - intervals[k+1])

37

38 spline = w1 * bspline(x, k, order - 1, intervals) +

39 w2 * bspline(x, k+1, order - 1, intervals)

40

41 return(spline)

42}

43

44find_intervals = function(knots, degree)

45 {

46

47 K = length(knots)

48

49 intervals = vector(mode = ’double’, length = 2*degree + K)

50

51 # support of knots

52 intervals[(degree+1):(degree+K)] = knots

53

54 # extreme

55 intervals[1:degree] = min(knots)

56 intervals[(degree+K+1):(2*degree+K)] = max(knots)

57

58 return(intervals)

59}

Listing 3: R code to obtain B-splines

##  S13 State level summary

Figure S27: State summary for California. (A) Weekly COVID-19 attributable
deaths reported by the CDC stacked by age groups (color). (B) Posterior mean
of the weekly COVID-19 attributable deaths stacked by age, estimated with a
regularised B-splines projected GP prior. Three weeks in the observation
period are selected and designated by a shape. The black line shows the all-
age weekly COVID-19 attributable deaths reported by JHU, that were used to
adjust the CDC data for reporting delays (3.4). (C) Posterior median (black
line) and 95% credible interval (ribbon) of the age-specific contribution to
COVID-19 attributable deaths in three weeks. Figure S28: State summary for
Florida. (A) Weekly COVID-19 attributable deaths reported by the CDC stacked
by age groups (color). (B) Posterior mean of the weekly COVID-19 attributable
deaths stacked by age, estimated with a regularised B-splines projected GP
prior. Three weeks in the observation period are selected and designated by a
shape. The black line shows the all-age weekly COVID-19 attributable deaths
reported by JHU, that were used to adjust the CDC data for reporting delays
(3.4). (C) Posterior median (black line) and 95% credible interval (ribbon) of
the age-specific contribution to COVID-19 attributable deaths in three weeks.
Figure S29: State summary for Georgia. (A) Weekly COVID-19 attributable deaths
reported by the CDC stacked by age groups (color). (B) Posterior mean of the
weekly COVID-19 attributable deaths stacked by age, estimated with a
regularised B-splines projected GP prior. Three weeks in the observation
period are selected and designated by a shape. The black line shows the all-
age weekly COVID-19 attributable deaths reported by JHU, that were used to
adjust the CDC data for reporting delays (3.4). (C) Posterior median (black
line) and 95% credible interval (ribbon) of the age-specific contribution to
COVID-19 attributable deaths in three weeks. Figure S30: State summary for
Illinois. (A) Weekly COVID-19 attributable deaths reported by the CDC stacked
by age groups (color). (B) Posterior mean of the weekly COVID-19 attributable
deaths stacked by age, estimated with a regularised B-splines projected GP
prior. Three weeks in the observation period are selected and designated by a
shape. The black line shows the all-age weekly COVID-19 attributable deaths
reported by JHU, that were used to adjust the CDC data for reporting delays
(3.4). (C) Posterior median (black line) and 95% credible interval (ribbon) of
the age-specific contribution to COVID-19 attributable deaths in three weeks.
Figure S31: State summary for Michigan. (A) Weekly COVID-19 attributable
deaths reported by the CDC stacked by age groups (color). (B) Posterior mean
of the weekly COVID-19 attributable deaths stacked by age, estimated with a
regularised B-splines projected GP prior. Three weeks in the observation
period are selected and designated by a shape. The black line shows the all-
age weekly COVID-19 attributable deaths reported by JHU, that were used to
adjust the CDC data for reporting delays (3.4). (C) Posterior median (black
line) and 95% credible interval (ribbon) of the age-specific contribution to
COVID-19 attributable deaths in three weeks. Figure S32: State summary for