A Bayesian spline-augmented piecewise exponential model with spatial frailty for under-five mortality in Nigeria

Introduction The Cox proportional hazards (PH) model is widely used in time-to-event research, but its validity depends on the PH assumption, which can be violated in child mortality studies where hazards vary with age. Piecewise exponential models (PEMs) relax this assumption by partitioning follow-up time into intervals. However, standard formulations impose constant hazards within each interval and typically ignore spatial heterogeneity, limiting their usefulness for public health analyses in which geographical variation in risk is important. Existing extensions address non-proportional hazards or spatial dependence separately, but rarely combine a smooth baseline hazard, interval-level temporal heterogeneity, and structured spatial frailty within a single, computationally tractable Bayesian framework. This study proposes and applies such a framework to under-five mortality (U5M) in Nigeria. Methods We formulated a Bayesian spline-augmented PEM incorporating a cubic-spline-smoothed baseline hazard, interval-specific Gaussian random effects, and a spatial frailty term with an intrinsic conditional autoregressive (ICAR) prior, within the latent Gaussian modeling framework of integrated nested Laplace approximation (INLA), with survival reformulated as a Poisson likelihood. The model was applied to U5M data for 103,439 children from the 2024 Nigeria Demographic and Health Survey (NDHS). Five nested model specifications (basic PEM, spline-augmented, spline plus interval effects, spline plus spatial effects, and the full spline-interval-spatial model) were compared across three sample sizes using DIC, WAIC, conditional predictive ordinates, posterior predictive checks, and calibration metrics. Results The global test of the PH assumption was significant (χ 2  = 898.66, p  < 2 x 10 -16 ), with all individual covariates significant at p < 0.05. Across all three sample sizes, the full spline-interval-spatial model achieved the lowest DIC and WAIC, the lowest CPO, and the highest posterior predictive correlation and calibration R 2 among the five candidates, with the largest single reduction in DIC (approximately 17,884 units at n = 104,557) attributable to the addition of interval-specific random effects. In the fitted model, twin birth (HR = 2.82; 95% CrI: 2.62–3.03), breastfeeding status (HR = 0.45; 95% CrI: 0.42–0.49), and term delivery (HR = 0.47; 95% CrI: 0.41–0.53) had the largest effects on U5M, with additional protective effects for longer birth intervals and higher maternal education, and elevated risk in the North West (HR = 1.39) and North East (HR = 1.43) relative to the North Central region. Posterior spatial frailty estimates showed positive residual clustering concentrated in the North West and North East zones. Discussion The results indicate that jointly modeling temporal and spatial heterogeneity yields the best-fitting and best-calibrated specification, although spatial frailty adds comparatively little once interval-level temporal effects are included, suggesting that residual heterogeneity in this setting is predominantly temporal rather than spatial. The Poisson-INLA formulation provides a computationally efficient alternative to MCMC-based spatial survival models, making it well suited to large-scale demographic surveys. From a public health perspective, the identified biological, maternal, and socioeconomic determinants point to the need for integrated interventions, while the persistence of unexplained spatial clustering after covariate adjustment indicates structural or contextual vulnerabilities in the northern zones in Nigeria.