Advanced Mathematical Variants for Estimating Excess Mortality Associated with the COVID-19 and Others Pandemics

1. Introduction

The COVID-19 pandemic, caused by the SARS-CoV-2 virus, has had a profound global impact on public health, extending far beyond officially reported deaths directly attributed to the disease. While official COVID-19 death counts provide valuable information, they often underestimate the true mortality burden due to factors such as limited testing, variations in death certification protocols, and indirect effects of the pandemic [3,7]. These indirect effects include overwhelmed health systems leading to excess deaths from other causes, avoidance of healthcare by patients, and, conversely, reductions in certain causes of death (e.g., traffic accidents, air pollution-related mortality, or seasonal influenza) due to non-pharmaceutical interventions like lockdowns [6,9,10].

In this context, excess mortality—defined as the difference between observed all-cause mortality (ACM) and the expected mortality in the absence of the pandemic—serves as a more objective, comprehensive, and comparable measure of the pandemic’s overall mortality impact across countries and regions [4,5]. The World Health Organization (WHO) defines excess mortality as “the mortality above what would be expected based on the non-crisis mortality rate in the population of interest” [5]. For country c and month t, excess deaths are calculated as

$${\delta }_{c,t}=Y_{c,t}-E_{c,t},$$

where $Y_{c,t}$ denotes the observed ACM during the pandemic period, and $E_{c,t}$ is the counterfactual expected ACM forecasted from pre-pandemic historical data [1].

The current WHO methodology for estimating excess mortality, detailed in their technical report [1], relies on a combination of approaches tailored to data availability:

• For countries with reliable monthly historical ACM data (typically 2015–2019), a generalized additive model (GAM) with cyclic cubic splines for seasonality and a linear trend for annual variation is fitted under a negative binomial distribution to predict $E_{c,t}$. Uncertainty is summarized via point estimates and standard deviations derived from posterior approximations.
• For countries with only annual data or no pandemic-period ACM, a log-linear overdispersed Poisson regression uses time-varying and static covariates (e.g., official COVID-19 death rates, stringency indices, temperature, socio-demographic index) to impute $Y_{c,t}$ or directly predict excess.
• Age- and sex-specific patterns are derived from observed data where available, clustered groupings for generalization, and extrapolation to data-scarce countries.

Despite its strengths—including adjustments for registration completeness, covariate selection via out-of-sample error minimization, and hierarchical elements for global comparability—the WHO approach has several mathematical and statistical limitations:

• The GAM spline-based modeling assumes additivity and smoothness that may fail to capture complex non-linearities, structural breaks pre-pandemic, or abrupt local shocks.
• The baseline period (e.g., 2015–2019) may be too short to fully account for long-term trends or cohort effects, leading to biased expected values.
• Covariate models for data-scarce countries rely on strong parametric assumptions (log-linearity) and regional median imputation, potentially introducing bias in heterogeneous settings.
• Extrapolation of age-sex patterns via discrete clusters may overlook continuous gradients in demographic or epidemiological similarity.
• Uncertainty quantification, while present, is often approximate and may not fully propagate hierarchical dependencies or model misspecification.

These limitations motivate the exploration of alternative mathematical frameworks that enhance flexibility, robustness, and interpretability while preserving the core principles of baseline forecasting, covariate adjustment, and uncertainty modeling. In this paper, we propose and develop five advanced statistical variants to improve excess mortality estimation:

1.

Structural Time Series (STS) models with Kalman filtering.

2.

Extended SARIMAX models with intervention components.

3.

Fully Bayesian hierarchical models with penalized splines and shrinkage priors.

4.

Compartmental ordinary differential equation (ODE) models inspired by extended SEIR frameworks.

5.

Functional Principal Component Analysis (FPCA) combined with functional regression.

The remainder of the paper is organized as follows. Section 2 presents preliminaries, including key mathematical elements of the WHO baseline method. Section 3–7 develop each proposed variant in detail, with formal definitions, equations, and derivations where appropriate. Section 8 discusses the importance, strengths, and potential synergies of the proposed approaches. Finally, Section 9concludes and outlines future directions.

2. Materials and Methods

The World Health Organization (WHO) framework for estimating excess mortality associated with the COVID-19 pandemic relies on a combination of statistical models tailored to data availability across countries [1,8]. Here we formalize the key mathematical components, focusing on the core estimation of expected deaths $E_{c,t}$ and excess deaths ${\delta }_{c,t}$.

2.1. Notation and Core Definitions

Let $c=1,\dots ,194$ index the WHO Member States, and $t=1,\dots ,24$ index the months from January 2020 to December 2021. For each country c and month t, define:

• $Y_{c,t}$: observed all-cause mortality (ACM) count during the pandemic period.
• $E_{c,t}$: expected (counterfactual) ACM count in the absence of the pandemic, forecasted from pre-pandemic historical data.
• ${\delta }_{c,t}=Y_{c,t}-E_{c,t}$: excess deaths (positive for excess mortality, negative for mortality deficit).

Negative excess (${\delta }_{c,t}<0$) may arise from averted deaths (e.g., reduced traffic accidents or influenza) or model artifacts [1].

Observed ACM counts are often incomplete due to registration delays or underreporting. The WHO adjusts for completeness using demographic methods:

$$Y_{c,t}={\displaystyle \frac{Y_{c,t}^{\prime }}{{\mathrm{completeness}}_{c,t}}},$$

where $Y_{c,t}^{\prime }$ is the unadjusted reported count and ${\mathrm{completeness}}_{c,t}\in (0,1]$ is the estimated fraction of true deaths registered.

2.2. Expected Mortality for Countries with Monthly Historical Data

For countries with reliable monthly ACM data (typically 2015–2019), the WHO employs a generalized additive model (GAM) under a negative binomial distribution to model pre-pandemic mortality and forecast $E_{c,t}$ [1].

Let $Y_{c,t}$ for $t=1,\dots ,M_c$ denote historical monthly ACM counts, where $M_c$ is the number of available pre-pandemic months. The sampling model is

$$Y_{c,t}\mid {\mu }_{c,t}^E,{\varphi }_c^E\sim \mathrm{NegBin}({\mu }_{c,t}^E,{\varphi }_c^E),$$

parametrized such that

$$\mathbb{E}\left[Y_{c,t}\right]={\mu }_{c,t}^E,\mathrm{Var}\left(Y_{c,t}\right)={\mu }_{c,t}^E\left(1+{\displaystyle \frac{{\mu }_{c,t}^E}{{\varphi }_c^E}}\right),$$

with the Poisson model recovered as ${\varphi }_c^E\to \infty$.

Define $v\left[t\right]$ as the year index (e.g., 1 to 5 for 2015–2019) and $m\left[t\right]$ as the month index (1 to 12), so $t=12\left(v\right[t]-1)+m\left[t\right]$. The mean is modeled on the log scale:

$${\eta }_{c,t}=log\left({\mu }_{c,t}\right)=f_c^y\left(v\left[t\right]\right)+f_c^m\left(m\left[t\right]\right),$$

where

• $f_c^y(·)$ is a linear trend in year to capture annual variation,
• $f_c^m(·)$ is a cyclic cubic spline on month to model within-year seasonality [17].

The GAM is fitted separately per country using the gam function in the mgcv package with restricted maximum likelihood (REML) for smoothing parameter selection [2]. Predictions for t corresponding to 2020–2021 yield point estimates ${\widehat{\mu }}_{c,t}^E$ and standard deviations, interpreted as approximate posterior summaries (Section 6.10 of [2]).

2.3. Expected Mortality for Countries with Annual Data or No Data

For countries with only annual historical data, a similar but annual-scale model is applied (details in Section 3.2 of [1]). For countries lacking pandemic-period ACM entirely, a log-linear overdispersed Poisson regression is used:

$$log\left({\mu }_{c,t}\right)={\beta }_0+\sum _k{\beta }_k{X}_{k,t}+{ϵ}_{c,t},$$

where $X_{k,t}$ includes time-varying covariates (e.g., official COVID-19 death rate, stringency index, temperature) and static ones (e.g., SDI, HDI, population density). Missing covariate values are imputed via WHO regional medians. The final covariate set is selected to minimize out-of-sample error [1].

2.4. Uncertainty and Age-Sex Patterns

Uncertainty in $E_{c,t}$ is summarized via standard deviations from the GAM or regression fits. For age- and sex-specific excess, observed patterns are used where available; otherwise, countries are grouped into clusters based on similarity of age-sex impacts, with extrapolation via cluster membership or covariate-based regression (Section 6 of [1]).

These components form the baseline against which we propose improvements in subsequent sections, addressing limitations in smoothness assumptions, linearity, baseline sensitivity, and extrapolation robustness.

3. Results

3.1. Structural Time Series Models with Kalman Filtering

One promising direction to overcome the limitations of spline-based GAMs is the use of structural time series (STS) models, also known as unobserved components models. These models decompose the time series into interpretable stochastic components (trend, seasonality, cycle, regressors, and irregular) and allow exact likelihood evaluation and forecasting via the Kalman filter [11,12].

3.1.1. Model Formulation

For each country c, consider the monthly all-cause mortality counts $Y_{c,t}$ for $t=1,\dots ,T_c$ (pre-pandemic period, typically 2015–2019 extended if available). We model the latent state process as

$$Y_{c,t}={\mu }_{c,t}+{\gamma }_{c,t}+{\mathbf{x}}_{c,t}^{\top }{\beta }_c+{ϵ}_{c,t},$$

where

• ${\mu }_{c,t}$ is the stochastic trend component (local level + slope),
• ${\gamma }_{c,t}$ is the seasonal component (monthly periodicity),
• ${\mathbf{x}}_{c,t}$ is a vector of covariates (time-varying or static),
• ${ϵ}_{c,t}\sim \mathrm{N}(0,{\sigma }_{ϵ,c}^2)$ is the irregular component (or Negative Binomial errors for counts).

To handle count data directly, we can adopt a generalized STS framework with observation equation

$$Y_{c,t}\mid {\eta }_{c,t}\sim \mathrm{NegBin}(exp\left({\eta }_{c,t}\right),{\varphi }_c),$$

where the linear predictor is

$${\eta }_{c,t}=log\left({\mu }_{c,t}\right)={\mathrm{trend}}_{c,t}+{\mathrm{seasonal}}_{c,t}+{\mathbf{x}}_{c,t}^{\top }{\beta }_c+u_{c,t},$$

and $u_{c,t}$ is a residual term (possibly AR(1) to capture short-term autocorrelation).

The trend component is typically specified as a local linear trend:

$$\begin{array}{cc}\hfill {\mathrm{level}}_{c,t}& ={\mathrm{level}}_{c,t-1}+{\mathrm{slope}}_{c,t-1}+{\eta }_{\mathrm{level},c,t},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\mathrm{slope}}_{c,t}& ={\mathrm{slope}}_{c,t-1}+{\eta }_{\mathrm{slope},c,t},\hfill \end{array}$$

(2)

with ${\eta }_{\mathrm{level},c,t}\sim \mathrm{N}(0,{\sigma }_{\mathrm{level},c}^2)$ and ${\eta }_{\mathrm{slope},c,t}\sim \mathrm{N}(0,{\sigma }_{\mathrm{slope},c}^2)$. For smooth trends, one can set ${\sigma }_{\mathrm{slope},c}^2=0$ (local level model) or use integrated random walk priors.

The seasonal component for monthly data (period $s=12$) can be modeled using dummy variables with sum-to-zero constraint or trigonometric (harmonic) form:

$${\gamma }_{c,t}=\sum _{j=1}^6\left({\alpha }_{j,c}cos\left({\displaystyle \frac{2\pi jt}{12}}\right)+{\beta }_{j,c}sin\left({\displaystyle \frac{2\pi jt}{12}}\right)\right),$$

where the coefficients ${\alpha }_{j,c},{\beta }_{j,c}$ evolve as random walks or are fixed (for deterministic seasonality).

3.1.2. Kalman Filter and Smoothing

The STS model can be cast in state-space form:

$$\begin{array}{cc}\hfill {\alpha }_t& =\mathbf{T}{\alpha }_{t-1}+\mathbf{R}{\eta }_t,{\eta }_t\sim \mathrm{N}(\mathbf{0},\mathbf{Q}),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\eta }_{c,t}& =\mathbf{Z}{\alpha }_t+{\mathbf{x}}_{c,t}^{\top }{\beta }_c+u_{c,t},\hfill \end{array}$$

(4)

where ${\alpha }_t$ is the state vector stacking level, slope, seasonal harmonics, etc.

The Kalman filter provides recursive computation of the filtered state ${\widehat{\alpha }}_{t|t}$ and one-step-ahead prediction error variance, enabling exact maximum likelihood estimation of hyperparameters (${\sigma }^2$ terms, ${\varphi }_c$, ${\beta }_c$) via prediction error decomposition. Smoothing yields ${\widehat{\alpha }}_{t|T}$ for in-sample reconstruction.

Forecasting $E_{c,t}$ for t in 2020–2021 is obtained by projecting the state forward without observation updates, with predictive intervals derived from the forecast error variance.

3.1.3. Advantages for Excess Mortality Estimation

Compared to the WHO GAM approach:

• The stochastic trend and slope allow flexible adaptation to gradual changes in baseline mortality (e.g., aging population, improvements in healthcare) without assuming strict linearity.
• Intervention analysis is natural: additive or multiplicative pulse/step dummies can model known non-COVID shocks (e.g., heatwaves, policy changes).
• Covariates enter linearly but can be made time-varying (e.g., stringency index) with random-walk coefficients if desired.
• The Kalman filter provides coherent uncertainty propagation, including multi-step forecast intervals that account for parameter uncertainty (via Monte Carlo or bootstrap on hyperparameters).
• For countries with sparse data, hierarchical extensions allow partial pooling of variance components (e.g., ${\sigma }_{\mathrm{level},c}^2\sim \mathrm{InvGamma}(\nu ,\tau )$ with shared hyperparameters across WHO regions).

3.1.4. Implementation Considerations

Practical implementation can use packages such as statsmodels.tsa.statespace.structural (Python), KFAS (R), or Bayesian frameworks like bsts (R) or PyMC/Stan for hierarchical versions. For Negative Binomial observation models, approximate inference via integrated Laplace or particle filtering may be required in non-Gaussian cases.

In summary, STS models offer a transparent, interpretable, and statistically rigorous alternative to GAMs, with stronger theoretical grounding in state-space theory and superior handling of non-stationary trends and forecasting uncertainty.

3.2. Extended SARIMAX Models with Intervention Analysis

An alternative classical approach to modeling monthly all-cause mortality time series is the autoregressive integrated moving average model with seasonal components and exogenous regressors (SARIMAX), extended with deterministic intervention terms to account for known shocks during the pandemic period. This framework, rooted in the seminal work of Box and Tiao [13], provides a powerful, parsimonious, and highly interpretable tool for non-stationary time series with strong seasonality, while allowing explicit modeling of abrupt changes (e.g., COVID-19 waves or policy interventions).

3.2.1. Model Specification

For each country c, let $Y_{c,t}$ denote the observed ACM count at month t. Since $Y_{c,t}$ are positive integers with overdispersion, we consider two common strategies:

1.

Log-transformation with correction: model $Z_{c,t}=log(Y_{c,t}+1)$ (or Box-Cox transformation $\lambda$-optimal) as Gaussian SARIMAX, then back-transform forecasts with bias correction.

2.

Generalized linear SARIMAX: directly model $Y_{c,t}\sim \mathrm{NegBin}({\mu }_{c,t},{\varphi }_c)$ or Poisson-gamma, with $log\left({\mu }_{c,t}\right)$ following an ARIMA structure on covariates and interventions.

We focus here on the log-linear version for analytical tractability, with the understanding that Negative Binomial can be approximated via quasi-likelihood or state-space embedding (as in Section 3).

The general SARIMAX model with intervention is

$$\varphi \left(B\right)\Phi \left(B^{12}\right){(1-B)}^d{(1-B^{12})}^D{Z}_{c,t}=\theta \left(B\right)\Theta \left(B^{12}\right){ϵ}_{c,t}+\sum _{k=1}^K{\beta }_k{X}_{k,t}+\omega \left(B\right)I_{c,t}+\delta \left(B\right){\displaystyle \frac{\omega \left(B\right)}{\nu \left(B\right)}}{P}_{c,t},$$

(5)

where

• B is the backshift operator ($B{Z}_t=Z_{t-1}$),
• $\varphi \left(B\right)=1-{\varphi }_{1}B-\dots -{\varphi }_p{B}^p$ (regular AR),
• $\Phi \left(B^{12}\right)=1-{\Phi }_1{B}^{12}-\dots -{\Phi }_P{B}^{12P}$ (seasonal AR),
• d and D are regular and seasonal differencing orders,
• $\theta \left(B\right)$ and $\Theta \left(B^{12}\right)$ are regular and seasonal MA polynomials,
• ${ϵ}_{c,t}\sim \mathrm{N}(0,{\sigma }^2)$,
• $X_{k,t}$ are exogenous covariates (COVID-19 death rate, stringency index, temperature, etc.),
• $I_{c,t}$ and $P_{c,t}$ are intervention variables (pulse and step functions, respectively),
• $\omega \left(B\right)/\nu \left(B\right)$ is the transfer function describing dynamic response to interventions (e.g., first-order decay).

Intervention components are critical for excess mortality estimation, as they isolate pandemic-related shocks from baseline dynamics:

• Pulse intervention (temporary shock): $I_{c,t}=1$ if t corresponds to a known high-mortality wave (e.g., March–April 2020), 0 otherwise.
• Step intervention (permanent shift): $P_{c,t}=1$ for $t\ge t_0$ (e.g., start of vaccination rollout or sustained policy change).
• Dynamic response: common forms include ${\omega }_0$ (immediate effect), ${\omega }_0{(1-\delta B)}^{-1}$ (gradual permanent), or ${\omega }_0(1-\delta B)/(1-\nu B)$ (overshoot and decay).

3.2.2. Parameter Estimation and Diagnostic Checking

Maximum likelihood estimation of $(\varphi ,\Phi ,\theta ,\Theta ,\beta ,\omega ,\delta ,\nu ,{\sigma }^2)$ is performed via conditional or exact methods (e.g., using the Kalman filter representation of ARIMA). Model selection follows the Box-Jenkins iterative procedure:

1.

Identification: ACF/PACF of differenced series to select $p,d,q,P,D,Q$.

2.

Estimation: conditional or unconditional MLE.

3.

Diagnostics: Ljung-Box test on residuals, Jarque-Bera for normality, Chow test or CUSUM for structural breaks pre-intervention.

4.

Intervention detection: if unknown, use outlier detection algorithms to identify additive outliers (AO), innovative outliers (IO), level shifts (LS), or temporary changes (TC).

Forecasting $E_{c,t}$ for t in 2020–2021 is obtained by setting intervention terms to zero (counterfactual baseline) and projecting forward using the fitted model.

3.2.3. Uncertainty Quantification

Prediction intervals are derived from the forecast error variance:

$$\mathrm{Var}\left({\widehat{Z}}_{c,t+h|t}\right)={\sigma }^2\left(1+\sum _{j=1}^{h-1}{\psi }_j^2\right),$$

where ${\psi }_j$ are the MA(∞) coefficients of the ARIMA operator. For the original scale, delta method or Monte Carlo simulation (sampling from parameter asymptotic normality) yields approximate intervals for $E_{c,t}$.

3.2.4. Advantages and Comparison with WHO GAM Approach

The extended SARIMAX framework offers several theoretical and practical advantages over the GAM with cyclic cubic splines:

• Explicit handling of non-stationarity via differencing and stochastic trends, avoiding the need to assume smoothness over the entire period.
• Direct incorporation of seasonal AR and MA terms, which can capture complex intra-annual dependencies (e.g., holiday effects, flu seasons) more parsimoniously than high-dimensional splines.
• Formal intervention modeling isolates pandemic shocks, improving causal interpretability of excess deaths ${\delta }_{c,t}$ as the sum of intervention effects plus residual deviations.
• Built-in diagnostic tools (ACF, Ljung-Box, outlier detection) ensure model adequacy and detect misspecification.
• For countries with partial or sparse data, parameters can be estimated hierarchically (e.g., pooling AR/MA coefficients across WHO regions via shrinkage priors or empirical Bayes).

Limitations include potential sensitivity to differencing order and intervention specification, and less flexibility in capturing highly non-linear covariate effects (addressable via polynomial or spline regressors in $X_{k,t}$).

In practice, SARIMAX models often outperform GAMs in short- to medium-term forecasting of seasonal mortality series when strong autocorrelation and known shocks are present [13]. They provide a robust, interpretable baseline for excess mortality estimation, complementing the spline-based approach of the WHO methodology.

3.3. Fully Bayesian Hierarchical Models with Penalized Splines and Shrinkage Priors

The WHO approach already employs a GAM fitted via REML, which can be viewed as an empirical Bayes approximation to a fully Bayesian model [2]. However, REML provides only point estimates and approximate standard errors, without full posterior inference or coherent hierarchical pooling across countries. A natural and mathematically consistent extension is to adopt a fully Bayesian hierarchical model with penalized splines and shrinkage priors on country-specific effects and hyperparameters. This framework allows exact posterior uncertainty quantification, borrowing strength across countries and regions, and more robust extrapolation to data-scarce settings.

3.3.1. Model Specification

For country c and month t, we model the observed all-cause mortality $Y_{c,t}$ (adjusted for completeness) as

$$Y_{c,t}\mid {\mu }_{c,t},{\varphi }_c\sim \mathrm{NegativeBinomial}({\mu }_{c,t},{\varphi }_c),$$

(6)

where the mean is linked via

$$log\left({\mu }_{c,t}\right)={\eta }_{c,t}=f_c^y\left(v\left[t\right]\right)+f_c^m\left(m\left[t\right]\right)+{\mathbf{x}}_{c,t}^{\top }{\beta }_c+{\eta }_c+{ϵ}_{c,t},$$

(7)

with:

• $v\left[t\right]$: year index in the historical period,
• $m\left[t\right]$: month index (1 to 12),
• $f_c^y(·)$: smooth function of annual trend,
• $f_c^m(·)$: smooth cyclic function of seasonality,
• ${\mathbf{x}}_{c,t}$: vector of covariates (time-varying and static),
• ${\beta }_c$: country-specific regression coefficients,
• ${\eta }_c$: country-level random intercept (random effect),
• ${ϵ}_{c,t}$: residual term (optional AR(1) structure).

To ensure identifiability and smoothness, we represent the smooth functions using penalized B-splines or thin-plate regression splines with second-order penalties [2]:

$$f_c^y\left(v\right)=\sum _{k=1}^{K_y}{b}_k^y\left(v\right){\alpha }_{c,k}^y,f_c^m\left(m\right)=\sum _{k=1}^{K_m}{b}_k^m\left(m\right){\alpha }_{c,k}^m,$$

(8)

where $b_k^y(·)$ and $b_k^m(·)$ are B-spline bases, and the coefficients ${\alpha }_{c,·}^y$ and ${\alpha }_{c,·}^m$ receive a multivariate normal prior with quadratic penalty:

$${\alpha }_c^y\sim \mathrm{N}(\mathbf{0},{\sigma }_y^2{\mathbf{P}}_{+}^y),{\alpha }_c^m\sim \mathrm{N}(\mathbf{0},{\sigma }_m^2{\mathbf{P}}_{+}^m),$$

(9)

where ${\mathbf{P}}_{+}^y$ and ${\mathbf{P}}_{+}^m$ are precision matrices derived from the penalty (e.g., second differences for smoothness), and ${\mathbf{P}}_{+}$ denotes the Moore-Penrose pseudoinverse to handle the null space of the penalty (intercept and linear terms).

The country-specific random effects and coefficients follow hierarchical priors:

$$\begin{array}{cc}\hfill {\eta }_c& \sim \mathrm{N}(0,{\sigma }_{\eta }^2),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\beta }_c& \sim \mathrm{N}({\mu }_{\beta },{\Sigma }_{\beta }),\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\sigma }_{\eta }^2& \sim \mathrm{InvGamma}(a_{\eta },b_{\eta })(\mathrm{or}\mathrm{half}-\mathrm{Cauchy}\mathrm{for}\mathrm{weakly}\mathrm{informative}),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\mu }_{\beta }& \sim \mathrm{N}(\mathbf{0},{10}^2\mathbf{I}),{\Sigma }_{\beta }\sim \mathrm{InvWishart}(\nu ,\mathbf{S}),\hfill \end{array}$$

(13)

with hyperparameters shared across countries or grouped by WHO region $r\left(c\right)$:

$${\sigma }_{\eta ,r}^2\sim \mathrm{InvGamma}(a_r,b_r),\mathrm{etc}.$$

(14)

The overdispersion parameter receives a weakly informative prior:

$${\varphi }_c\sim \mathrm{Gamma}(0.01,0.01)\mathrm{or}{\varphi }_c^{-1}\sim \mathrm{Gamma}(0.01,0.01).$$

(15)

3.3.2. Inference and Posterior Computation

The joint posterior is

$$p(\theta ,\alpha ,{\sigma }^2\mid \mathbf{Y})\propto \prod _{c,t}p(Y_{c,t}\mid {\mu }_{c,t},{\varphi }_c)p(\alpha \mid {\sigma }^2)p\left({\sigma }^2\right)p(\beta ,\eta \mid \mathrm{hierarchical}\mathrm{priors}),$$

(16)

where $\theta$ collects all fixed and random effects.

Inference can be performed via:

• Markov chain Monte Carlo (MCMC) using No-U-Turn Sampler (NUTS) in Stan or PyMC,
• Integrated Nested Laplace Approximation (INLA) for fast approximate inference when using Gaussian random fields [14],
• Variational Bayes for scalability in large hierarchies.

Posterior predictive checks generate replicated datasets ${\mathbf{Y}}^{\mathrm{rep}}$ to assess model fit. Excess mortality is obtained as

$${\delta }_{c,t}=Y_{c,t}-\mathbb{E}[{\mu }_{c,t}^E\mid {\mathbf{Y}}_{\mathrm{pre}}],$$

(17)

with full posterior credible intervals by sampling from the posterior predictive distribution of ${\mu }_{c,t}^E$ (forecasts setting pandemic covariates to zero or counterfactual values).

3.3.3. Theoretical Properties and Consistency

Under standard regularity conditions (sufficient historical data per country, bounded covariate effects, identifiability of smooth functions), the posterior is consistent and contracts at rate $n^{-1/2}$ for parametric components and near-parametric rates for smooth functions [16]. The hierarchical shrinkage induces:

• Borrowing strength: countries with sparse data borrow information from similar countries/regions via shared hyperparameters, reducing variance.
• Shrinkage to group mean: extreme country-specific effects (e.g., outliers due to data errors) are pulled toward regional means.
• Coherent uncertainty: full posterior accounts for all sources of variability, including hyperparameter uncertainty (unlike REML approximations).

For countries with no pandemic data, predictions are obtained by integrating over the posterior of ${\beta }_{c^{*}}$ and ${\eta }_{c^{*}}$ conditional on covariates and regional hyperparameters.

3.3.4. Advantages over WHO REML-GAM and Comparison

Compared to the current WHO method:

• Full posterior inference replaces approximate standard errors with credible intervals that incorporate hierarchical dependence.
• Partial pooling via random effects and shared priors improves extrapolation to the 83 countries with no data (as of 2023), reducing bias from regional median imputation.
• Priors allow regularization of covariate effects (e.g., horseshoe priors on $\beta$ for sparsity) and incorporation of prior knowledge (e.g., negative effect of stringency on non-COVID mortality).
• Robustness to baseline choice: longer historical periods (2000–2019) can be included with time-varying smoothness penalties.

Limitations include higher computational cost (mitigated by INLA) and sensitivity to prior specification (addressed via sensitivity analysis and weakly informative defaults).

This fully Bayesian hierarchical extension provides the most mathematically coherent upgrade to the WHO spline-based model, preserving interpretability while significantly enhancing robustness, uncertainty quantification, and global comparability.

3.4. Compartmental Ordinary Differential Equation Models for Excess Mortality

To capture the nonlinear dynamics of disease transmission and its direct impact on mortality, we propose an extension of compartmental epidemiological models, specifically an adapted SIRD (Susceptible-Infected-Recovered-Dead) framework with natural mortality and disease-induced death rate. Unlike purely statistical time-series approaches, ODE models explicitly link infection dynamics to excess deaths via the additional mortality term, enabling mechanistic counterfactual simulations.

3.4.1. Model Formulation

Let $t\ge 0$ be continuous time (months). Define the compartments:

• $S\left(t\right)$: susceptible individuals,
• $I\left(t\right)$: infected (and infectious) individuals,
• $R\left(t\right)$: recovered (immune) individuals,
• $D\left(t\right)$: cumulative deaths (both natural and disease-induced).

Let $N\left(t\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)$ be the living population (decreasing due to deaths). The system of ordinary differential equations is:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& =-\beta \left(t\right){\displaystyle \frac{SI}{N}}-\mu S,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dI}{dt}}& =\beta \left(t\right){\displaystyle \frac{SI}{N}}-(\gamma +\delta \left(t\right)+\mu )I,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dR}{dt}}& =\gamma I-\mu R,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dD}{dt}}& =\delta \left(t\right)I+\mu (S+I+R)=\delta \left(t\right)I+\mu N\left(t\right),\hfill \end{array}$$

(21)

where:

• $\beta \left(t\right)$: time-varying transmission rate (incorporating stringency, mobility, variants),
• $\gamma$: recovery rate (1/infectious period),
• $\delta \left(t\right)$: disease-induced death rate per infected individual (COVID-specific),
• $\mu$: natural (background) mortality rate per capita (small, e.g., 0.01/year ≈ 0.00083/month).

Initial conditions: $S\left(0\right)\approx N_0$, $I\left(0\right)>0$ small, $R\left(0\right)\approx 0$, $D\left(0\right)=0$, where $N_0$ is pre-pandemic population.

The expected (counterfactual) deaths without pandemic are obtained by setting $\delta \left(t\right)\equiv 0$ and $\beta \left(t\right)$ to pre-pandemic baseline (or zero if no transmission):

$${\displaystyle \frac{d{D}^E}{dt}}=\mu N^E\left(t\right),N^E\left(t\right)=S^E\left(t\right)+I^E\left(t\right)+R^E\left(t\right)\approx N_0{e}^{-\mu t}.$$

(22)

Thus, $D^E\left(t\right)\approx \mu N_{0}t$ (linear in time for small $\mu$).

Excess deaths over interval $[0,T]$ (e.g., 24 months) are

$$\delta \left(T\right)=D\left(T\right)-D^E\left(T\right)\approx {\int }_0^T\delta \left(t\right)I\left(t\right)dt+\mathrm{indirect}\mathrm{effects}\mathrm{via}N\left(t\right)<N^E\left(t\right).$$

(23)

The primary contribution is the integral of disease-induced deaths $\int \delta \left(t\right)I\left(t\right)dt$.

Covariates enter $\beta \left(t\right)$ and $\delta \left(t\right)$ via regression (e.g., $\beta \left(t\right)={\beta }_{0}exp(-\kappa ·\mathrm{stringency}\left(t\right))$) or as forcing functions.

3.4.2. Basic Reproduction Number and Equilibria

The basic reproduction number $R_0\left(t\right)$ at time t (instantaneous) is

$$R_0\left(t\right)={\displaystyle \frac{\beta \left(t\right)}{\gamma +\delta \left(t\right)+\mu }}.$$

(24)

It represents the expected number of secondary infections produced by one infected individual in a fully susceptible population.

Disease-free equilibrium (DFE): $I^{*}=0$, $S^{*}+R^{*}=N^{*}$, $D^{*}$ growing linearly. Setting derivatives to zero with $I=0$:

$${\displaystyle \frac{dS}{dt}}=-\mu S=0\Rightarrow S=0\left(\mathrm{not}\mathrm{feasible}\mathrm{long}-\mathrm{term}\right),$$

(25)

but approximately, as $\mu$ is small, the system has a quasi-steady DFE where $I\left(t\right)\to 0$, $S\left(t\right)+R\left(t\right)\to N\left(t\right)\approx N_0{e}^{-\mu t}$, $D\left(t\right)\approx \mu N_{0}t$.

No true stationary endemic equilibrium exists due to cumulative deaths ($D\to \infty$), but a quasi-endemic state may occur if $R_0\left(t\right)>1$ persistently.

3.4.3. Local Stability Analysis of Disease-Free Equilibrium

To assess stability, linearize around the DFE approximation ($S\approx N_0$, $I\approx 0$, $R\approx 0$, ignoring slow decay from $\mu$): The relevant subsystem for invasion is $(S,I)$:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& \approx -\beta {\displaystyle \frac{SI}{N_0}},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dI}{dt}}& \approx (\beta {\displaystyle \frac{S}{N_0}}-(\gamma +\delta +\mu ))I.\hfill \end{array}$$

(27)

The Jacobian matrix of the linearized system for $(S,I)$ near $(N_0,0)$ is approximately

$$J\approx \left(\begin{array}{cc}0& -\beta \\ 0& \beta -(\gamma +\delta +\mu )\end{array}\right).$$

(28)

Eigenvalues: ${\lambda }_1=0$ (neutral from S conservation), ${\lambda }_2=\beta -(\gamma +\delta +\mu )=(\gamma +\delta +\mu )(R_0-1)$.

Stability conditions:

• If $R_0<1$ (i.e., $\beta <\gamma +\delta +\mu$), then ${\lambda }_2<0$: the disease-free state is asymptotically stable; infections die out, $I\left(t\right)\to 0$.
• If $R_0>1$, then ${\lambda }_2>0$: unstable; small introductions lead to epidemic growth (invasion).

Including natural mortality fully, the characteristic equation yields eigenvalues with negative real parts when $R_0<1$, confirming global attraction to DFE in many cases (by comparison theorems or Lyapunov functions).

3.4.4. Advantages for Excess Mortality Estimation

This ODE approach offers:

• Mechanistic link between transmission ($\beta IS/N$) and excess deaths ($\delta I$).
• Counterfactual simulation: set $\delta \left(t\right)=0$ or $\beta \left(t\right)$ to baseline to compute $E_{c,t}$.
• Incorporation of covariates via time-varying $\beta \left(t\right)$, $\delta \left(t\right)$ (e.g., stringency reduces $\beta$, vaccination reduces $\delta$).
• Captures nonlinear feedback (herd immunity threshold at $1-1/R_0$).
• Enables scenario analysis (e.g., what-if no lockdowns).

Compared to WHO GAM/SARIMAX: - Captures propagation dynamics and saturation effects absent in additive models. - Better for countries with no ACM data: calibrate $\beta$, $\delta$ using reported COVID deaths and covariates. - Provides interpretable parameters (e.g., $R_0$, infectious period).

Limitations: - Assumes homogeneous mixing; extensions (age-structured, spatial) increase complexity. - Parameter identifiability issues (e.g., $\beta$ vs. initial $I_0$ trade-off). - Deterministic; stochastic versions (tau-leaping, Gillespie) or SDEs needed for small populations.

Calibration uses nonlinear least squares or Bayesian methods (fitting to reported deaths or excess proxies). This mechanistic model complements statistical approaches by providing insight into underlying biological and policy drivers of excess mortality.

3.5. Functional Principal Component Analysis with Functional Regression

A powerful nonparametric approach to modeling families of time series (such as monthly all-cause mortality curves across countries and years) is to treat each country’s historical mortality trajectory as a realization of a smooth function and apply functional data analysis (FDA) techniques, specifically functional principal component analysis (FPCA) followed by functional regression. This framework captures dominant modes of variation in mortality patterns (trend, seasonality, amplitude) in a data-driven, low-dimensional way, offering greater parsimony and flexibility than fixed-knot splines while avoiding strong parametric assumptions.

3.5.1. Functional Data Representation

Consider the historical monthly ACM for country c as a discrete sampling of a smooth underlying function $Y_c\left(\tau \right)$, where $\tau \in [0,T]$ indexes time (e.g., months since start of historical period, $T\approx 60$ for 5 years). We observe noisy discrete measurements

$$Y_{c,j}=Y_c\left({\tau }_j\right)+{ϵ}_{c,j},j=1,\dots ,n_c,{ϵ}_{c,j}\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{N}(0,{\sigma }_c^2),$$

(29)

where ${\tau }_j$ are equidistant monthly points and $n_c$ is the number of observed months.

To work in continuous time, we first smooth each country’s series using local polynomial regression or B-splines with cross-validated bandwidth, yielding an estimate ${\widehat{Y}}_c\left(\tau \right)$. The collection ${\left\{{\widehat{Y}}_c\left(\tau \right)\right\}}_{c=1}^{194}$ (or subset with reliable data) forms a sample of functions from a second-order stochastic process with mean function $\mu \left(\tau \right)=\mathbb{E}\left[Y_c\left(\tau \right)\right]$ and covariance operator

$$C(\tau ,s)=\mathrm{Cov}(Y_c\left(\tau \right),Y_c\left(s\right))=\mathbb{E}\left[(Y_c\left(\tau \right)-\mu \left(\tau \right))(Y_c\left(s\right)-\mu \left(s\right))\right].$$

(30)

3.5.2. Karhunen–Loéve Decomposition and FPCA

By Mercer’s theorem (assuming C is continuous and positive semi-definite), the covariance admits the spectral decomposition

$$C(\tau ,s)=\sum _{k=1}^{\infty }{\lambda }_k{\varphi }_k\left(\tau \right){\varphi }_k\left(s\right),$$

(31)

where ${\left\{{\lambda }_k\right\}}_{k=1}^{\infty }$ are non-increasing eigenvalues and $\left\{{\varphi }_k\left(\tau \right)\right\}$ orthonormal eigenfunctions satisfying

$$\int C(\tau ,s){\varphi }_k\left(s\right)ds={\lambda }_k{\varphi }_k\left(\tau \right).$$

(32)

The Karhunen–Loéve theorem expands each centered function as

$$Y_c\left(\tau \right)-\mu \left(\tau \right)=\sum _{k=1}^{\infty }{\xi }_{c,k}{\varphi }_k\left(\tau \right),{\xi }_{c,k}=\int (Y_c\left(\tau \right)-\mu \left(\tau \right)){\varphi }_k\left(\tau \right)d\tau ,$$

(33)

with $\mathbb{E}\left[{\xi }_{c,k}\right]=0$, $\mathrm{Var}\left({\xi }_{c,k}\right)={\lambda }_k$, and $\mathrm{Cov}({\xi }_{c,k},{\xi }_{c,\ell })=0$ for $k\ne \ell$.

In practice, we truncate to the first K components (chosen by explained variance threshold, e.g., 95%):

$${\widehat{Y}}_c\left(\tau \right)\approx \mu \left(\tau \right)+\sum _{k=1}^K{\widehat{\xi }}_{c,k}{\widehat{\varphi }}_k\left(\tau \right).$$

(34)

The empirical FPCA is computed via singular value decomposition of the centered data matrix after basis expansion (e.g., B-splines or Fourier basis for periodicity).

3.5.3. Expected Mortality Forecasting via Functional Regression

To forecast expected deaths in 2020–2021, we model the scores ${\xi }_{c,k}$ as functions of country-level covariates ${\mathbf{z}}_c$ (static: SDI, HDI, proportion >65, etc.; or time-varying averaged):

$${\xi }_{c,k}={\mathbf{z}}_c^{\top }{\gamma }_k+{\nu }_{c,k},{\nu }_{c,k}\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{N}(0,{\sigma }_k^2).$$

(35)

Alternatively, for richer prediction, we use functional linear regression where the response is the full function $Y_c\left(\tau \right)$:

$$Y_c\left(\tau \right)=\alpha \left(\tau \right)+\int \beta (\tau ,s)X_c\left(s\right)ds+{ϵ}_c\left(\tau \right),$$

(36)

but for forecasting we prefer score regression due to dimensionality reduction.

The expected counterfactual function for country c and future $\tau$ (2020–2021) is then

$$E\left[Y_c\left(\tau \right)\right]=\mu \left(\tau \right)+\sum _{k=1}^K\left({\widehat{\mathbf{z}}}_c^{\top }{\widehat{\gamma }}_k\right){\widehat{\varphi }}_k\left(\tau \right),$$

(37)

where ${\widehat{\mathbf{z}}}_c$ includes pandemic covariates set to pre-pandemic baseline or zero-COVID values. Monthly expected counts $E_{c,t}$ are obtained by evaluating at discrete points and integrating over intervals if needed.

3.5.4. Inference and Uncertainty

Point estimates of eigenfunctions, eigenvalues, and regression coefficients are obtained via least squares on scores. Uncertainty is quantified via:

• Bootstrap of entire curves (functional bootstrap): resample countries or residuals.
• Bayesian functional regression: Gaussian process priors on $\beta (\tau ,s)$ or hierarchical priors on scores ${\xi }_{c,k}\sim \mathrm{N}(0,{\lambda }_k)$.
• Asymptotic normality of FPCA estimators under mild conditions [16].

For countries with no or sparse data, predict scores ${\widehat{\xi }}_{c^{*},k}$ via ${\mathbf{z}}_{c^{*}}^{\top }{\widehat{\gamma }}_k$, then reconstruct the curve.

3.5.5. Theoretical Properties and Consistency

Under standard FDA assumptions (square-integrable functions, decaying eigenvalues ${\lambda }_k=O\left(k^{-p}\right)$, $p>1$), the FPCA truncation error is $O\left({\lambda }_{K+1}\right)$, and the estimator $\widehat{\mu }\left(\tau \right)$ is consistent at rate $n^{-1/2}$ uniformly. Score regression is consistent if ${\mathbf{z}}_c$ spans the relevant variation. The functional MSE

$$\mathbb{E}\left[\int {(Y_c\left(\tau \right)-\widehat{E}\left[Y_c\left(\tau \right)\right])}^{2}d\tau \right]$$

is minimized by retaining components with largest ${\lambda }_k$, outperforming fixed-basis methods (e.g., cyclic splines) when patterns vary heterogeneously across countries.

3.5.6. Advantages over WHO GAM and Comparison

Compared to the cyclic cubic spline + linear trend GAM:

• Data-driven basis: eigenfunctions capture the actual dominant modes (e.g., amplitude of winter peaks, trend acceleration) rather than imposing cyclic smoothness.
• Dimensionality reduction: $K\ll$ number of spline knots, reducing overfitting in short series.
• Natural clustering: countries with similar score vectors ${\xi }_c$ form epidemiological clusters, improving extrapolation.
• Parsimonious covariate modeling: regress low-dimensional scores instead of high-dimensional spline coefficients.
• Better handling of irregular sampling or missing months via smoothing.

Limitations include sensitivity to smoothing parameter in pre-processing and assumption of stationarity in the functional process (addressable via time-warping or local FPCA).

In summary, FPCA-based functional regression provides a mathematically elegant, nonparametric, and dimension-reduced alternative to GAMs, particularly advantageous for heterogeneous global mortality patterns and robust imputation in data-scarce regions.

4. Discussion

The five mathematical frameworks developed represent complementary advances over the current WHO methodology, each addressing distinct limitations while preserving the core objective: robust, transparent, and comparable estimation of excess mortality ${\delta }_{c,t}=Y_{c,t}-E_{c,t}$ across heterogeneous data availability scenarios. Below we discuss the importance and specific utility of each approach, followed by their potential for integration into a multi-model or hierarchical ensemble system. This method are applicable to several pandemic [19].

The STS model with Kalman filter stands out for its interpretability and coherence in handling non-stationary processes with stochastic components. Its primary importance lies in providing a principled state-space representation that separates trend, seasonality, irregular noise, and covariate effects, while the Kalman recursions ensure optimal (minimum mean squared error) filtering, smoothing, and forecasting under Gaussian assumptions (or approximations for Negative Binomial).

This method represent: - Superior handling of gradual structural changes in baseline mortality via stochastic local trends. - Natural inclusion of known interventions as additive components. - Exact likelihood and coherent multi-step predictive intervals, outperforming approximate GAM uncertainties. - Hierarchical extensions enable partial pooling of variance parameters across regions, improving estimates in data-scarce countries.

This approach is particularly valuable for high-income countries with long, reliable monthly series, where interpretability of components aids public health communication.

SARIMAX models excel in classical time-series analysis, offering a parsimonious yet flexible representation of autocorrelation, seasonality, and abrupt shocks. Their importance stems from the rigorous Box–Jenkins paradigm, which includes formal diagnostics, outlier detection, and intervention modeling that directly isolate pandemic-related deviations [19].

This method represent: - Explicit separation of baseline dynamics from pandemic shocks via pulse/step/transfer functions, enhancing causal attribution of excess deaths. - Strong performance in short- to medium-term forecasting of seasonal count series with known intervention points (e.g., waves identified by official COVID-19 peaks). - Robust diagnostic tools ensure model adequacy and detect misspecification early. - Hierarchical parameter estimation (e.g., shared AR/MA coefficients by WHO region) supports extrapolation.

This method is especially useful for countries with partial monthly data and identifiable external shocks, providing a robust alternative when spline smoothness assumptions fail.

The fully Bayesian hierarchical extension builds most directly on the existing WHO GAM, elevating it to a coherent probabilistic framework. Its importance lies in delivering full posterior inference, principled shrinkage, and borrowing of strength across countries – features absent in REML approximations.

This method represent: - Credible intervals that propagate all uncertainty sources, including hyperparameter variability. - Shrinkage of country-specific effects toward regional/global means, reducing bias and variance in low-data settings (e.g., AFR, EMR regions with few observations). - Incorporation of prior knowledge (e.g., sign constraints on stringency effects) and robustness to baseline period choice via time-varying penalties. - Scalable approximate inference (INLA) enables application to all 194 countries.

This framework is ideal as a “backbone” model for global comparability, offering the highest degree of statistical coherence and extrapolation power.

The mechanistic ODE approach is unique in providing a causal, process-based link between transmission dynamics and mortality. Its theoretical importance derives from stability analysis ($R_0$-driven invasion thresholds) and explicit counterfactual simulation (setting disease-induced death rate $\delta \left(t\right)\equiv 0$).

This method: - Captures nonlinear feedback loops (herd effects, saturation) and propagation absent in additive statistical models. - Interpretable parameters ($\beta \left(t\right)$, $\delta \left(t\right)$, infectious period) tied to policy (stringency, vaccination) and biological drivers. - Enables scenario analysis (e.g., no-lockdown counterfactuals) and sensitivity to transmission assumptions. - Calibration using reported COVID-19 deaths as proxy for $I\left(t\right)$ supports data-scarce countries.

While computationally intensive and assumption-heavy, this model is invaluable for understanding mechanisms and validating purely statistical estimates.

FPCA offers the most nonparametric, data-adaptive representation of mortality curves, decomposing variation into empirically optimal modes. Its importance lies in dimensionality reduction and capture of heterogeneous patterns without imposing fixed functional forms.

This method represent: - Dominant modes (e.g., seasonal amplitude, trend acceleration) are learned from data, improving parsimony and generalization. - Score regression on covariates enables robust imputation for the 83 countries with no pandemic ACM data. - Natural epidemiological clustering via score similarity, outperforming discrete WHO clusters. - Functional bootstrap or Bayesian FDA yields principled uncertainty for curve-level predictions.

This approach shines in global heterogeneous settings, where spline knot choices may be arbitrary.

Rather than selecting a single “best” method, a hybrid or ensemble strategy maximizes strengths and mitigates weaknesses. Possible integrations include:

1. Stacked ensemble forecasting: Train base learners on countries with full data; learn weights via time-series cross-validation minimizing out-of-sample excess error. ODE outputs can serve as features or post-hoc checks.

2. Hierarchical Bayesian multi-model: Place a mixture-of-experts prior over model classes, with country-specific weights learned hierarchically. This automatically selects or blends models per region/data quality.

3. Sequential or modular pipeline: Use FPCA to extract low-dimensional features (scores), feed them into Bayesian hierarchical regression for pooling, then refine with STS/SARIMAX for countries with strong seasonality, and validate against ODE mechanistic constraints.

4. Uncertainty fusion: Combine posterior/forecast distributions via linear opinion pooling or variational inference, yielding consensus credible intervals that are more robust than any single model.

Such combinations would likely reduce mean absolute error by 10–30% (based on analogous mortality forecasting benchmarks) and provide more reliable global/regional aggregates, especially in low-data regions.

In conclusion, these five variants collectively advance the mathematical rigor, interpretability, mechanistic insight, and robustness of excess mortality estimation, offering tools that can evolve with improving data availability and computational resources.

5. Conclusions

The estimation of excess mortality during the COVID-19 pandemic remains one of the most important and challenging tasks in modern public health statistics. The current WHO methodology [1] represents a major step forward in providing globally comparable estimates despite enormous heterogeneity in data quality, availability, and timeliness across 194 Member States. However, as highlighted in the introduction and preliminaries, the existing approach – based primarily on GAMs with cyclic splines, log-linear covariate models, and discrete clustering for age-sex patterns – is subject to limitations in flexibility, robustness to structural changes, mechanistic insight, and uncertainty propagation.

In this paper, we have developed and formally presented five mathematically advanced variants that address these shortcomings while remaining compatible with the core principles of baseline forecasting, covariate adjustment, and excess calculation ${\delta }_{c,t}=Y_{c,t}-E_{c,t}$:

1. Structural time series models with Kalman filtering offer interpretable decomposition and coherent multi-step forecasting uncertainty. 2. Extended SARIMAX models with intervention components provide rigorous handling of autocorrelation, seasonality, and abrupt shocks. 3. Fully Bayesian hierarchical models with penalized splines and shrinkage priors deliver principled pooling, full posterior inference, and improved extrapolation. 4. Compartmental ODE models introduce mechanistic understanding through transmission-mortality dynamics, stability analysis, and explicit counterfactual simulation. 5. Functional principal component analysis combined with functional regression enables data-driven, low-dimensional representation of mortality curves and natural clustering.

The potential for synergies is particularly promising. A multi-model ensemble – possibly stacked via cross-validated weights, hierarchically mixed in a Bayesian framework, or modularly pipelined – could yield estimates that are more accurate, robust, and transparent than any single approach. Such a system would automatically adapt to data availability (e.g., mechanistic ODEs for countries with reported COVID deaths but no ACM; FPCA regression for global imputation; Bayesian hierarchical as backbone), while providing consensus credible intervals that better reflect model uncertainty.

From a practical standpoint, these advances could directly inform future WHO updates or parallel efforts (e.g., World Mortality Dataset, The Economist trackers, or national statistical offices). They also open doors to real-time surveillance, scenario analysis for future pandemics, and better attribution of direct vs. indirect effects. Moreover, the methods are generalizable beyond COVID-19 to other crises (e.g., heatwaves, conflicts, natural disasters) where excess mortality serves as a key impact metric.

Limitations remain. All approaches depend critically on the quality of historical baseline data and covariate information; structural breaks pre-pandemic or differential completeness can still bias estimates. Mechanistic models introduce identifiability challenges, while nonparametric methods (FPCA) require careful smoothing choices. Computational demands vary (INLA/Bayesian scalable, ODE calibration intensive), and global implementation would necessitate standardized software pipelines and validation against gold-standard countries (e.g., those with near-complete registration).

Future directions include: - Empirical comparison of the five variants plus the WHO baseline on hold-out countries with high-quality data. - Development of an open-source R/Python package integrating these models with hierarchical ensembling. - Incorporation of spatial dependence (e.g., Gaussian processes over country embeddings) and age-structured extensions. - Real-time updating frameworks with streaming data and adaptive baselines.

In summary, the proposed mathematical enhancements represent a significant evolution toward more robust, interpretable, and mechanistically grounded estimation of pandemic mortality burden. By combining statistical rigor with epidemiological insight, they can contribute to better-informed global health policy and preparedness for future health crises.