Sub-national levels and trends in contraceptive prevalence, unmet need, and demand for family planning in Nigeria with survey uncertainty

Background Ambitious global goals have been established to provide universal access to affordable modern contraceptive methods. To measure progress toward such goals in populous countries like Nigeria, it’s essential to characterize the current levels and trends of family planning (FP) indicators such as unmet need and modern contraceptive prevalence rates (mCPR). Moreover, the substantial heterogeneity across Nigeria and scale of programmatic implementation requires a sub-national resolution of these FP indicators. The aim of this study is to estimate the levels and trends of FP indicators at a subnational scale in Nigeria utilizing all available data and accounting for survey design and uncertainty. Methods We utilized all available cross-sectional survey data from Nigeria including the Demographic and Health Surveys, Multiple Indicator Cluster Surveys, National Nutrition and Health Surveys, and Performance, Monitoring, and Accountability 2020. We developed a hierarchical Bayesian model that incorporates all of the individual level data from each survey instrument, accounts for survey uncertainty, leverages spatio-temporal smoothing, and produces probabilistic estimates with uncertainty intervals. Results We estimate that overall rates and trends of mCPR and unmet need have remained low in Nigeria: the average annual rate of change for mCPR by state is 0.5% (0.4%,0.6%) from 2012-2017. Unmet need by age-parity demographic groups varied significantly across Nigeria; parous women express much higher rates of unmet need than nulliparous women. Conclusions Understanding the estimates and trends of FP indicators at a subnational resolution in Nigeria is integral to inform programmatic decision-making. We identify age-parity-state subgroups with large rates of unmet need. We also find conflicting trends by survey instrument across a number of states. Our model-based estimates highlight these inconsistencies, attempt to reconcile the direct survey estimates, and provide uncertainty intervals to enable interpretation of model and survey estimates for decision-making.


Input Data
In this section we describe the thirteen household surveys that were analyzed for the article 'Sub-national levels and trends in contraceptive prevalence, unmet need, and demand for family planning in Nigeria with survey uncertainty.'

Demographic and Health Surveys
We processed the four available Demographic and Health Surveys (DHS) for Nigeria. The 1990 DHS selected 34 households from each of the 299 enumeration areas (EAs), which were stratified by urban/rural status and selected via probability proportional to size (pps), and was designed to produce national and regional (4 regions at that time) estimates [1]. A stratified two-stage sampling design was used for the 2003 DHS to provide estimates of health indicators at the national level, the six health regions, and by urban/rural status. Households were select from a complete listing of households within the 365 selected EAs [2]. Similarly, the 2008 and 2013 DHS employed stratified multi-stage sampling designs consisting of 888 EAs and 904 EAs, respectively [3,4] and were intended to be representative at the national, regional, and state level. An additional survey was conducted in 1999, however the DHS Program was not centrally involved with the study and does not distribute the collected data.
To calculate unmet need we relied on the recoded definitions used by the United Nations Population Division (UNPD), which applied the revised DHS unmet need definition [5] to previous surveys. To compute the revised unmet need in our analysis, we relied on the UNPD script available at https: //github.com/PhilUeff/PDU_FP-Indicators.
To ensure state-level values reflect current geographical divisions, we mapped the GPS coordinates associated with the sampled clusters from each survey onto the 2013 Nigeria DHS shapefile and used their corresponding state names. Each DHS survey contained a few clusters without recorded GPS coordinates (Table 1). For these, we assigned the clusters to the states that they were recorded as belonging to at the time. The assignment is ambiguous for one DHS 1990 cluster because the corresponding recorded state has split between the time of the survey and 2013; this is noted in Table 1. Table 2 compares national all-women estimates for contraceptive prevalence rate (CPR), modern contraceptive preva-lence rate (mCPR), total unmet need (UnN), and demand satisfied by modern methods (DSM), as computed by the DHS survey official reports, the UN recode, and our analysis. While not presented, traditional contraceptive prevalence rate can be computed as CPR -mCPR.
Contraceptive prevalences match among all sources. The differences in unmet need between the DHS official reports and the UN computed values prior to 2012 reflect the revised unmet need calculation. Dashes indicate estimates that were not reported by the DHS.  [6]. A twostage sampling design was used for the 2011 MICS, which was designed to provide estimates and the state and national level. Forty EAs were selected with equal probability within each state and all 1,480 EAs were included in field work [7]. The 2016-17 MICS was designed to provide estimates at the national and state level and relied on a two-stage sampling design. Sixty EAs were selected from each state, except for Lagos and Kano which were over sampled at 120 EAs. Of the 2,340 sampled EAs only 2,239 were included in field work because 101 EAs from Borno, Yobe, and Adamawa states were excluded due to security concerns [8].
The UNPD has recoded the 2007 and 2011 MICS surveys, but at the time of this writing, had not yet recoded the 2016 survey. We used the UNPD recode scripts from https:// github.com/PhilUeff/PDU_FP-Indicators to process the 2007 and 2011 MICS survey data, then computed the survey estimates as with the DHS surveys. When we used the UN code on the data, we found a small difference in the 2011 estimated unmet need from the value reported by the UNPD. For the 2016 MICS survey, we checked that the variables were encoded in the same form as the 2011 MICS survey. The unavailable estimates from the MICS official reports and 2016 UNPD recode are indicated with dashes. Finally, certain questions needed to calculate unmet need were not asked in the 2007 MICS survey. As a result, unmet need and demand satisfied for this survey are marked with 'NA'. The MICS official reports only present estimates over currently-married/in-union women. Table 3 displays the estimates reported by MICS, calculated by the UN, and calculated via our direct analysis of the microdata.

Performance Monitoring and Accountability 2020
We used four rounds of the Performance Monitoring and Accountability 2020 (PMA2020) surveys conducted in Nigeria. We obtained the data for the 2014-2016 PMA2020 rounds [9,10,11] from IPUMS [12], which had recoded the PMA2020 surveys. We computed survey estimates using this data. The 2017 PMA2020 survey [13] had not been recoded at the time of this writing by IPUMS, so this survey was recoded manually. The 2014 and 2015 rounds were designed to provide state level estimates for Kaduna and Lagos states and relied on a two-stage sampling design. The 2016 and 2017 rounds were designed to provide national estimates relying on three-stage sampling design. At the first stage seven states were selected and then within each stage two-stage sampling was implemented. Post stratification was used to generate national estimates [14]. Table 4 provides the point estimates reported by PMA2020 and those calcuated in our analysis of the microdata.

National Nutrion and Health Surveys
We analyzed the microdata from two National Nutrition and Health Surveys (NNHS) from 2014 and 2015. Both surveys were conducted using the Standardized Monitoring and Assessment of Relief and Transition (SMART) methods and relied on two-stage cluster and designed to provide estimates at the national and state level [15,16]. Unfortunately, the NNHS does not include the variables required to calculate unmet need nor parity, so the NNHS data was only included in the estimates of mCPR for all women. Additionally, the survey reports highlight that in 9 Local Governmental Areas (LGAs) were excluded for security reasons and thus results for Borno are not representative of the whole state.

Modeling family planning indicators
Family planning indicators in Nigeria are estimated from multiple complex surveys, sometimes in the same year and can be quite noisy between surveys at the state level. In an effort to understand the underlying population rates from which these survey observations were drawn, we are adapting a previously developed two-step space-time smoothing approach that acknowledge complex sampling designs to model the family planning indicators.
The first step of our approach requires estimating the state (i), year (t), and survey (s) specific direct estimates of proportions and corresponding variance via the Hajek estimator [17] p its = j x j,its w j,its / j w j,its where x j,its is the binary indicator for modern contraception, traditional contraception, unment need, or demand satisfied for woman j, in sampled in area i, time t, and survey s, and w j,its is her corresponding sampling weight.
In the second step we use a three-stage Bayesian hierarchichal model. In the first stage we rely on the working likelihood based on the asymptotic distribution, Y its ∼ N η its , V its,DES  where Y its = logit[ p its ] and V its,DES is the design-based variance of Y its . This pseudo-likelihood for a binary variable was initially developed to estimate zipcode-level smoking prevalence in Washington state [18] relying on a single survey instrument. This model was then adapted for the complex indicator of under-five child mortality [19] and extended the classic inseperable spatio-time models of [20,21] to include surveyspecific random effects.
At the second stage we assume where µ is a shared interecept, the temporally structured random effects α are assigned second order random walk priors [22,21], the spatial random effects θ are assigned either the Besag, York, and Mollie (BYM) [23] or the scaled BYM [24], the space-time interactions δ are assinged the 'type II' temporally structured interaction [21], and the independent random effects γ, ν, φ, and ψ are assigned independent mean-zero Normal distributions. Note, if θ ∼BYM then θ i = U i + V i where U is assigned an intrinsic conditional autoregression prior (ICAR) [25] and V are assigned an independent mean-zero Normal prior. At the third stage of the model we assume a default diffuse prior for µ and Gamma(1, 5e − 5) priors were assigned to, σ −2 α , σ −2 θ (or to both inverse variances for standard BYM), σ −2 γ , σ −2 δ , σ −2 ν , σ −2 φ , and σ −2 ψ . Finally, to generate our estimates of the underlying rates we adopt the approach of Mercer (2015) and draw k=1,...1000 posterior samples for all parameters and calculate The posterior median is used as the estimate and the uncertainty intervals are defined by the 2.5% and 97.5% quantiles.

Model Selection
The twelve possible models considered for mCPR, traditional contraceptive prevalence, unmet need, and demand satisfied are described in Table 5. Each model was fit for mCPR, traditional contraceptive prevalence, unment need and demand satisfied for all women and the four age and parity combinations. For each model fit we calculated the sum of the log conditional predictive ordinate (LCPO) [26], the deviance information criteria (DIC) [27], and the Watanabe-Akaike information criterion (WAIC) [28]. The results for mCPR and unment need for all women are shown in Table 6 and Table 7, respectively with with the lowest DIC and WAIC and highest LCPO in bold. Table 8 shows the selected model for all outcomes and ageparity groups. Table 5: Random effects models considered for time t, state i, and survey s. All models contain an intercept µ, temporally independent γ t ∼ N(0, σ 2 γ ), and temporally structured effect α t ∼ RW2(σ 2 α ). The survey and survey interaction random effects were assigned mean-zero Normal priors. In models designated with 'a' θ i ∼ BY M(σ 2 θ ) and in models designated 'b' θ i ∼ BY M2(σ 2 θ ).
Model outcome

Decomposition of Variance
To assess the relative importance of the individual random effects in explaining the variance in the observed data we calculated the percent of the total variance described by each random   Table 9 shows the percent of the variance described by each random effect where the total variance for model 6 is and for outcomes and age-parity groups using model 4, the σ 2 φ term is removed. Demand Satisfied 4b

Trends in family planning indicators in PMA2020 states
The observed data and smoothed estimates for mCPR, traditional contraceptive prevalence, unment need, and demand satistified in the seven states (Anambra, Kaduna, Kano, Lagos, Nasarawa, Rivers, and Taraba) are shown in Figures 1 through 7.
6. Maps of family planning indicators for all age-parity subgroups Figures 8,9,10, and 11 show the 2017 state-level estimates or the four age-parity subgroups for mCPR, unmet need, traditional contraceptive prevalence, and demand satisfied, respectively.

Small area estimation at the second administrative level
The health policies are often enacted at the second adminstrative level, called local government areas (LGAs) in Nige-ria, and ideally we would provide estimates at the LGA-level. Unfortunately, MICS, DHS, PMA2020, and NNHS do not inlcude the LGA locations for the sampled clusters. However, the DHS have provided jittered [29] GPS locations which can be assigned LGA membership. In this section we investigate and compare the precision attained at the state-and LGA-level. We apply the methods described in Section 2, removing the surveyrelated random effects, to the DHS data from 1990. 2003, 2008, and 2013.
The left panel of Figure 12 display the state-level estimates and 95% credible intervals for mCPR and the right panel displays the same, but at the LGA-level. From the figure it is clear that in many of the LGAs the uncertainty is intolerably large. This is not a surprising result, as many LGAs have few or no sampled clusters. By contrast differences between states and within states over time can be reliably estimated; see Figure 12 in the Appendix and Figure 4 in the main article, respectively. With currently available data it is not possible to generate substate level estimates for FP indicators with reasonable precision. FP estimates at the LGA rely exclusively on DHS data and are thus not timely nor precise. Table 9: Percent of total variance described by each random effect, where the interpretation is σ 2 θ for spatial, σ 2 γ for independent temporal, σ 2 α for structured (RW2) temporal, σ 2 δ for space-time interaction, σ 2 ν for surveys, σ 2 φ for survey-space interactions, and σ 2 ψ for survey-time interaction.