The common sources of data for cause-specific mortality include vital registration systems, sample registration systems, nationally representative household surveys and sentinel Demographic Surveillance Sites (DSS) for epidemiological studies. With an exception of a few countries, such as South Africa, reliable and functioning vital registration systems have been presented a challenge in supporting attribution of causes of child death in many low-middle income countries, particularly in sub- Saharan Africa [29, 30].

The main source of data for researchers to guide policy makers in a developing country such as Nigeria is the National DHS conducted by the Data Measure program. The United States Agency for International Development (USAID) has provided the technical assistance and funding to conduct surveys in several developing countries, thereby promoting global understanding of public health. The DHS program collects survey data nationally on a variety of socio-demographic and health related issues. The survey collected information about the background of the respondents, specifically collected information on fertility levels, marriage, fertility preference, awareness and use of family planning methods, child mortality and child nutrition. Detailed information and procedures about the data collection, and questionnaires have been published elsewhere by [31].

The 2013 NDHS survey conducted by the DHS measure used a multi-stage cluster design consisting of 40,320 households in 904 clusters with 372 in urban areas and 532 in

Fig. (1). Map of Nigeria showing 37 districts (36 states and the Federal Capital Territory (FCT), Abuja).

rural areas. The survey successfully interviewed 38,948 women occupied in 38, 520 households nested in 886 clusters. This yielded a household response rate for women of 99%. Data extracted from the 2013 NDHS for the present study are: the number of children born between 2008 and 2013, the number of children alive and counts of child deaths at the time of the survey, the proportion of poorest and poor households, the number of cases (children) experiencing diarrhoea two weeks prior to the survey, the number of households using solid cooking fuels such as, coal, charcoal, fire wood, cow dung and agricultural crop residues.

For the purpose of the present study, Fig. (1) shows the geographical map of Nigeria showing 36 states (districts) and the Federal Capital Territory, Abuja. Nigeria comprises of six geopolitical regions; North-East, North-West, North-Central, South-East, South-South, and South-West which are sub-divided into 36 administrative states and the Federal Capital Territory (FCT). The population groupings within the geopolitical regions and states are relatively homogeneous. Also, the people's cultural beliefs such as the demographic characteristics, arid environment factors and socio-cultural structures are considered similar within the geopolitical zones and states.

In disease mapping, the first step is the removal of the effect of the confounding factors on the risk estimate in the study population through distribution standardization. Standardization of mortality rates or disease incidence is a basic tool in both demography [32] and epidemiology [33, 34]. The most frequently used method in epidemiology is the traditional method for estimating the relative risk is the internal standardization method, which calculates the expected disease counts as functions of the observed numbers of cases. However, such models are regarded as incoherent and not generative probabilistically according to [35], because the observed count appears on both sides of the equation.

Consider the death counts, Y_k aggregated data over a state (district), say, 37 (k = 1;::: ; 37) states, where the mother, k resides in Nigeria. In this study, the SMR is calculated as

(1)

and

In equation (1), Y_k is a random variable representing the number of observed cases (under-five deaths) in each k^th state(district) and n_k represents the number of children at risk (number of under-five children in each state, k). In addition SMR_k is calculated as the ratio of observed number of child death cases to the expected number of cases in the k^th state, representing the risk of each k^th area (state). Whenever the value of SMR is greater (lower) than one (1), it indicates that the area (state) k has a higher (lower) risk than the average disease risk of the whole region. For example, for SMR_k, it can be said that area k has a 25% higher risk of the disease (childhood morbidity). These quantities SMR are plotted as a crude map. This estimator is unbiased, and is frequently used by epidemiologists. However, this estimate is based only on a sample size of one and hence it is not really statistically useful because it is a saturated model. Some of the advantages and disadvantages of a crude map of the SMR have been highlighted by Lawson et al. [15].

In recent years, the attempts to map incidence and mortality from diseases such as cancer have been explored. Such maps usually display either relative rates in each district or province, as measured by a SMR or similar index. The standard models are detailed in the literature on the empirical methods and its applications can be found in [36, 37]. Clayton [38] has earlier provided approach for estimating the SMR using Bayesian methodology. According to [26] the mean of the estimator, _k and its variance, which will be large if the expected number of incidence cases is small. This is one of the disadvantages of using the SMR. Other disadvantages are discussed by [39], showing that the SMR is based on a ratio estimator.

Fig. (2). Descriptive summary maps: (a) observed child death counts, (b) crude standardized mortality ratio SMR (c) proportion of poor households (d) smoothed relative risk RR map from independent Poisson-Gamma distribution.

From Fig. (2a), (the raw death counts map), a cluster of child mortality was occurred and concentrated in the northern part of Nigeria. With reference to Fig. (2c), most northern states are darker than the southern regions. This indicates that there are more economically deprived households in northern Nigeria than in the southern regions. Reference to Fig. (2b) , in (SMR), shows few states with an unusual low or (high) mortality incidence, while such districts (states) are bordered /surrounded by relatively high counts. For instance, the observed (isolated) low incidence of mortality in Borno State, which is surrounded by states with relatively high death counts, is an evidence of geographical disparities with no clear patterns. Another scenario can be seen in a state like Zamfara (with a high mortality count (130)), which shares boundaries with states in the north -western region with relative high mortality risk states such as Kebbi (109.55) and Sokoto (124.77). These scenarios may better be handled by a spatial random effect model or the BYM model, which exhibit borrowing strength that is regions close to each other share similarities in their prevalence.

Fig. (1a) Displays the map of observed mortality counts across the 37 states(districts) in Nigeria, Fig. (2b) is the map of SMRs of the child mortality and corresponding Table 1 with SMRs that vary widely around their mean, 0.920, (standard deviation: 0.306). Although, the evidence of observed lower SMRs were recorded in the southern regions of the country, the geographic variation in child mortality with clustering of high mortality observed in the northern states with a relatively low mortality prevalence (isolated area) of Borno state is apparent. No clear spatial pattern emerges from the map.

Fig. (1b) also depicts the empirical estimation of SMR of child mortality. The geographical patterns of child mortality distribution are similar for Fig. (1a), raw mortality count map and the crude SMR 2(b). The clusters of high mortality or concentrated mortality found in the northern regions can be attributed to partly unobserved heterogeneity and environmental factors. The smooth SMR map (d) shows evidence of localized spatial smoothness and neighbouring states exhibited similar patterns of mortality risk, while other neighbouring regions far apart showed different (local) risks patherns. In Figs. (2b and 2d), regions coloured black show the SMRs with RR greater than one indicates significantly excess(higher) mortality . The light coloured regions are states signify low prevalence (RR less than one) of child mortality, while the grey coloured regions are not significant.

The map of Relative Risk (RR mean) from the independent Poisson model is shown in Fig. ( 2d). The map depicts that the model could not capture the geographical variation in the spatial pattern of the actual mortality count data.

For instance, discontinuities can be seen in some states with clustering of high mortality rates in the north when compared with the observed counts map. Katsina State recorded the highest expected counts and posterior mean RR, but the Poisson model (no random effect) would classify Katsina lower than the actual mortality level. Another scenario in the spatial disparities was also observed in Ekiti State in south-west Nigeria with small expected counts, but the state was an elevated high risk. This dispersion can be attributed to small population size. The independent Poisson sometimes under-estimates the mortality risk such as Balyesa and Lagos, perhaps these could result from a high expected value (denominator).

A careful inspection of the expected counts in Fig. (2) reveals that higher child mortality risk were detected in some states of the north regions resulting from empirical computation of SMR i.e. Kano, due to large expected count of 197.827 (divisor)(refer to Table 2). Four other states were considered with expected counts of 46.569, 51.153, 61.420 and 131.64 corresponding to approximate percentile values of 10th, 25th, 50th and 90th of the expected counts respectively. It is worth mentioning that the choice of unusually low relative risk values (SMR) (10^th percentile of the expected counts) would establih an epidemiological importance. Interestingly, some districts (states) had unusual low death rates surrounded by neighbouring states with a fairly high mortality risk. The expected mortality rates of the four states correspond to the percentiles (10^{th -} 90^th percentile): Abia (46.57), Osun (51.12), Plateau Sate (61.42) and Jigawa (131.64). For example, Plateau state (61.42) had relative low prevalence but is surrounded by states had relatively high mortality. In such cases, the mortality rates in those neighbouring states may have a substantial influence on the smoothing effects on states that share borders. This scenario can be handled by the spatial Conditional Auto-Regressive (CAR) model.

To conclude this section, in comparing the smooth maps of the SMR map and the PG map, it shows that there was no clear difference in the smooth risk maps from both estimates. The empirical approach makes epidemiological sense and provides better understanding of mortality prevalence across the regions in Nigeria. These maps are primarily used as a tool for identifying regions with unusually (low) high risk area, so that further attention can be given to these priority districts (states).

Mapping mortality rates or disease incidence could provide important information in many epidemiological studies for resource allocation and disease management. To estimate and map crude mortality rates, particularly rare disease aggregated at the administrative unit or regional level can be statistically challenging if the high variability of population sizes over a small area is not taken into account. To mitigate the problem, an exploratory data analysis was carried out by mapping the standardized mortality ratio as suggested by [32]. The following four models are explored to capture the effects of spatial dependence and overdispersion in the data

Model 1: The Poisson-Gamma model is sometimes used to model the relative risk of the number of child mortality in a district (state). The relative risk combines with the Poisson likelihood function for the death counts and Gamma prior distribution to yield a Gamma posterior distribution for the relative risk [36].

Let y_i and E_i; i = 1;::: ; n, denote the observed and expected number of death cases in district (state) i. We assume the death count that y_i ~ Poisson(Eϑ), where ϑ is the unknown relative risk and Poisson mean µ_i is modeled as

(2)

We assume that ϑ = Gamma(a,b) for i = 1::: 37 in our study n=37 districts (states) in Nigeria. By combining the likelihood and the prior distribution, the posterior mean or the relative is obtained as

where ; represents a weighted average that indicates how much the posterior mean shrunk towards the individual expectation, E_i as explained in [37]. One advantage of the Poisson-gamma model is that it provides a simplified way to accommodate over-dispersion in the model. A drawback is that this Poisson-gamma model does not permit the inclusion of covariate(s) [36, 38].

Model 2: Clayton and Kaldor [16] first proposed a Poisson log normal model that combines the relative risk and a normally distributed random variable. The model includes area-specific random effects or spatially unstructured random effects, v_i and ϑ is the overall level of the relative risk .

Using equation (1) above, µ_i = Eϑ, the log normal model for the relative risk becomes

(3)

where the linear link function η = log(ϑ) = X'β+v_i

v_i is the spatially unstructured random effects that were modeled as using the Gaussian prior distribution with a zero mean and the variance, i.e.v_i ~ N(0, ) where represents specific area variance. X is a vector of covariates (such as proportion of poor households, unimproved source of drinking water, unprotected toilet, children having diarrhoea, the proportion of mothers using solid fuels (coal, wood, agricultural residues cow dung etc) as cooking method . Thus, the relative risk provides a more flexible alternative to the independent Poisson model, as stated in [39].

Model 3: The conditional autoregressive (CAR) model has been widely used for the analysis of spatial data in different areas, such as demography, geography and epidemiology. This model was introduced by [40] as a spatial methodology to estimate disease risk, which assumed a spatial dependence with neighboring regions. The u_i is the spatially structured (correlated) random effects were modeled using the conditional autoregressive prior distribution as suggested by [40].

Using equation (1) above, u_i, the CAR model for the relative risk becomes

(4)

and the linear link function becomes η = log(ϑ) = X'β+u_i

where u_i where area i ~ j are adjacent (neighbours), and w_ij = 1 and zero if they are not. X is a vector of covariates (such as proportion of poor households, unimproved source of drinking water, the proportion of households using unprotected toilets, the number of children having diarrhoea, proportion of mothers using solid fuels (coal, wood, agricultural residues, cow dung etc) as cooking methods.

Model 4: Besag, York and Mollie (BYM) model was first introduced by [16] and later extended by [40]. BYM model is then split into two spatial random and heterogeneity components and it is formulated through the following equation. The death count assumes, y_i ~ Poisson(ϑ) the log relative risk is modeled through equation (1) above, µ_i = Eϑ, and the BYM model for the relative risk becomes

(5)

and the linear link function becomes η = log(ϑ) = X'β + u_i + v_i

where v_i and µ_i are unstructured and structured spatial random effects respectively. They are model as v_i ~ N(0, ) and where area i ~ j are adjacent (neighbours), w_ij = 1 and zero if they are not. X is a vector of covariates as stated above . The ϑ reflects the amount of extra Poisson variation in the data and represents specific area variance as stated in [39]. The precision parameters τ_u² and τ_v² control the variability of u and v respectively. The parameter estimation was executed via the Bayesian Markov Chain Monte Carlo Convergence of the MCMC, which was reached at 15000 iteration after a burn-in period of 5,000 samples and the thinning was done at every 90th element of the chain. The statistical inference is based on full Bayesian framework and prior distributions were specified for the model parameters. The posterior estimates are used to explain the model results of the UH, CAR and the BYM model which are presented in Table 4.

The model performance was investigated via Deviance Information Criterion (DIC) which is due to Spiegelhalter et al. [41] given as

(6)

where D is the posterior mean of the deviance and is the vector of model parameters. pD is the number of effective parameters in the model that penalizes its complexity. DIC takes into account both the model fit (summarized by D) and model complexity (captured by PD) when comparing models. Therefore, the model having the smaller value of DIC is the most preferred one as it achieves a more optimal combination of fit and parsimony.

The parameter estimation was done using Bayesian Markov Chain Monte Carlo via Gibbs Sampling. The convergence of the MCMC was achieved at 15,000 iterations after a burn-in period of 5,000 samples and thinning of every 90th element of the chain. The hyper-prior prior distributions assumed for the precision parameters, τ_u², τ_v² and τ_β² are Gamma distributions as τ_u², ~ Г(0.05,0.005), τ_v² ~ Г(0.05,0.005) and τ_β², ~ Г(0.05,0.005) respectively. The coefficients of the covariates of the regression model are assumed to be normally distributed given as, β ~ N(0,0.005). All model analyses were carried out in WinBUGS after [41] and data manipulation was done in R programming [42]

Table 2 presents the number of child deaths, total births, expected deaths, and relative frequency distribution. The study involved 31482 children born between 2008 and 2013, out of which 2886 children died before reaching the age of five. Zamfara recorded the highest child mortality and relative frequency of 229 (7.96) and the second highest occurred in Kano, 221(7.66). Both states are found in the north-western region of Nigeria. The lowest under-five mortality was recorded in Osun state of 21(0.73).

Table 3 presents the estimates of the parameters and goodness of fit for the hierarchical models discussed in the previous section. The non-spatial method (P-Gamma model) does not account for autocorrelation in the residuals, although they appear to perform reasonably well overall.

Although the CAR model and BYM model each provides important information about clustering of the childhood mortality relative risk pattern, one would recommend that the BYM is the best fitted model for Nigerian child mortality data, since it yielded the lowest value of the DIC = 285:310 and with a lower pD= 25.04. The CAR model had DIC= 286.40 and pD=(24.16) as the goodness of measure and it competes closely with the BYM model. However, the BYM model is the most preferred one due to its robustness and at the same time one can evaluate the proportional of variation that can be attributed to spatial dependence (clustering) and the variation due to random heterogeneity effect structure of the mortality prevalence.

Table 4 presents the posterior statistics of the fitted hierarchical models. It can be observed that the posterior mean of P-G model is 0.923: 95% CI (0.826, 1.030), which is approximately the same as the mean of the SMR of 0.920 and standard deviation, 0.306. The overall population parameters, a = 10:310; (6:232; 15:970) and b = 11:200(6:680; 17:350) from the Poisson Gamma model. The Poisson -log normal (PLN) model yielded a precision variance of, τ_v²= 41.76 with a standard deviation of 0.168. This indicates that the relative risk of child mortality at any given state is similar (less heterogeneous) to that of its neighbours. The CAR model's precision variance, τ_u² =14:34; (5,006 to 35.47) and standard deviation of 0.291, which indicates that the geographic patterns of under five mortality exhibits more of clustering across the selected administrative units (states)in Nigeria. The precision variance parameter of the BYM model has CAR precision variance, = 56.98; 95% CI (6.104, 339.0) and σ_u = 0.291. In other words, the small value of standard deviation,σ_u = 0.291 of spatial structured random effects, which means that the neighbours are not independent. The spatial heterogeneity component of variation in the BYM model has precision variance, = 330.60, 95%CI (23.84, 2101) and σ_v = 0.099. From the BYM model analysis, one can deduce the proportion of the variation that is due to clustering as = 69.06% and the proportion of variability attributed to the heterogeneity random effect is 1 -α = 30.93%.

The results revealed further that the geographic patterns of the under-five mortality prevalence in Nigeria exhibit more clustering than the spatial heterogeneity variation, as evidenced from the estimates. The geographic pattern of variation of the under-five mortality can be attributed to clustering from the exposure to local environmental factors, underlying ecological indices or severity of poverty index at the community level.

Furthermore, the risk factors are presented along with posterior statistics in Table 4. The results revealed that the estimated intercept, relative risk effect of the models are: PLN β ₀ = -0.137; 95%CI (-0.209, -0.075), CAR: β ₀ = -0.137, (-0.182, -0.092), and BYM model: β ₀ = -0:138, 95% CI (-0.200, -0.080). These risk effects are significantly different from zero and negative. These models (CAR and BYM) consolidate the result of the UH model that indicates the overall child mortality risk. A negative coefficient intercept indicates a decreasing relative risk of childhood mortality by keeping the (fixed covariates) determinant factors of under-five mortality constant. The household poverty variables are significant and positive for all the models (UH, CAR and the BYM) with these parameter estimates UH: 1.653, 95% CI (0.773, 2.491), CAR: 2.088 95%CI (1.088, 3.165), BYM: 2.003, 95%CI (1.101, 3.006). The results showed that the household poverty level would increase the relative mortality risks among the children who belong to the most economically deprived households. Other covariates in the model were not significant for the childhood mortality. However, the children who suffered from diarrhoea and who used unhygienic toilets/ sanitation had a higher tendency to die before reaching the age of five (i.e. positive association with the under- five child mortality), although they were not significant in this case. Children from mothers who used solid fuels (such as charcoal, coal, wood or agricultural residues) for cooking and drank from unprotected water are negatively insignificant.

Table 5 presents the results of the conditional autoregressive (CAR) model with the classification of the states according to the relative risk (RR) value of childhood mortality prevalence and significance probability (RR > 1) for UH model. The geographical variation in the relative risk values range from 0.438 to 1.910. The relative risk above 1(RR >1.000) indicates that the under-five mortality prevalence are higher in those states. The lowest estimated risk value occurred in Osun state: 0.0.572 (0.438, 0.714) and highest risk was recorded in Zamfara:1.680 (1.479 to 1.910). In the risk map displayed in Fig. (3), the geographical variation ranges from 0.420 to 1.922. Out of the 37 districts, the BYM model classified six (6) states as having a high relative risk of under- five mortality (RR >1.000). The relative risk value for the mortality ranges from the lowest Osun state: 0.566(0.42, 0.73) to the highest risk in Zamfara: 1.696 (1.491, 1.992). These showed that six (6) states had a relative risk significantly above 1 Table 6.

The probability risk map displayed in Fig. (4) represents smooth map of mortality for the BYM model and the states with relative risk value greater than 1. This is considered as an indication of a lower prevalence of under-five child mortality detected in the south western states and a high prevalence of mortality was found in the northern regions of the country.

Fig. (3). Relative Risk of childhood mortality prevalence and significance probability (RR > 1.000) for the CAR model based on 2013 Nigeria DHS.

Fig. (4). Relative Risk of childhood mortality prevalence and significance probability (RR > 1.000) for the BYM model based on 2013 Nigeria DHS.