Multilevel analysis of predictors of multiple indicators of childhood vaccination in Nigeria

Background Substantial inequalities exist in childhood vaccination coverage levels. To increase vaccine uptake, factors that predict vaccination coverage in children should be identified and addressed. Methods Using data from the 2018 Nigeria Demographic and Health Survey and geospatial data sets, we fitted Bayesian multilevel binomial and multinomial logistic regression models to analyse independent predictors of three vaccination outcomes: receipt of the first dose of Pentavalent vaccine (containing diphtheria-tetanus-pertussis, Hemophilus influenzae type B and Hepatitis B vaccines) (PENTA1) (n = 6059) and receipt of the third dose having received the first (PENTA3/1) (n = 3937) in children aged 12-23 months, and receipt of measles vaccine (MV) (n = 11839) among children aged 12-35 months. Results Factors associated with vaccination were broadly similar for documented versus recall evidence of vaccination. Based on any evidence of vaccination, we found that health card/document ownership, receipt of vitamin A and maternal educational level were significantly associated with each outcome. Although the coverage of each vaccine dose was higher in urban than rural areas, urban residence was not significant in multivariable analyses that included travel time. Indicators relating to socio-economic status, as well as ethnic group, skilled birth attendance, lower travel time to the nearest health facility and problems seeking health care were significantly associated with both PENTA1 and MV. Maternal religion was related to PENTA1 and PENTA3/1 and maternal age related to MV and PENTA3/1; other significant variables were associated with one outcome each. Substantial residual community level variances in different strata were observed in the fitted models for each outcome. Conclusion Our analysis has highlighted socio-demographic and health care access factors that affect not only beginning but completing the vaccination series in Nigeria. Other factors not measured by the DHS such as health service quality and community attitudes should also be investigated and addressed to tackle inequities in coverage. The study used data from the nationally representative cross-sectional 2018 Nigeria Demographic and Health Survey (NDHS) which was implemented by the National Population Commission (NPC) with technical assistance provided by Inner City Fund (ICF) International through The Measure DHS Program [23]. Data collection was conducted from 14 August to 29 December 2018 with pre-test conducted from 30 April to 20 May 2018. The survey utilised a stratified two-stage sampling approach. In the first stage, 1389 enumeration areas (EAs) or clusters were selected as the primary sampling units. At the second stage, 40427 households were selected. Stratification was achieved by separating each administrative level 1 area (i.e., the 36 states and the Federal Capital Territory) into urban and rural areas, and samples were selected independently within each stratum. Detailed information on the description of the methods employed in the study is available elsewhere [23]. The primary outcome variables/indicators in this study are receipt of PENTA1 vaccine (n = 6059) and receipt of PENTA3 having received PENTA1 (PENTA3/1—the converse of dropout, n = 3937) among children aged 12–23 months, and receipt of MV (n = 11839) among children aged 12–35 months. For each of the three indicators, we assessed the binomial outcome: any evidence of vaccination versus no evidence of vaccination. For each of PENTA1 and MV, we additionally assessed the multinomial outcome: no evidence of vaccination, card invalid/history evidence of vaccination and card valid vaccination, to assess potential variation in associations that could be caused by misclassifying vaccination status when a verbal history of vaccination is accepted [24] (see Table A in S1 File). Valid and invalid vaccine doses were defined according to WHO guidance [25]. The study considered covariates at the child-, household- and community-levels. The selection of these covariates was informed by literature on the predictors of vaccination coverage or of child health outcomes in general, expert opinion, and availability in the 2018 NDHS or other sources, as detailed in Table A in S1 File [6,11,17,21,26–35]. The study however excluded some pre-selected DHS covariates due to missingness or multicollinearity (see supplementary materials). These variables include preceding birth interval, antenatal and postnatal care, maternal receipt of tetanus toxoid vaccination (these variables had >15% missing cases for MV and were excluded from the analyses to make the results comparable across all the three indicators), maternal decision-making (whether mother decides health care, visits, and purchases)—and region of residence. Also, among similar covariates (e.g., mother’s occupation and employment status), one was selected for inclusion in our model based on the literature and expert opinion. The geospatial covariates retained in this study include travel time to health facility, enhanced vegetation Index (EVI) and livestock density index. Tertiles of the distribution of these covariate were used to allow similar number of observations for each tertile. We present further description about these covariates and their relevance in Fig A and Table B in S1 File [18–22,29,36,37]. Other geospatial covariates considered included distance to conflict locations, maximum number of conflicts, night light intensity, annual aridity index, maximum temperature, annual precipitation/rainfall, proximity to national borders, water, protected areas, and slope [18–20] were excluded from final models due to multicollinearity (as was EVI) or non-significance after adjusting for DHS variables. Cross-tabulations and single-level logistic regression analysis. We tabulated each outcome against each of the selected covariates separately to explore relationships and used Chi-squared tests to determine the significance of the associations. We then fitted frequentist single level simple logistic regression models to obtain the corresponding crude odds ratios (cORs) and associated 95% confidence intervals (CI). These results were later compared with results from the multiple multilevel analyses to determine changes in statistical significance and direction of effects. Multiple multilevel binomial regression analysis (any evidence of vaccination) and interaction effects. We fitted Bayesian multilevel [38,39] binomial regression models to estimate adjusted odds ratios (aORs) and corresponding 95% credible intervals, accounting for the hierarchical structure of the data (child/household, community, and stratum levels) and, intrinsically, the survey design (clustering and stratification) through the last two hierarchies (see Fig B in S1 File). A detailed description of the model is included in the supplementary information (S1 File). We investigated whether child/household covariate effects could be modified by the geospatial covariates by introducing interaction terms between both sets of covariates. To incorporate the interaction terms, we first fitted the main effects model using both DHS and geospatial covariates and then introduced the interactions between selected DHS and geospatial covariates sequentially, retaining only those that were significant in the final model. Multiple multilevel multinomial regression analysis (vaccination according to source of evidence). The study also employed a Bayesian multinomial multivariable multilevel modelling approach to estimate the adjusted relative risk (aRR) and associated 95% credible intervals for covariates significantly associated with PENTA1 and MV using the multinomial outcomes defined previously. No interaction terms were considered for the multinomial analyses due to model complexities and non-convergence challenges. We computed summary measures [30,40] of the amount of residual variation attributable to the hierarchies in the binomial models. These included the variance partitioning coefficient (VPC) which measures residual variation between clusters/communities in different strata; the median odds ratio (MOR) which quantifies residual community level variation in the likelihood of vaccination on the odds ratio scale, and the percent change in variance (PCV) which measures change in residual variation due to the inclusion of covariates in the models [41–45]. Also, although prediction was not the main goal, we evaluated the discriminatory or predictive power of the fitted models using the area under the receiver operating characteristic (AUROC) curve (see supplementary materials for details). All analyses were implemented in Stata version 16 [46], MLwiN version 3.05 [47], and the R programming language version 4.0.3 [48]. Additionally, we used the runmlwin [49] program to run the MLwiN multilevel modelling software from within Stata. We utilized MCMC algorithms with a burn-in length of 1000, a monitoring chain length of 60000, and thinning of 20. Convergence of the MCMC chains was assessed via visual inspection of the trace and autocorrelation plots of the parameters. Ethical approval was obtained from the Nigeria National Health Research Ethics Committee and the ICF Institutional Review Board for the main NDHS [23], and from the Ethics and Research Governance, University of Southampton, United Kingdom. Written informed consent was obtained from all study respondents. However, the data were analysed anonymously in the present study. All methods were performed in accordance with the relevant guidelines and regulations.