Introduction {#Sec1}
============

The invention of and advancements in the childhood vaccination has revolutionized human healthcare and facilitated to averting morbidity and mortality from many vaccine-preventable infectious diseases^[@CR1]^. Immunization not only has helped preventing life-threatening diseases but also benefitted mankind by increasing their life expectancy (life span has increased by 15 to 25 years since the commencement of vaccines) and in improving quality of life^[@CR2],[@CR3]^. However, these successes cannot mask the failures in shape of millions of deaths of the unvaccinated children annually, across the globe^[@CR4]^. Estimated 1.5 million children died of vaccine-preventable diseases in each of 2013 and 2017 respectively which construes that the world could have avoided an additional huge number of deaths a year had the universal immunization coverage achieved^[@CR5],[@CR6]^.

Despite the fact that overall vaccination coverage in 2018 remained around 86%^[@CR6]^, there remained a wide disparity between different World Health Organization (WHO) regions, e.g., the Americas, Europe and Western Pacific maintained over 90% immunization coverage of DTP3 containing vaccine whilst these rates were low for most of the countries in Africa and some in Asia^[@CR7]^. Low levels of immunization coverage in the low-and middle-income countries pose a serious challenge to achieving universal vaccination coverage goals^[@CR8],[@CR9]^. For any child to develop adequate immunity against the particular disease, adherence to complete vaccination schedule is of critical importance, as, incomplete vaccination leads to partial immunity and the disease risk persists^[@CR10],[@CR11]^. There have been fatal outbreaks of vaccine-preventable infectious diseases in the developing countries which reflect the presence in large numbers of non-and under-vaccinated children^[@CR12],[@CR13]^. As a matter of fact, the bottommost ten countries with lowest vaccination coverage belonged to either low or lower-middle income group^[@CR7]^. Approximately 19.4 million eligible children did not receive DTP3 dosage worldwide---around 40% (8 million) of these children lived in war-afflicted areas^[@CR6]^.

In Pakistan, government patronage for vaccination coverage through the Expanded Programme on Immunization (EPI) has been there for decades yet the low coverage areas have been stubbornly in existence across the country^[@CR14]^. This has increased vulnerability to vaccine preventable diseases and not surprisingly Pakistan has been amongst countries with highest mortality and morbidity rates in children under 5 years of age^[@CR15],[@CR16]^. Under-5 deaths make nearly half of all the deaths in Pakistan, which is around 8 to 10% of all deaths in the developed countries^[@CR17]^. Pakistan has the lowest vaccination coverage statistics amongst the sub-continent countries (Bangladesh, Sri Lanka and India)^[@CR18]^. These neighbouring countries had also overcome their respective indigenous polio transmission since 2011 through their massive scale high quality vaccination campaigns^[@CR18]--[@CR20]^, unlike Pakistan which is still struggling with Polio endemicity (along with Afghanistan and Nigeria)^[@CR21]^.

Statistics from different sources in Pakistan have reported varying percentages of vaccination coverage, ranging from 40 to 80% with wide variations between provinces as well as between male and female children^[@CR22]^. The higher percentages for vaccination coverage are generally reported by the government while some researchers and non-governmental agencies have reported low vaccination coverages across Pakistan over the time^[@CR23]--[@CR25]^. This variation may be attributed to the methodology adopted in collection of data (recall/record base data, estimation of the numerators/denominators), over/miss-reporting of the collected data, training/skill of the field teams and some other socio-geographical factors (parents' knowledge, religious beliefs and accessibility to health facility).

EPI is virtually the exclusive provider of immunization services in Pakistan, conducting around 97% of the total immunization activities in the country. EPI in Pakistan vaccinates more than 5 million children under 1 year of age annually, in order to provide them protection against eight vaccine preventable infectious diseases through its routine immunization services^[@CR26]^. Vaccination coverage with third dose of diphtheria, pertussis and tetanus (DTP3) is a commonly used indicator for the routine immunization services' performance of various countries; whereas, Bacillus Calmette Guerin (BCG), third dose of oral polio vaccine (OPV3), pentavalent3 (containing DTP3) and Measles1 are considered to be yardstick of completed vaccination schedule^[@CR27]^. Therefore, this study did not discuss rest of the vaccines which are included in EPI Pakistan routine immunization schedule but are further than the criteria of labelling a child as fully immunized (e.g. for *Haemophilus Influenzae*, Meningitis, Rota virus, Pneumonia and 2nd dose of Measles).

Epidemiologists have been interested in the geographical clustering of diseases^[@CR28]^ and other public health problems such as communities with non/low vaccination^[@CR29]^. Low vaccination coverage areas tend to exist in clusters and spatial scan statistics are an important tool to detect them for their subsequent targeting by the concerned authorities^[@CR29]^. There are only few studies available on vaccination coverage discussing country-wide data and even scarce studies using modern statistical methods (e.g., spatial analysis) in Pakistan. The available data on vaccination coverage in Pakistan is limited and on many occasions faulty and over-estimated^[@CR30]^. The available studies on vaccination coverage either discuss secondary data from surveys across Pakistan or are about one or few districts and predominantly about single antigen^[@CR25],[@CR31]--[@CR33]^. Our study has tried to contribute in these lacking areas of spatial analysis on vaccination coverage across Pakistan.

Results {#Sec2}
=======

Total number of children vaccinated in each year (2012 to 2016) and for each vaccination type are shown in Table [1](#Tab1){ref-type="table"}.Table 1Number of vaccinated children for each EPI vaccination type in Pakistan, 2012--2016.Vaccination type20122013201420152016BCG6,002,9836,093,2166,087,5305,825,4056,217,308OPV35,307,8555,365,1505,340,1325,243,0345,370,649Pentavalent35,255,2365,378,4425,338,4995,213,1645,363,685Measles15,428,9665,602,6425,354,0795,185,4885,444,831

Descriptive analysis {#Sec3}
--------------------

The descriptive statistics were indicated in Table [2](#Tab2){ref-type="table"}. In this table, high standard deviation and significant difference between minimum and maximum values indicate widely spread out vaccination coverage status among Pakistan districts. The coverage rates with highest mean values for all the vaccine types were found in 2012, while the lowest mean vaccination coverage rates were seen in 2015, e.g., mean BCG vaccination coverage rate was 70.89 in 2012, which was highest in all the study years, whereas, in 2015 it was 61.79 and lowest of all the study years.Table 2Summary of descriptive statistics of EPI vaccination coverage rates (per 100) by type at the district level in Pakistan, 2012--2016.YearVaccine typeMeanStd. deviationMinimumQuantilesMaximum2550752012BCG70.8915.9424.5362.5273.7780.6498.02OPV 362.3814.7317.5656.6167.1171.2193.62Pentavalent 362.0914.5817.5656.5865.8170.5193.62Measles 162.1716.1615.8054.1666.5574.7587.932013BCG69.3415.188.0263.3573.0778.9497.69OPV 360.2915.703.9855.3165.5370.5296.23Pentavalent 360.5515.295.3455.3165.5670.6396.23Measles 161.4617.296.1052.4067.1672.9495.342014BCG66.8315.8110.9859.9371.4676.7097.87OPV 357.4616.251.2751.3162.7870.0892.70Pentavalent 357.4716.183.8251.3462.7870.1092.70Measles 157.1816.825.6847.4361.7670.4680.512015BCG61.7916.103.5554.4766.1273.2395.14OPV 355.8316.111.7748.4760.9767.4885.79Pentavalent 355.6315.941.7748.4760.9766.9985.79Measles 154.4817.572.5142.5259.7868.3384.482016BCG67.1414.3718.1957.4871.0276.7097.47OPV 357.6114.565.8349.3261.4268.1187.60Pentavalent 357.5714.544.8249.3261.4267.7787.60Measles 157.5516.337.5047.2861.7269.9292.39

Spatiotemporal and spatial clustering analysis {#Sec4}
----------------------------------------------

The sequence plot for vaccination rates per 100 children for each type of vaccination were shown in the Fig. [1](#Fig1){ref-type="fig"} for national level as well as for all the provinces in Pakistan. The figure indicated that the rates (percentage) for BCG vaccination coverage were highest amongst all types of vaccination and there was a decline in 2015 in the rates of all types of vaccination coverage at the district level in Pakistan.Figure 1Sequence plot for EPI vaccination rates by type at national and provincial level in Pakistan, 2012--2016.

Figures [2](#Fig2){ref-type="fig"} exhibited vaccination coverage rates per 100 children for all EPI vaccination types in Pakistan at the district level. Vaccination coverage rates for BCG were higher than the rest of the vaccination types and mostly lied in high to very high categories. The rates for OPV3 and pentavalent3 were similar to one another, as, according to the EPI vaccination schedule, both these types are administered at the same time. The vaccination coverage rates for BCG, OPV3 and pentavalent3 were found to be relatively higher in Azad Kashmir (north), Punjab (east) and most districts of Sindh (southeast), as compared to the rest of the country. The coverage rates for measles1 were found low to very low in most districts of Balochistan, Khyber Pakhtunkhwa (KPK) and Federally Administered Tribal Areas (FATA), especially in 2015.Figure 2EPI vaccination coverage rates by type per 100 children at the district level in Pakistan, 2012--2016.

Table [3](#Tab3){ref-type="table"} described the primary clusters for EPI low vaccination coverage by type in each of the study year at the district level in Pakistan. The yearly information on primary clusters including; radii of clusters, target population, vaccination rates per 100,000 children and the *P* values showing statistical significance were highlights of this table.Table 3Summary of primary high-risk clusters' information of Expanded Program of Immunization (EPI) annual vaccination coverage by type at the district level in Pakistan, 2012--2016.Type of vaccineYearPrimary cluster radius (km)Primary cluster target populationNumber of vaccinated children in primary clusterVaccination rate/100,000RR*P* valueBCG2012123.33229,703126,33754,886.30.72 \< 0.00012013023,2421,8658,029.60.11 \< 0.00012014123.33241,469115,95348,051.70.65 \< 0.0001201576.7770,46413,30018,887.40.27 \< 0.0001201699.71112,51243,86338,904.40.54 \< 0.0001OPV3201287.22165,17071,09942,956.90.64 \< 0.00012013023,2429253,982.50.6 \< 0.0001201487.22173,97755,60231,980.60.49 \< 0.0001201591.3873,71010,72414,558.60.23 \< 0.0001201699.71112,51224,49121,722.40.35 \< 0.0001Pentavalent 3201287.22165,17071,57743,245.70.65 \< 0.00012013023,2421,2405,338.70.8 \< 0.0001201487.22173,97756,20832,329.20.5 \< 0.0001201591.3873,71010,72614,561.30.24 \< 0.0001201699.71112,51224,01121,296.60.34 \< 0.0001Measles 12012199.41179,75482,84745,993.70.67 \< 0.00012013153.6596,70428,27929,262.30.42 \< 0.00012014180.86865,089420,78048,672.40.75 \< 0.00012015195.57261,95194,19235,981.80.58 \< 0.00012016191.941,454,135731,13150,175.40.77 \< 0.0001

The size of primary clusters for measles1 low coverage were much larger than those for the other vaccination types included in the study, e.g., in 2014 and 2016 measles 1 clusters comprised of 22 and 29 districts, respectively; while for OPV3 and pentavalent3 vaccination types, there were seven districts per cluster in each of the corresponding study year. The measles1 clusters were also inconsistent and comprised of districts from different provinces for different years in the study.

Figure [3](#Fig3){ref-type="fig"} represented the primary, secondary and no cluster areas of low vaccination coverage clusters for all vaccination types in each year of study at the district level.Figure 3EPI low vaccination coverage by type clusters at district level in Pakistan, 2012--2016.

The districts detected in the primary clusters in each year of the study and for each type of vaccinations were as follows.

Table [4](#Tab4){ref-type="table"} depicts province wise details of the most significant clusters at the district level in Pakistan during the study period. Districts in low vaccination coverage primary clusters for BCG, in 2012, (total 10 districts, 5 each in FATA and KPK); in 2013, (1 district in FATA); in 2014, (total 10 districts, 5 each in FATA and KPK); in 2015, (total 5 districts, 3 in FATA, 1 each in KPK and Balochistan); in 2016, (total 7 districts, 4 in FATA, 2 in KPK and 1 in Balochistan).Table 4Province-wise distribution of primary high-risk clusters of Expanded Program of Immunization (EPI) annual vaccination coverage by type at the district level in Pakistan, 2012--2016.Type of vaccineYearTotal number of districts in primary clustersNumber of districts in FATANumber of districts in KPKNumber of districts in BalochistanNumber of districts in SindhNumber of districts in PunjabNumber of districts in GB/federal capitalNumber of districts in AJKBCG20121055----------201311------------20141055----------20155311--------20167421--------OPV32012743----------201311------------2014743----------20156411--------20167421--------Pentavalent 32012743----------201311------------2014743----------20156411--------20167421--------Measles 1201213--13--------201392--7--------2014221110----1----201518----18--------201629720----11--

Districts in low vaccination coverage primary clusters for OPV3 and pentavalent3 were alike for each of the corresponding years, i.e., in 2012, (total 7 districts, 4 in FATA and 3 in KPK); in 2013, (1 in FATA); in 2014, (total 7 districts, 4 in FATA and 3 in KPK); in 2015, (total 6 districts, 4 in FATA and 1 each in KPK and Balochistan); in 2016, (total 7 districts, 4 in FATA, 2 in KPK and 1 in Balochistan).

Districts in low vaccination coverage primary clusters for Measles 1: in 2012; (all 13 districts in Balochistan); in 2013 (total 9 districts, 7 in Balochistan and 2 in FATA) in 2014; (total 22 districts, 11 in FATA, 10 in KPK and 1 in Punjab); in 2015; (all 18 districts in Balochistan); in 2016; (total 29 districts, 20 in KPK, 7 in FATA and 1 each in Punjab and Federal capital).

Discussion {#Sec5}
==========

This study explored EPI immunization coverage of the children vaccinated against vaccine preventable diseases, including Childhood Tuberculosis (BCG), Poliomyelitis (OPV3), Diphtheria, Pertussis, Neonatal Tetanus, *Haemophilus Influenzae* type B, Hepatitis B (pentavalent3 which contains DTP3) and Measles (Measles 1) over a period of 5 years (2012 to 2016) at the district level in Pakistan. During the study period, 30,226,442 children were vaccinated with BCG, 26,626,820 children were vaccinated with OPV3, 26,548,979 were vaccinated with pentavalent3 and 27,016,006 were vaccinated with measles1 vaccine throughout the country, respectively.

The spatial cluster analysis showed that most of the low vaccination districts for BCG, OPV3 and pentavalent3 were from FATA and KPK; while for measles1 vaccination, most of the detected districts belonged to Balochistan Province. There was a consistent pattern of distribution of these clusters, especially for OPV3 and pentavalent3 and in some years of study for BCG as well. However, measles1 clusters were inconsistently distributed in different years of study and over larger areas, especially those found in Balochistan. Only two districts from Punjab were detected for measles1 primary clusters one each in 2014 and 2016 respectively. However, none of the districts from Punjab, Sindh, Gilgit Baltistan and Azad Kashmir were seen in any of the primary spatial clusters of low vaccination coverage during the study period.

These findings were consistent with country wide surveys in Pakistan. A policy brief by "Research and Development Solutions" in 2012 stated that routine coverage rates for fully immunized children were significantly low in Pakistan, especially in FATA, KPK, Balochistan and Sindh, moreover the study also agreed that one in every five children remains unimmunized at the national level while in the rural population this rate is even higher with two out of every three children being non-immunized^[@CR17]^.

Herd immunity is defined as the resistance a population develops against the spread of an infectious disease resulted from adequately high proportion of individuals being immune to that disease---especially through vaccination. There are different vaccination levels required to achieve herd immunity for various diseases, e.g. the herd immunity thresholds for Vaccine-Preventable Diseases are; Tuberculosis (BCG) not known, Polio (OPV) 50--95%, Diphtheria 83--85%, Tetanus N/A, Pertussis 80--94% (DTP) and Measles 92--94%^[@CR34]^. Our study found that, at the national level, average vaccination rates (all types of vaccination, BCG, OPV3, pentavalent3, measles1) per 100 children were 60.60, 67.12, 58.53, 58.47 and 58.29 respectively. Whereas, according to previous surveys, average vaccination rate for each type of vaccination varied significantly in assessments done by different bodies. For example, the latest Multiple Indicator Cluster Surveys (MICS) of Punjab (63.8% male, 64.4% female)^[@CR35]^, Sindh (43.1% male, 43.4% female)^[@CR36]^, KPK (55.5% male, 55.5% female)^[@CR37]^ and Balochistan (3.5% male, 4.8% female)^[@CR38]^, and Pakistan Social and Living Standards Measurement Survey (PSLM 2014--2015) Punjab (70% male, 70% female), Sindh (45% male, 45% female), KPK (60% male, 56% female) and Balochistan (30% male, 25% female) clearly depict marked variation in the vaccination coverage statistics across different provinces in Pakistan. Punjab is the most populous province harbouring more than 50% of total population in Pakistan, followed by Sindh which is the second largest. Despite that our study showed better vaccination coverage statistics for both these provinces as compared to FATA, Balochistan and KPK, yet the number of non/under-vaccinated children are still very high owing to the size and density of their population; suggesting that the decision makers should not lose focus on them. Also some antigens showed similar statistics as in our study, e.g. average vaccination rate per 100 children for BCG was 88, 85 and 89, according to PSLM 10--11, 12--13 and 14--15 respectively^[@CR39]--[@CR41]^; and 85 and 88 according to PDHS 12--13 and 17--18 respectively^[@CR42],[@CR43]^. BCG coverage was the highest among all the vaccination types and measles1 coverage was the lowest; similar trends were observed in some of the preceding assessments. However, there were varying trends in some of the districts and provinces in different years of this study.

Another study showed marked variation from our results with only 31% fully vaccinated children in Pakistan which is alarming and at the same time necessitate ponder over^[@CR24]^. This huge discrepancy in the reported vaccination coverage status through different sources creates doubt on the health of data. Reporting bias in administrative data can be one of these factors as WHO and United Nations Children Fund (UNICEF) have already raised questions on the quality of administrative data reported by various member countries^[@CR44]^. The numerator errors/bias (errors/bias in reporting number of doses administered by health staff either inadvertently or purposely to magnify their performance) has been accounted for in this study to a fair extent but there are other issues like wrong denominator estimations due to varying population growth rates in the districts and numerator/denominator mismatches due to migrations in and out of the districts and provinces which might have affected the results to some extent. Majority of the population in Pakistan live in the rural areas with urbanization trends increasing with the passage of time. Migration within the country is common, e.g. one third of these migrants moved from rural to the urban settings^[@CR45]^. According to Pakistan Labour Force Survey 2014--2015, rural to urban migration accounted for 22.5% of the total flow of migrants within and between districts and out of these migrations, 43.8% in Balochistan, 25.1% in Punjab, 21.1% in KPK and 13.9% in Sindh^[@CR46]^. These migration trends may explain the low vaccination coverage due to more missed children in Balochistan and higher vaccination trends in Sindh^[@CR23],[@CR47]^. Also, the heterogeneities in coverage of different districts and provinces is a key factor in non-achievement of desired periodical targets and eventually of herd immunity at the subnational and national levels. Spatial heterogeneity in vaccination coverage impacts the health status of the residents as it can delay disease elimination, not only in the developing but also in the developed countries where nation-wide vaccination coverage rates are usually higher^[@CR48],[@CR49]^.

All the national surveys (e.g., PSLM, Pakistan Demographic and Health Survey (PDHS) conducted in the past used similar measures and methods to estimate immunization coverage rates and yet produced varying results^[@CR40]--[@CR43]^. This was because of relying on recall of mother's memory (as done in PDHS) or any one of the household member's memory (as done in PSLM) for the estimation. On the other hand, assessment on the basis of child's vaccination card record has resulted in underestimation because of poor record keeping practices by the household members^[@CR50]^. In fact, both these methods (record based and recall based) cannot be relied upon independently, as they provided inaccurate estimation of the immunization rates^[@CR7]^. This was confirmed by a study done in one district of Sindh where the accuracy of BCG vaccination coverage was assessed using mothers' recall, vaccination card record and the record maintained by health staff (Vaccinator) at the healthcare facility. Mothers' recall was moderately accurate, while vaccination card record and health staff record were very less accurate^[@CR51]^. The vaccination coverage record maintained at respective healthcare facilities in Pakistan is often deficient as it usually reflects aggregated number of the vaccinated children only and lacks demographic information (e.g. name, age, sex and home address) of individual children. The vaccination process, however is carried out through a well versed healthcare system administered by trained staff with monitoring and evaluation being a part and parcel of it having varied efficiencies and capacities across the districts in Pakistan^[@CR52]^. Hence in the current scenario it is the best available data to perform country wide analysis and do the forecasting for future vaccination requirements.

Besides Balochistan, very low immunization coverage statistics were also observed in FATA and KPK in our study, which is not a surprise. A few districts from Balochistan and KPK while most districts from FATA have been conflict struck with compromised law and order situations which resulted in humanitarian crises; millions of temporarily displaced persons, lack of basic equipment and professional staff in healthcare facilities, and poor access to health services. And this all added up to already lower immunization coverage status in these provinces^[@CR21],[@CR53],[@CR54]^. In 2016 only, an estimated 20% of households in FATA did not vaccinate their children for any routine EPI vaccine^[@CR21]^.

Formal registrations of the new-borns in Pakistan has not been widely practiced, but the trends have changed over the years, especially after the establishment of National Database Authority Pakistan (NADRA) in 2000. The registrations of the children increased between the ages of 2 to 4 as the birth certificates during the school admissions are mandatory, but there is much room of improvement in the current practices of child registration at birth^[@CR42]^, which can contribute in bettering the vaccination coverage in the country.

The performance of Expanded Program on Immunization in Pakistan has not been up to the mark, despite the fact that, achievement of universal immunization in the country remained a priority national agenda since long; yet, the success achieved is below par, and much work is still to be done^[@CR55],[@CR56]^. There has been huge amount of funding into the Polio Eradication Initiative of EPI program in Pakistan; yet, few years back the coverage worsened in some areas due to poor management^[@CR17]^. It is believed by the public health experts that Polio eradication in Pakistan has mainly been driven by international pressure, rather than by the demands of local communities or health staff or management or by the political leadership. The recent surveys have shown improvements in EPI working and handling the situation in collaboration with its supporting organizations like WHO and UNICEF^[@CR21]^, but there is still long journey to travel for the realisation of our health goals.

Our study had certain strengths to its credit. Firstly, this is the first ever study in Pakistan to discuss country-wide vaccination coverage data collected through EPI Pakistan, to our knowledge. Secondly, this is the first study to detect the spatiotemporal variation in EPI vaccination coverage in Pakistan at the district level and explored the low vaccination coverage districts for all the antigens used in routine immunization program in the country. Thirdly, the methods used in this study are simple yet robust and are easy to understand for the readers as well as decision makers. Lastly, our study provided scientific evidence on low vaccination coverage clusters which transcended districts and provincial administrative boundaries so that the policy makers can make decisions, targeting clusters with low vaccination rates.

There were a few limitations in this study as well. Firstly, the data in the study was provided by EPI National cell and may be biased as it was collected by their own staff throughout the country. Secondly, the demographic details of individual vaccinated children (e.g., age, sex, home address etc.) were not available and hence not used in our study. Thirdly, the health staff involved in vaccination activities in different districts had different levels of skills, expertise and motivation and this may have been masked in the finalized data. Fourthly, the retrospective extrapolation of the population at a constant rate of 2.4% across all the districts is a limitation as the population growth in all the districts do not follow a uniform pattern, yet it was the best option as the growth rates for individual districts were not available. Lastly, the data did not reflect vaccination coverage in urban and rural populations separately.

Conclusion {#Sec6}
==========

Our study detected low vaccination coverage clusters for each type of EPI vaccination in Pakistan at the district level, and has helped us understand and visualize immunization trends in Pakistan over 5 years. The desired success in getting rid of the burden of vaccine preventable diseases and improving child health cannot be achieved without attaining homogeneous levels of vaccination coverage throughout the country. The low vaccination coverage districts in FATA, KPK and Balochistan are a major concern; however, at the same time Punjab and Sindh with relatively higher vaccination coverage should not be considered risk free. The efficiency of EPI at national, provincial and district level needs to be enhanced. Decision makers, e.g., political leadership and the public health experts need to design in-depth research to probe the low-vaccination coverage in detail and also should amend the existing policies and implementation strategies to better the current situation.

Methods {#Sec7}
=======

Data collection {#Sec8}
---------------

The district level, vaccination coverage data for each antigen included in the routine immunization program in Pakistan from 2012 to 2016 were collected from the Federal EPI Islamabad Pakistan. Immunization services are delivered to the community by more than 7,000 EPI centres throughout the country with each centre serving an average population of 30,000 individuals approximately. Immunization services are delivered to the community through fixed sites as well as by mobile or outreach approaches^[@CR26]^. The data on vaccination and its coverage is reported from sub-district level to district level, and then to the respective provinces and eventually reaches the Federal EPI via vaccine Logistics Management and Information System (vLMIS)^[@CR26]^.

The data used in this study were from 142 of 146 districts in Pakistan. Two of the districts which were not included in our study (Frontier Bannu and Frontier Lakki Marwat) were from Federally Administered Tribal Areas (FATA). The non-inclusion of the data from these two districts was due to the fact that rest of the districts in FATA (other than these two; Frontier Bannu and Frontier Lakki Marwat) remained conflict struck during most period of the study^[@CR31]^, and the public sector healthcare facilities from these two districts were designated as make-shift vaccination hub catering most of the vaccination load from neighbouring districts. As a result, in trying to compensate for the missed non-vaccinated children the achieved vaccination targets in these districts reached even 150% or sometimes around 200%, which could not be reflected in the data analysis and therefore were excluded. Besides these two districts, district named Panjpai as shown in the map (in Balochistan Province) is an aberration in the map being used for current analysis and no such district exist in Pakistan; also the data from Indian held Kashmir was neither available nor included in our study.

Population data for each of the districts in Pakistan were downloaded from Pakistan Bureau of Statistics' official website depicting latest statistics from census completed in 2017^[@CR57]^. As the last population census in Pakistan was done in 1998 and between these two censuses (1998 to 2017), the national average annual growth rate was estimated to be 2.4%. The 2017 district wise population was retrospectively extrapolated, depreciating by 2.4% annually and gave us the population for each of the preceding year in the study respectively (2016, 2015, 2014, 2013, and 2012). The provinces and administrative units in Pakistan are shown in the supplemental file [S1](#MOESM1){ref-type="media"} for the reference purpose (for further demographic information on provinces in Pakistan see under the heading, "[Methods](#Sec7){ref-type="sec"}" and sub-heading, "Study area" in the following link: <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6025434/)59>.

Data analysis {#Sec9}
-------------

### Data organization, descriptive analysis and sequence plots {#Sec10}

The raw vaccination coverage data by type, obtained from Federal EPI Pakistan were organized in Microsoft Excel as per the requirements of spatial analysis. EPI Pakistan calculates the target population of the eligible children for vaccination on the basis of WHO criteria and the average national population growth rate of the country^[@CR58]^. Children less than 1 year of age are estimated to be 3.5% of the total population which was taken as target population by EPI for the vaccination coverage in Pakistan^[@CR26]^. The sub-district level health staff on some occasions has been found over-reporting the vaccination coverage data to meet their performance indicators^[@CR7]^ and the in-migrations between the districts were also not accounted for in the 3.5% target population. These two factors were considered to produce bias in the reported numbers of vaccinated children. To control this bias, we increased the target population by 1%; i.e., the target population in our study was set to be at 4.5% rather than 3.5%. The explanation for this adjustment is that, variation in the vaccination coverage in Pakistan reported in the literature has been between 40 and 80% (mean of which is 60% i.e. maximum reported vaccination coverage is reduced by 20%). The average difference between vaccination rates of BCG for 3.5% and 4.5% target population respectively for all the districts for each study year was calculated to be (19 to 21%) \~ 20%. Thus, the bias in our data was minimized by increasing the target population by 1% (3.5% to 4.5%) and the administrative data rates were reduced on the average by 20% for all types of vaccination included in this study.

The annual vaccination rate (VR) per 100 children, separately for each year and for each type of vaccination (BCG, OPV3, Penta3 and Mealsles1) were calculated at the district level. The descriptive statistics for vaccination rate of each type of vaccination were calculated for all 5 years of study separately (2012 to 2016). All the calculated summary statistics included mean values, standard deviation, minimum and maximum values, and quartiles of 25, 50 and 75 respectively.

Sequence plots for each type of vaccination coverage were plotted with years on the X-axis and vaccination rates by type on the Y-axis at the provincial levels and displayed in figure. The descriptive statistics and sequence plots were performed using statistical package SPSS 13.0^[@CR59]^.

Spatiotemporal analysis {#Sec11}
-----------------------

GeoDa 1.8 software, a GIS tool used for creating and utilizing maps, discovering and sharing geographic information, compiling geographic data and analysing mapped information, was applied to create maps of spatiotemporal pattern of vaccination coverage^[@CR60]^.

For spatial cluster analysis, the already organized data from Microsoft Excel sheet containing annual vaccination coverage data by type and district wise annual population and coordinates of centroids (latitudes and longitudes) were then run in the SaTScan 9.4.4^[@CR61]^. SaTScan is applied to detect spatial clusters of diseases and for distribution of resources like interventions against any specific diseases^[@CR28]^. SaTScan has also been used to detect the spatial heterogeneity in the vaccination coverage in the field of epidemiology.

In current study, for the detection of low-rates of vaccination coverage by type, we performed purely spatial scanning analysis using the discrete Poisson model for the districts' population at risk. SaTScan has the provision of performing Poisson based modelling for the number of events in a fixed interval of time in any geographical setting having non-random and unique coordinates (assigned by researcher). The underlying spatial heterogeneities are also adjusted for the population while running this discrete Poisson model in SaTScan. To avoid overlapping of primary clusters by the secondary ones, selection of "No geographic overlap" was done as a default setting. Initially, we performed scan statistics using different settings of maximum radii for detection of low coverage clusters (e.g. 100 km, 200 km, 300 km and 400 km), and 5%, 10% and 15% for maximum population at risk and ultimately the final parameters were selected having cluster radii limit of 200 km and the limit for spatial cluster size to be for maximum 10% of the population at risk. Another spatial cluster study had adopted the same settings criteria using SaTScan software^[@CR62]^. By using these finalized settings, primary and secondary clusters of low vaccination coverage by type were identified. Primary clusters were the most statistically significant clusters with low vaccination coverage; while, the secondary clusters were additional statistically significant clusters with low vaccination coverage by type in each year from 2012 to 2016.

The spatial patterns of the low vaccination coverage by type exhibited by their respective vaccination coverage rates (percentage) and all the derived information from SaTScan (information on primary and secondary clusters of low vaccination coverage by type, 2012--2016) were then displayed via maps using GeoDa 1.8 software.

Ethical clearance {#Sec12}
-----------------

For this study, vaccination coverage by type was obtained from Federal EPI Cell via letter no. {F.No EPI/Gen/1-2016} dated 9th Feb 2018. An exemption approval letter was obtained from National Bioethical Committee Pakistan via letter number No.4-87/NBC-279-Exempt./17/1139. The research protocol was approved by Institutional Review Board of Xi'an Jiaotong University, China.

Supplementary information
=========================

 {#Sec13}

Supplementary file1 (PDF 308 kb)

BCG

:   Bacillus Calmette Guerin

DTP

:   Diphtheria tetanus pertussis

EPI

:   Expanded program on immunization

FATA

:   Federally Administered Tribal Areas

GIS

:   Geographical Information System

KPK

:   Khyber Pakhtunkhwa

NADRA

:   National Database Authority Pakistan

OPV

:   Oral polio vaccine

PDHS

:   Pakistan Demographic and Health Survey

PSLM

:   Pakistan social and living standards measurement survey

UNICEF

:   United Nations Children Fund

VR

:   Vaccine rate

vLMIS

:   Vaccine Logistics Management Information System

WHO

:   World Health Organization

**Publisher\'s note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information
=========================

is available for this paper at 10.1038/s41598-020-67839-0.

We are highly thankful to Federal EPI Pakistan for providing us the EPI Vaccination Coverage data for this research. G.Z. was funded by The National T&S Major Project of China (Approval No.: 2018ZX10721202); M.F.U. and S.Z. were funded by Chinese Scholarship Council (CSC) for their PhD studies and Masters studies, respectively at Xi'an Jiaotong University, China; and X.Q. was funded by China Postdoctoral Science Foundation (Approval No.: 2015M570842) and National Natural Science Foundation of China (Approval No.: 81502892).

M.F.U., S.Z., X.Q. and G.Z. designed the study. M.F.U. and S.Z. acquired the data, processed the data, implemented all statistical analysis and drafted the manuscript. G.Z. and X.Q. supervised the research and revised the study protocol. W.H. provided statistical assistance. All the authors read and approved the final manuscript.

The Federal EPI Cell, provided vaccination coverage data and took undertaking from the authors to use these data exclusively for this study. As, we are not authorized to share it, so, the data is unavailable to access.

The authors declare no competing interests.
