Formulation of a spatiotemporal model for the analysis of neonatal mortality amidst SDG interventions: The case of Uganda

TL;DR

This study creates a spatiotemporal model to analyze neonatal mortality in Uganda during SDG interventions.

Contribution

The paper introduces a Bayesian dynamic linear model for spatiotemporal analysis of neonatal mortality.

Findings

The model was applied to simulated data over 14 years to assess neonatal mortality trends.

Bayesian techniques like Kalman filtering and Gibbs sampling provided reliable forecasts.

The model helps evaluate the impact of SDG interventions on neonatal mortality rates.

Abstract

This study aimed to formulate a dynamic linear model within a Bayesian framework to conduct a spatiotemporal analysis of neonatal mortality in Uganda during SDG interventions. This study formulated a model based on appropriate health-related covariates while considering the spatial and temporal dimensions of the data whose variable of interest (dependent variable) was a quantitative variable measuring the monthly rates of neonatal mortality (number of newborns dying within their first 28 days of life) at the district level. Through Markov chain Monte Carlo (MCMC) simulations, the applicability of the model could be assessed using simulated data covering 14 years, starting in January 2010, to evaluate the situation before and after the implementation of interventions to achieve the SDGs targets. Using a Bayesian approach through the Kalman filtering technique, the parameters of the formulated model were estimated. This study used the same technique through Gibbs sampling to extract meaningful information from the simulated data and provide reliable forecasts for the rates of neonatal mortality.

Figures

Methodological flowchart summarizing the data design, Bayesian model specification & parameter estimation, and validation process.

Model diagnostics using trace and density plots for posterior distributions showing MCMC convergence of model parameters: (in Fig 2(2A)–2(2G) ((2A)) Trace and density plots for beta_cost, beta_healthfacdensity, beta_hcareaccess, beta_hcenteraccess, beta_maternaleduc and beta_hholdinc. (Set 1).((2B)) Trace and density plots for beta_regioncentral, beta_regionwestern, beta_regioneastern, beta_regionnorthern, beta_seasonDry and beta_seasonWet. (Set 2). ((2C)) Trace and density plots for beta_SDGintroAfter, beta_SDGintroBefore, beta_RMNCAH, beta_NHSDP, beta_ENAP, and beta_SMGL. (Set 3). ((2D)) Trace and density plots for beta_UNNSC, beta_MPDSR, beta_QUINH, beta_RMNCAH_NHSDP, beta_RMNCAH_SMGL and beta_RMNCAH_UNNSC. (Set 4). ((2E)) Trace and density plots for beta_RMNCAH_MPDSR, beta_RMNCAH_QUINH, beta_NHSDP_SMGL, beta_NHSDP_UNNSC, beta_NHSDP_MPDSR and beta_NHSDP_QUINH. (Set 5). ((2F)) Trace and density plots for beta_ENAP_SMGL, beta_ENAP_UNNSC, beta_ENAP_MPDSR, beta_ENAP_QUINH, beta_SMGL_UNNSC and beta_SMGL_MPDSR. (Set 6). ((2G)) Trace and density plots for beta_SMGL_QUINH, beta_UNNSC_MPDSR, beta_UNNSC_QUINH and beta_MPDSR_QUINH. (Set 7).

Equations

Full text

Introduction

The burden of neonatal mortality has continued to increase in most countries [1]. A healthy economy with a minimal mortality rate is desirable for every country worldwide. Neonatal (newborn) mortality is defined for this study as death within the first 28 days of life per 1000 live births as established by the World Health Organization [2]. As noted by [1], reducing neonatal mortality is an essential part of SDG 3, Section 3.2.2: countries should aim to reduce neonatal mortality to at least 12 per 1,000 live births and under five mortality to at least 25 per 1,000 live births by 2030 and achieving this requires an understanding of the levels of and trends in neonatal mortality. According to a study by [3], 5.2 million children died before reaching their fifth birthday in 2019, with almost half of those deaths, 2.4 million occurring in the first month of life, despite the countries’ efforts to reduce death rates, for instance, UN member states’ interventions in terms of the MDGs and SDGs.

In Uganda and globally, most studies on neonatal mortality, such as [4–6], have focused on risk factors while ignoring progress and forecasting the situation given that the targeted 2030 to achieve the SDGs is fast approaching. Moreover, these studies did not consider the influence of time or space. This study aimed to formulate a dynamic linear model to analyze neonatal mortality using a case study of Uganda to simultaneously investigate its persistent patterns over time and space and illuminate any unusual patterns. As mentioned by [7], the majority of time series models include well-known models such as autoregressive (AR) and moving average (MA) models. The autoregressive integrated moving average (ARIMA) and autoregressive moving average (ARMA) are helpful for handling stationary data; otherwise, they become limited. This is further supported by [8–12], as they suggest that state space models offer a very rich class of models that have several advantages. For example, they do not require stationarity, which eliminates the need to transform the data since data transformation leads to the loss of some important components in the data, which at times leads to less accurate results. To formulate a dynamic linear model, this study simulated data for health care-related factors and health care policies that were put in place as a way of achieving SDG 3.2. The health policies used are shown in Table 1.

Table 1: Policies in place to reduce neonatal mortality in Uganda.

Materials and methods

Study design and overview

In this study, a simulation-based design was used to illustrate the application of the formulated Bayesian dynamic linear model (BDLM) in evaluating spatiotemporal trends in neonatal mortality and the impact of health policy interventions aligned with the Sustainable Development Goals (SDGs) framework. The simulated data represent monthly neonatal mortality rates and related health indicators in Uganda from January 2010 to December 2023, reflecting the influence of key health care initiatives implemented during the SDGs period. The rationale for choosing this period was to cover the time before and after the implementation of interventions to achieve the SDGs and to assess their impact.

Methodological framework

To make the methodology presentation easier to follow, this study constructed a visual diagram, as shown in Fig 1.

Methodological flowchart summarizing the data design, Bayesian model specification & parameter estimation, and validation process.

Study variables

The dependent variable was neonatal status, which was measured as the monthly number of liveborn children who died within their first 28 days of life per district.

The independent variables (covariates) were as follows:

Cost was measured as the monthly amount of money the district receives from the Ministry of Finance through the Ministry of Health for the implementation of the seven policies obtained by dividing the annual amount by 12.Year was measured as the calendar year for which the data were recorded.Month was measured as the month of the calendar year for which the data were recorded.District was measured as the name of the district from which the data were recorded.Region was measured as the main administrative region in which the district was located since these data are reported monthly from districts but not regions. Therefore, districts were regrouped to obtain regions during data cleaning and editing. The regions were coded as 1 for Central, 2 for Western, 3 for Eastern and 4 for Northern.HealthfacDensity was measured as the number of health centers in the district.HcareAccess was measured as the number of government health care providers in the district.HcenterAccess was measured as the distance in meters to the nearest health center from the farthest household in the district.MaternalEduc was measured as the average number of antenatal visits by women in the district per month.HholdInc was measured as the average monthly household income in the district.Seasons were measured as the changes in weather conditions in the country captured as (1 Dry, 0 Wet).SDGintro was measured as the time of the start of the SDGs during the observation period was captured as (1 After, 0 Before).HcarepolicyInfo was measured as the time of the start of the 7 health care policies shown in Table 1 during the observation period captured as (1 After, 0 Before).Observ was measured as the observation time running from 1 in January 2010 to December 2023.

Data simulation procedures and availability

A synthetic time series was generated through simulations to emulate realistic variations in neonatal mortality associated with health care system performance and policy interventions. A simulated dataset comprising 1,008 observations, incorporating health system indicators such as antenatal care coverage, facility access, and skilled birth attendance, along with temporal components that capture secular trends and seasonal fluctuations, was used. This reflects Uganda’s neonatal mortality context before and during SDG implementation, using baseline mortality of 27 deaths per 1,000 live births based on a report by the Ministry of Health, Uganda, as reported in [15]. We also based on Sample data from the DHIS2 (District Health Information System 2) since this was reported as a better reporting tool by [17] and the UN Inter-agency Group for Child Mortality Estimation (UN IGME) [18]. The simulation was performed in R (version 4.3.2) to ensure reproducibility and control over model parameters, thereby allowing transparent evaluation of methodological performance under known data-generating conditions. The simulated dataset, along with the R script used for data generation, model estimation, validation, and figure production, is publicly accessible, as it is simulated.

Bayesian inference and prior specification.

Bayesian inference was performed using Markov chain Monte Carlo (MCMC) methods implemented in JAGS (Just Another Gibbs Sampler). The following prior distributions were specified:

Regression coefficients:

$$[eqn]$$

for the continuous variables and

$$[eqn]$$

for the categorical variables. Observation variance:

$$[eqn]$$

State variance:

$$[eqn]$$

Initial state:

$$[eqn]$$

These priors are weakly informative, chosen for their ability to identify the model without influencing posterior estimates. The scale of the priors was informed by prior neonatal mortality studies and time-series modeling literature. Four independent MCMC chains were run, each with 1,000 burn-in iterations followed by 5,000 sampling iterations. Convergence was assessed through trace plots, posterior density plots, and the Gelman–Rubin diagnostic, with $[eqn]$ indicating satisfactory convergence. Posterior means and 95% credible intervals (CrIs) were reported for all the model parameters.

Model implementation

Nature of and rationale for the model.

The study opted to formulate a dynamic linear model to monitor the impact of the SDG interventions on neonatal mortality, as dynamic linear models constitute a flexible class of models that can effectively address typical features in policy monitoring and large datasets.

For instance, this study aimed to construct a dynamic linear model to assess the situation after the onset of the SDGs, having previously evaluated the situation. To achieve this objective, this study used a general Gaussian dynamic linear model, as presented in Eq 1:

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ represents the observation for the $[eqn]$ state at time t.

$[eqn]$ denotes a function of the x covariates for the $[eqn]$ state at time t specified as

$$[eqn]$$

$[eqn]$ is a vector of time-carrying parameters

$[eqn]$ is a transition matrix depicting how the model parameters evolve over time.

$[eqn]$ and $[eqn]$ are the observation and evolution errors, respectively.

On the basis of (1), the model is formulated as follows (3):

$$[eqn]$$ $$[eqn]$$

Where:

$[eqn]$ represents the $[eqn]$ district neonatal mortality rate at time t.

$[eqn]$ are the m covariates, including the maternal health care policies and interventions, in which case a code 0 is used for Before and 1 for After the introduction of a given policy/intervention during the observation time for a given policy or intervention.

$[eqn]$ is the interaction covariate between the $[eqn]$ covariate, including the maternal health care policies and interventions, with $[eqn]$ as their respective coefficients,

$[eqn]$ denotes the time-varying intercept.

$[eqn]$ are the time-varying coefficients associated with each of the m covariates, including maternal health care policies and interventions.

$[eqn]$ and $[eqn]$ are the observation and evolution error terms, respectively, where $[eqn]$ and $[eqn]$ .

Model parameter estimation.

A Bayesian approach using the Kalman filtering technique was used for estimating the study model parameters because the Bayesian approach focuses on how a particular part of the data depends on the other data parts since this study had to compare the neonatal mortality situation in the country before and after the introduction of the SDGs, implying that past values had a considerable influence on future values in this study.

Considering that Bayesian inference has its roots in Bayes rule, as noted in the studies by [13] and [14], to obtain the final parameter estimates, this study made use of the theorem through the posterior density function as expressed in (4);

$$[eqn]$$

The final parameter estimates are obtained using a quadratic mean loss function (squared error loss) method, which involves obtaining a full posterior distribution of parameters and then calculating their average as the final parameter estimate.

Rewriting the formulated model by the study presented in (3), which is expressed in the form of two equations, observation and state equations, we obtain (5);

$$[eqn]$$ $$[eqn]$$

where

$[eqn]$ is the observed data $[eqn]$ is a p-dimensional design vector of covariates and the interaction terms.

$[eqn]$ is a p-dimensional time-varying parameter vector (for both covariates and the interaction terms).

$[eqn]$ and $[eqn]$ are the observation and evolution error terms, respectively, where $[eqn]$ and $[eqn]$ .

Since the formulated model is a multivariate Gaussian state space model, based on the standard results of a multivariate Gaussian distribution pertaining to its marginal and conditional distributions, the random vector of parameters $[eqn]$ also has a Gaussian distribution, and the same applies to the other respective marginal and conditional distributions.

Therefore, since all the relevant distributions are Gaussian, and bearing in mind that the process of Bayesian inference involves passing from a prior distribution, $[eqn]$ , to a posterior distribution, $[eqn]$ , we estimated the model parameters through the Kalman filtering technique to obtain the means and variances to finally derive the parameter estimates (posterior means) through the following procedure:

Obtaining the prior distribution (based on the evolution equation reflected in the second line of (3)).

Considering the fact that a dynamic linear model is a Gaussian time space model, we used a normal prior distribution for a p-dimensional state vector by first obtaining the initial values where $[eqn]$ , and we then attained one step ahead of the prediction of the parameter vector distribution $[eqn]$ by first obtaining the mean as given in (6)

$$[eqn]$$

Then, we obtained the predicted covariance in (7) denoted as $[eqn]$

$$[eqn]$$

Therefore, $[eqn]$ which is the prior distribution.

Since our prior has been derived on the basis of a normal distribution, we can write its probability density function as given in (8):

$$[eqn]$$

iiObtaining the likelihood distribution function (based on the observation equation reflected in the first line of (3)).

We start this step by obtaining the predicted (one step) distribution function of $[eqn]$ (the observed value): $[eqn]$

Let $[eqn]$

Nevertheless, we first obtained the mean, $[eqn]$

$$[eqn]$$

We then obtained the variance, $[eqn]$

$$[eqn]$$

Therefore, $[eqn]$ , which is the observed value likelihood distribution.

Since the observed value variable follows a normal distribution, we can write its likelihood function probability density function as follows;

$$[eqn]$$

iiiObtaining the posterior distribution

From (4), an expression derived from the Bayes theorem, we can obtain the posterior probability density function by multiplying (11) by (8) and hence obtaining the probability density function of the posterior distribution function as given in (12);

$$[eqn]$$

Since both the prior and likelihood distributions are Gaussian, the resulting posterior is also Gaussian. If we let the posterior $[eqn]$ , from the first property of a Gaussian distribution pertaining to the probability density function, we can obtain the values of $[eqn]$ by simplifying the right-hand side of (12);

$$[eqn]$$

From (13), since the posterior $[eqn]$ , it implies that;

$$[eqn]$$

Therefore, $[eqn]$

From this expression, we used a quadratic mean loss function (squared error loss), which involves computing the full posterior distribution of the parameters and averaging them to obtain the final parameter estimate.

Model validation.

To validate the model, we conducted residual tests on autocorrelation, heteroscedasticity, and normality using the Durbin–Watson test, the Breusch Pagan test, and the Shapiro–Wilk test, respectively. The following properties were checked for:

Autocorrelation: The residuals should not be autocorrelated.Homoscedasciticity: The residuals should have a constant variance.Normality: The residuals are expected to be normally distributed.

To verify the above model properties, tests were carried out using the following hypotheses:

$[eqn]$ : Residuals are not autocorrelated vs. $[eqn]$ : Residuals are autocorrelated. $[eqn]$ : Residuals have constant variance Vs $[eqn]$ : Residuals do not have constant variance. $[eqn]$ : Residuals are normally distributed $[eqn]$ vs. $[eqn]$ : Residuals are not normally distributed.

Test statistics, critical values, and the decision rule.

While testing for autocorrelation of residuals, the Durbin–Watson test was used. To obtain its test statistic, the following formula was used:

$$[eqn]$$

where

$[eqn]$ = the residuals from the regression model.T = the number of observations.

With respect to the critical value and the decision rule, we assume that the DW (value) ranges between 0 and 4; since a DW value approaching 2 indicates no autocorrelation, a value approaching 0 indicates positive autocorrelation, and a value approaching 4 indicates negative autocorrelation, taking into consideration the probability value of the test statistic in comparison to the level of significance.

The Breusch Pagan test was used to assess the homoscedasticity of residuals, whereby to obtain its test statistic, the following formula was used:

$$[eqn]$$

where

n = number of observations

$[eqn]$ = R-squared value obtained by regressing the squared residuals from the original regression on the independent variables, their squares, and their cross products. This regression is known as auxiliary regression.

With respect to the critical value and the decision rule, because the test statistic of Breusch Pagan follows a chi-square distribution with degrees of freedom equal to the number of independent variables in the auxiliary regression, the test statistic value was compared with the critical (tabulated) value to ascertain whether to reject the null hypothesis basically considering the significance level and the reported probability value of the test statistic.

The Shapiro–Wilk test was used to determine if the residuals were normally distributed, whereby the following formula was used to obtain its test statistic, denoted as W:

$$[eqn]$$

where

$[eqn]$ = ordered residuals from smallest to largest.

$[eqn]$ = mean of the residuals.

$[eqn]$ = constants derived from the expected values of the order statistics of a standard normal distribution computed using statistical software.

n = number of observations.

The critical value and the decision rule were determined on the basis that W ranges between 0 and 1; values of W close to 1 suggest that the residuals are approximately normally distributed, whereas lower values indicate deviations from normality.

Furthermore, The performance of the Bayesian dynamic linear model formulated in this study has also been evaluated in our separate, independently published study, which compares it with mixed-effects interrupted time-series models using real-world neonatal mortality data from Uganda (2015–2023). Whereas the present work emphasizes model formulation and Monte Carlo simulation-based validation, the other study focuses on empirical comparison [16].

Model sensitivity analysis

To assess model fitness, a residual analysis was performed. Having estimated posterior means and credible intervals for the model parameters, trace and means plots were constructed, which enabled convergence to be assessed. To evaluate predictive performance, the model was extended to generate forecasts of neonatal mortality for the first 10 months after December 2023, aligned with the national health planning cycles. In addition to using trace plots as a way of assessing model convergence visually and numerically, convergence was assessed using the Gelman–Rubin diagnostic ( $[eqn]$ ), with values less than 1.1 indicating satisfactory convergence.

Results and discussion

Model validation

The following conclusions have been made regarding model validation on the basis of the results in Table 2:

Table 2: Models validation tests results.

While testing for autocorrelation of residuals with the Durbin-Watson test, though the Durbin-Watson value (DW = 0.07253) ranges between 0 and 2, showing no autocorrelation of the residuals, its p-value ≤ 0.05, indicating rejection of the null hypothesis, which implies residuals are autocorrelated.

While testing for homoscedasticity with the Breusch-Pagan test, because the p-value (0.3755) is greater than 0.05, we fail to reject the null and conclude that the residuals have constant variance.

While testing whether the residuals are normally distributed with the Shapiro–Wilk test, based on the fact that Shapiro–Wilk’s test value W ranges between 0 and 1, since W (W = 0.9963) is close to or equal to 1, this suggests that the residuals are normally distributed but its p-value is less than 0.05 which implies rejection of the null hypothesis and hence concluding that residuals aren’t normally distributed.

The model fails to fulfill the assumptions of non-autocorrelation and normality of residuals based on the probability values, but this does not invalidate our model, but instead gives justification for its use based on the fact that unlike classical linear regression, Bayesian inference does not rely on independence or normality assumptions for valid parameter estimation. Instead, such diagnostics highlight temporal dependence and unobserved heterogeneity in the data, which motivate the use of a Bayesian framework capable of explicitly modeling complex data-generating processes and propagating uncertainty through the posterior distribution [19–22]

Findings from [16] show that the Bayesian Dynamic Linear Model outshone Mixed Effects Interrupted Time Series Model as a preferred model for Neonatal mortality analysis amidst SDGs onset impact evolution

Model convergence

As shown in Fig 2 depicted in Fig 2(2A)–2(2G), the trace and density plots indicate convergence across almost all parameters except one MaternalEduc, measuring the number of women doing their first antenatal visit per month in the district, whose Gelman–Rubin diagnostic ( $[eqn]$ ) statistic (1.37) is greater than the threshold value of 1.1, as further reported by the Gelman–Rubin diagnostic ( $[eqn]$ ) results in Table 3. This implies that the MCMC algorithm successfully captured the true posterior distribution of our model parameters, meaning that our parameter estimates, uncertainties, and forecasts are trustworthy. We have also included in the Supporting Information a file containing the graphical illustration of the convergence test obtained from Table 3.

Table 3: Gelman–Rubin diagnostic (R^) and 95% upper credible intervals for model parameters.

Model diagnostics using trace and density plots for posterior distributions showing MCMC convergence of model parameters: (in Fig 2(2A)–2(2G) ((2A)) Trace and density plots for beta_cost, beta_healthfacdensity, beta_hcareaccess, beta_hcenteraccess, beta_maternaleduc and beta_hholdinc. (Set 1).((2B)) Trace and density plots for beta_regioncentral, beta_regionwestern, beta_regioneastern, beta_regionnorthern, beta_seasonDry and beta_seasonWet. (Set 2). ((2C)) Trace and density plots for beta_SDGintroAfter, beta_SDGintroBefore, beta_RMNCAH, beta_NHSDP, beta_ENAP, and beta_SMGL. (Set 3). ((2D)) Trace and density plots for beta_UNNSC, beta_MPDSR, beta_QUINH, beta_RMNCAH_NHSDP, beta_RMNCAH_SMGL and beta_RMNCAH_UNNSC. (Set 4). ((2E)) Trace and density plots for beta_RMNCAH_MPDSR, beta_RMNCAH_QUINH, beta_NHSDP_SMGL, beta_NHSDP_UNNSC, beta_NHSDP_MPDSR and beta_NHSDP_QUINH. (Set 5). ((2F)) Trace and density plots for beta_ENAP_SMGL, beta_ENAP_UNNSC, beta_ENAP_MPDSR, beta_ENAP_QUINH, beta_SMGL_UNNSC and beta_SMGL_MPDSR. (Set 6). ((2G)) Trace and density plots for beta_SMGL_QUINH, beta_UNNSC_MPDSR, beta_UNNSC_QUINH and beta_MPDSR_QUINH. (Set 7).

Parameter estimates

From the means (parameter estimates) and quantiles of the parameters reported in Tables 4 and 5, most of the variables have a positive influence on the dependent variable (neonatal mortality) with the exemption of health center access and household income as reported under means whereas regarding the quantiles, most of the variables significantly influence the outcome variables since their extreme quantiles are nonzero. On the basis of the simulated data, none of the health care policies were found to interact with one another in terms of the influence of neonatal mortality.

Table 4: Empirical mean and standard deviation for each variable plus the standard error of the mean for model parameters estimated using the Bayesian dynamic linear model (BDLM).

Table 5: Quantiles for each variable.

Forecasts

Using the formulated study model, we were able to make a 10-month prediction of Uganda’s neonatal mortality rates, as shown in Table 6. The study predicted that the neonatal mortality rate of Uganda would stand at 7.7 deaths on average in the first 10 months after December 2023.

Table 6: Monthly forecasts of neonatal mortality for Uganda for January–October 2024 using the formulated Bayesian dynamic linear model.

Conclusion

A dynamic linear model that can be used to monitor health care policies individually and assess the possibility of interactions between any two policies such that the wastage of resources on duplicated policies is avoided was developed in this study. It can be used to generally assess the health care situation in a country.

Using the formulated model, the study was able to make a 10-month prediction of neonatal mortality, hence demonstrating the applicability of the model.

Supporting information

S1 FileR code script. The R script used for data processing and analysis is provided as a supplementary file named the S1 File.R.(R)

S2 FileSimulated dataset. An Ms. Excel-simulated dataset used for analysis is provided as a supplementary file named S2 File.xlsx.(XLSX)

S6 FileGelman’s Rubin’s diagnosis. This is the Rubin diagnosis result. This file contains graphical plots representing the Gelman–Rubin diagnostic ( $[eqn]$ ) and 95% upper credible intervals for the model parameters.(PDF)