Caesarean section delivery & determining factors in northeastern State, India: A Bayesian spatio-temporal analysis of areal data

Data

This study used cross-sectional data from three consecutive rounds of National Family Health Survey (NFHS 2005-06, 2015-16 and 2019-21)^7,8, and age-referenced shapefile. The NFHS was conducted from 2005-06, 2015-16, and 2019-21, under the stewardship of the Ministry of Health and Family Welfare, Government of India.

This study considered all mothers who delivered in the last five years of the respective survey in the northeastern States of Arunachal Pradesh, Assam, Manipur, Meghalaya, Mizoram, Nagaland, Sikkim, and Tripura. We used NFHS-3, NFHS-4, and NFHS-5 datasets for information on the number of CS deliveries from 2001 to 2005, 2011 to 2015, and 2016 to 2019, respectively. The unit of analysis is the most recent birth in the last five years preceding the survey.

Outcome variable

The outcome variable is the reported number of CS deliveries in each district of the northeastern States of India.

Predictor variables

Based on the review of existing literature, we identified the following predictor variables and maternal characteristics: maternal age (15-29 and 30-49 yr), educational level (at most secondary level education, higher education), mothers’ BMI (<30 kg/m², ≥30 kg/m²), alcohol consumption (yes and no), use of smokeless tobacco (yes and no); child characteristics: child’s birth order (<3, ≥3); household characteristics: religion (Hindu and others), social groups (unreserved and others), wealth quintile of the household (poor/middle, and rich); and community characteristics: place of residence (rural and urban). Additionally, the study also considered the type of healthcare facility for delivery (private and public) and the distance to the health facility (no problem and problem)^3,4,7.

Statistical analysis

We estimated the district-wise prevalence of CS deliveries in 15-49 yr old women from 2011 to 2019, across 104 districts in the northeastern States. The svyset command in STATA was used to take into account the complex survey sampling design. Further, we computed the global Moran’s I statistics (0.473, P<0.001) to test for spatial autocorrelation.

We used SIR to compare the observed number of occurrences of an event in a population relative to what might be the expected number of occurrences.

We applied Poisson log-linear model, represented as follows: $y_{kt}\sim Poisson\left(E_{kt}{\theta }_{kt}\right)$ $fork\left(area\right)=1,\dots .,104,t\left(time\right)=2011,\dots ,2019$

Risk of CS in the area k with time t is denoted by, $\ln \left({\theta }_{kt}\right)={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+{\beta }_4{x}_4$

where, $y_{kt}$ is the number of CS in local authority area k during year t, ${\theta }_{kt}$ is the RR of CS in local authority area k during year t and $E_{kt}$ is the expected count. In addition, $x_1$ , $x_2$ , $x_3,\text{and}{x}_4$ denote the proportion of mothers with higher education, between ages 30-49 yr, and those with obesity and within the richest wealth quintile, respectively.

We employed Bayesian spatio-temporal modelling using MCMC simulation that deals with spatial and spatio-temporal data accounting for both space and time variability in CS risk modelling⁸. In MCMC simulation, we proposed the general form: $\begin{array}{c}{y}_{kt}\sim Poisson\left(E_{kt}{\theta }_{kt}\right)\\ fork=1,\dots .,104,t=2011,\dots ,2019 \end{array}$ $\ln \left({\theta }_{kt}\right)={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+{\beta }_4{x}_4+{\psi }_{kt}$

where, $\text{log}\left({\theta }_{kt}\right)$ is expressed as a sum of several components, including spatial and temporal structures that take into account the fact that neighbouring areas and consecutive times may have more similar risks. The spatio-temporal interaction can also be included in dif ferent areas, which may have different time trends, but these may be more similar in neighbouring areas⁸.

Here, ${\psi }_{kt}$ is the spatio-temporal random effect for local authority area k and time-period t.

Our aim here is to quantify the evolution of the spatial pattern in CS risk over time, so we used the spatially autocorrelated first-order autoregressive process model given by, ${\psi }_t={\rho }_T{\psi }_{t-1}+{\epsilon }_t$

where, ${\psi }_t=({\psi }_{1t},\dots \dots \mathrm{..},{\psi }_{kt})$ denotes the vector of random effects for all areal units (local authority) at time t., and ${\epsilon }_t$ = ( ${\epsilon }_{1t}$ , ……, ${\epsilon }_{kt}$ ) the vector of errors. Temporal autocorrelation is controlled by the mean function ${\rho }_T{\psi }_{t-1}$ , while spatial autocorrelation is controlled by the covariance structure of ${\epsilon }_t$ . This model was used by Rushworth et al (2014)⁹, and is available in R as a package called CARBayesST.

For assessing the convergence of the Markov chains, we drew a trace plot of the samples for each regression parameters $\beta ={\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4$ . We provided estimated posterior median, 95% credible intervals, Geweke diagnostic value for all three chains in supplementary table I.

We computed estimated RR for four predictor variables by constructing a matrix of MCMC samples for the regression parameters ( ${\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4$ ) from all chains using the formula given in equation (3) for $x_1$ variable. Similarly, we computed for $x_2,$ $x_3$ and $x_4$ . $RR(x_1,\xi )=\frac{Riskofc\text{-}sectionif{x}_{1}increasedby\psi }{Riskofc\text{-}sectiongiventhecurrentvalueof{x}_1}$ $\begin{array}{c}=\frac{\exp ({\beta }_0+{\beta }_1(x_{1k}+\xi )+{\beta }_2{x}_{2kt}+{\beta }_3{x}_{3kt}+{\beta }_4{x}_{4kt}+{\psi }_{kt})}{\exp ({\beta }_0+{\beta }_1{x}_k+{\beta }_2{x}_{2kt}+{\beta }_3{x}_{3kt}+{\beta }_4{x}_{4kt}+{\psi }_{kt})}\\ =\exp \left({\beta }_i\xi \right)\end{array}$

We used STATA version 18 (Stata Corp, College Station, TX, USA) and R version 4.3.2 to fulfil the study objectives. Throughout the analysis, an α level of 0.05 was adopted to indicate statistical significance.

Trend analysis

The trends of CS rate in northeastern States and India from 2001-2019 are displayed in figure 1. The overall CS rates in the northeastern States increased by nearly 15 per cent over the last two decades. It increased from 5.4 per cent in 2001 to 19.9 per cent in 2019. At the Indian level, the CS rates increased from 6.8 per cent in 2001 to 23.7 per cent in 2019. Nearly all northeastern States observed an increasing trend in CS rates from 2001 to 2019.

Close

Fig. 1.

Trends of CS rate in northeastern States and India from 2001 to 2019.

Export to PPT

Local indicator spatial autocorrelation (LISA)

Based on the results of the LISA analysis, we observed a substantial increase in the number of high-high LISA clusters from 2011 to 2019, i.e., districts with a high prevalence of CS rate were surrounded by districts where higher CS rates were observed. In 2011 and 2012, Imphal West and Thoubal districts of Manipur had high-high LISA clusters (Fig. 2A and B). Further, in 2013, three districts of Manipur and Lakhimpur of Assam were in high-high LISA category (Fig. 2C). In 2014, the south district of Sikkim and Dhalai district of Tripura were also in high-high LISA category (Fig. 2D). Further, in 2015, south Tripura and Gomati districts of Tripura, Imphal West district of Manipur, south district of Sikkim, and Dibrugarh district of Assam exhibited high-high LISA clustering (Fig. 2E).

Close

Fig. 2.

LISA of CS delivery rate in northeastern States from 2011 to 2019. Maps were reproduced with permission from Georeference data received from DHS website.

Export to PPT

High-high LISA clusters in south and east districts of Sikkim, Sepahijala and south Tripura districts of Tripura, Imphal West district of Manipur, and Kamrup and Shivsagar districts of Assam were observed in 2016 (Fig. 2F). Similar trends of high-high LISA categories were observed in 2017 (Fig. 2G). In 2018, three districts (north, west, and south) of Sikkim, two districts of Assam (Golaghat and Shiv Sagar), Imphal West, Thoubal districts of Manipur, and Sepahijala and South Tripura district of Tripura were identified as belonging to the high-high LISA category (Fig. 2H). Moreover, three districts of Sikkim (north, south and west district), Sepahijala, West Tripura, and Gomati were from high-high LISA category in 2019 (Fig. 2I).

Figure 3 showed that a higher RR value was found in the upper-western districts of Arunachal Pradesh, lower and western part of Assam, and north district in Sikkim from 2015 to 2019 (a detailed interpretation of Fig. 3 is provided in supplementary material).

Close

Fig. 3.

RR for CS delivery rate in northeastern States from 2011 to 2019. Maps were reproduced with permission from Georeference data received from DHS website.

Export to PPT

Spatio-temporal analysis

Figure 4 indicated that the chains appear to have converged, as there is no change in the mean or variance of the samples between the three chains. To draw inferences from the MCMC model on the effect of the predictor variables on CS deliveries, the posterior median RR was estimated for a fixed increase in each of the covariates. Table shows that the estimated posterior median RR for four predictor variables was above 1, which indicated that the variables were significantly related to CS delivery risk.

Close

Fig. 4.

Trace plots of the regression parameters using MCMC samples from each chain.

Export to PPT

To estimate the average temporal trend, we first estimated the average (mean) risk across the k=104 districts for each year and MCMC sample, which yielded the posterior distribution of these spatial averages for each year. Then, we computed the posterior median RR surface ${\stackrel{}{\theta }}_{kt}$ and 95% credible interval for the spatially averaged risk. Figure 5 displays that the estimated spatio-temporal trend shows a clear upward trend in CS delivery risk over the nine-year study period.

Close

Fig. 5.

Posterior median relative risk surface ${\stackrel{}{\theta }}_{kt}$ and 95% credible interval for the spatio-temporal trend.

Export to PPT