Bayesian Dynamic Generalized Additive Model for Mortality during COVID-19 Pandemic

The COVID-19 pandemic has had profound impacts across the globe. Researches focus on various aspects such as health disparities linked to racial and socio-economic factors and the healthcare system’s adaptation in terms of testing, contact tracing, and vaccine rollouts (Kretzschmar et al.,, 2020; Peretti-Watel et al.,, 2020; Abedi et al.,, 2021). A key area of investigation is excess mortality, which provides an overarching view of the pandemic’s impact on human health. This includes factors such as government lockdown measures and disruptions to non-COVID healthcare services (Wang et al.,, 2022; Msemburi et al.,, 2023). While excess mortality offers a general perspective, examining the consequences of specific causes of death as a result of these factors during the pandemic is crucial for developing more targeted future mitigation strategies. For example, the pandemic has contributed to increases in deaths from chronic conditions as observed by Shiels et al., (2022). There also has been a rise in accidental deaths, homicides or suicides (Pell et al.,, 2020; Mitchell and Li,, 2021; Dmetrichuk et al.,, 2022). Contributing to this research field, we analyze Italian monthly death counts recorded in 21 Italian regions from 2015 to 2020, categorized according to the International Classification of Diseases, 10th Revision (ICD-10) (WHO,, 1992), to understand the effects on specific human mortality during the pandemic better.

We apply the generalized additive model (GAM) (Hastie and Tibshirani,, 1987; Hastie,, 2017) to study non-linear relationships between response, cause-specific mortality rate, and continuous covariates, including lockdown intensity level and age, in terms of smooth functions. GAMs have demonstrated considerable potential to study COVID-19 mortality (Koum Besson et al.,, 2021; Zhu et al.,, 2021; Basellini and Camarda,, 2022). In the more general context of spatial-temporal data modeling, GAMs are well-suited for analyzing changes over time and across different geographic regions, crucial in understanding the differentiating consequences. One way of employing GAMs to model the change of spatial pattern over time is via a tensor product smoother (Wood,, 2017; Feng,, 2022). Alternatively, Clark and Wells, (2023) introduce dynamic GAMs where extra dynamic spatial random effects are incorporated into the mean of response variable, offering a solution to forecasting discrete time series while estimating relevant nonlinear predictor associations that conventional generalized linear models (GLM) plus spatial temporal random effects are not able to account for (Nazia et al.,, 2022; Yin et al.,, 2023). In our study, we follow Clark and Wells, (2023)’s approach and assume that the mortality rate of specific cause in a region for a certain month combines a GAM regression component and random effects. The random effects are further assumed to be correlated as opposed to being independent, facilitating statistical inference on dependence structure between geographical regions and various causes of death (Rao and Molina,, 2015; Smith et al.,, 2015). We adopt a Bayesian approach and assign suitable priors on model parameters to better infer correlation structure from the data.

The model framework is a special case of non-Gaussian state space models where the observation equation is Poisson distribution. Inference with state space models includes filtering to estimate the current state given past observations and smoothing to estimate past states given all observations up to the current time (Petris et al.,, 2009; Durbin and Koopman,, 2012). Filtering and smoothing is exact in linear Gaussian state space models where Kalman Filter can be applied. When the model is non-Gaussian, particle filter approximates the states using a set of weighted particles (Doucet et al.,, 2001). The posterior samples in Clark and Wells, (2023) are drawn either in the Gibbs sampling software JAGS (Plummer et al.,, 2003) or with the Hamiltonian Monte Carlo in Stan (Carpenter et al.,, 2017). However, the approximation for the target deteriorates as the dimension increases (Beskos et al.,, 2014; Van Leeuwen et al.,, 2019). Due to the high dimensionality inherent in the data and model, these sampling algorithms are impractical. To address the issue, we use a variational inference algorithm for fast approximation (Jordan et al.,, 1999; Blei et al.,, 2017); more specifically, the joint variational density of fixed effect coefficients and random effects takes the Gaussian form. Various methods exist for optimizing the mean and variance of this approximating Gaussian distribution. For instance, Titsias and Lázaro-Gredilla, (2014) propose to parameterize the density in terms of its mean and a lower triangular scale matrix whereas Tan and Nott, (2018) incorporate sparsity in the Cholesky factors of the precision matrix. Both developed stochastic gradient methods for optimization. Our approach instead involves Newton’s method and fixed point iteration as in Arridge et al., (2018) to achieve faster convergence in fewer iterations. This aligns with recent advancements with recent development in Gaussian variational approximations for high-dimensional state space models (Quiroz et al.,, 2023), but distincts in that we do not define structure of the proposed variational approximation thanks to our chosen optimization technique.

The rest of the paper is organized as follows. In Section 2, we formulate the model and specify the priors imposed on parameters. In Section 3, we demonstrate how to derive the ELBO in this model setup and present the variational algorithm for posterior inferences. We then apply the model and the algorithm to the Italian monthly mortality data in Section 4. Finally, Section 5 gives some concluding remarks and points to future work.

Let $Y_{n,t}$ be the mortality count of instance $n$ at time $t$, $n=1,\dots,N,t=1,\dots,T$. For each $n$, covariates such as ordinal $a_{n}$ for age categorical $g_{n}$ for gender, $l_{n}$ for $L$ geographical locations and $k_{n}$ for $K$ different causes of death are available. Additionally, we are interested in the lockdown effect on outcome variable $Y_{n,t}$, therefore we include a stringency index that quantifies the intensity of government restriction policies for each geographical location over time and we denote it by $r_{n,t}$, which depends on $n$ through $l_{n}$. We assume that $Y_{n,t}$ is Poisson random variable with rate equal to

where $\epsilon_{n,t}$ is the offset, $\mathbf{x}^{\bm{\beta}}_{n}\bm{\beta}$ are parametric terms which include gender effect. $f^{r}(r_{n,t})$ and $f^{a}(a_{n})$ are natural cubic splines for government intervention effect and age effect respectively and they take the form

where $J$ is the number of knots, $u^{r}_{j},u^{a}_{j}$ are natural cubic spline bases, $\beta^{r}_{j}$ and $\beta^{a}_{j}$ are coefficients to be estimated. The smoothness penalization terms associated with the bases are defined as

Here $S^{r}_{1},S^{r}_{1},S^{a}_{1},S^{a}_{2}$ contain known coefficients. $\bm{\beta}^{r}=\left(\beta^{r}_{1},\dots,\beta^{r}_{J}\right)^{\prime}$ and $\bm{\beta}^{a}=\left(\beta^{a}_{1},\dots,\beta^{a}_{J}\right)^{\prime}$ are vectors of spline coefficients. From a Bayesian perspective, this is equivalent to imposing the following multivariate normal priors on $\bm{\beta}_{r}$ and $\bm{\beta}_{a}$

$\lambda^{r}_{1},\lambda^{r}_{2},\lambda^{a}_{1}$ and $\lambda^{a}_{2}$ control the level of roughness.

As for the two-way interaction terms, $f^{k,r}(k,r)$ models interactions between cause of death and stringency index in a non-parametric manner. $f^{k,a}(k,a)$ accounts for interactions between causes of death and age while $f^{g,a}(g,a)$ captures interactions between gender and age. The three interactions terms are constructed in the following way. For each level $k$ of causes of death or each $g$ of gender, a unique spline is specified (Wood,, 2017). We formulate $f^{k,r}(k,r)$ as an example, $f^{k,a}(k,a)$ and $f^{g,a}(g,a)$ are modeled in a similar way. For each $k$ other than the baseline, $f^{k,r}(k,r)$ is assumed to be natural cubic spline such that

The $J(K-1)$ dimensional vector $\bm{\beta}^{k,r}=\left(\beta^{k,r}_{2,1},\dots,\beta^{k,r}_{2,J},\beta^{k,r}_{3,1}\dots,\beta^{k,r}_{K,J}\right)^{\prime}$ is jointly penalised by $S^{k,r}_{1}$ and $S^{k,r}_{2}$ with smoothing parameters $\lambda^{k,r}_{1}$ and $\lambda^{k,r}_{2}$ and the penalisation term is $\lambda^{k,r}_{1}(\bm{\beta}^{k,r})^{\prime}S^{r}_{1}\bm{\beta}^{k,r}+\lambda^{k,r}_{2}(\bm{\beta}^{k,r})^{\prime}S^{r}_{2}\bm{\beta}^{k,r},$ which, in terms of Bayesian prior, translates to

The last term $z^{*}_{n,t}$ embeds the spatial-temporal structure in the data; in fact, $z^{*}_{n,t}$ depends on $n$ through $l_{n}$ and $k_{n}$, therefore in total, we have $LK$ time series of length $T$. $z^{*}_{n,t}$ can be modeled with great flexibility. We make the following latent state assumption. Let $\mathbf{z}^{*}_{t}=(z^{*}_{1,t},\dots,z^{*}_{N,t})^{\prime}$ and $\mathbf{z}_{t}=(z_{1,t},\dots,z_{LK,t})^{\prime}$ such that

where $\mathbf{1}$ stands for a vector whose elements are all equal to 1. The dimension of $\mathbf{1}$ is determined by the number of age groups and gender. The assumption implies that residual mortality rates $z^{*}_{n,t}$ of the same causes of death in the same region are identical regardless of age and gender, which is reasonable since we have already taken into account age and gender effect in GAM component. We further assume that marginally the latent state $\mathbf{z}_{t}\sim\mathcal{N}\left(\bm{\mu},\Sigma\right)$ and $P\left(\mathbf{z}_{t}-\bm{\mu}\right)$ follows a simple autoregressive (AR) process such that

where $P$ is an upper triangular matrix such that $P^{\prime}P=\Omega=\Sigma^{-1}$, i.e. P are the Cholesky factors of the precision matrix $\Omega$. $\Phi$ contains multivariate autoregressive coefficients. For simplicity, we require $\Phi$ to be diagonal with diagonal entries $\bm{\phi}=c(\phi_{1},\dots,\phi_{LK})^{\prime}$ such that $-1<\phi_{1},\dots,\phi_{LK}<1$. The pre-multiplying of $\mathbf{z}_{t}$ by $P$ is to ensure that marginally $\mathbf{z}_{t}\sim\mathcal{N}\left(\mathbf{0},\Sigma\right),t=0,\dots,T$ where $\Sigma$ incorporates dependence structure of geographical locations as well as mortality causes. Reorganizing $\mathbf{z}_{t}$ as matrices $Z_{t}$ of size $L\times K$, we assume that

which is a matrix normal distribution equivalent to

where $\bigotimes$ stands for the Kronecker product between $\Sigma^{k}$ and $\Sigma^{l}$, the covariance matrices of the causes of death dependence and regional dependence respectively. We impose the following priors on the model parameters. Firstly, $\bm{\mu}$ is assigned a normal prior $\mathcal{N}\left(\mathbf{0},\sigma^{2}_{\mu}I_{LK}\right)$. $\Omega^{k}=\left(\Sigma^{k}\right)^{-1}$ is the precision matrix that underlines the conditional independence of mortality causes. For simplicity, we stick to the usual Wishart distribution

The same reasoning applies to the precision matrix $\Omega^{l}=\left(\Sigma^{l}\right)^{-1}$ and it follows that

We complete the model with a hierarchy of prior specification on the linear regression coefficients, on penalization parameters that control smoothness of splines and on the autoregressive coefficients, on the cause specific mean. Denoting by $\mathbf{x}_{n}$ the vector containing $\mathbf{x}^{\bm{\beta}}_{n}$ and all bases in GAMs for $n$, the hierarchical model can be summarized as

To make posterior inference, one may devise suitable Markov Chain Monte Carlo (MCMC) algorithms to obtain posterior samples. However, the algorithm may fail to deliver desirable convergent output within reasonable time when it is difficult to explore the geometry of the target distribution due to the sheer dimension of the data. Therefore we resort to variational inference approach for fast approximation. The target posterior distribution is

and the goal is to find the optimal $q^{*}\left(\bm{\beta}^{*},\bm{\lambda},\mathbf{z}_{0:T},\bm{\mu},\bm{\phi},\Omega^{k},\Omega^{l}\right)$ from a pre-specified family of distributions such that the Kullback–Leibler (KL) divergence, defined as

with expectation taken with respect to $q\left(\cdot\right)$, is minimized. This is equivalent to maximizing the evidence lower bound (ELBO) given by

We further assume that the family of candidate approximation $q(\cdot)$ can be factorized as

where $q\left(\bm{\beta}^{*},\mathbf{z}_{0:T}\right)$ is multivariate normal distribution with mean $\mathbf{m}$ and covariance matrix $M$. This is essentially the variational Gaussian approximation (VGA) that has been widely implemented in literature and its theoretical properties when applied to Poisson data are studied by Arridge et al., (2018). $M$ can be full or it can be sparse block diagonal matrix so that there is no correlation between $\bm{\beta}^{*}$ and $\mathbf{z}_{0:T}$, which means that $q\left(\bm{\beta}^{*},\mathbf{z}_{0:T}\right)$ can be further factorized as the product of two multivariate normal densities, $q\left(\bm{\beta}^{*}\right)$ and $q\left(\mathbf{z}_{0:T}\right)$. The sparse version of $M$ greatly reduces algorithm complexity. However, in our application, since the dimension of $\mathbf{z}_{0:T}$ is overwhelmingly larger than $\bm{\beta}^{*}$, the block diagonal assumption does not produce much gain; therefore we stick to the full matrix characterisation of $M$. Furthermore, we choose the following variational densities implied by mean field approximation: for each element in $\bm{\lambda}$, it is a point mass at $\lambda^{q,r}_{1},\lambda^{q,r}_{2},\lambda^{q,a}_{1},\lambda^{q,a}_{2},\lambda^{q,k,r}_{1},\lambda^{q,k,r}_{2},\lambda^{q,k,a}_{1},\lambda^{q,k,a}_{2},\lambda^{q,g,a}_{1},\lambda^{q,g,a}_{2}$. Together they are $\bm{\lambda}^{q}$. The Dirac measure choice as a variational density avoids evaluating expectation with respect to otherwise non-trivial variational densities. For $q\left(\bm{\mu}\right)$, we assume independent normal approximation densities $q\left(\bm{\mu}\right)=\mathcal{N}\left(\bm{\mu}^{q},\text{diag}\left[\left(\bm{\sigma}^{q}\right)^{2}\right]\right)$. For $q\left(\bm{\phi}\right)$, even though it is analytically possible to compute the expectation when assuming that $(\phi_{1}+1)/2\overset{i.i,d}{\sim}\text{Beta}(\alpha^{q,\phi}_{1},\beta^{q,\phi}_{1}),\dots,(\phi_{LK}+1)/2\overset{i.i,d}{\sim}\text{Beta}(\alpha^{q,\phi}_{LK},\beta^{q,\phi}_{LK})$, the optimization scheme we adapt diverges as only the means of beta variational distributions is identifiable, therefore we use Dirac measures on $q\left(\bm{\phi}\right)$ so they are $\bm{\phi}^{q}$. Finally, we set $q(\Omega^{k})=\text{Wishart}(\delta^{q,k},D^{q,k})$ and $q(\Omega^{l})=\text{Wishart}(\delta^{q,l},D^{q,l})$.