Analysis of the spatial distribution of under-5 mortality rate in local areas of South Africa

Abstract

Policymakers and healthcare service managers demand reliable, accurate and disaggregated information about child deaths at the sub-national level to plan and monitor healthcare service delivery and health outcomes. In support of this demand, this research aimed at providing reliable local municipality estimates of the under-5 mortality rate (U5MR) in South Africa. The paper used a small area estimation approach to improve the precision of local municipality estimates of U5MR by linking data from the 2016 Community Survey (CS) and the 2011 Population Census (PC). The diagnostic measures and validation of the reliability of the results showed that the local municipality estimates of U5MR produced by small area estimation are more efficient and precise than direct estimates of U5MR based only on the CS data. Further, accurate and cost-effective local municipality estimates of U5MR were produced without the need for more resources through combining the available data sources. This was achievable since the research did not require a separate survey for this purpose. The results can be used to monitor U5MR at the local level in South Africa since they link directly with the Sustainable Development Goals (SDGs).

Keywords

Under-5 mortality small area estimation Fay-Herriot model community survey population census

1. Introduction

Mortality indicators, especially child mortality rate, is a representative demographics index and has been used as one of the main indicators of quality of life for a particular country or region [1]. Global reductions in under-5 mortality rate (U5MR), the probability that a child will die before reaching the age of five, has long served as an important proxy measure for overall improvements in population health and progress towards achievement of development goals [2].

There has been sustainable global progress in reducing child mortality since 1990 [2]. The global U5MR has declined from 19.6 million in 1950 to 5.4 million in 2017 [2]. Similar to the rest of the world, national official data indicates that South Africa has made substantial progress in reducing U5MR [3]. The number of child deaths has reduced substantially during the past decade, especially for U5MR in South Africa and which subsequently was 37–40 per 1 000 live births in 2015 [3]. One of the main goals of South Africa is reducing U5MR. Goal 3.2 aims to reduce U5MR to less than 25 per 1 000 live births by 2030 [3].

Nevertheless, the progress in reducing U5MR at the sub-national level (e.g., district, municipality, ward) remains unclear. In other words, there has been very limited effort in obtaining sub-national U5MR estimates in South Africa. This is unfortunate, as these estimates provide more localized spatial distribution [4]. For example, regarding economic history, one can learn new information from the spatial distribution of mortality rates in the sub-national administrative levels [4]. Sub-national estimates could be used as a basis for identifying marginalized households and communities to ensure the equitable and inclusive provision of assistance to the ones most in need. Furthermore, the U5MR in local municipalities is also an essential tool for monitoring progress towards the SDGs to ensure that certain districts or municipalities are not altogether left behind.

The aim of this paper is twofold. First, the primary aim is to account for inequality in U5MR between local municipalities in South Africa. In other words, we report the estimated U5MR in South Africa using a small area estimation method, underpinned by an approach that does not receive much attention in the mortality literature. Second, this paper analyzes geographical (spatial) variations of U5MR in local municipalities using spatial mapping. It is vital for decision-makers at all levels to better understand the spatial heterogeneity in child survival.

From the above aims, it is possible to summarize the following questions: (i) How does the inequality in U5MR vary across local municipalities? (ii) How is precision improved when model-based U5MR is used instead of direct survey estimates? (iii) What policy implications can be proposed based on this study? To explore these questions, this research relied on the 2016 Community Survey (CS) and the 2011 Population Census (PC) of South Africa.

In a nutshell, this study is a contribution to limited literature on U5MR in South Africa where no significant studies are as yet available to understand the mortality inequality between local municipalities using the above two South African data sets. In order to obtain local municipality estimates of U5MR, this research applies the Fay-Herriot (FH) and spatial FH models [5]. The advantage of using these models is that they only require summary data (area-level aggregate data), not unit-level data that might not be available for the researcher due to confidentiality concerns [6].

The remainder of the paper is structured as follows. Section 2 describes the data sources and model specifications. In Section 3, the results are provided, while Section 4 discusses the findings. The conclusions are then presented in Section 5.

2. Materials and methods

2.1 Data sources

South Africa has nine provinces which in turn, are divided into district municipalities, then local municipalities. We applied the FH and spatial FH models to derive local municipality estimates of U5MR using the available birth history data from the 2016 CS and the 2011 PC. The following variables were required to perform this analysis.

The target variable was drawn from the 2016 CS data conducted by Statistics South Africa (Stats SA). The main goal of this survey was to provide indicators such as population count, an estimate of household count, fertility, mortality, migration, employment, unemployment, the extent of poverty in households. The CS data was accessed from Stat SA website: http://www.statssa.gov. za/. There are $m=214$ sampled local municipalities in the 2016 CS. The variable of interest is U5MR at the local municipality level. In addition to the 2016 CS data, auxiliary aggregate data were taken from the 2011 South African Census.

The auxiliary variables were categorical (i.e., proportions of individuals in each category). Around 40 covariates were selected from the 2011 PC to consider for the modelling. These were gender (male, female); race (Black African, Coloured, Indian or Asian, and White); age (0–14, 15–24, 25–34, 35–44, 45–54, 55–64 and 65–120); employment status (employed, unemployed, not economically active, head employed, head unemployed and head not economically active); employment sector (formal sector, informal sector and private household); marital status (married, living together like married partners, never married, widower/widow, separated, divorced); education level (no schooling, some primary, completed secondary, some secondary, grade 12, tertiary and other); urban area, farm area, Income in South African Rand (No income, R 1–R 76 800, R 76 801–R 614 400, and R 614 401 or more). Reference categories were removed from the analysis since the proportions of individuals in the categories of the auxiliary variable summed up to one [7, 8]. The 2011 South African PC data were accessed from the DataFirst website: https://www.datafirst.uct. ac.za/dataportal/index.php/catalog/485.

Out of these datasets, suitable covariates were chosen using correlation analysis followed by step-wise regression analysis [9]. First, the correlation analysis was run to determine, for example, if one of the covariates (i.e., unemployed category) has a reasonably good correlation with U5MR. Such assessment was replicated for all target and auxiliary variables. Finally, the following variables were identified for further analysis: Black African, no schooling, some primary attained, grade 12, married, divorced, farm area, head unemployed, formal sector, no income and R 614 401 or more.

2.2 Methodology

The FH model has been widely used in small area estimation for a variety of reasons: (i) it uses area level summary data and not unit-level data that might be unavailable to the analyst, because of confidentiality concerns [6]; (ii) its simplicity [10]; and (iii) its ability to provide design-consistent estimators [10].

Consider a finite population of size $N$ is partitioned into $m$ domains of sizes $N_{1},\ldots,N_{m}$ , where $i=1,\ldots,m$ refers to domain $i$ . For this research, each domain represented one local municipality for the purpose of the analysis. In addition, each household or observation $d$ , was assigned to one local municipality, $d=1,\ldots,N_{i}$ . Since U5MR figures were not available for the entire population of South Africa, a sample was drawn from this population using a complex sampling design with sample sizes $n_{i}$ for the local municipalities $i=1,\ldots,m$ . Suppose $y_{i}$ represents the direct survey estimator of $\theta_{i}$ for local municipality $i$ . The aim was to estimate the local municipality estimates for U5MR $\theta_{i}=\bm{x}^{\prime}_{i}\bm{\beta}+v_{i},i=1,\ldots,m$ , where the random effects $v_{i}$ .

2.2.1 The FH area-level model

The widely-used FH model [5] consists of two levels as follows:

Level 1:
Sampling Model:

$\displaystyle y_{i}|\theta_{i}\lx@stackrel{{\scriptstyle\text{iid}}}{{\sim}}N[% \theta_{i},D_{i}],i=1,\ldots,m.$

$D_{i}$ represents the sampling variance of the direct estimator to be estimated from the sample survey.
Level 2:
Linking Model:

$\displaystyle\theta_{i}\lx@stackrel{{\scriptstyle\text{ind}}}{{\sim}}N[\bm{x}^% {\prime}_{i}\bm{\beta},\sigma^{2}_{v}],i=1,\ldots,m.$

where $\sigma^{2}_{v}$ denotes the variance of the random effect, $v_{i}$ . Level 1 and Level 2 models are based on the sample survey (in this case, the 2016 CS) and the borrowed strength component, respectively [11]. The standard FH model links the sampling model to the population model as follows:

$\displaystyle\begin{cases}y_{i}=\bm{x}^{\prime}_{i}\bm{\beta}+v_{i}+e_{i},∼{}i% =1,\ldots,m\\ v_{i}\lx@stackrel{{\scriptstyle\text{iid}}}{{\sim}}N(0,\sigma^{2}_{v}),\\ e_{i}\lx@stackrel{{\scriptstyle\text{ind}}}{{\sim}}N(0,D_{i})\end{cases}$ (1)

where $y_{i}$ is the direct survey estimate for small area $i$ , $D_{i}$ is the known sampling variance of the direct survey estimator in the area $i$ , $v_{i}$ is the area specific random effect that can be used to relate $\theta_{i}$ to $\bm{x}^{\prime}_{i}$ , $\bm{x}^{\prime}_{i}$ is a $k\times 1$ area-specific covariate vector, $\bm{\beta}=(\beta_{1},\ldots,\beta_{k})^{\prime}$ is a $k\times 1$ vector of regression coefficients, and $e_{i}$ is the sampling error. The random effect $v_{i}$ and the sampling error $e_{i}$ are assumed to be independently distributed random variables.

When the unknown parameters in Model Eq. (1) are replaced by their estimators, then the Empirical Best Linear Unbiased Predictor (EBLUP FH here after) for each small area $i$ is:

$\displaystyle\hat{\theta}^{\rm EB}_{i}=y_{i}-\hat{\gamma}_{i}(y_{i}-\bm{x}^{% \prime}_{i}\hat{\bm{\beta}}),i=1,\ldots,m$ (2)

where $\hat{\bm{\beta}}=(\bm{X}^{\prime}\hat{\bm{\Sigma}}^{-1}\bm{X})^{-1}\bm{X}^{% \prime}\hat{\bm{\Sigma}}^{-1}y$ is the weighted least square (WLS) estimator of the vector of $\bm{\beta}$ , $\hat{\gamma}_{i}=\frac{D_{i}}{D_{i}+\hat{\sigma}^{2}_{v}},i=1,\ldots,m$ is the shrinkage factor and $\hat{\bm{\Sigma}}=\operatorname{diag}(\hat{\sigma}^{2}_{v}+D_{1},\ldots,\hat{% \sigma}^{2}_{v}+D_{m})$ .
2.2.2 The spatial FH model

The FH model given in Eq. (1) assumes that the neighbouring small area estimates are spatially uncorrelated. In many applications, however, the random effects between the neighbouring small areas are correlated [6]. The proposed extensions of the FH model by [12] allow spatial association instead of independence between neighbouring small areas following a Simultaneous Autoregressive (SAR) process. This process is defined as [13]:

$\displaystyle v=\rho Wv+u\implies v=(I-\rho W)^{-1}u,$ (3)

where $\rho$ is the spatial autoregressive coefficient, $\bm{W}=\{w_{ij}\}_{m\times m}$ is a spatial proximity matrix that indicates whether the small areas are neighbours or not $(i,j=1,\ldots,m)$ and $\bm{u}=col_{1\leqslant i\leqslant m}(u_{i})_{(m\times 1)}$ is an uncorrelated vector of random effects with zero mean and constant variance $\sigma^{2}_{u}$ . The random effects, $\bm{v}$ has mean zero and covariance matrix $G$ given by $G=\sigma^{2}_{u}((I-\rho W)^{\prime}(I-\rho W))^{-1}$ . Incorporating this substitution into Eq. (1) the spatially-correlated random effects takes the form:

$\displaystyle\hat{\theta}=X^{\prime}\beta+(I-\rho W)^{-1}u+e.$ (4)

The spatial EBLUP (EBLUP SFH here after) estimator of $\theta_{i}$ is given by

$\displaystyle\hat{\theta}^{\text{SEB}}_{i}=x_{i}\hat{\beta}+b^{\prime}_{i}\hat% {G}\hat{\bm{\Sigma}}^{-1}(y_{i}-x^{\prime}_{i}\hat{\beta}),$ (5)

where $b^{\prime}_{i}$ is the $1\times m$ vector $(0,\ldots,0,1,0,\ldots,0)$ with 1 in the ith position, $\hat{\bm{\beta}}=(\bm{X}^{\prime}\hat{\bm{\Sigma}}^{-1}\bm{X})^{-1}\bm{X}^{% \prime}\hat{\bm{\Sigma}}^{-1}y$ is the Generalized Least Square (GLS) estimator of $\beta$ with $\hat{\bm{\Sigma}}=\hat{G}+\bm{D}$ , $\hat{G}$ is obtained by replacing $\sigma^{2}_{u}$ and $\rho$ with their estimates and $\bm{D}=\operatorname{diag}(D_{1},\ldots,D_{m})$ . The vector of regression coefficients $\bm{\beta}$ , $\sigma^{2}_{u}$ and $\rho$ can be estimated iteratively via ML or REML methods using the Nelder-Mead algorithm and the scoring algorithm [6]. For this research, the REML approach was used in the application of both the FH and spatial FH models.

The most common practical problem in small area estimation is measuring the variability associated with the EBLUP FH and EBLUP SFH. The variability under these estimators was measured using mean squared error (MSE). For example, for the REML estimator, an unbiased analytical estimator of the MSE is:

$\displaystyle\widehat{\textit{MSE}}(\hat{\theta}^{\rm EB}_{i})=g_{1i}(\hat{% \sigma}^{2}_{v})+g_{2i}(\hat{\sigma}^{2}_{v})+2g_{3i}(\hat{\sigma}^{2}_{v}).$ (6)

Details about the specification of $g$ components, the estimation of parameters and uncertainty of EBLUP FH and EBLUP SFH can be found in [14, 15, 16, 17].

2.2.3 The log-transformed FH and spatial FH models

Since the relationship between U5MR and the auxiliary variables is non-linear on the original scale, the log-transformed FH and spatial FH models are considered in this application. Following [18, 19], the direct estimator and its variance can be transformed as follows:

$\displaystyle y^{*}_{i}=\log(y_{i})$ $\displaystyle\operatorname{var}(y^{*}_{i})=(y_{i})^{-2}\operatorname{var}(y_{i})$ $\displaystyle x^{*}_{i}=\log(x_{i})$

where $*$ refers the transformed scale.

The EBLUP FH estimator in the transformed scale can be obtained as follows [18, 19]:

$\displaystyle\hat{\theta}^{\text{EB*}}_{i}=y^{*}_{i}-\hat{\gamma}^{*}_{i}(y^{*% }_{i}-{\bm{x}^{*\prime}_{i}}\hat{\bm{\beta}}^{*}),i=1,\ldots,m$ (7)

where $\hat{\bm{\beta}}^{*}={(\bm{X}^{*\prime}\hat{\bm{\Sigma}}^{*-1}\bm{X}^{*})^{-1}% \bm{X}^{*\prime}\hat{\bm{\Sigma}}^{*-1}\bm{y}^{*}}$ is the GLS estimator of the vector of $\bm{\beta}^{*}$ , $\hat{\gamma}^{*}_{i}=\frac{D^{*}_{i}}{D^{*}_{i}+\hat{\sigma}^{2}_{v}}$ and $D^{*}_{i}=\operatorname{var}(y^{*}_{i})$ is an estimate for the sampling error variance. The point estimates in the log-transformed scale need to be a beck-transformed to the original scale for the sake of an appropriate interpretation of the results [18, 19]. Several bias-corrected back-transformations are proposed in the literature. For example, [18] proposed the back transformation of $\hat{\theta}^{\text{EB*}}_{i}$ along with its MSE estimate:

$\displaystyle\hat{\theta}^{\rm EB,\text{Final}}_{i}=\exp\{\hat{\theta}^{\text{% EB*}}_{i}+0.5\widehat{\text{MSE}}(\hat{\theta}^{\text{EB*}}_{i})\},$ $\displaystyle\widehat{\text{MSE}}(\hat{\theta}^{\rm EB,\text{Final}}_{i})=\exp% \{\hat{\theta}^{\rm EB^{*}}_{i}\}^{2}\widehat{\text{MSE}}(\hat{\theta}^{\text{% EB*}}_{i})$ (8)

where $\widehat{\text{MSE}}(\hat{\theta}^{\text{EB*}}_{i})$ is the MSE estimate in the log-transformed scale. This transformation is referred to as crude back-transformation [6]. In addition, [20] developed a bias-corrected small area predictor for the FH model.

In a similar manner to the FH model, the spatial FH model was fitted by using $y^{*}_{i}$ as dependent variable and the log-transformed sampling error variance $D^{*}_{i}$ . To the best knowledge of the author, there is no bias corrected back transformation available for the spatial FH model. For this reason, the point estimates were converted back to the original scale through a direct back-transformation using the exponential function. According to [20], the bias from the variability of parameter estimates vanishes when $m$ is large. In this application, we considered 214 local municipalities; thus correcting the bias from the variability of parameter estimates may not be important. However, the bias from the non-linearity of the back transformation cannot be ignored, even when the number of small areas $m$ is large [20]. This is the limitation of using direct back-transformation in the spatial FH model.

In addition to the point estimates, the prediction intervals were computed using the following equations:

$\displaystyle\text{CI}^{\text{EB,Final}}_{i}=\hat{\theta}^{\rm EB,\text{Final}% }_{i}\mp Z_{\alpha/2}\sqrt{\widehat{\text{MSE}}(\hat{\theta}^{\rm EB,\text{% Final}}_{i})}$ $\displaystyle\text{CI}^{\text{SEB,Final}}_{i}=\hat{\theta}^{\rm SEB,\text{% Final}}_{i}\mp Z_{\alpha/2}\sqrt{\widehat{\text{MSE}}(\hat{\theta}^{\rm SEB,% \text{Final}}_{i})},$ (9)

where $Z_{\alpha/2}$ denotes the $100(1-\alpha/2)\%$ percentile of the standard normal distribution. These confidence intervals were not accurate enough for a small to moderate number of small areas [21, 22]. However, the coverage probabilities of these confidence intervals were close to the nominal level $0.95$ when the number of small areas, $m$ increased to 60 [21, 22]. Note that in this study, the number of small areas, $m=214$ , which is large.

The goodness-of-fit statistics such as log-likelihood, the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) were used to choose the best models (see [6]).

2.2.4 Generalized variance function

In most cases, the sampling variance, $D^{*}_{i}$ , was assumed to be known. This assumption may not be realistic in some applications. This known sampling variance assumption can be relaxed by estimating $D^{*}_{i}$ from the 2016 CS data and then smoothing the estimated variances to get a more stable estimate of $D^{*}_{i}$ [6]. The generalized variance function (GVF) approach of [23] was applied to obtain smoothed estimates of the sampling variances. For this research, following [8], the bias corrected GVF was considered by fitting the following simple linear regression model:

$\displaystyle\operatorname{var}(y^{*}_{i})=b_{0}+b_{1}y^{*}_{i}+\epsilon_{i}$ (10)

where $\operatorname{var}(y^{*}_{i})$ is the dependent variable, $y^{*}_{i}$ as independent variable, $b_{0}$ and $b_{1}$ are the least square estimates and $\epsilon_{i}\lx@stackrel{{\scriptstyle\text{iid}}}{{\sim}}N(0,\sigma^{2}),i=1,% \ldots,m$ . Then the GVF motivated from the above model was given by:

$\displaystyle\widehat{\text{GVF}}=\exp\left(\frac{\hat{\sigma}^{2}}{2}\right)% \exp(\hat{b}_{0}+\hat{b}_{1}y^{*}_{i})$ (11)

the factor $(\exp(\frac{\hat{\sigma}^{2}}{2})=2.193)$ is the bias-correction term in the log-linear analysis for the U5MR estimate. Underestimation of the true variance will occur when we ignore the correction term in the GVF method [8]. Figure 1 shows the estimate for the sampling variance ( $D^{*}_{i}$ ) and predicted variances $(\hat{\textit{GVF}})$ versus the direct survey estimates for U5MR. From this plot, it is easy to observe that the GVF smooths out the unreliable and noisy estimated variance [8, 6].

Table 1

Model comparison based on the log-likelihood, AIC and BIC values

Models	loglike	AIC	BIC
Log-transformed FH model	$-$ 136.48	296.95	337.34
Log-transformed spatial FH model	$-$ 217.44	460.87	504.63

Figure 1.

The dispersion plots for GVF fit.

3. Results

The response variable is the direct survey estimates of U5MR. The FH and spatial FH models were fitted using U5MR as a response variable and the covariates were chosen above as auxiliary variables. Table 1 presents the log-likelihood, AIC and BIC values of the log-transformed FH and spatial FH models. The goodness-of-fit statistics favoured the log-transformed FH model instead of the spatial FH model.

The regression parameters and their corresponding p-values are presented in Table 2. The $\hat{\bm{\beta}}^{*}$ coefficients and p-values for each category were considered in interpreting the results in Table 2. For instance, education level (some primary education), employment status (head unemployed) and income level (no income) were found to be significant predictors of U5MR, with positive $\hat{\bm{\beta}}^{*}$ coefficients. Meanwhile, employment sector (formal) had a negative $(\hat{\beta}^{*}_{9}=-1.68)$ meaning for every one unit increase in employment sector ‘formal’, the target variable decreased by the $\hat{\beta}^{*}_{9}$ coefficient value (i.e. $-$ 1.68). We also observed that the more household with lower income $(\hat{\beta}^{*}_{10}>0)$ the greater was the U5MR and conversely, the more households with high income $(\hat{\beta}^{*}_{11}<0)$ the smaller is the target variable. Notably, Black African, no schooling, married and divorced had no influence on U5MR.

Table 2
Estimates of the log-transformed FH model parameters

Auxiliary variables	$\hat{\beta}^{*}$	Estimate	Std. Error	$t$ -value	$p$ value
$x^{*}_{1}$ (Black African)	$\hat{\beta}^{*}_{1}$	0.06	0.16	0.36	0.72
$x^{*}_{2}$ (No schooling)	$\hat{\beta}^{*}_{2}$	0.05	0.16	0.28	0.78
$x^{*}_{3}$ (Some primary)	$\hat{\beta}^{*}_{3}$	0.67	0.31	2.16	0.03
$x^{*}_{4}$ (Some secondary)	$\hat{\beta}^{*}_{4}$	$-$ 3.21	0.85	$-$ 3.79	0.00
$x^{*}_{5}$ (Grade 12)	$\hat{\beta}^{*}_{5}$	$-$ 1.10	0.48	$-$ 2.28	0.02
$x^{*}_{6}$ (Married)	$\hat{\beta}^{*}_{6}$	$-$ 0.67	0.65	$-$ 1.03	0.30
$x^{*}_{7}$ (Divorced)	$\hat{\beta}^{*}_{7}$	0.39	0.20	1.95	0.05
$x^{*}_{8}$ (Head unemployed)	$\hat{\beta}^{*}_{8}$	0.62	0.25	2.48	0.01
$x^{*}_{9}$ (Formal sector)	$\hat{\beta}^{*}_{9}$	$-$ 1.68	0.82	$-$ 2.05	0.04
$x^{*}_{10}$ (No income)	$\hat{\beta}^{*}_{10}$	0.88	0.34	2.57	0.01
$x^{*}_{11}$ (R 614,401 or more)	$\hat{\beta}^{*}_{11}$	$-$ 1.00	0.46	$-$ 2.19	0.03

The direct, EBLUP FH and EBLUP SFH estimates of U5MR among the South African local municipalities are summarized in Table 3 using percentiles and mean. The median value was equal to 10.12; therefore 50% of the local municipalities (about 107) showed an estimated U5MR greater or equal to 10.12. Furthermore, the local municipality estimates of U5MR generated by the EBLUP FH approach ranges from 2.53 to 33.28 with an average of 11.66. The range of EBLUP SFH estimates appeared to be slightly larger (2.23 vs. 33.34) than the EBLUP FH (2.53 vs. 33.28). This is more likely due to the additional spatial information used in the spatial FH model [24].

Table 3

Summary of direct and model based estimates of U5MR using the FH and spatial FH models

Estimates	Min	P25	Median	Mean	P75	Max
Direct	1.52	7.00	9.92	11.45	14.53	33.33
EBLUP FH	2.53	7.47	10.12	11.66	14.50	33.28
EBLUP SFH	2.23	7.47	10.03	11.62	14.38	33.34

Figure 2.

Bias diagnostics plots with $y=x$ line (red) and regression line (blue) for (a) the log-transformed FH model and (b) the log-transformed spatial FH model.

3.1 Validation of the reliability of the small area estimates

The research considered bias diagnostics, the goodness of fit diagnostics, coefficients of variation (CV) and root MSE to validate the reliability of the small area estimates.

3.1.1 Bias diagnostics

Bias diagnostics were used to examine if the model-based estimates were less extreme when compared to the direct survey estimates [25, 6]. If the model-based estimate is unbiased, then a plot of this estimate against the direct survey estimate would show the scatter plot evenly centered around the unit line $(y=x)$ [25, 6]. Figure 2 shows that the bias scatter plot of the direct survey estimates against the model-based estimates and studied the divergence of the fitted regression line (red) from $y=x$ by testing the intercept $=0$ and slope $=1$ (blue). This figure shows that the slope of model-based estimates of U5MR regression on the true values were close to one (i.e., $\hat{\beta}_{\text{EBLUP FH}}=0.991$ and $\hat{\beta}_{\text{EBLUP SFH}}=0.989$ ), indicating that all the plots include all estimates without outlying estimates. In addition, the model based estimates are approximately unbiased since they lie along the line $y=x$ for most of the local municipalities [25].

3.1.2 Goodness of fit diagnostic

The goodness-of-fit diagnostic proposed by [26] also considered testing whether the direct survey estimates and model-based estimates are statistically different. The hypotheses of interest can be stated as follows:

$H_{0}$ : The direct survey estimates and the model-based estimates are statistically equivalent vs $H_{a}$ : The direct survey estimates and the model-based estimates are statistically different.

The Wald statistic for the goodness-of-fit diagnostic for the log-transformed FH model is given as follows

$\displaystyle W^{*}=\sum_{i}\left[\frac{(y^{*}_{i}-\hat{\theta}^{\text{EB*}}_{% i})^{2}}{\widehat{\text{GVF}}\!+\!\widehat{\text{MSE}}(\hat{\theta}^{\text{EB*% }}_{i})}\right]\!\sim\!\chi^{2}(m),$ (12)

where $\widehat{\text{MSE}}(\hat{\theta}^{\text{EB*}}_{i})$ is the unbiased MSE estimate of $\theta_{i}$ based on the log-transformed FH model. The value of the Wald statistic follows a chi square distribution with $m$ degrees of freedom. The value of the Wald statistics are 16.77 and 31.52 for U5MR estimate for the FH and spatial FH models. These values are less than $\chi^{2}_{0.05{,}214}\approx 256$ indicating that there is no statistical difference between the direct survey and model-based estimates.

Table 4

Summary of the CVs of direct and model-based estimates along with the PIs

				Percent improvement (PI)
	Direct	EBLUP FH	EBLUP SFH	EBLUP FH	EBLUP SFH
Min	0.47	0.47	0.46	0.01	$-$ 0.53
P25	3.5	3.44	3.42	1.58	1.45
Median	6.26	5.92	5.9	4.71	4.97
Mean	7.71	6.52	6.51	7.20	7.83
P75	9.69	8.67	8.44	9.55	10.81
Max	75.13	28.78	33.43	61.7	55.62

Figure 3.

The coefficient of variation of direct, EBLUP FH and EBLUP SFH estimates for each local municipality in South Africa.

3.1.3 Coefficients of variation (CV)

To compare the improved precision of the EBLUP estimates with directs survey estimates, the percent CV was computed. The CV measures the sampling variability as a percentage estimate [25]. The CVs of direct survey estimator were given as $CV(y^{*}_{i})=\sqrt{\widehat{\text{GVF}}}/y^{*}_{i}$ , where $y^{*}_{i}$ is the ith local municipality level direct survey estimates with estimated MSE of $\sqrt{\widehat{\text{GVF}}}$ [27]. The CVs of the EBLUP FH and EBLUP SFH estimators are given as $CV(\hat{\theta}^{\text{EB*}}_{i})=\sqrt{\widehat{\text{MSE}}(\hat{\theta}^{% \text{EB*}}_{i})}/\hat{\theta}^{\text{EB*}}_{i}$ and $CV(\hat{\theta}^{\text{SEB*}}_{i})=\sqrt{\widehat{\text{MSE}}(\hat{\theta}^{% \text{SEB*}}_{i})}/\hat{\theta}^{\text{SEB*}}_{i}$ respectively. In addition to CVs, the percent of improvement (PI) was also calculated as follows: $PI_{i}=100[CV(\widehat{\text{GVF}})-CV(\widehat{\text{MSE}}(\hat{\theta}^{% \text{EB*}}_{i}))]/CV(\widehat{\text{GVF}})$ . Table 4 shows the summaries of the CVs of direct and model-based estimators along with the PIs. The PI for the EBLUP FH varied from 0.01% to 61.70%. The average PI was 7.20%. Moreover, Fig. 3 presents the CVs of model-based estimates and CVs of the corresponding direct estimators for U5MR. The CVs of model-based estimators were smaller than the CVs of the corresponding direct estimators. The gain in precision of the model-based estimators tended to be larger on the right-hand side of the plot (local municipalities with small sample sizes). This fact was observed by [25, 6, 28] in the application of FH and spatial FH models.

3.1.4 Root MSE

The MSE of the model-based estimates and the MSE of the direct survey estimates should be comparable if the fitted model fits the data reasonably well [29, 6]. For this reason, this research examined the two MSEs (i.e., the MSE of the model based estimator and the MSE of the direct estimator). Figure 4 plots the estimated model-based root MSEs of the EBLUP FH and EBLUP SFH and the root MSE of the direct estimators of U5MR. It can be easily observed that the root MSE of model-based estimators tended to be smaller than the root MSE of the corresponding direct survey estimators. This is the case since minimizing the model-based MSE is one of the properties of the class of the model-based unbiased linear predictors [29, 6].

Figure 4.

The root MSEs of direct and model based estimates for each local municipality in South Africa.

Figure 5.

The plots of residuals versus estimates for the log-transformed FH model (top) and the log-transformed spatial FH model (bottom).

3.2 Diagnostics measures to check model assumptions

Checking model assumptions is crucial since the random effects in the FH model are identically independently distributed normal with a mean of zero [6]. In addition, the optimality properties of the model-based estimates depend on the extent to which distributional assumptions are true [28]. For this research, standard diagnostics such as residual vs estimates (both direct and model based estimates) were used for this purpose. Figure 5 presents the standardized residuals versus model-based estimates (left) and standardized residuals versus direct estimates (right). The direct takeaway from these plots is that both graphs have smoother shapes because they do not include anomalous observations [8].

Figure 6 presents the distribution of standardized residuals (left), histograms of standardized residuals (centre) and standard normal $q$ - $q$ plots of standardized residuals (right). The left plots show that standardized residuals appear to be randomly distributed around zero. The centre plots show that histograms of these residuals appear to be symmetric. The right plots show that the q-q plots of these residuals are well-distributed along a straight line. These plots provide evidence in support of the normality assumption of the standardized residuals. The Shapiro-Wilk normality test was also used to assess the normality assumption of the standardized residuals. In this study, the $p$ -values of Shapiro-Wilk normality test on the standardized residuals are 0.816 and 0.981 for the log-transformed FH and spatial FH models, which suggests non-rejection of the hypothesis of normality.

Figure 6.

Distribution of the standardized residuals (left), histogram of the standardized residuals along with it normal fit (centre) and q-q plots for the standardized residuals (right) for the log-transformed FH model (top plots) and spatial FH model (bottom plots).

Finally, Fig. 7 shows the spatial mapping of local municipality-wise U5MR produced using the EBLUP FH (top left) and EBLUP SFH (top right) estimates along with their corresponding CVs under the log-transformed FH (bottom left) and spatial FH (bottom right) models.

Figure 7.

Map of U5MR in South Africa for the log-transformed FH model (left) and the log-transformed spatial FH model (right).

4. Discussion

This article provides the local municipality estimates of U5MR in South Africa using the 2016 CS and the 2011 PC. In particular, since U5MR depends on the binary outcome coded as zero for dead and one for alive, neither linear mixed models nor unit-level model-based estimation approaches can be used [8]. “Instead, area-level models provide an easy-to-apply solution” [8]. Therefore, for this research, log-transformed FH and spatial FH models were applied to obtain local municipality estimates of U5MR in South Africa. These models have not received much attention in the mortality literature.

The regression analysis results indicated that auxiliary variables such as education level (some primary education, some secondary education and attainment of Grade 12), employment status (head unemployed), income level (no income and R 614 401 or more) and employment sector (formal) were strong predictors of U5MR. The results of this study are discussed here in light of related literature.

In the South African context, and indeed in other parts of the world, U5MR is linked to socio-economic factors such as poverty and/or lack of income, food insecurity, unemployment, education, the global food price increase and so on. Income (lack of income) increased U5MR, while income (R 614 401 or more) decreased U5MR. Several studies have demonstrated the relationship between mortality and poverty/lack of income. For example, researchers such as [30, 31, 32] have argued that there is a positive relationship between poverty and mortality rates (particularly child mortality). Furthermore, [33] emphasized in his essay that the survival conditions of poor households are worse when compared to non-poor households. In another study, [34] argued that persons living in poverty have higher mortality rates than those out of poverty. In other words, child mortality rates are low in high income earners compared to low income earners [35]. Education level (some primary education, some secondary education and attainment of grade 12) and employment status (head of household unemployed) were also strong predictors of U5MR. Mothers with relatively low levels of education are less likely to be employed, which prevents them from investing in better health outcomes for their children [36]. In addition, lack of income and/or poverty is the principal driver of South Africa’s food insecurity [37], and food insecurity is an important determining factor of child mortality [38]. The food prices increases due to economic crises result in death in famine situations, which in turn leads to an increase in childhood mortality [39].

4.1 Which local municipality is most at risk?

The local municipality-wise U5MR produced by the EBLUP FH approach in South Africa ranges from a minimum value of 2.53 observed for the Drakenstein local municipality in the Western Cape province to a maximum value of 33.28 observed for the Khai-Ma local municipality in the Northern Cape province.

The lowest ( $<$ 4%) and highest ( $>$ 30%) estimates of U5MR, along with their 95% confidence interval (CI) are presented in Table 5. Along with Fig. 7, Table 5 shows the South African local municipalities where the U5MR was highest ( $>$ 30 deaths per 1000 live births) are those situated in Mpofana with 30.41 (95% CI $=$ 29.15, 31.68), Laingsburg with 30.71 (95% CI $=$ 29.46, 31.96), Maquassi Hills with 32.15 (95% CI $=$ 31.00, 33.29), Tswaing with 32.47 (95% CI $=$ 31.35, 33.59) and Khai-Ma with 33.28 (95% CI $=$ 32.22, 34.35). On the other hand, local municipalities where the U5MR was smallest ( $<$ 4 deaths per 1000 live births) were those situated in the Drakenstein with 2.53 (95% CI $=$ 1.29, 3.76) in Western Cape, Midvaal with 2.96 (95% CI $=$ 1.55, 4.38) in Gauteng, Cederberg with 3.45 (95% CI $=$ 1.83, 5.06) also in the Western Cape and Mossel Bay with 3.29 (95% CI $=$ 1.76, 4.81) in the Western Cape. This pattern is not surprising, since by poverty and inequality standards, the Western Cape along with Gauteng and the Free State are the richest provinces, while Limpopo, Eastern Cape, and KwaZulu-Natal are the poorest provinces in South Africa [40]. Furthermore, the Western Cape is the richest province in terms of average annual consumption expenditure [41].

Table 5
Direct, EBLUP FH and EBLUP SFH estimates of U5MR (per 1000 live births), CVs of direct, EBLUP FH and EBLUP SFH estimators ( $\times$ 100) and 95% confidence intervals of EBLUP FH and EBLUP SFH estimates for the South African local municipalities with the lowest U5MR ( $<$ 4 deaths per 1000 live births) and highest U5MR ( $>$ 30 deaths per 1000 live births). For example, the lower and upper limits of the confidence interval for $\theta_{i}$ based on the EBLUP FH estimate is obtained as: $(L_{i},U_{i})$ , where $L_{i}=\hat{\theta}^{\rm EB,\text{Final}}_{i}-Z_{\alpha/2}\sqrt{\hat{\rm MSE}(% \hat{\theta}^{\rm EB,\text{Final}}_{i})}$ and $U_{i}=\hat{\theta}^{\rm EB,\text{Final}}_{i}+Z_{\alpha/2}\sqrt{\hat{\rm MSE}(% \hat{\theta}^{\rm EB,\text{Final}}_{i})}$

							Confidence interval (CI)
							EBLUP FH		EBLUP SFH
Municipality	Direct	CV	EBLUP FH	CV	EBLUP SFH	CV	$L_{i}$	$U_{i}$	$L_{i}$	$U_{i}$
	U5MR ( $<$ 4 deaths per 1000 live births)
Drakenstein	1.52	75.13	2.53	28.78	2.23	33.43	1.29	3.76	0.00	4.79
Midvaal	1.90	46.97	2.96	23.83	3.20	20.84	1.55	4.38	0.70	5.69
Cederberg	2.19	37.52	3.45	20.47	2.83	25.74	1.83	5.06	0.27	5.39
Mossel Bay	2.40	33.05	3.29	21.07	2.77	24.24	1.76	4.81	0.26	5.28
	U5MR ( $>$ 30 deaths per 1000 live births)
Mpofana	30.47	0.62	30.41	0.62	30.41	0.62	29.15	31.68	28.41	32.41
Laingsburg	30.77	0.61	30.71	0.60	30.75	0.60	29.46	31.96	28.75	32.75
Maquassi Hills	32.18	0.52	32.15	0.52	32.11	0.52	31.00	33.29	30.12	34.11
Tswaing	32.50	0.51	32.47	0.51	32.48	0.51	31.35	33.59	30.48	34.47
Khâi-Ma	33.33	0.47	33.28	0.47	33.34	0.46	32.22	34.35	31.34	35.33

The reliability of the model-based estimates of U5MR was examined using bias diagnostics, the goodness of fit diagnostics, CVs and root MSEs. These measures indicated that the model-based estimates are superior to the direct survey estimates, since the model-based estimates borrow strength from auxiliary variables from the 2011 PC [29]. Moreover, the standard diagnostics such as residual vs estimates (both direct survey- and model-based) were used to asses the models assumptions. The empirical results from the model diagnostics have shown that the assumptions are satisfactorily met.

4.2 U5MR maps

Figure 7 (top) maps the EBLUP FH (left) and EBLUP SFH (right) estimates of U5MR. These maps give a visual illustration of estimated local municipality-wise U5MR in South Africa. They provide an important source of information on the distribution of local municipality estimates of U5MR with yellow colour corresponding to a better situation (lower estimate for the U5MR), while red corresponds with a higher estimate for the U5MR. Figure 7 (bottom) maps the corresponding CVs of the EBLUP FH (left) and EBLUP SFH (right) estimates (shown in yellow) indicate local municipalities with lower-values CVs. For example, in the northern part of South Africa (North West and Northern Cape provinces) there are many local municipalities with a high U5MR (shown in red). The values of CV are low (shown in yellow) in this part of South Africa. The Western Cape and Gauteng, in contrast, have many local municipalities with low U5MR (shown in yellow) with high corresponding CVs (shown in red). These maps are currently used to identify the policy priorities for reducing U5MR, for instance to allocate funds to local jurisdictions.

5. Conclusion and policy implication

Estimation of mortality-related variables at the sub-national scale is a challenging statistical problem because sample sizes are too small and simple parametric models may only poorly capture the relationships between these variables [42]. At the same time, there is a growing demand for accurate, reliable and cost-effective U5MR among policy- and decision-makers, as well as private and public sector administrators using data from different sources. As a consequence, many area- and unit-level small area models have been proposed in the literature [25]. In this context, the FH and spatial FH models have the potential to make a real impact by providing smoothed and stabilized estimates [4].

Having good small area estimates is important because most policies will generally be based on this type of information. The results indicate that root MSEs and CVs of the direct survey estimates are generally larger than model-based estimates. This is the case because the direct estimates of U5MR in local municipalities with a sparse sample in a 2016 CS do not provide adequate precision [42, 6]. From this, we can conclude that a model-based estimate of U5MR is more reliable than a direct survey estimate.

Both local municipality estimates, as well as the spatial map of U5MR have great potential in assessing mortality issues, as they give useful information in identifying municipalities with higher child mortality rates for targeted interventions. Most importantly, the local municipality level U5MR statistics are useful for (i) relevant Departments and Ministries of South Africa to implement policy interventions for targeted local municipalities with high U5MR, as well as for ensuring effective administrative and financial decisions (e.g., budget allocation) (ii) international organizations for their strategic planning and policy research.

Furthermore, accurate and cost-effective U5MR statistics at the local municipality level were produced without the need for more resources/surveys through combining the available data sources (i.e., the 2016 CS and the 2011 PC).

Footnotes

Acknowledgments

The author would like to acknowledge Statistics South Africa for the 2016 CS data and Data First for the 2011 Census data. In addition, the author thanks the editor and two anonymous referees for their valuable comments and suggestions; these led to a considerable improvement in the article. All errors are mine.

References

Anggreyani

Kurnia

. Small area estimation for estimating the number of infant mortality in West Java, Indonesia. AIP Conference Proceedings. 2016; 1707(1): 080002.

Burstein

Henry

Collison

Marczak

Sligar

Watson

, et al. Mapping 123 million neonatal, infant and child deaths between 2000 and 2017. Nature. 2019; 574(7778): 353-358.

Bamford

Mckerrow

Barron

Aung

. Child mortality in South Africa: Fewer deaths, but better data are needed. South African Medical Journal. 2018; 108(3): 25-32.

Sugasawa

Kubokawa

Ogasawara

. Empirical uncertain bayes methods in area-level models. Scandinavian Journal of Statistics. 2017; 44(3): 684-706.

Fay

Herriot

. Estimates of income for small places: An application of James-Stein procedure to census data. Journal of the American Statistical Association. 1979; 74: 269-277.

Rao

JNK

Molina

. Small area estimation. 2nd Edition John Wiley and Sons, New York. 2015.

Benavent

Morales

. Multivariate Fay-Herriot models for small area estimation. Computational Statistics & Data Analysis. 2016; 94: 372-390.

Esteban

Morales

Perez

Santamaria

. Small area estimation of poverty proportions under area-level time models. Computational Statistics & Data Analysis. 2012; 56(10): 2840-2855.

Chandra

Salvati

Sud

. Disaggregate-level estimates of indebtedness in the state of Uttar Pradesh in India: An application of small-area estimation technique. Journal of Applied Statistics. 2011; 38(11): 2413-2432. Available from: http://ro.uow.edu.au/cssmwp/69.

10.

Lahiri

. An adjusted maximum likelihood method for solving small area estimation problems. Journal of Multivariate Analysis. 2010; 101(4): 882-892.

11.

Chatterjee

Lahiri

. Parametric bootstrap approximation to the distribution of EBLUP and related prediction intervals in linear mixed models. Annals of Statistics. 2008; 36(3): 1221-1245.

12.

Salvati

. Small area estimation by spatial models: The spatial empirical best linear unbiased prediction (Spatial EBLUP). Working Paper 2004/03, Firenze: Universita degli Studi di Firenze.

13.

Anselin

. Spatial econometrics: Methods and models. Journal of the American Statistical Association. 1988; 85(411): 905-907.

14.

Datta

Lahiri

. A unified measure of uncertainty of estimated best linear unbiased predictors in small area estimation problems. Statistica Sinica. 2000; 10: 613-627.

15.

Chandra

Salvati

Chambers

. Small area estimation for spatially correlated populations – a comparison of direct and indirect model-based methods. Statistics in Transition. 8(2): 331-350.

16.

Petrucci

Pratesi

Salvati

. Geographic information in small area estimation: Small area models and spatially correlated random area effects. Statistics in Transition. 2005; 7(3): 609-623.

17.

Pratesi

Salvati

. Small area estimation in the presence of correlated random area effects. Journal of Official Statistics. 2009; 25(1): 37-53.

18.

Neves

Silva

DBDN

Correa

. Small domain estimation for the Brazilian service sector survey. Estadistica. 2013; 65: 13-17.

19.

Kreutzmann

Marek

Salvati

Schmid

. Estimating regional wealth in Germany: How different are East and West really? Discussion Paper, Deutsche Bundesbank. 2019.

20.

Chandra

Aditya

Kumar

. Small-area estimation under a log-transformed area-level model. Journal of Statistical Theory and Practice. 2018; 12(3): 497-505.

21.

Diao

Smith

Datta

Maiti

Opsomer

. Accurate confidence interval estimation of small area parameters under the fay-herriot model. Scandinavian Journal of Statistics. 2014; 41(2): 497-515.

22.

Shiferaw

Galpin

. A corrected confidence interval for a small area parameter through the weighted estimator under the basic area level model. Journal of the Iranian Statistical Society. 2019; 18(1): 17-51.

23.

Wolter

. Introduction to variance estimation, 2nd Edition. Springer. 2007.

24.

Pratesi

Salvati

. Small area estimation: the EBLUP estimator based on spatially correlated random area effects. Statistical Methods and Applications. 2008; 17(1): 113-141.

25.

Chandra

Salvati

Chambers

. Small area estimation under a spatially non-linear model. Computational Statistics & Data Analysis. 2018; 126: 19-38.

26.

Brown

Chambers

Heady

Heasman

. Evaluation of small area estimation methods – An application to unemployment estimates from the UK LFS. Proceedings of Statistics Canada Symposium 2001: Achieving Data Quality in a Statistical Agency: A Methodological Perspective. 2001.

27.

Luzi

Solari

Rocci

. A study of small area estimation for italian structural business statistics. Journal of Official Statistics. 2018; 34(2): 543-555.

28.

Molina

Marhuenda

. sae: An R package for small area estimation. R Journal. 2015; 7(1): 81-98.

29.

Morales

Santamaria

. Small area estimation under unit-level temporal linear mixed models. Journal of Statistical Computation and Simulation. 2019; 89(9): 1592-1620.

30.

Fritzell

Kangas

Hertzman

Blomgren

Hiilamo

. Cross-temporal and cross-national poverty and mortality rates among developed countries. Journal of Environmental and Public Health. 2013; 2013: 915490-915490.

31.

Engels

. The condition of the working class in England. Collected Works of Marx and Engels. 1845; 4: 295-596.

32.

Rowntree

. Poverty: A study of town life. Macmillan, London, UK. 1901.

33.

Malthus

. An essay on the principle of population. 1798.

34.

Marmot

. The influence of income on health: Views of an epidemiologist. Health Affairs. 2002; 21(2): 31-46.

35.

Campbell

, De Pee

Sun

Kraemer

Thornelyman

Moenchpfanner

, et al. Relationship of household food insecurity to neonatal, infant and under-five child mortality among families in rural Indonesia. Food and Nutrition Bulletin. 2009; 30(2): 112-119.

36.

Green

Peters

. Region of birth and child mortality among black migrants to south africa: Is there a foreign-born advantage? Migration Letters. 2016; 13(3): 359-376.

37.

Misselhorn

Hendriks

. A systematic review of sub-national food insecurity research in South Africa: Missed opportunities for policy insights. PLOS ONE. 2017; 12(8).

38.

Oluwatoyin

. Impact of food availability on child mortality: A cross country comparative analysis. Agricultural Economics-zemedelska Ekonomika. 2017; 63(6): 283-297.

39.

Christian

. Impact of the economic crisis and increase in food prices on child mortality: Exploring nutritional pathways. Journal of Nutrition. 2010; 140(1).

40.

Chitiga-Mabugu

Ngepah

Nhemachena

Motala

Mathebula

Mupela

. State of poverty and its manifestation in the nine provinces of South Africa. Human Sciences Research Council. 2014.

41.

Stats

. An analysis of household expenditure and income data using the LCS 2014/2015. Statistical Release P0310. 2017.

42.

Pratesi

. Analysis of poverty data by small area estimation. John Wiley and Sons. 2015; 1-436.

Analysis of the spatial distribution of under-5 mortality rate in local areas of South Africa

Abstract

Keywords

1. Introduction

2. Materials and methods

2.1 Data sources

2.2 Methodology

2.2.1 The FH area-level model

Table 2 Estimates of the log-transformed FH model parameters

3.1.1 Bias diagnostics

3.1.2 Goodness of fit diagnostic

3.1.4 Root MSE

4.1 Which local municipality is most at risk?

5. Conclusion and policy implication

Footnotes

Acknowledgments

References

Table 2
Estimates of the log-transformed FH model parameters