Hierarchical Bayesian autoregressive smooth transition time series models

1 Introduction

Timely public health responses to infectious disease outbreaks critically depend on accurate and robust forecasting of disease burden. Time series models have long served as essential tools for capturing temporal patterns in infectious disease data, including seasonality, trend shifts, and outbreak dynamics [1, 2]. Classical linear autoregressive (AR) models are widely used due to their interpretability, computational simplicity, and well-established theoretical properties [3–6]. These models are particularly effective for stationary or near-stationary processes, where the mean, variance, and autocorrelation structure remain stable over time [7, 8].

However, infectious disease data often exhibit non-linear dynamics, including abrupt regime shifts, asymmetric transmission patterns, and time-varying covariate effects [8, 9]. Extensions of AR models using higher-order terms or interactions may partially address these features, but at the cost of increased model complexity and risk of overfitting [10, 11]. Consequently, non-linear time series models have gained attention for their ability to capture the complex evolution of disease incidence [12–14]. Prominent approaches include threshold autoregressive (TAR) models [15, 16], Markov-switching autoregressive (MSAR) models [17, 18], and smooth transition autoregressive (STAR) models [19–21].

STAR models are particularly attractive for infectious disease applications due to their flexibility in modeling gradual, continuous transitions between regimes. Logistic STAR (LSTAR) models are suitable for asymmetric or threshold-like transitions, whereas exponential STAR (ESTAR) models capture smooth, continuous exponential shifts. This unifying structure nests linear AR and TAR models as special cases, facilitating robust modeling of epidemic waves, intervention effects, and non-linear dynamics [19, 22–24].

Extensions of non-linear time series modeling also include ARCH and GARCH frameworks [25, 26], as well as EGARCH and threshold heteroscedastic models [24, 27, 28], which account for volatility clustering and asymmetric uncertainty. These developments acknowledge that epidemiological shocks—such as surges in cases—often induce greater uncertainty than declines. Further refinements, including Double Smooth Transition GARCH (DST-GARCH) models [29, 30], introduce multiple transition estimation and adaptive MCMC strategies for improved computational efficiency and identifiability.

Despite these advances, applications of STAR models in infectious disease forecasting remain limited. Standard methods often assume linearity within regimes and lack hierarchical structures to jointly model multiple related populations [31–34]. Maximum likelihood estimation may struggle in these settings, particularly when data are sparse, unevenly reported, or exhibit latent regime shifts [29, 35, 36].

Bayesian hierarchical models offer a natural solution, enabling information pooling across regions, handling latent regime changes, and quantifying parameter uncertainty [37–39]. Simulation-based inference techniques, such as Markov Chain Monte Carlo (MCMC), facilitate estimation in complex non-linear structures and allow the incorporation of prior epidemiological knowledge [39–41].

This study addresses these gaps by utilizing a Bayesian Hierarchical Smooth Transition Autoregressive (BH-STAR) model for infectious disease forecasting. The modeling approach accommodates non-linear transitions, asymmetric effects, and country-level heterogeneity, leveraging hierarchical Bayesian structuring for information sharing and MCMC for computational tractability. In particular, we compare the BH-STAR models against two baseline competitors: a standard autoregressive (AR) model and a heterogeneous autoregressive (HAR) model. The AR model represents the simplest linear alternative, while the HAR model incorporates multiple time scales and has been widely used in infectious disease forecasting [36]. Non-linear alternatives such as threshold autoregressive (TAR) and Markov-switching autoregressive (MSAR) models are conceptually related, but STAR models offer the advantage of smooth transitions, which better reflect gradual changes in epidemic dynamics (e.g., waning immunity or phased interventions). A systematic comparison with TAR/MSAR is left for future work. The models are compared through simulation studies and empirical applications to African COVID-19 data, benchmarking performance against conventional AR models.

2 Materials and methods

2.1 COVID-19 data description and study settings

We obtained monthly counts of confirmed COVID-19 cases and deaths from 53 African countries. The data span March 2020 to December 2022 and were sourced from publicly available COVID-19 data at the Johns Hopkins Coronavirus Resource Center (https://coronavirus.jhu.edu/about/how-to-use-our-data). This period captures multiple epidemic waves, varying reporting capacities, and heterogeneous policy responses across the continent.

2.1.1 Spatial distribution of COVID-19 burden

Figure 1 presents the spatial distribution of COVID-19 burden across Africa. Panel A shows cumulative confirmed cases, revealing substantial heterogeneity: South Africa, Morocco, Tunisia, Egypt, and Ethiopia reported the highest case counts (deep red shading), while many Central and West African countries experienced lower cumulative numbers (yellow to light orange). Panel B displays reported deaths, with South Africa accounting for a disproportionate share of mortality, followed by Tunisia, Egypt, and Morocco.

Panel C maps incidence rates per 100,000 population, which normalizes for population size. Notably, some smaller populations (e.g., Cabo Verde, Mauritius, Botswana) exhibit higher incidence rates than their raw case counts would suggest. Panel D presents the case fatality rate (CFR) per 1,000 cases, where Sudan, Liberia, and several Central African countries show elevated CFRs, potentially reflecting differences in healthcare access, reporting quality, or age structure.

Each metric offers distinct advantages for understanding epidemic dynamics. Cumulative confirmed cases (Panel A) provide an intuitive measure of the absolute burden on healthcare systems and are easily communicated to policymakers, though they are strongly influenced by population size and testing capacity. Reported deaths (Panel B) are less susceptible to testing variability than cases, as severe outcomes are more consistently ascertained, but they lag behind infections and depend on accurate cause-of-death certification. Incidence rates per 100,000 population (Panel C) enable fair comparisons across countries with different population sizes, revealing hidden hotspots where transmission intensity is high despite modest absolute case counts. Case fatality rate per 1,000 cases (Panel D) reflects disease severity and healthcare system performance, though it should be interpreted cautiously because it is sensitive to under-ascertainment of mild infections and differences in age structure. Together, these four metrics provide a multi-faceted view of the pandemic, balancing absolute burden, population-adjusted risk, mortality impact, and healthcare system resilience.

2.1.2 Temporal dynamics of cases and deaths

Figures 2, 3 display monthly time series of observed COVID-19 cases and deaths, respectively, for selected countries (Kenya, Nigeria, Senegal, South Africa, Tanzania). These plots reveal distinct epidemic waves with varying timing and amplitude across countries. South Africa experienced a sharp peak in mid-2020 followed by subsequent waves, while Kenya and Nigeria showed more gradual increases. Tanzania reported consistently low counts throughout the period, though this may reflect under-reporting. The temporal heterogeneity across countries motivates the use of hierarchical models that can borrow strength across regions while accommodating country-specific dynamics.

3 Hierarchical smooth transition time series modeling and evaluation

3.1 Hierarchical smooth transition time series modeling

Let yit∈ℝ+ denote confirmed COVID-19 cases for country i at time t. We extend the classical autoregressive (AR) framework to a hierarchical smooth transition autoregressive (STAR) model structure to accommodate non-linear dynamics and capture regime-switching behaviors inherent in infectious disease progression. This modeling approach allows for heterogeneity across countries while leveraging shared information through a Bayesian hierarchical framework.

3.2 Classical AR specification

The basic AR(p) specification for each country is modeled dynamically. As shown in Equation 1, the base ne transmission structure is defined as

yit=ϕ0i+∑j=1pϕjiyi,t-j+ϵit, ϵit~N(0,σi2),(1)

where ϕ_*ji*∈ℝ are country-specific autoregressive coefficients. Here, (e_*it*) denotes the country specific error term for country (i) at time (t), and (σ_*i*) denotes the country specific error standard deviation. Although effective for short-term stationarity, this model assumes linearity and constant parameters, which are often violated during epidemic waves or policy shifts.

3.3 Smooth transition autoregressive (STAR) model

The AR(p) generalizes to a STAR model where regime shifts occur smoothly over time via a bounded transition function G_*it*∈[0, 1]. This structural extension is formalized in Equation 2:

yit=xit⊤[ϕi(1)(1-Git)+ϕi(2)Git]+ϵit,(2)

where xit=(1,yi,t-1,…,yi,t-p)⊤ and ϕi(k) are regime-specific AR coefficients (k = 1, 2).

The threshold parameter locates the midpoint of the transition. Two transition functions are considered in this study: the logistic transition function and the exponential transition function. Their mathematical specifications are given in Equations (3) and (4), respectively.

• Logistic:

Git(L)=[1+exp(-γi(sit-θi))]-1,(3)

• Exponential

Git(E)=1-exp(-γi(sit-θi)2),(4)

In Equations (3) and (4), (s_*it*) denotes the transition variable for country (i) at time (t), typically taken as a lagged value of the response variable. Often the lag 1 is used to capture the most recent epidemic intensity and has been commonly used in epidemiological STAR applications [42]. The parameter (θ_*i*) is the country specific threshold parameter that determines the midpoint of the transition, while (γ_*i*>0) is the smoothness parameter controlling the rate at which the transition occurs. Larger values of (γ_*i*) produce a more abrupt regime change (For example, following a sudden public health intervention), while smaller values imply a more gradual and smoother transition.

3.4 Bayesian hierarchical specification

We complete the model by specifying prior distributions for all parameters and adopting a hierarchical structure that allows partial pooling of information across countries. Let ϕi(k)=(ϕ0i(k),ϕ1i(k),…,ϕpi(k))⊤ denote the regime-specific autoregressive coefficients for country i, and let γ_*i* and θ_*i* be the smoothness and threshold parameters in the transition function G_*it* (Eq. 3 or 4).

At the first level of the hierarchy, we assume

ϕi(k)~N(μϕ(k),Σϕ(k)), k=1,2,log(γi)~N(μγ,τγ2), θi~N(μθ,τθ2),log(σi)~N(μσ,τσ2).

The hyperparameters μϕ(k), μ_γ, μ_θ, μ_σ, and the variance components Σϕ(k), τγ2, τθ2, τσ2 are assigned weakly informative priors:

μϕ(k)~N(0,10I), μγ~N(0,2), μθ~N(0,2), μσ~N(0,1), Σϕ(k)~Inverse-Wishart(p+2,I),τγ,τθ,τσ~Half-Cauchy(0,2.5).

To improve sampling efficiency, we adopt a non-centered parameterization for the country-level random effects [43]. The full posterior is proportional to

∏i=153∏t=1Tip(yit∣ϕi(1),ϕi(2),γi,θi,σi,xit)×p(ϕi(k)∣μϕ(k),Σϕ(k)) ×p(γi∣μγ,τγ)×p(θi∣μθ,τθ) ×p(σi∣μσ,τσ)×p(hyperparameters).

Shrinkage is quantified by the posterior variance of the hyperparameters; countries with sparse data borrow strength from the pooled estimates.

3.5 Simulation study

A simulation study was conducted to evaluate the finite-sample properties of the proposed Smooth Transition Autoregressive (STAR) models under controlled non-linear data-generating mechanisms. The standard STAR data-generating process was specified in its reparameterized form as

yit=α0+ϕ~1yi,t-1+ϕ~2yi,t-1G(yi,t-1;γ,θ)+βi+ϵit,

where G(·) denotes a smooth transition function (logistic or exponential), θ represents the threshold parameter, β_*i* captures country-level random effects, and ϵ_*it* is an idiosyncratic error term. The composite parameters are defined as ϕ~1=α1+ϕ1 and ϕ~2=α2+ϕ2, ensuring identification of the effective autoregressive and non-linear regime-dependent effects. For simulation purposes, the underlying structural components were calibrated such that α₁ = 0 and α₂ = 0, implying ϕ~1=ϕ1 and ϕ~2=ϕ2, thereby preserving interpretability of the autoregressive parameter of interest without affecting inferential conclusions.

3.5.1 Data-generating process

Data were generated from a hierarchical autoregressive process with random effects and a non-linear regime-switching component governed by the smooth transition mechanism. The true data-generating parameters were held fixed across all simulation runs, including the autoregressive coefficient of primary interest (ϕ₁ = 0.6), baseline intercept (α₀ = 0.5), regime-dependent parameter (ϕ₂ = −0.4), and variance components for the observation error (σ_ϵ = 0.5) and random effects (σ_*b* = 0.7), while the transition location was fixed at θ = 0. The smoothness parameter γ was systematically varied to represent low (γ = 1), moderate (γ = 10), and fast transition (γ = 50) regimes, thereby controlling the degree of regime separation and non-linear smoothness in the data-generating process.

For the exponential DGP, the transition function was specified as:

GESTAR(yi,t-1;γ,θ)=1-exp(-γ(yi,t-1-θ)2).

For the logistic DGP, the transition function was specified as:

GLSTAR(yi,t-1;γ,θ)=11+exp(-γ(yi,t-1-θ)).

3.5.2 Simulation design

For each simulated dataset, competing models including HLSTAR, HESTAR, AR, and HAR specifications were estimated using Bayesian Markov chain Monte Carlo methods. The transition smoothness parameter γ was treated as known and fixed within each model to ensure that comparisons focused on structural specification. For each replication, we generated time series of length T = 120 (months) for N = 5 countries. This relatively small number of countries and a moderate time series length reflect typical data availability in the African COVID-19 context, where many countries have limited reporting periods.

We fitted the BH-LSTAR and BH-ESTAR models as described in Section 3 and compared them with a pooled AR(2) model. Model performance was assessed over 100 Monte Carlo replications and 2,000 MCMC iterations using three metrics:

• Bias: The systematic deviation of the posterior mean estimator of ϕ₁ from its true value (0.6), calculated as Bias=1R∑r=1R(ϕ^1(r)-0.6).
• Mean Absolute Error (MAE): The average absolute difference between observed and predicted values, measuring overall predictive accuracy.
• Root Mean Squared Error (RMSE): The square root of the average squared prediction error, penalizing large deviations more heavily.

Table 1 summarizes the simulation scenarios considered in this study. The results are presented in the following sections. First, both STAR models (HESTAR and HLSTAR) consistently outperform the linear AR model in terms of MAE and RMSE across all scenarios, confirming the value of allowing non-linear regime-switching dynamics. Second, and more critically, the results reveal distinct asymmetries in the consequences of model misspecification.

Table 1. Simulation results for BH-STAR models under exponential and logistic data-generating processes.

Values are averages over 100 Monte Carlo replications, with time series of length T = 120 months for N = 5 countries. Bias is reported for the autoregressive parameter ϕ₁ (true value = 0.6). MAE and RMSE are based on one-step-ahead predictions. The smoothness parameter γ is treated as known within each scenario.

Under the exponential DGP, the correctly specified HESTAR model yields near-zero bias across all values of γ (bias ≈ −0.057, −0.041, −0.037), with consistently low MAE (≈ 0.398) and RMSE (≈ 0.501). The HLSTAR model, although misspecified under exponential dynamics, exhibits only moderate and stable bias (bias = 0.142, −0.049, −0.171 for γ = 1, 10, 50 without explosive behavior. This indicates that the logistic transition function can approximate the symmetric exponential transition reasonably well.

Under the logistic DGP, the correctly specified HLSTAR model shows almost no bias in all smoothness regimes (bias ≈ −0.018, −0.165, −0.148) and stable predictive performance (MAE ≈ 0.397, RMSE ≈ 0.500). In contrast, the misspecified HESTAR model exhibits rapidly increasing bias as γ increases (bias = −0.278, −1.662, −4.525), reflecting a fundamental incompatibility between the exponential transition function and the asymmetric regime-switching mechanism of the true DGP.

This asymmetry arises from fundamental differences in the functional geometry of the transition functions. The exponential transition function is symmetric and smooth around zero, allowing the logistic specification to approximate it reasonably well even under misspecification. Consequently, the induced parameter distortion remains limited and does not amplify substantially with increasing γ.

By contrast, the logistic transition function introduces directional asymmetry, which the exponential form cannot reproduce. When HESTAR is fitted to logistic data, the model fails to capture this asymmetry and compensates by reallocating non-linear dynamics into the linear autoregressive component. This induces strong parameter confounding and leads to increasingly severe bias, particularly as γ grows and the transition function approaches a near-threshold regime γ = 50, where the logistic function approximates a step function.

Importantly, these discrepancies in parameter bias are not reflected in predictive performance. Across both DGPs, MAE and RMSE remain remarkably stable. For instance, under the logistic DGP at γ = 50, HESTAR exhibits extreme bias (–4.525) yet maintains MAE ≈ 0.402 and RMSE ≈ 0.507, which are comparable to the correctly specified HLSTAR (MAE ≈ 0.396, RMSE ≈ 0.500). Similarly, under the exponential DGP, the predictive performance of HLSTAR remains close to that of the correctly specified HESTAR model (MAE ≈ 0.401–0.406 vs. 0.398; RMSE ≈ 0.502–0.508 vs. 0.498).

This pattern is expected given the simulation design, which used only five countries and 120 data points per replication. With limited data, predictive metrics average out individual deviations, whereas parameter bias captures structural misalignment that does not necessarily translate into large forecast errors, especially when the outcome is on a transformed scale (e.g., log-cases). In the data-generating process, the non-linear component enters through the transition function, which governs how the effective autoregressive dynamics evolve across regimes. When the transition function is correctly specified, the model successfully recovers the structural decomposition of the conditional mean. When misspecified, the model can still achieve similar predictive accuracy by adjusting other parameters (e.g., autoregressive coefficients or intercepts) to compensate, even if the individual transition parameters are severely biased.

These findings suggest an asymmetry in the consequences of model misspecification:

• Misspecifying exponential dynamics with a logistic transition (HLSTAR under an exponential DGP) results in moderate, stable bias.
• Misspecifying logistic dynamics with an exponential transition (HESTAR under a logistic DGP) leads to severe, rapidly increasing bias, especially as γ increases.

Thus, the impact of misspecification is inherently directional, depending on the extent to which one transition function can approximate the other. For COVID-19 data in Africa, where epidemic waves often exhibit asymmetric rise-and-fall patterns (rapid increases followed by slower declines), the logistic transition function may be more appropriate. This is supported by the empirical results in Section 3.6, where the LSTAR model produced narrower predictive intervals and more precise parameter estimates than its exponential counterpart.

3.6 Empirical application: COVID-19 forecasting in Africa

The Bayesian Hierarchical Smooth Transition Autoregressive (BH-STAR) model was applied to monthly confirmed COVID-19 cases across 53 African countries. The training period spanned March 2020–September 2022, with October–December 2022 used for validation and out-of-sample forecasting. Convergence of the MCMC chains was assessed using the potential scale reduction factor R^ [44] and the effective sample size (ESS).

A log transformation with a small offset (log(y+0.01)) was applied to case counts to stabilize variance and handle near-zero values. This transformation also helps to normalize the error distribution and reduce the influence of extreme values. The logarithm of the lagged case counts served as the transition variable s_*it* with lag d = 1, capturing the immediate past incidence as the driver of regime shifts.

Reporting gaps were present in several countries, with sporadic missing months. Because our model is likelihood-based and we use MCMC, missing values in the response variable y_*it* were treated as latent variables and imputed during sampling. Specifically, at each MCMC iteration, missing y_*it* were drawn from their conditional posterior distribution given the observed data and current parameters. This approach accounts for the uncertainty due to incomplete records and avoids ad-hoc imputation. Countries with consistently sparse reporting (e.g., fewer than 12 months of data) still contribute information through the hierarchical structure, which pulls their estimates toward the continental averages.

3.6.1 Prior specification for smoothing parameters and scale parameters in Bayesian STAR models

In a Bayesian STAR framework, the choice of prior for the smoothing parameter γ should reflect both theoretical constraints and empirical regularities observed in applications. Since γ must be strictly positive and is often weakly identified with a relatively flat or multi-modal likelihood, priors play a stabilizing role in estimation. The smoothing parameter is known to be weakly identified with a relatively flat likelihood, often leading to imprecise estimation and sensitivity to initialization [42, 45]. In a Bayesian context, this motivates the use of weakly informative priors to regularize estimation [46].

A log-normal prior, specified as log(γ)~N(μ,σ2), is commonly preferred because it enforces positivity while allowing flexible, right-skewed uncertainty; typical choices such as μ = 0 and σ∈[1, 2] center the parameter around moderate transition speeds without overly restricting the tail. Alternatively, a Gamma prior (e.g., γ~Gamma(2, 1)) provides mild regularization toward moderate values, while bounded Uniform priors are less recommended due to sensitivity to the upper bound and potential over-weighting of extreme regime sharpness.

In the context of COVID-19 data in Africa—characterized by noise, reporting inconsistencies, and structural breaks—moderate values of γ are generally more plausible, as very large values may spuriously capture noise as abrupt regime shifts. A practical and defensible specification is therefore log(γ)~N(0,1.52), which balances flexibility and regularization. To further enhance robustness, especially in multi-country analyses, a hierarchical prior structure can be adopted, allowing country-specific γ_*i* to be partially pooled via log(γi)~N(μγ,σγ2), thereby improving stability and borrowing strength across countries. The prior distribution for the hyperparameter (μ_γ), representing the population mean of the country specific smoothness parameters, is specified as [μ_γ~Normal(0, 2)].

For scale parameters such as the residual standard deviation σ, the between country standard deviation β_sigma, and the hierarchical variance σ_γ, we used half-normal priors rather than the conventional half-Cauchy. While half-Cauchy priors offer mild shrinkage and good behavior near zero [47], their heavy tails can occasionally produce unrealistically large posterior draws when the likelihood is flat. Thus we used the half-normal prior, σ~Half-Normal(0, τ), where τ controls the scale. The half-normal is a special case of the half-t with ν = ∞; it provides a finite prior mean and variance, and yields slightly stronger regularization toward zero compared to the half-Cauchy. For weakly identified scale parameters, choosing τ = 1 or τ = 2 offers a sensible weakly informative prior that avoids extreme values while still allowing the data to dominate if sufficient information exists.

For the BHESTAR and BHLSTAR models, we therefore specified:

σ~Half-Normal(0,1), βsigma~Half-Normal(0,1), σγ~Half-Normal(0,1).

This choice retains the desirable property of penalizing very large standard deviations without the occasional pathological behavior of the Cauchy tail.

3.6.2 Convergence assessment and posterior predictive checks

Convergence of the MCMC chains was assessed using trace plots. Figures 4, 5 display trace plots for the logistic STAR (LSTAR) and exponential STAR (ESTAR) models, respectively, for all global parameters over 4,000 iterations with four chains (red, green, blue, black).

For the LSTAR model (Figure 4), the chains for α₀, α₁, α₂, ϕ₁, ϕ₂, θ, σ, β_σ, μ_γ, σ_γ show good mixing. All four chains overlap substantially after the warmup period, with no evidence of parallel behavior, drift, or divergent transitions. The red, green, blue, and black chains explore the same posterior region, indicating that the sampler has converged to the target distribution.

For the ESTAR model (Figure 5), the trace plots for most parameters also show acceptable mixing. However, mild parallel tendencies are visible for γ_μ and γ_*sigma* in the first 1,000 iterations, where the chains appear to explore slightly different regions before converging. By iteration 2,000, all chains have stabilized and overlap reasonably well. This behavior suggests that the hierarchical smoothing parameters may be weakly identified in the exponential specification, but the chains ultimately reach a common distribution.

For both models, the Gelman–Rubin R^ statistic was below 1.01 for all parameters, and effective sample sizes exceeded 200 for all quantities of interest, confirming that the chains have converged and that posterior inferences are reliable.

3.6.2.1 Posterior predictive checks

Posterior predictive checks were performed to assess how well the fitted BH-STAR models reproduce the observed data. Figures 6 and 7 present posterior predictive checks for the ESTAR and LSTAR models, respectively, for the five selected countries (Kenya, Nigeria, Senegal, South Africa, and Tanzania). The training period (March 2020–September 2022) is used for model fitting, and the test period (October–December 2022) evaluates out-of-sample forecasts. For all five countries, the observed values lie within the credible intervals, indicating adequate model fit.

For the ESTAR model (Figure 6), the posterior predictive medians generally track the observed values well during the training period. In Kenya, the 95% credible intervals comfortably contain the observed counts, indicating adequate model fit. Senegal shows similar performance, with predictions closely matching observations. Tanzania, which reported consistently low counts, has narrow credible intervals that capture the observed values. During the test period (September—December 2022), the model maintains reasonable predictive accuracy, though the credible intervals widen slightly, reflecting increased uncertainty when forecasting beyond the observed data.

For the LSTAR model (Figure 7), the posterior predictive medians are largely consistent with the ESTAR results. However, the LSTAR model exhibits narrower credible intervals for most countries, particularly during the test period. In Nigeria, the 95% credible intervals are tighter while still encompassing the observed values, suggesting that the logistic transition function provides more precise forecasts for these data. Kenya and Senegal also show stable predictions with good coverage. South Africa (partial data shown) indicates that the model captures the general level of cases, though longer test series would be needed for a full assessment.

Overall, both models demonstrate satisfactory posterior predictive performance, with observed values falling within the 95% credible intervals for the majority of months. The LSTAR model appears to offer slightly tighter uncertainty bounds, which may be advantageous for public health decision-making where precise forecasts are valued. Formal out-of-sample metrics (MAE, RMSE, MdAPE) are reported in Section 3.6.6.

3.6.3 Posterior parameter summaries

Posterior summaries for key model parameters are presented in Table 2. The AR(1) coefficient in regime 1 (ϕ⁽¹⁾) has a posterior mean of 0.65 (95% credible interval: 0.42, 0.87), indicating strong positive autocorrelation during high-incidence periods. In contrast, the AR(1) coefficient in regime 2 [ϕ⁽²⁾] is much smaller (mean 0.12, 95% CI: −0.05, 0.31), suggesting that case counts are nearly uncorrelated or weakly mean-reverting during low-transmission phases. The transition smoothness parameter γ has a posterior mean of 6.87 (95% CI: 3.21, 11.02), indicating moderately rapid transitions between regimes. The threshold parameter θ (mean 1.48, 95% CI: 0.78, 2.12) locates the midpoint of the transition on the log-scale. The error standard deviation σ is estimated at 0.42 (95% CI: 0.28, 0.59). These results are consistent with epidemic dynamics: strong persistence during waves, weaker dependence during lulls, and smooth but relatively quick regime shifts.

Table 2. Posterior summaries for selected BH-STAR parameters across 52 countries. Posterior summaries include the posterior mean, standard deviation (SD), and 95% credible interval.

3.6.4 Posterior predictive performance

The posterior predictive distribution was used to generate three-month-ahead forecasts for each country, covering the testing period October 2022 to December 2022. Figures 8, 9 present the forecast medians and 95% predictive intervals alongside the observed monthly case counts for selected countries (Kenya, Nigeria, South Africa, Tanzania) under the LSTAR and ESTAR models.

For both models, all observed values during the testing window fall within the 95% predictive intervals, indicating that the proposed framework produces well-calibrated probabilistic forecasts. In some instances, such as Kenya in November 2022, the observed count exceeds the forecast median substantially; however, the credible intervals remain sufficiently wide to accommodate this deviation, reflecting the model's honest quantification of uncertainty in the presence of volatile epidemic dynamics. The ability to maintain coverage despite data sparsity and reporting irregularities confirms the practical usability of the BH-STAR models for real-time outbreak monitoring.

3.6.5 Regime identification

Posterior summaries of the transition function G_*it*, which represents the probability of being in the high-incidence regime, reveal clear regime-switching patterns corresponding to epidemic waves. Figures 10, 11 present the estimated regime transition probabilities by month for the LSTAR and ESTAR models, respectively, for Kenya, Nigeria, South Africa, and Tanzania.

For the LSTAR model (Figure 10), South Africa exhibits a sharp rise in transition probability from approximately 0.45 in March 2020 to near 1.00 by May 2020, remaining at or above 0.95 throughout the pandemic. This reflects the country's sustained periods of high incidence and multiple distinct waves. Nigeria follows a similar but slightly delayed pattern, with probabilities exceeding 0.95 since May 2020. Kenya shows a more gradual increase, reaching 0.75 in March 2020 and climbing to 1.00 by May 2020, remaining elevated thereafter. Tanzania presents a unique pattern: the probability fluctuates around 0.50 for most of 2020, with modest increases to approximately 0.95 in later months, consistent with its reported low case counts and potential under-reporting.

For the ESTAR model (Figure 11), the transition probabilities exhibit a different temporal structure. South Africa rose from 0.00 in March 2020 to 0.95 by May 2020, maintaining high probabilities (0.95–0.99) throughout. Nigeria and Kenya follow similar trajectories. Notably, Tanzania shows near-zero probability in March 2020, rising sharply to 0.95 by May 2020 and remaining elevated thereafter—a pattern that differs from the LSTAR model's more moderate estimates for Tanzania. This discrepancy may reflect the exponential transition function's tendency to assign more extreme probabilities (close to 0 or 1) compared to the logistic specification.

Across both models, peaks in G_*it* align with observed epidemic waves. For South Africa, the sustained high probability corresponds to the country's severe and prolonged first and second waves (mid-2020 and early 2021). Nigeria and Kenya show similarly persistent high probabilities, reflecting their ongoing transmission throughout the study period. Tanzania's lower and more variable probabilities under the LSTAR model (oscillating around 0.50) suggest genuine uncertainty about its regime status, likely due to sparse reporting. Under the ESTAR model, Tanzania's probabilities are more extreme, indicating less flexibility in capturing intermediate states.

These results demonstrate that the BH-STAR framework successfully identifies dynamic transitions in disease activity, with the logistic specification offering more nuanced probability estimates for countries with ambiguous or sparse data. The regime probabilities can serve as an operational tool for public health agencies, providing early warning of increasing transmission when G_*it* crosses a predefined threshold (e.g., 0.5 or 0.75).

3.6.6 Forecast accuracy evaluation

Table 3 presents the aggregated forecast accuracy metrics for the 52 countries over the three month validation period (October to December 2022). To comprehensively evaluate forecasting performance, both scale dependent measures, namely Mean Absolute Error (MAE) and Root Mean Squared Error (RMSE), and percentage based relative measures, namely Mean Absolute Percentage Error (MAPE) and Median Absolute Percentage Error (MdAPE), were employed. The scale dependent metrics quantify the absolute magnitude of forecast errors, whereas the percentage based metrics facilitate standardized comparisons across countries with differing disease burdens and population sizes.

Table 3. Mean validation forecast accuracy across the 52 countries, October–December 2022.

The nonlinear Bayesian Hierarchical Smooth Transition Autoregressive (BH STAR) models consistently outperformed both the linear AR and HAR benchmarks across all four evaluation metrics. Specifically, the logistic specification (HLSTAR) achieved the highest overall forecasting accuracy, with an MAE of 0.58, an RMSE of 0.79, a MAPE of 14.2%, and an MdAPE of 12.8%. The exponential specification (HESTAR) produced similarly strong performance, with an MAE of 0.61, an RMSE of 0.82, a MAPE of 14.8%, and an MdAPE of 13.1%. The HAR model showed moderate improvements over the linear AR model (MAE: 0.91 vs. 1.22; RMSE: 1.12 vs. 1.58; MAPE: 21.4% vs. 27.5%; MdAPE: 18.9% vs. 24.1%). Overall, both BH STAR models substantially reduced forecast errors relative to the linear benchmarks, with country specific reductions in MAE and RMSE ranging from approximately 30% to 50% and overall improvements exceeding 50% compared with the aggregated AR benchmark.

Notably, the inclusion of MdAPE provides a robust complement to MAPE, as the median is less sensitive to extreme country specific forecast errors. The consistently higher MAPE values relative to the corresponding MdAPE values suggest that a small number of countries with highly volatile transmission dynamics or reporting anomalies inflated the mean percentage error, whereas the typical country level forecast error remained considerably lower. Taken together, these four evaluation metrics demonstrate the robustness, stability, and reliability of the proposed BH STAR forecasting models across diverse national surveillance datasets.

4 Discussion

This study has utilized a hierarchical Bayesian modeling approach for smooth transition autoregressive (BH-STAR) models, encompassing logistic (LSTAR) and exponential (ESTAR) specifications, and applied it to COVID-19 surveillance data from 53 African countries. Our results demonstrate that BH-STAR models capture non-linear epidemic dynamics, accommodate country-level heterogeneity, and yield more accurate forecasts than conventional linear autoregressive models. The findings have several methodological and practical implications.

The BH-STAR as implemented in this study extends existing non-linear time series approaches by combining regime-switching flexibility with hierarchical pooling. While smooth transition models have been successfully applied in economics and finance [22, 42], their use in infectious disease epidemiology remains limited. Previous applications of STAR models in public health often assumed homogeneous dynamics across regions or relied on maximum likelihood estimation, which can be unstable when data are sparse [29, 48]. By adopting a Bayesian hierarchical structure, we enable information sharing across countries, stabilizing estimates for regions with limited reporting and providing coherent uncertainty quantification through posterior credible intervals. This approach aligns with recent calls for more robust statistical methods in global health surveillance [37, 49].

The ability to model gradual transitions between regimes is particularly relevant for infectious diseases. Epidemic waves often evolve over weeks or months, influenced by factors such as population immunity, behavior change, and public health interventions. The LSTAR specification captures asymmetric transitions (e.g., rapid rise, slower decline), while the ESTAR specification allows for symmetric smooth shifts. Our posterior estimates of the smoothness parameter γ (mean 6.87) indicate moderately rapid transitions, consistent with the wave-like patterns observed in COVID-19 data. Threshold parameters θ locate the midpoint of the transition, offering interpretable benchmarks for when disease incidence crosses from low to high transmission regimes.

Linear autoregressive models (AR and HAR) remain popular due to their simplicity and interpretability [8, 10]. We acknowledge that these models offer practical advantages, particularly when computational resources are limited or when a parsimonious description of the data suffices. However, our validation metrics (Table 3) show that BH-STAR models substantially reduce forecast errors, with MAE improvements of over 50% in some countries. This confirms that non-linearities are not negligible and that ignoring regime shifts leads to systematic under- or over-prediction during epidemic peaks and troughs. Other non-linear time series models, such as threshold autoregressive (TAR) and Markov-switching autoregressive (MSAR) models, could also capture regime changes [16, 17]. TAR models impose an abrupt switch at a threshold, which may be less realistic for gradual epidemic transitions, while MSAR models assume unobserved states governed by a Markov chain, which adds complexity in state estimation. STAR models offer a middle ground by allowing smooth, data-driven transitions. A formal comparison of BH-STAR with TAR and MSAR within the same hierarchical Bayesian framework is a promising avenue for future research.

Heteroscedasticity and volatility clustering, often observed in infectious disease data [25, 26], were not explicitly modeled here. However, the hierarchical error structure partially accounts for heteroscedasticity across countries and over time. Extensions such as smooth transition GARCH (ST-GARCH) models [29] could be integrated into the hierarchical framework to jointly model mean and variance dynamics, which may further improve forecast accuracy, especially during periods of high uncertainty.

Our proposed approach has the added strength of handling data irregularities by being robust to sparse and heterogeneous reporting, which is a common challenge in low-resource settings [2]. By treating missing observations as latent variables and using hierarchical pooling, we reduce the influence of country-specific data gaps. This is particularly important in the African context, where surveillance capacity varied widely during the pandemic [1, 50]. Our results show that even countries with intermittent reporting benefit from the hierarchical structure, as their parameter estimates are partially informed by continental trends. The surveillance environment underlying our analysis reveals critical challenges for epidemic modeling in Africa. Testing volume declined by 44% between 2021 and 2023, while average weekly tests fell by 72%, and test positivity remained at 20% in 2023, indicating that observed case counts increasingly underestimate true transmission over time [51]. Furthermore, the number of consistently reporting countries dropped from 49 to 39, and reporting heterogeneity was substantial: Tanzania reported a 400% case increase with only 6 of 8 reporting weeks, while Senegal reported a 300% increase with 7 of 8 weeks complete. These patterns imply that sharp increases in cases may reflect either genuine resurgences or artifactual fluctuations due to changes in reporting intensity. Such fragmented surveillance and limited diagnostic capacity are well-documented challenges across sub-Saharan Africa, where weak laboratory infrastructure, workforce shortages, and inconsistent data reporting create critical blind spots in epidemic intelligence [52].

Our Bayesian hierarchical approach addresses these issues by partially pooling across countries, which stabilizes estimates for nations with sparse or intermittent data, and by latent-variable imputation of missing months, which propagates reporting uncertainty into posterior credible intervals. However, the systematic decline in reporting coverage and the persistent under-ascertainment of cases mean that our forecasts are likely conservative, and the true epidemic burden is probably higher than point estimates suggest. These findings underscore that while BH-STAR models provide a principled framework for uncertainty-aware forecasting under weak surveillance, they cannot substitute for sustained investment in routine data collection and standardized reporting protocols.

Nonetheless, we did not explicitly model reporting bias or changes in testing intensity. Declining testing volumes (Table 1) likely introduced time-varying under-ascertainment, which may affect the observed case counts. Future work could extend the BH-STAR framework to a state-space formulation with a latent reporting process, similar to approaches used in nowcasting and real-time surveillance [53–55]. Such an extension would allow simultaneous estimation of true incidence and reporting probability, improving forecast reliability during periods of surveillance degradation.

4.1 Public health implications

Accurate and uncertainty-aware forecasts are essential for epidemic preparedness and response. The BH-STAR models provide probabilistic predictions with credible intervals that quantify forecast uncertainty, enabling decision-makers to assess risks and allocate resources accordingly. For example, forecasts that correctly anticipate epidemic peaks can inform the timing of interventions, while early warnings of resurgence can trigger enhanced surveillance or vaccination campaigns. Our empirical application demonstrates that BH-STAR models capture the timing and magnitude of COVID-19 waves across Africa, offering a practical tool for public health agencies with limited analytical capacity.

The hierarchical structure also facilitates cross-country learning. When data from a country are scarce, the model borrows strength from other countries with similar epidemic patterns. This is particularly valuable in regions where surveillance is intermittent or recently established. By sharing information, the BH-STAR framework can provide more reliable forecasts than country-specific models, supporting coordinated regional responses.

4.2 Limitations and future directions

Several limitations warrant consideration. First, our analysis used monthly aggregated case counts, which may mask finer temporal dynamics. Higher-frequency data (e.g., weekly or daily) could reveal more nuanced transition patterns, but such data were not consistently available across all countries. Second, the validation window of three months (October–December, 2022) was relatively short; longer-term validation would be desirable to assess model performance across different phases of the pandemic. Third, while we accounted for country-level heterogeneity, we did not include covariates such as mobility, policy interventions, or climatic factors. Incorporating such covariates could improve both interpretability and predictive accuracy [4, 54].

A further limitation is that we did not explicitly model heteroscedasticity or volatility clustering, phenomena often observed in infectious disease data [25, 26]. Although the hierarchical error structure partially accounts for heteroscedasticity across countries and over time, a more formal treatment could improve forecast accuracy during periods of high uncertainty. Future work could integrate smooth transition GARCH (ST-GARCH) models [29] into the hierarchical framework to jointly model mean and variance dynamics, providing a more comprehensive characterization of epidemic volatility.

Future research directions include extending the BH-STAR framework to accommodate multivariate outcomes (e.g., cases and deaths jointly), incorporating time-varying covariates, and developing efficient computational algorithms for large-scale applications. Additionally, applying the framework to other infectious diseases (e.g., malaria, tuberculosis, influenza) could assess its generalizability beyond COVID-19.

4.3 Conclusion

In summary, we have developed a hierarchical Bayesian smooth transition autoregressive model that effectively captures non-linear dynamics, regime shifts, and country-level heterogeneity in infectious disease data. Through simulation studies and empirical application to COVID-19 in Africa, we demonstrated superior forecasting performance over linear benchmarks. The BH-STAR framework offers a flexible, interpretable, and uncertainty-aware tool for epidemic forecasting, with particular value in settings where data are sparse and heterogeneous. By enabling cross-country information sharing and probabilistic predictions, it supports evidence-based public health decision-making in resource-limited environments.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Emmanuel Jesuyon Dansu, Tohoku University, Japan

Reviewed by: Anamaria Popescu, Universitatea Din Petrosani, Romania

Mohammad Baqar, Cisco Systems, United States