Predictive spatial risk model of poliovirus to aid prioritization and hasten eradication in Nigeria
 Alexander M UpfillBrown^{1}Email author,
 Hil M Lyons^{1},
 Muhammad A Pate^{2},
 Faisal Shuaib^{3, 4, 5},
 Shahzad Baig^{4, 6},
 Hao Hu^{1},
 Philip A Eckhoff^{1} and
 Guillaume ChabotCouture^{1}
DOI: 10.1186/174170151292
© UpfillBrown et al.; licensee BioMed Central Ltd. 2014
Received: 4 March 2014
Accepted: 9 May 2014
Published: 4 June 2014
Abstract
Background
One of the challenges facing the Global Polio Eradication Initiative is efficiently directing limited resources, such as specially trained personnel, community outreach activities, and satellite vaccinator tracking, to the most atrisk areas to maximize the impact of interventions. A validated predictive model of wild poliovirus circulation would greatly inform prioritization efforts by accurately forecasting areas at greatest risk, thus enabling the greatest effect of program interventions.
Methods
Using Nigerian acute flaccid paralysis surveillance data from 20042013, we developed a spatial hierarchical Poisson hurdle model fitted within a Bayesian framework to study historical polio caseload patterns and forecast future circulation of type 1 and 3 wild poliovirus within districts in Nigeria. A Bayesian temporal smoothing model was applied to address data sparsity underlying estimates of covariates at the district level.
Results
We find that calculated vaccinederived population immunity is significantly negatively associated with the probability and number of wild poliovirus case(s) within a district. Recent case information is significantly positively associated with probability of a case, but not the number of cases. We used lagged indicators and coefficients from the fitted models to forecast reported cases in the subsequent sixmonth periods. Over the past three years, the average predictive ability is 86 ± 2% and 85 ± 4% for wild poliovirus type 1 and 3, respectively. Interestingly, the predictive accuracy of historical transmission patterns alone is equivalent (86 ± 2% and 84 ± 4% for type 1 and 3, respectively). We calculate uncertainty in risk ranking to inform assessments of changes in rank between time periods.
Conclusions
The model developed in this study successfully predicts districts at risk for future wild poliovirus cases in Nigeria. The highest predicted district risk was 12.8 WPV1 cases in 2006, while the lowest district risk was 0.001 WPV1 cases in 2013. Model results have been used to direct the allocation of many different interventions, including political and religious advocacy visits. This modeling approach could be applied to other vaccine preventable diseases for use in other control and elimination programs.
Keywords
Polio eradication Spatial epidemiology Risk modeling Disease mappingBackground
Since the commitment by the World Health Assembly to eradicate poliovirus in 1988, the disease burden has dramatically declined by more than 99% to 223 cases reported in 2012 [1, 2]. Yet, in the face of growing financial and political investments, polio remains endemic in Nigeria, Pakistan, and Afghanistan and has been repeatedly exported to other previously poliofree countries—leading the 65th World Health Assembly to declare polio eradication a “programmatic emergency for global public health” in 2012 [3].
Though substantial, the resources of the Global Polio Eradication Initiative (GPEI), including vaccines, specially trained personnel, and social mobilization campaigns, are limited and must be targeted to highrisk areas within endemic countries in order to maximize impact [4]. The GPEI supports the surveillance and survey collection designed to monitor program performance in these areas [5]. However, converting these data sources into operationally useful and scientifically valid measures of future risk can be challenging.
Furthermore, before 2010, wild poliovirus (WPV) epidemics affected large portions of northern Nigeria, and 80 to 90 different districts were reporting cases within a sixmonth period. Thus, the need for prioritization was low, as cases were scattered across the whole of northern Nigeria and outbreak response was the focus of the problem. The lull in reported WPV cases in 2010 coupled with the small number of WPV “infected” districts scattered across northern Nigeria in 2011 suggested that the epidemic in Nigeria had become more focused; a much smaller number of districts were contributing to WPV transmission in Nigeria. Increased focus by the program was needed to anticipate which highrisk districts were at highest risk for future cases in order to prevent further transmission. Programmatically, the need to identify the districts at highest risk was intensified by a planned surge in technical and administrative capacity, supported by the World Health Organization (WHO).
Currently, one systematic method for risk assessment at the district level has been described in the literature: the WHO’s method of regional polio risk assessment, which uses a weighted linear combination of available indicators [6]. This approach is currently used by all WHO regional offices for supplementary immunization activity (SIA) campaign planning and outbreak risk assessment. As Lowther and colleagues recognize [6], this approach has limitations: these weights are based on expert opinion, not statistical modeling, and its historical predictive accuracy has yet to be demonstrated.
Previous work has demonstrated the accuracy of WPV outbreak prediction using hierarchical statistical modeling in the African region at the country/province level [7]. This model incorporated measures of connectedness between areas and population immunity and demonstrated good predictive accuracy. However, in endemic areas, operations are carried out at the district administrative level or below; at this spatial scale, underlying causal factors, such as migration, are poorly measured.
To better support modeling of spatial heterogeneity and correlation, recent work focusing on other infectious diseases has increasingly applied Bayesian spatial modeling methods [8]. These methods have been used to map infection prevalence in unmeasured areas for malaria [9, 10], schistosomiasis [11–15], soiltransmitted helminths [16–18], and filariasis [19].
Only recently in infectious disease research have these methods been used to aid in forecasting disease prevalence or incidence in a future time period [20, 21]. Though spatial models of infectious disease describe reported case counts, these models rarely apply zeroinflated Poisson and Poisson hurdle models to adjust for excess zeros [22]; however, these models, with and without a spatial component, have been used to analyze other ecological and public health processes [23–26]. Rarely have zeroinflated models been used to forecast future case counts.
We present a predictive spatial hierarchical Poisson hurdle model for serotype 1 and serotype 3 WPV (WPV1 and WPV3) transmission. This model is currently incorporated in districtlevel prioritization planning in Nigeria. Using this modeling framework, we identify the most important risk factors of the presence of WPV and the historical number of WPV cases, and also highlight areas in which these indicators are weakly informative. Due to the sparsity of the data underlying the estimates of districtlevel covariates, we employ a twopart methodology in which these estimates are smoothed using a temporal hierarchical model prior to their use as inputs to the spatial Poisson hurdle models. We examine this model’s historical ability to forecast districts that will report one or more WPV cases six months in the future.
Methods
Description of the data
Individually reported dose histories, vaccinespecific campaign exposures, a probability model for campaign participation, and vaccinespecific OPV efficacies based on AFP surveillance data [33] were used to calculate the probability that an individual is protected against paralysis from each WPV serotype (Additional file 1, Eq. S1) [34]. Individual campaign exposure was defined as the number of SIA campaigns occurring between birth and the date of paralysis onset, based on the historical SIA calendar in the district.
Both polio AFP and NPAFP cases were linked to districts in Nigeria and binned into sixmonth periods: December through May or June through November. In the following, the former period is referred to as the first half of a year and the latter as the second half of a year; for example, the first half of 2008 refers to the period from December 2007 through May 2008.
Districtlevel OPVinduced population immunity was taken to be the mean calculated OPVinduced immunity of all NPAFP cases in children under five yearsold occurring in a district within a sixmonth period. District zerodose fraction was computed as the proportion of NPAFP cases for children 0 to 59 months old for which there were reported zero OPV doses in a district within a given sixmonth period.
Recent caseload is equal to the total number of WPV1 or WPV3 cases in a district in the previous six months. Recent neighboring caseload, also serospecific, is the total number of WPV1 or WPV3 cases in bordering districts in the previous six months. We used AFP and campaign data from June 2004 through May 2013.
Statistical methods
WPV1 and WPV3 caseloads were modeled using a Poisson hurdle model: a twopart model consisting of a Bernoulli component that models the probability of reporting one or more WPV cases and a truncated Poisson portion that models the number of WPV cases in infected districts (districts reporting one or more WPV cases; see Additional file 1, Eq. S8, S9, S10) [35, 36]. For each component, we included spatially and nonspatially dependent random effects to account for unobserved spatially correlated risk factors [37, 38]. Together, these terms enable both local and global borrowing of information, respectively. A bivariate conditional autoregressive (CAR) prior was used to model spatial random effects, and a bivariate normal prior was used to model nonspatial random effects, as we expect areas with greater rates of WPV transmission to report larger case counts [26].
A Bayesian estimation approach was utilized, and models were fit using WinBUGS 1.4.3 [39] and the R2WinBUGS and CODA packages [40, 41] in R [42]. Diffuse priors and hyperpriors were used for the regression coefficients and the random effect precisions. The convergence of Markov chain Monte Carlo (MCMC) models was checked through visual inspection of chains and the GelmanRubin potential scale reduction factor [43, 44]. The selected model was determined using an iterative subtractive approach to arrive at the most parsimonious model with the lowest deviance information criterion (DIC), a well described method for quantifying the quality of a fit [25, 45] (see Additional file 1 for full model details).
The accuracy of model predictions was estimated using prospective sampling: predictions were compared to future data excluded from the initial analysis [46]. The selected model was applied historically using data smoothed over a specified time period; the historical caseload in an upcoming window was predicted using the model coefficients estimated using data from prior time periods and covariates from the most recent time period^{b}. The predictive power of the model was assessed using receiver operator characteristic (ROC) analysis using the “pROC” package in R [47]. The empirical area under the ROC curve (AUC) was used to quantify the accuracy of model predictions of infected districts (districts reporting >1 WPV1 case) [48]. As a rule of thumb, an AUC of 0.5 to 0.7 indicates marginal predictive power, an AUC of 0.7 to 0.9 indicates moderate predictive power, and an AUC above 0.9 indicates strong predictive power [49].
Uncertainties in the predicted district caseloads and district ranks were estimated using 1,000 samples from the posterior distribution of parameters. District rank was determined by ordering districts according to predicted WPV1 risk with 1 being most at risk and 774 being least at risk. As districts are prioritized based on rank above a certain threshold, we used samples from the posterior distribution of district ranks to estimate the probability of a district belonging to the N districts of highest rank.
Results
Influence of data smoothing model
Data smoothing resulted in substantial changes in the estimated regression coefficients; the impact of the smoothing model is highlighted in Figure 3. The smoothing model shrinks extreme district estimates from the first half of 2013 towards the mean (Figure 3A, B). We find the calculated population immunity to be low in southern districts (see Figure 2 for map of Nigeria); see Discussion for possible explanations of this phenomenon. The reduction of rapid fluctuations in values over time is highlighted in Figure 3C and D.
The smoothing effectively captures major trends in the data while reducing the amount of variability between adjacent time points: fewer than 1.5% of changes in the smoothed calculated immunity between adjacent time periods are greater than 0.1 (none are greater than 0.15), compared to 49.7% of changes greater than 0.1 for the empirically calculated immunity. Additionally, fewer than 3% of relative immunity decreases between adjacent time periods were greater than 20% in the smoothed immunity model, while 43% of immunity decreases were greater than 20% for empirical population immunity estimates.
Factors associated with WPV1 transmission
Covariate estimates for models using data through May 2013
Component  Variable  Full  Selected  Null 

Poisson  Intercept  11.21  11.51  12.68 
(11.88, 10.51)  (11.83, 11.17)  (13.03, 12.37)  
Population Immunity  2.50  1.80  
(3.32, 1.66)  (2.20, 1.44)  
ZeroDose Fraction  0.72  
(1.47, 0.02)  
Density  0.02  
(0.07, 0.11)  
Sqrt Recent Cases  0.05  
(0.11, 0.02)  
Sqrt Neighboring Recent Cases  0.02  
(0.02, 0.05)  
Bernoulli  Intercept  0.91  1.11  3.60 
(1.64, 0.16)  (1.42, 0.80)  (3.80, 3.43)  
Population Immunity  4.86  4.74  
(5.69, 4.01)  (5.27, 4.21)  
ZeroDose Fraction  0.19  
(1.12, 0.75)  
Density  0.02  
(0.11, 0.08)  
Sqrt Recent Cases  0.19  0.19  
(0.07, 0.31)  (0.07, 0.32)  
Sqrt Neighboring Recent Cases  0.17  0.17  
(0.10, 0.23)  (0.10, 0.23)  
DIC  9263.6  9266.8  10121.3  
N  13932  13932  13932  
Variance  Bernoulli CAR^{a}  1.51  1.50  1.48 
Components  (1.06, 2.08)  (1.04, 2.04)  (1.08, 1.98)  
Bernoulli Ind^{b}  0.10  0.10  0.09  
(0.06, 0.17)  (0.05, 0.18)  (0.05, 0.15)  
Poisson CAR^{c}  0.49  0.46  0.67  
(0.24, 0.86)  (0.22, 0.79)  (0.39, 1.05)  
Poisson Ind^{d}  0.14  0.13  0.15  
(0.08, 0.23)  (0.07, 0.22)  (0.08, 0.23)  
Covariance CAR^{e}  0.10  0.07  0.42  
(0.23, 0.46)  (0.23, 0.39)  (0.12, 0.77)  
Covariance Ind^{f}  0.01  0.01  0.01  
(0.05, 0.04)  (0.07, 0.04)  (0.06, 0.04) 
In models of WPV3 in Nigeria, only the OPVbased population immunity was significantly associated with the number of cases in a district in the selected model, while population immunity and recent cases in neighboring districts were associated with the presence of WPV3 cases in a district (see Additional file 1: Table S2).
A null model without covariates is included in Table 1 as a reference; this model only makes use of population and historical outcome information, through spatial and nonspatial random effects. Interestingly, the random effects variance estimates are only marginally higher in this model than in the selected model, though the covariance between spatial random effects in the Poisson and binomial portions is significantly greater than zero. Fixed effect covariates therefore account for a fairly small proportion of the model variance captured by the random effects alone in the null model. This suggests that the timevarying covariates contain additional information not captured by random effects alone.
Predicting WPV1 circulation in districts in Nigeria
Coefficient estimates from the selected model applied historically
Component  Variable  2012  2011  2010  2009  2008 

Poisson  Intercept  11.42  11.50  11.60  11.80  12.20 
(11.76, 11.09)  (11.87, 11.17)  (11.96, 11.24)  (12.22, 11.39)  (12.79, 11.69)  
Population  1.95  1.80  1.60  1.62  1.28  
Immunity  (2.40, 1.51)  (2.30, 1.31)  (2.09, 1.10)  (2.14, 1.11)  (1.92, 0.61)  
Binomial  Intercept  0.81  0.54  0.75  1.43  1.89 
(1.14, 0.46)  (0.91, 0.17)  (1.16, 0.35)  (1.81, 1.02)  (2.42, 1.36)  
Population  5.24  5.69  5.03  3.42  3.17  
Immunity  (5.82, 4.65)  (6.34, 5.07)  (5.71, 4.36)  (4.18, 2.71)  (4.19, 2.20)  
Sqrt Recent  0.19  0.18  0.18  0.13  0.18  
Cases  (0.06, 0.32)  (0.05, 0.30)  (0.05, 0.31)  (0.01, 0.26)  (0.04, 0.33)  
Sqrt Neighboring  0.17  0.18  0.16  0.06  0.02  
Recent Cases  (0.10, 0.24)  (0.11, 0.25)  (0.09, 0.23)  (0.01, 0.13)  (0.07, 0.10)  
DIC  8732.4  8306.9  8107.3  7798.5  6200.7  
N  12384  10836  9288  7740  6192  
Variance  Bernoulli CAR^{a}  1.47  1.49  1.41  1.59  1.81 
Components  (1.03, 2.03)  (1.02, 2.06)  (0.97, 1.98)  (1.08, 2.23)  (1.26, 2.52)  
Bernoulli Ind^{b}  0.11  0.11  0.11  0.11  0.12  
(0.06, 0.19)  (0.06, 0.19)  (0.06, 0.19)  (0.06, 0.19)  (0.06, 0.22)  
Poisson CAR^{c}  0.47  0.49  0.50  0.51  0.61  
(0.22, 0.80)  (0.24, 0.85)  (0.24, 0.86)  (0.25, 0.85)  (0.30, 1.02)  
Poisson Ind^{d}  0.14  0.14  0.14  0.14  0.18  
(0.07, 0.22)  (0.07, 0.24)  (0.08, 0.23)  (0.08, 0.23)  (0.10, 0.28)  
Covariance CAR^{e}  0.04  0.04  0.06  0.26  0.36  
(0.28, 0.34)  (0.31, 0.41)  (0.28, 0.38)  (0.10, 0.63)  (0.06, 0.76)  
Covariance Ind^{f}  0.01  0.01  0.01  0.02  0.03  
(0.07, 0.04)  (0.07, 0.04)  (0.07, 0.04)  (0.07, 0.04)  (0.10, 0.04) 
Interestingly, the selected WPV1 model did not perform better on average than a null model, suggesting that an aggregate measure of spatially smoothed historical risk is equally adept at predicting future infected districts as a model that considers population immunity, recent cases, and recent cases in neighboring districts. The mean AUC of the null WPV1 model over the last three years was 0.86 (SE=0.01). Similarly, the predictive accuracy of the null WPV3 model was comparable to that of the selected WPV3 model; the mean AUC over the last three years was 0.84 (SE=0.02).
In 2009, both the selected and null WPV1 models were very poorly predictive. This result is largely due to a flareup of WPV1 in southern Nigeria in 2009, in districts that had not previously reported a case (see Additional File 1: Figure S3E, F). The relatively poor predictive power of the model in the first half of 2010 (AUC=0.71) is associated with an uncharacteristically small number of reported WPV1 cases (two cases).As such, a useful quantity for policy makers is the expected proportion of cases contained within a set of the highest risk districts (Figure 6B). Historically in Nigeria, 85 to 90% of future cases have appeared in the 200 highest risk districts, and 55% of future cases have appeared in the 100 highest risk districts under the selected model.
Uncertainty in prioritization according to predicted caseload
Discussion
The smoothing model handles missing data and smooths over unrealistic changes in indicators, such as rapid sixmonth fluctuations in the underfive immunity estimated from sparse NPAFP samples. In the absence of any vaccination, the underfive OPVinduced population immunity should decrease by roughly 89% in a sixmonth period (see Figure 1B)^{c}. Though thresholds based on this or similar considerations are not explicit in our smoothing model, the smoothed variations between time periods are more realistic than those in empirical data. Smoothing models could be improved by limiting changes between time points to demographically constrained values.
Despite its limitations [50], the OPVbased population immunity, calculated using efficacies derived from casecontrol population studies [29, 33], is strongly associated with the presence and case count of WPV1 at the district level in Nigeria. It is the only covariate significantly associated with the number of cases in a district.
It is expected that the recent caseload is significantly associated only with the expected presence of WPV1, as recent circulation has a mixed effect on the future case count. While a large number of cases within a time period suggests potentially higher transmission (temporally and spatially) from a large infectious reservoir, WPV1 circulation will also boost natural immunity in an area, providing additional protection and reducing the expected caseload in the ensuing time period [50]. Indeed, we find that both recent caseload within a district and in surrounding districts are significantly associated with the probability of a WPV1 case, but not the number of WPV1 cases.
Although increased population density is thought to increase poliovirus transmission potential [51], we find no association in either portion of the model. Due to data limitations, population density was calculated at the district level, which may be a poor representation of the experienced population density on more functional scales of transmission, such as the population density immediately surrounding the household reporting a WPV case. A more finegrained analysis is needed to better understand the impact of population density of WPV transmission.
The ability to detect an association between zerodose fraction, an indicator of possible clustering of vulnerability within a district, and the presence and number of WPV1 cases was compromised by a strong correlation between population immunity and zerodose fraction. The collinearity between these indicators resulted in a negative association between zerodose fraction and the number of WPV1 cases in a district, while in a model without population immunity, zerodose fraction has a strong positive association with the number of cases in a district.
Although the historical predictive accuracy is high, covariates in the selected model do not substantially improve the forecasting ability of the model over the null model. In a null model, the random effects capture the mean historical frequency and number of WPV1 cases within a district. This result suggests that although dynamic indicators are associated with WPV1 transmission, historical transmission patterns are stronger predictors of future transmission than available model covariates.
The decrease in performance in 2009 is caused by a flareup in WPV1 cases in southern districts with no history of WPV1 transmission in our dataset; because the spatial random effects alone are strong predictors of future caseloads, this modeling approach places very little risk in areas with no historical cases in the dataset. The selected model performed as poorly as the null model during this time period, indicating that factors other than calculated OPVderived population immunity and local transmission contributed to this outbreak.
The outputs from the predictive models are well suited for use in prioritization by public health organizations. In Nigeria, model forecasts are used to inform subnational SIA planning and allocation of specialized personnel. Prioritization analysis played a critical role in the distribution of technical and administrative field staff—supported by both WHO and the Nigerian government—to the highest risk Local Government Areas across the north during the capacity surge in June 2012. WHO alone grew an initial staff of 744 to over 2,900, an increase of nearly 300%. In addition to supporting SIA planning, monitoring, and evaluation, these personnel have supported householdbased micro planning, intensified AFP surveillance activities, and helped strengthen routine immunization activities. Outreach from the federal government to state governors and district chairmen to increase local political buyin was also heavily increased around this time; prioritizations based on the predictive risk model results were also used to direct this surge in advocacy. These efforts can be partly credited for the absence of cases in 2013 from the northwest of Nigeria, the area of focus for political engagement in the latter half of 2012. F Since June 2013, model outputs have been combined with results and input from the National Primary Health Care Development Agency (NPHCDA), Centers for Disease Control and Prevention (CDC), and WHO to categorize a subset of districts as “high risk” and “highest risk”. This categorization has been used to direct a number of interventions across partner organizations in Nigeria. Management support teams, comprising high level personnel from NPHCDA and other GPEI partner agencies, are deployed to prioritized districts seven to ten days before a SIA campaign to address challenges limiting vaccination coverage. Prioritized districts are selected for advocacy visits targeted at local political, traditional, and religious leaders. Additionally, supplemental logistics funds are provided to prioritized districts to enhance the ability of teams to vaccinate hardtoreach and scattered settlements. The tracking of vaccinators using GPSenabled smart phones has been targeted to prioritized districts. Often, these districts receive the first implementation of an intervention planned to eventually deploy across northern Nigeria.
Outside of national program planning, district prioritization categorizations are closely monitored by state governors and district chairmen. Because of additional support, a prioritized district reporting a WPV case is held more accountable than a nonprioritized district.
The quantification of uncertainty in risk rank can also be a useful tool for policy makers. As changes in prioritization are an additional strain on resources (due to required reallocation of people and materials), the quantification of uncertainty in rank can enable objective decisions regarding changes to the prioritization of an area. If a district is newly ranked in the top 100 or 200 highest risk areas, policy makers may only want to prioritize if there is a high level of certainty that the district belongs in the highest risk group.
One aspect not considered in the model was the migration and transport structures connecting nonneighboring districts. With data to inform such structures, we may be better able to predict infection in naive southern districts. WPV1 case information was only available beginning in 2004 for this study. More historical data, with instances of transmission in these areas, could also have improved the forecasting accuracy of the model during these time periods. Seasonal variations and trends, such as those used in predictive models for meningitis in France and Mali [20, 21], could be incorporated, though the sixmonth time scale of model predictions required due to data sparsity may be too large for this technique to be useful.
Additional covariates or more representative data are needed to more fully understand WPV transmission in Nigeria. In the selected model, a substantial amount of residual spatial variation remains, as is demonstrated in Figure 4. This variation may be due to other factors such as poverty, malnutrition, sanitation, and level of health services, which influence WPV transmission potential and population vaccine efficacy [52, 53]. In addition, the inclusion of a number of operational factors could greatly improve the model. Indicators capturing district management performance, training quality, vaccinator selection, population accessibility (the presence of hardtoreach areas within districts, which may be seasonal), and noncompliance are currently missing from the model. Such indicators are an important part of regular program operations. Furthermore, these indicators could be more dynamic than current smoothed indicators, thus improving the responsiveness of model outputs to shortterm changes in performance.
Based on the magnitude of districtlevel random effects, we find that in large portions of northern Nigeria covariates underestimate the probability of the presence of WPV1 in a district, while in southern Nigeria, these covariates overestimate the probability of WPV1 presence. The latter could occur partly because we assume that vaccine efficacy is the same throughout the country, which results in relatively low population immunity estimates in southern districts (Figure 3). This observation is a direct result of the differential number of SIA campaigns executed in northern and southern Nigeria historically: over the last few years, more than twice as many campaigns have been carried out in northern states than in southern states. There is evidence and theory to suggest that vaccine efficacy is not uniform across Nigeria; it may actually be higher in southern Nigeria [33, 54], possibly due to a lower burden of nonpolio enteroviruses, which are known to interfere with OPV efficacy [55, 56]. OPVderived population immunity estimates based on geographically sensitive vaccine efficacies could resolve this anomaly and strengthen the measured association between calculated population immunity and WPV transmission.
There are several possible sources of statistical bias in our analysis. Though we treated NPAFP as a random sample of the population in an area, NPAFP may be underrepresentative of subpopulations in a district, especially higher risk subpopulations in areas with worse sanitation and less access to health services. This NPAFP bias is also likely to vary by location. Though the annual rate of NPAFP (8/100,000 under fifteen) is higher than the minimum WHO guideline (2/100,000 under fifteen), small sample sizes limit the temporal and spatial resolution of this analysis.
Another limitation of this analysis is the recall bias that arises from the oral history collected from the mothers of reported AFP cases. It is possible that they may over or underreport the number of OPV doses received by the child. It is also possible that the literacy level of the caregivers may influence the accuracy of reported doses. Our analysis does not make allowances for the impact of insecurity on the observed AFP detection rates in some states, such as Borno, Yobe and Kano; the effects may be more meaningful in data collected since December 2012.
There likely exist important heterogeneities in WPV transmission and associated factors below the district level: SIA activities are often planned and carried out at the subdistrict (ward) level in Nigeria. The mean population immunity may be poorly representative of atrisk populations, which will attenuate estimated relationships in the model. Certain populations, such as nomadic peoples [57], may be missed by routine AFP surveillance; in this case, there may be additional WPV cases not included in the analysis. Poorly performing wards may persist in districts with low average risk, as each ward has a focal person responsible for vaccinator selection and training. The population diversity within a district, which often include both rural and urban environments, suggests that the ward level is a more representative level for analysis, although we are severely constrained by a lack of historical data and sparse sample sizes at the subdistrict level.
Conclusion
The model developed in this study successfully predicts districts at risk for future wild poliovirus cases in Nigeria. Furthermore, the smoothing model handles missing data and smooths over unrealistic changes in important covariates. By quantifying uncertainty in risk ranking, we minimize the likelihood that decision makers respond to stochastic fluctuations in predicted risk over time. Predictive model outputs are well suited for use by public health organizations: in Nigeria, model forecasts are used to inform subnational SIA planning and to direct diverse interventions.
In order to more fully understand WPV transmission in Nigeria, additional covariates or more representative data are necessary. Heterogeneity in indicators or outcomes below the district level were not considered by this analysis due to data limitations, though there is evidence that this may be important. In addition, better metrics of connectedness between districts—informed by road networks or other transportation structures–could also have improved our model.
In addition to its usefulness for polio eradication through maximizing the impact of political and financial resources, this modeling approach could be extended to other vaccine preventable diseases (VPDs), such as measles. With adequate surveillance data, this type of predictive model provides a principled way of prioritizing administrative areas in future control or elimination efforts.
Endnotes
^{a} This is compatible with a relatively low life expectancy, especially in northern Nigeria. The 2008 Demographic and Health Survey (DHS) in Nigeria estimates an underfive proportion of 0.171 [58]. In this analysis, the exact fraction is not particularly crucial; as a single value is used for all districts, the relative populations of the districts are preserved.
^{b} For example, to make predictions about the second half of 2008 (June to November 2008), we smoothed indicators through May 2008, and estimated a model with data from July 2004 through May 2008. We then used population immunity and caseload information from the first half of 2008 (December 2007 through May 2008) to make predictions and rank districts according to model outputs.
^{c} Waning of individual immunity has been documented [59, 60]; however, this should not be relevant on a sixmonth time scale.
Abbreviations
 AFP:

acute flaccid paralysis
 AUC:

area under the curve
 CDC:

Centers for Disease Control and Prevention
 CAR:

conditional autoregressive
 CI:

credible interval
 DHS:

Demographic and Health Survey
 DIC:

deviance information criterion
 GPEI:

Global Polio Eradication Initiative
 MCMC:

markov chain monte carlo
 NPAFP:

nonpolio acute flaccid paralysis
 NPHCDA:

National Primary Heath Care Development Agency
 OPV:

oral polio vaccine
 ROC:

receiver operator characteristic
 SIA:

supplementary immunization activity
 VPD:

vaccine preventable disease
 WHO:

World Health Organization
 WPV:

wild poliovirus
 WPV1:

wild poliovirus type 1
 WPV3:

wild polio virus type 3.
Declarations
Acknowledgements
The authors would like to thank Bill and Melinda Gates for their active support of this work and their sponsorship through the Global Good Fund. We would also like to thank Kevin McCarthy for helpful comments on the manuscript. Productive discussions with colleagues at the Institute for Disease Modeling are likewise greatly appreciated. Finally, we would like to thank all those in Nigeria involved in AFP surveillance and laboratory testing.
Authors’ Affiliations
References
 Aylward B, Yamada T: The Polio endgame. N Engl J Med. 2011, 364: 22732275. 10.1056/NEJMp1104329.View ArticlePubMedGoogle Scholar
 World Health Organization: Polio eradication & endgame strategic plan 20132018. [Online]. Available from: http://www.polioeradication.org, Geneva.
 World Health Assembly: Poliomyelitis: Intensification of the Global Eradication Initiative. 2012, Geneva: World Health OrganizationGoogle Scholar
 National Primary Health Care Development Agency: 2012 Nigeria polio eradication emergency plan. 2012, [Online]. Available from: http://www.polioeradication.org, Abuja.Google Scholar
 Grassly NC: The final stages of the global eradication of poliomyelitis. Philos Trans R Soc London Ser B. 2013, 368: 2012014010.1098/rstb.2012.0140. doi:10.1098/rstb.2012.0140.View ArticleGoogle Scholar
 Lowther SA, Roesel S, O’Connor P, Landaverde M, Oblapenko G, Deshevoi S, Ajay G, Buff A, Safwat H, Salla M, Tangermann R, Khetsuriani N, Martin R, Wassilak S: World health organization regional assessments of the risks of poliovirus outbreaks. Risk Anal. 2013, 33: 664679. 10.1111/risa.12032. doi:10.1111/risa.12032.View ArticlePubMedGoogle Scholar
 O’Reilly KM, Chauvin C, Aylward RB, Maher C, Okiror S, Wolff C, Nshmirimana D, Donnelly CA, Grassly NC: A statistical model of the international spread of wild poliovirus in Africa used to predict and prevent outbreaks. PLoS Med. 2011, 8, doi:10.1371/journal.pmed.1001109.Google Scholar
 Best N, Richardson S, Thomson A: A comparison of Bayesian spatial models for disease mapping. Stat Methods Med Res. 2005, 14: 3559. 10.1191/0962280205sm388oa. doi:10.1191/0962280205sm388oa.View ArticlePubMedGoogle Scholar
 Brooker S, Clements ACA, Hotez PJ, Hay SI, Tatem AJ, Bundy DAP, Snow RW: The codistribution of Plasmodium falciparum and hookworm among African schoolchildren. Malar J. 2006, 5: 9910.1186/14752875599. doi:10.1186/14752875599.View ArticlePubMedPubMed CentralGoogle Scholar
 Hay SI, Guerra CA, Gething PW, Patil AP, Tatem AJ, Noor AM, Kabaria CW, Manh BH, Elyazar IRF, Brooker S, Smith DL, Moyeed RA, Snow RW: A world malaria map: Plasmodium falciparum endemicity in 2007. PLoS Med. 2009, 6: 1000048doi:10.1371/journal.pmed.1000048.View ArticleGoogle Scholar
 Clements ACA, Garba A, Sacko M, Touré S, Dembelé R, Landouré A, BosqueOliva E, Gabrielli AF, Fenwick A: Mapping the probability of schistosomiasis and associated uncertainty, West Africa. Emerg Infect Dis. 2008, 14: 16291632. 10.3201/eid1410.080366. doi:10.3201/eid1410.080366.View ArticlePubMedPubMed CentralGoogle Scholar
 Clements ACA, Lwambo NJS, Blair L, Nyandindi U, Kaatano G, Kinung’hi S, Webster JP, Fenwick A, Brooker S: Bayesian spatial analysis and disease mapping: tools to enhance planning and implementation of a schistosomiasis control programme in Tanzania. Trop Med Int Heal. 2006, 11: 490503. 10.1111/j.13653156.2006.01594.x. doi:10.1111/j.13653156.2006.01594.x.View ArticleGoogle Scholar
 Clements ACA, Brooker S, Nyandindi U, Fenwick A, Blair L: Bayesian spatial analysis of a national urinary schistosomiasis questionnaire to assist geographic targeting of schistosomiasis control in Tanzania, East Africa. Int J Parasitol. 2008, 38: 401415. 10.1016/j.ijpara.2007.08.001. doi:10.1016/j.ijpara.2007.08.001.View ArticlePubMedPubMed CentralGoogle Scholar
 Vounatsou P, Raso G, Tanner M, N’goran EK, Utzinger J: Bayesian geostatistical modelling for mapping schistosomiasis transmission. Parasitology. 2009, 136: 16951705. 10.1017/S003118200900599X. doi:10.1017/S003118200900599X.View ArticlePubMedGoogle Scholar
 Wang XH, Zhou XN, Vounatsou P, Chen Z, Utzinger J, Yang K, Steinmann P, Wu XH: Bayesian spatiotemporal modeling of Schistosoma japonicum prevalence data in the absence of a diagnostic ‘gold’ standard. PLoS Negl Trop Dis. 2008, 2: 25010.1371/journal.pntd.0000250. doi:10.1371/journal.pntd.0000250.View ArticleGoogle Scholar
 Pullan RL, Gething PW, Smith JL, Mwandawiro CS, Sturrock HJW, Gitonga CW, Hay SI, Brooker S: Spatial modelling of soiltransmitted helminth infections in Kenya: a disease control planning tool. PLoS Negl Trop Dis. 2011, 5: 95810.1371/journal.pntd.0000958. doi:10.1371/journal.pntd.0000958.View ArticleGoogle Scholar
 Raso G, Vounatsou P, Gosoniu L, Tanner M, N’Goran EK, Utzinger J: Risk factors and spatial patterns of hookworm infection among schoolchildren in a rural area of western Côte d’Ivoire. Int J Parasitol. 2006, 36: 201210. 10.1016/j.ijpara.2005.09.003. doi:10.1016/j.ijpara.2005.09.003.View ArticlePubMedGoogle Scholar
 Raso G, Vounatsou P, Singer BH, N’Goran EK, Tanner M, Utzinger J: An integrated approach for risk profiling and spatial prediction of Schistosoma mansonihookworm coinfection. Proc Natl Acad Sci. 2006, 103: 69346939. 10.1073/pnas.0601559103. doi:10.1073/pnas.0601559103.View ArticlePubMedPubMed CentralGoogle Scholar
 Diggle P, Thompson M, Christensen O, Rowlingson B, Obsomer V, Gardon J, Wanji S, Takougang I, Enyoun P, Kamgno J, Remme J, Boussinesq M, Molyneux D: Spatial modelling and the prediction of Loa loa risk: descision making under uncertainty. Ann Trop Med Parasitol. 2007, 101: 499509. 10.1179/136485913X13789813917463. doi:10.1179/136485907X229121.View ArticlePubMedGoogle Scholar
 Paul M, Held L: Predictive assessment of a nonlinear random effects model for multivariate time series of infectious disease counts. Stat Med. 2011, 30: 11181136. doi:10.1002/sim.4177.PubMedGoogle Scholar
 Stanton MC, Agier L, Taylor BM, Diggle PJ: Towards realtime spatiotemporal prediction of district level menigitis incidence in subSaharan Africa. J R Stat Soc A. 2013, [Online]. doi:10.1111/rssa.12033.Google Scholar
 Musal M, Aktekin T: Bayesian spatial modeling of HIV mortality via zeroinflated poisson models. Stat Med. 2013, 32: 267281. 10.1002/sim.5457. doi:10.1002/sim.5457.View ArticlePubMedGoogle Scholar
 Agarwal DK, Gelfand AE, CitronPousty S: Zeroinflated models with application to spatial count data. Environ Ecol Stat. 2002, 9: 341355. 10.1023/A:1020910605990.View ArticleGoogle Scholar
 Min Y, Agresti A: Random effect models for repeated measures of zeroinflated count data. Stat Model. 2005, 5: 119. 10.1191/1471082X05st084oa. doi:10.1191/1471082X05st084oa.View ArticleGoogle Scholar
 Neelon BH, O’Malley AJ, Normand SLT: A Bayesian model for repeated measures zeroinflated count data with application to outpatient psychiatric service use. Stat Modelling. 2010, 10: 421439. 10.1177/1471082X0901000404. doi:10.1177/1471082X0901000404.View ArticlePubMedPubMed CentralGoogle Scholar
 Neelon B, Ghosh P, Loebs PF: A spatial Poisson hurdle model for exploring geographic variation in emergency department visits. J R Stat Soc Ser A. 2013, 176: 389413. 10.1111/j.1467985X.2012.01039.x. doi:10.1111/j.1467985X.2012.01039.x.View ArticleGoogle Scholar
 Marx A, Glass JD, Sutter RW: Differential diagnosis of acute flaccid paralysis and its role in poliomyelitis surveillance. Epidemiol Rev. 2000, 22: 298316. 10.1093/oxfordjournals.epirev.a018041.View ArticlePubMedGoogle Scholar
 Jenkins HE, Aylward RB, Gasasira A, Donnelly CA, Mwanza M, Corander J, Garnier S, Chauvin C, Abanida E, Pate MA, Adu F, Baba M, Grassly NC: Implications of a circulating vaccinederived poliovirus in Nigeria. N Engl J Med. 2010, 362: 23602369. 10.1056/NEJMoa0910074. doi:10.1056/NEJMoa0910074.View ArticlePubMedGoogle Scholar
 Jenkins HE, Aylward RB, Gasasira A, Donnelly CA, Abanida EA, KoleoshoAdelekan T, Grassly NC: Effectiveness of immunization against paralytic poliomyelitis in Nigeria. N Engl J Med. 2008, 359: 16661674. 10.1056/NEJMoa0803259. doi:10.1056/NEJMoa0803259.View ArticlePubMedGoogle Scholar
 O’Reilly KM, Durry E, ul Islam O, Quddus A, Abid N, Mir TP, Tangermann RH, Aylward RB, Grassly NC: The effect of mass immunisation campaigns and new oral poliovirus vaccines on the incidence of poliomyelitis in Pakistan and Afghanistan, 2001–11: a retrospective analysis. Lancet. 2012, 380: 491498. 10.1016/S01406736(12)606485. doi:10.1016/S01406736(12)606485.View ArticlePubMedPubMed CentralGoogle Scholar
 National Population Commission of Nigeria: Report of Nigeria’s national population commission on the 2006 census. Popul Dev Rev. 2007, 33: 206210.Google Scholar
 GADM: GADM Database of Global Administrative Areas, Version 2.0. [Online]. Available from: http://www.gadm.org [Accessed 10th Apr 2013].
 Mangal TD, Aylward RB, Mwanza M, Gasasira A, Abanida E, Pate MA, Grassly NC: Key issues in the persistence of poliomyelitis in Nigeria: a casecontrol study. Lancet Glob Health. 2014, 2: 9097. 10.1016/S2214109X(13)701682.View ArticleGoogle Scholar
 Behrend M, Hu H, Nigmatulina KR, Eckhoff P: A quantitative suvery of the literature on poliovirus infection and immunity. Int J Infect Dis. 2014, 18: 413.View ArticlePubMedGoogle Scholar
 Mullahy J: Specification and testing of some modified count data models. J Econom. 1986, 33: 341365. 10.1016/03044076(86)900023. doi:10.1016/03044076(86)900023.View ArticleGoogle Scholar
 Pohllmeier W, Ulrich V: An econometric model of the twopart decisionmaking process in the demand for health care. J Hum Resour. 1995, 30: 339361. 10.2307/146123.View ArticleGoogle Scholar
 Besag J, York J, Mollié A: Bayesian image restoration, with two applications in spatial statistics. Ann Inst Statist Math. 1991, 43: 159. 10.1007/BF00116466.View ArticleGoogle Scholar
 Besag J, Kooperberg C: On conditional and intrinsic autoregressions. Biometrika. 1995, 82: 733746. doi:10.1093/biomet/82.4.733.Google Scholar
 Lunn DJ, Thomas A, Best N, Spiegelhalter D: WinBUGS – A Bayesian modelling framework: concepts, structure, and extensibility. Stat Comput. 2000, 10: 325337. 10.1023/A:1008929526011.View ArticleGoogle Scholar
 Sturtz S, Ligges U, Gelman A: R2WinBUGS: a package for running WinBUGS from R. J Stat Softw. 2005, 12: 116.View ArticleGoogle Scholar
 Plummer M, Best N, Cowles K, Vines K: CODA: convergence diagnosis and output analysis for MCMC. R News. 2006, 6: 711.Google Scholar
 R Core Team: R: A Language and Evironment for Statistical Computing. 2013, Vienna: R Foundation for Statistical ComputingGoogle Scholar
 Gelman A, Rubin DB: lnference from iterative simulation using multiple sequences. Stat Sci. 1992, 7: 457472. 10.1214/ss/1177011136.View ArticleGoogle Scholar
 Brooks SP, Roberts GO: Assessing convergence of Markov chain Monte Carlo algorithms. Stat Comput. 1998, 8: 319335. 10.1023/A:1008820505350.View ArticleGoogle Scholar
 Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A: Bayesian measures of model complexity and fit. J R Stat Soc Ser B. 2002, 64: 583639. 10.1111/14679868.00353. doi:10.1111/14679868.00353.View ArticleGoogle Scholar
 Fielding AH, Bell JF: A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ Conserv. 1997, 24: 3849. 10.1017/S0376892997000088. doi:10.1017/S0376892997000088.View ArticleGoogle Scholar
 Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, Müller M: pROC: an opensource package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics. 2011, 12: 7710.1186/147121051277.View ArticlePubMedPubMed CentralGoogle Scholar
 Brooker S, Hay SI, Bundy DAP: Tools from ecology: useful for evaluating infection risk models?. Trends Parasitol. 2002, 18: 7074. 10.1016/S14714922(01)022231.View ArticlePubMedPubMed CentralGoogle Scholar
 Greiner M, Pfeiffer D, Smith RD: Principles and practical application of the receiveroperating characteristic analysis for diagnostic tests. Prev Vet Med. 2000, 45: 2341. 10.1016/S01675877(00)00115X.View ArticlePubMedGoogle Scholar
 Thompson KM, Pallansch MA, Tebbens RJD, Wassilak SG, Cochi SL: Modeling population immunity to support efforts to end the transmission of live polioviruses. Risk Anal. 2013, 33: 647663. 10.1111/j.15396924.2012.01891.x. doi:10.1111/j.15396924.2012.01891.x.View ArticlePubMedGoogle Scholar
 Hu H, Nigmatulina K, Eckhoff P: The scaling of contact rates with population densitity for the infectious disease models. Math Biosci. 2013, 244: 125134. 10.1016/j.mbs.2013.04.013.View ArticlePubMedGoogle Scholar
 Grassly NC, Fraser C, Wenger J, Deshpande JM, Sutter RW, Heymann DL, Aylward RB: New strategies for the elimination of polio from India. Science. 2006, 314: 11501153. 10.1126/science.1130388. doi:10.1126/science.1130388.View ArticlePubMedGoogle Scholar
 Webster P: A poliofree world?. Lancet. 2005, 366: 359360. 10.1016/S01406736(05)670080.View ArticlePubMedGoogle Scholar
 Orenstein WA, Bernier RH, Hinman AR: Assessing vaccine efficacy in the field. Epidemiol Rev. 1988, 10: 212241.PubMedGoogle Scholar
 Triki H, Ould Mohamed Abdallah MV, Ben Aissa R, Bouratbine A, Ben Ali Kacem M, Bouraoui S, Koubaa C, Zouari S, Mohsni E, Crainic R, Dellagi K: Influence of host related factors on the antibody response to trivalent oral polio vaccine in Tunisian infants. Vaccine. 1997, 15: 11231129. 10.1016/S0264410X(97)000017.View ArticlePubMedGoogle Scholar
 Patriarca PA, Wright PF, John TJ: Factors affecting the immunogenicity of oral poliovirus vaccine in developing countries: review. Rev Infect Dis. 1991, 13: 926939. 10.1093/clinids/13.5.926.View ArticlePubMedGoogle Scholar
 Callaway E: Polio’s moving target. Nature. 2013, 496: 290292. 10.1038/496290a.View ArticlePubMedGoogle Scholar
 National Population Commission [Nigeria]: Nigeria Demographic and Health Survey 2008. 2009, Abuja, Nigeria: Technical report, National Population Commission and ICF MacroGoogle Scholar
 Grassly NC, Jafari H, Bahl S, Sethi R, Deshpande JM, Wolff C, Sutter RW, Aylward RB: Waning intestinal immunity after vaccination with oral poliovirus vaccines in India. J Infect Dis. 2012, 205: 15541561. 10.1093/infdis/jis241. doi:10.1093/infdis/jis241.View ArticlePubMedGoogle Scholar
 Mayer BT, Eisenberg JNS, Henry CJ, Gomes MGM, Ionides EL, Koopman JS: Successes and shortcomings of polio eradication: a transmission modeling analysis. Am J Epidemiol. 2013, 177: 12361245. 10.1093/aje/kws378. doi:10.1093/aje/kws378.View ArticlePubMedPubMed CentralGoogle Scholar
 The prepublication history for this paper can be accessed here:http://www.biomedcentral.com/17417015/12/92/prepub
Prepublication history
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.