Association between pertussis vaccination coverage and other sociodemographic factors and pertussis incidence using surveillance data

2.1. Pertussis cases

Pertussis case data was obtained from the Public Health Seattle and King County (PHSKC) Department of Communicable Diseases and Immunizations. Pertussis is a nationally notifiable disease and cases are reported to PHSKC by primary care providers and diagnostic laboratories within 24 h of detection. The clinical case definition of pertussis used was a cough illness lasting 2 weeks or more with at least 1 of the following: paroxysms of coughing or inspiratory “whoop”, posttussive vomiting, or apnea (with or without cyanosis) for infants up to 1 year of age (CDC, 2014). Suspected, probable, and confirmed pertussis cases of all ages diagnosed in King County between January 1, 2010, and December 31, 2017 were included in this analysis. A confirmed case is defined as a case of acute cough illness of any duration with isolation of B. pertussis from a clinical specimen or polymerase chain reaction (PCR) positive for B. pertussis. A probable case is defined as (in absence of a more likely diagnosis) an illness meeting the clinical criteria, or an illness with cough of any duration with at least one of the clinical symptoms and contact with a laboratory confirmed case (epidemiologic link). Suspected pertussis cases are cases with cough lasting ≥ 2 weeks with no other symptoms, or cough of any duration with one of the casedefining symptoms without lab confirmation or epidemiologic link, or an epidemiologic link with cough of any duration and no other symptoms and no lab confirmation, or PCR positive for B. pertussis but no documentation of cough or case-defining symptoms (CDC, 2014; Washington State Department of Health, 2016). Demographic data (date of birth, gender, home address), date of diagnosis, vaccination status, number of DTaP doses received, and date of last DTaP dose was available for all reported cases.

2.2. Geocoding and aggregating pertussis cases at school-district level

We used ArcGIS 10.1 for geocoding cases’ home addresses (ESRI, 2012) and aggregated pertussis cases at the school district level. When the street address was missing or incorrect, we used zipcode for geocoding. School district level population estimates for the census year 2010 were obtained from the National Historical Geographical Information System (Appendix A Table 1) (Manson et al., 2017). A map of King County school districts (N=18) is in Appendix B Fig 1.

2.3. Estimating vaccination coverage at school district level

We obtained DTaP immunization records for children 0–9 years old born or living in King County, Washington between January 1, 2008, and December 31, 2017, from the Washington State Immunization Information System (WA-IIS). WA-IIS is a lifetime registry that tracks immunization records for people of all ages in Washington State (Washington State Department of Health, 2015). Healthcare providers such as primary care physicians, hospitals, and healthcare plans voluntarily report patient immunizations to the WA-IIS. Additionally, birth certificates of children born in King County are loaded into the registry fortnightly. The study cohort was restricted to children born in King County after 2008 to ensure data completeness and accuracy. Ninety-nine percent of children aged 4 months - 5 years have 2 or more immunizations recorded in the WA-IIS (Washington State Department of Health, 2015). Using WA-IIS data, we created a retrospective cohort of 316,404 children aged 3 months to 9 years. Vaccine name and date of receipt for all pediatric vaccines recommended from birth through 9 years of age was available for each child. Demographic information included date of birth, sex, current residential address, residential address at birth, county of residence, and insurance information. Home addresses (or zip codes when home addresses were incorrect or not available) of WA-IIS participants were geocoded by the WA Department of Health staff. We spatially overlaid the geocoded home addresses of WA-IIS participants onto a shapefile of King County school districts to obtain school district of residence for each child. DTaP doses are recommended at ages 2, 4, and 6 months of age, with booster doses at ages 15–18 months and 4–6 years (Havers et al., 2020). We used participants’ dates of birth and dates of DTaP vaccination to calculate age-appropriate vaccination status for each DTaP dose for each child. We estimated annual DTaP vaccination coverage at the school district level as proportion of 19–35 month old children living in a school district with ≥4 DTaP doses in a given year from 2011 to 2017 (the 2008 birth cohort turned 35 months old in 2011). We assumed that vaccine coverage in 2010 was the same as in 2011. We obtained school district level sociodemographic factors, namely, the proportion of population in each school district that are foreign born, White, non-citizens, speak a language other than English at home, proportion of population over 16 years old that have not completed high school, proportion of households with more than 4 people living in them, and median income, from the 2010 US Census data (U.S. Census Bureau, 2010). This study was reviewed and approved by the Washington State Institutional Review Board and PHSKC Research Administration Review Committee.

2.4. Notation

Here we present some notation and assumptions common to the ecological vaccine model and endemic–epidemic model. Let $Y_{it}$ and $N_{it}$ be the number of cases and the total population in school-district $i$ at time $t$. Total population for King County summed over all school districts is given by $N$. Let $x_{it}$ be the time-varying vaccination coverage estimated as proportion of 19–35 month old children vaccinated against pertussis with ≥ 4 DTaP doses. Let ${\lambda }_{it}$ be the force of infection, i.e., risk of infection at time $t$ for an individual who was susceptible at time $t-1$. For our analysis, we assume a time step of four weeks or a month which is the approximate generation time for pertussis (Vynnycky and White, 2010). Assuming a constant hazard rate between time steps, the probability that a susceptible individual at time $t-1$ gets infected at time $t$ is given by ${\lambda }_{it}$. Assuming that time until infection is independent for all susceptible individuals (Halloran et al., 2010), the number of new cases in area $i$ at time $t$ can be modeled as:

Assuming ${\lambda }_{it}$ is small, $1-e^{-{\lambda }_{it}}\approx {\lambda }_{it}$. When the number of infections is small and the population is large, a Poisson distribution approximates the Binomial distribution. Thus, eq. (1) can be written as $Y_{it}\mid {\lambda }_{it}\sim Poisson\left(S_{i,t-1}{\lambda }_{it}\right)$. Let ${\mu }_{it}=S_{i,t-1}{\lambda }_{it}$, then a general form of eq. (1) is,

Both the endemic–epidemic and ecological vaccine models assume that the number of infections is negligible compared to the number of susceptibles. Under this assumption, the number of susceptibles at time $t,S_{it}$, can be approximated by the initial number of susceptibles, thus $S_{it}\sim N_{it}$. We assume that the population in school districts remains constant over time, i.e. $N_{it}\sim N_i$. Because this is a partially vaccinated population, the initial number of susceptibles is given by $\left(1-x_{it}\right)N_i$ for the endemic–epidemic model and $\left(1-\varphi x_{it}\right)N_i$ for the ecological vaccine model, where $x_{it}$ is the vaccine coverage in school district $i$ at time $t$ and $\varphi$ is the direct vaccine effectiveness on susceptibility (reduction in individual’s risk of disease after vaccination). Thus, depletion of susceptibles is not explicitly modeled in either model. In the context of the ecological vaccine model, CI stands for credible intervals as this model uses the Bayesian framework. For the endemic–epidemic models, CI denotes confidence intervals.

2.5. The ecological vaccine model

Infectious disease surveillance data is often available aggregated over space and time. Using ecological regression models for infectious diseases data can be problematic because they do not account for the dependent nature of infectious disease data and can introduce ecological bias. The main feature of the ecological vaccine model developed by Fisher and Wakefield (Fisher and Wakefield, 2020) is that it reduces ecological bias when using area-level infectious disease data and its primary goal is inference. The authors started with an individual-level infectious disease model with vaccination as a parameter to model how vaccine reduces an individual’s risk of infection, under two modes of vaccine action: leaky and all-or-none (Halloran et al., 2010). In the leaky model, the individual’s risk of infection is reduced by a constant proportion for all vaccination individuals while in the all-or-none model (or primary vaccine failure), vaccinated individuals are fully protected against infection but the vaccine fails to take in some individuals. These individual-level vaccine models are then aggregated to area level to give ecological vaccine models. Under certain simplifying assumptions, the authors constructed an ecologically consistent model from an individual level model accounting for vaccine coverage. The detailed derivation of the ecological vaccine model can be found here (Fisher and Wakefield, 2020). The impact of vaccination in this model is defined as the ability of the vaccine to reduce susceptibility against infection. This ecological vaccine model is fit using the Bayesian framework in R package <monospace>rstan</monospace> (Guo et al., 2020) and provides posterior estimates and corresponding posterior credible intervals of epidemiologically relevant parameters.

In a partially vaccinated population, let $\varphi$ be the reduction in a vaccine recipient’s risk of infection, which can be interpreted as the vaccine effectiveness after ≥ 4 DTaP vaccine doses (Fisher and Wakefield, 2020). Given that $x_{it}$ is the vaccine coverage in school district $i$ at time $t$, the number of susceptibles in school district $i$ can be written as $\left(1-\varphi x_{it}\right)N_i$. Using the Bayesian framework, we can incorporate prior knowledge about pertussis vaccine effectiveness into the ecological vaccine model. Randomized controlled trials have estimated the vaccine efficacy of the DTaP vaccine to be ~ 85% (95% CI: 81%, 87% ) (Fulton et al., 2016). We fit the ecological vaccine model as:

where ${\mu }_{it}$ is the total risk (epidemic plus endemic risk) and $v_{it}$ and ${\lambda }_i$ are the endemic and epidemic pertussis risk components. The school district specific random effects $a_i$ and $b_i$ are assumed to be independent; ${\omega }_t=\frac{2\pi t}{13}$; an informative beta prior Beta (10, 2.5) was used for $\varphi$ with median 0.78 and 90% of the mass between 0.66 and 0.99. We assumed normal priors with mean 0 and variance 5 for ${\alpha }_{AR}$ and ${\alpha }_{EN}$ and variance 10 for $\gamma$ and $\delta$. We assumed frequency-dependent transmission in the formulation of ${\lambda }_i$ (Keeling and Rohani, 2008). Hamiltonian Monte Carlo sampling via R package <monospace>rstan</monospace> was used to fit this model (Guo et al., 2020). We adapted code published previously for an ecological vaccine model for measles data to include time-varying vaccination coverage (Fisher and Wakefield, 2020). As a sensitivity analysis, we fit the same model with non-informative priors on $\varphi$, with median 0.50 and 90% of the mass is between 0.05 and 0.95, to check the influence of priors on the estimate of vaccine effectiveness.

We also computed school-district specific time-varying autoregressive components and their 95% credible intervals from the ecological vaccine model. This parameter may be interpreted as the effective reproductive number, $R_{eff}$, which is defined as the average number of new cases per infectious case in a partially vaccinated population (Fisher and Wakefield, 2020). The time-varying autoregressive parameters were calculated as:

Fitted values were calculated as:

where $Y_{i,t-1}$ was the observed number of cases in school $i$ and month $t-1$.

2.6. The endemic–epidemic model

The endemic–epidemic model is motivated by the Poisson branching process with immigration. Total pertussis incidence ${\mu }_{it}$ is split into two components: the endemic component with rate $v_{it}$ and the epidemic component with rate ${\lambda }_{it}{Y}_{it-1}$, such that ${\mu }_{it}=v_{it}+{\lambda }_{it}{Y}_{it-1}$, where $Y_{it-1}$ is the observed pertussis count in school district $i$ and month $t-1$ (Held et al., 2005). The epidemic component can be further decomposed into the autoregressive and neighborhood components. In a model with spatial data, the autoregressive component models cases arising from infected individuals from the same area, while the neighborhood component captures cases arising from infected individuals in neighboring areas. The endemic component represents the background number of cases or remaining cases not explained by these two components. In this analysis, the endemic component captures incidence in an area related to sociodemographic factors. Each component can be modeled with a log linear model with covariates and fixed or random effects (Meyer et al., 2016). The model can be fit in the R <monospace>surveillance</monospace> package (Höhle et al., 2011) and likelihood estimation is done using the quasi-Newton algorithm. When random effects are included, penalized and marginal log-likelihoods are maximized alternately until convergence (Meyer et al., 2016).

In the simple model (Model 1), the endemic component included the population as an offset modeled as the fraction of the population that live in school district $i$ denoted by $e_i$. To account for the temporal variance in incidence, the endemic component included an overall linear trend and a sinusoidal wave of frequency ${\omega }_t=\frac{2\pi t}{13}$. The endemic component is written as :

where $v_{it}$ is the endemic risk of pertussis, ${\alpha }_{EN}$ is the endemic intercept, assumed to be constant over the region, and $\beta$ is the parameter associated with temporal trend.

We included vaccination coverage in the epidemic component of the model as log proportion of children 19–35 months old that had fewer than 4 doses of DTaP vaccine because we are interested in the effect of pertussis vaccination on size and occurrence of pertussis epidemics. Effect of vaccination coverage on disease incidence has been modeled similarly for measles using the endemic–epidemic model before (Herzog et al., 2011). The epidemic component is written as:

where ${\lambda }_{it}$ is the epidemic pertussis risk, ${\alpha }_{AR}$ is the epidemic intercept, $x_{it}$ is the vaccination coverage in school district $i$ at time $t$, and ${\beta }_v$ is the parameter associated with under-vaccination. The proportion of 19–35 month old children under-vaccinated or susceptible in the school-districts is $1-x_{it}$. The intercept ${\alpha }_{AR}$ is assumed to be constant over all areas. In this simple endemic-epidemic model, we do not include a neighborhood component to be able to compare the ${\alpha }_{AR}$ and ${\alpha }_{EN}$ estimates with estimates from the ecological vaccine model. Overall, Model 1 can be written as:

where ${\mu }_{it}$ is the total pertussis risk, $e_i$ is the population fraction in school district $i$ used as a multiplicative offset in the endemic component, and $Y_{i,t-1}$ is the observed number of cases in school district $i$ at time $t-1$.

We fit a separate endemic–epidemic model (Model 2) with sociodemographic covariates in the endemic component and the epidemic component split into autoregressive and neighborhood components to measure spatio-temporal dependence. School districts that shared a boundary were defined as neighbors. A matrix of transmission weights, $w_{ji}$, which represent the flow of cases from school district $j$ to school district $i$ when $j\ne i$, was included in the neighborhood component. The model assumes that the epidemic can only arrive from adjacent areas. Thus, if two school districs share a boundary, the assigned weight was 1, otherwise 0. To reflect that people likely commute to densely populated metropolitan areas, we scaled the school district’s risk with respect to its population fraction, $e_i$. The neighborhood component can be written as:

where ${\zeta }_i$ is the neighborhood associated risk of pertussis, ${\alpha }_{NE}$ is the intercept associated with the neighborhood component, and ${\beta }_{Pop}$ is the parameter associated with the population fraction of school-district $e_i$.

Let ${\beta }_z$ be a vector of parameters associated with a vector of school-district level sociodemographic covariates denoted by ${\beta }_i$, namely the proportion of population in each school district that are foreign born, White, non-citizens, speak a language other than English at home, proportion of population over 16 years old that have not completed high school, proportion of households with more than 4 people living in them, and median income. These factors included in the model were selected a priori based on previous research. Household size was selected as it could potentially impact pertussis transmission (Magpantay and Rohani, 2015). Median income, proportion of residents who were White, foreign-born, non-citizens, non-native English speakers, and who did not complete high school were selected as proxies for socioeconomic status and have been previously associated with vaccination coverage (Hill et al., 2016, 2018; Sánchez-González et al., 2017; Larson et al., 2014; Varan et al., 2017). Combining the endemic, epidemic, and neighborhood components, we get Model 2 as:

where $Y_{j,t-1}$ is the number of cases in school-district $j$ at time $t-1$.

3.1. Descriptive analysis

Between 2010 and 2017, 1885 pertussis cases of all ages were reported in the 18 school districts in King County, WA. There was one large epidemic in 2012 with 894 cases and a smaller one in 2015 with 250 cases. The largest number of cases in a school district in a single month was 50 during the 2012 epidemic and occurred in Seattle school district, which also recorded the highest number of cases overall (n = 481). Because pertussis is a rare disease, many school districts recorded zero cases during several months. Tukwila school district recorded only 3 cases over the span of 10 years. Pertussis incidence per 100,000 by school district is in Fig. 1a.

Data from WA-IIS was used to estimate vaccine coverage at the school-district level which is displayed in Fig. 1b. Even though we used 19–35 month olds to estimate vaccination coverage with ≥ 4 DTaP doses, we assumed that the vaccine coverage for the entire population of King County is the same as the coverage estimated for this analysis. Vaccination coverage is higher in school districts in northern and northeastern King County, compared to school districts in the south. Vashon Island (in black) has the lowest vaccine coverage of all school districts. Within each school district, vaccine coverage appeared to increase with time between 2010 and 2017 (Table 1).

Pearson correlation coefficients showing correlations between DTaP vaccine coverage among 19–35 month old children and disease incidence calculated over a rolling 3-year window are in Table 2. With the exception of 2015, vaccine coverage with ≥ 4 doses of DTaP was negatively correlated with pertussis incidence. This association was strong and statistically significant in the years 2012 (R: −0.61, 95% CI: −0.20,−0.84), 2013 (R: −0.61, 95% CI: −0.19,−0.83) and 2014 (R ∶ −0.60, 95% CI: −0.19,−0.83). Thus, there is some indication that areas with higher vaccine coverage showed lower disease incidence, especially during a period of high incidence.

3.2. Ecological vaccine model

Summaries of posterior medians of fixed effects and 95% credible intervals from the ecological vaccine model are in Table 3. Using a strong prior for the vaccine effectiveness $\varphi$, we estimated it to be 83% (95% CI: 63%, 95%). Thus, the vaccine effectiveness associated with receiving ≥ 4 doses of DTaP compared to receiving < 4 DTaP doses is statistically significant. With a strong prior on $\varphi$, our results agree with estimates of efficacy of DTaP vaccine found in the literature (~85%) (Fulton et al., 2016). As a sensitivity analysis, we ran the same model using a uniform prior for $\varphi$ and vaccine effect was estimated to be 79% (95% CI: 33%, 96%) (Appendix A Table 2). The uniform prior resulted in a slightly lower estimate for $\varphi$, but credible intervals were wider. This vaccine effectiveness estimate does not differentiate between primary vaccine failure and failure due to leakiness. With a uniform prior, the estimate of epidemic intercept ${\alpha }_{AR}$ in this model was also smaller with wider credible intervals. Using a uniform prior did not change the endemic intercept ${\alpha }_{EN}$ by much.

The school-district-specific time-varying $R_{\text{eff}}$ are plotted in Appendix C Fig 1. All estimates are below 1 and vary slightly with time within school-districts. There was no apparent effect of DTaP vaccine coverage on $R_{eff}$. Due to small number of cases within school-districts, the credible intervals of $R_{eff}$ were quite wide.

Fig. 2 shows the observed number of cases, incidence per 100,000 people, and model fits for each school district obtained from the ecological vaccine model. District-specific estimates of $R_{eff}$ and endemic risk $v_{it}$ are given in each panel. The model seems to fit the data well especially for areas with a large number of cases such as Federal Way, Kent, Lake Washington, and Seattle. The correlation coefficients between average DTaP vaccination coverage (measured as total children 19–35 months old with ≥ 4 doses of DTaP between 2010 and 2017 divided by total children 19–35 months old for each school district) and $R_{\text{eff}}$ $(r=-0.05$; 95% CI: −0.50, 0.42) and average DTaP vaccine coverage and endemic risk $v_{it}(r=0.18$; 95% CI: −0.31, 0.59) are small and not statistically significant. Population density (measured as persons per square mile) is positively correlated with both $R_{eff}(r=0.24$; 95% CI: −0.25, 0.63) and endemic risk ($r=0.58$; 95% CI: 0.16, 0.82). Thus, $R_{eff}$ in this study was not statistically significantly correlated with either vaccination coverage or population density. Endemic pertussis risk was significantly correlated with population density.

Appendix B Fig 2(a) and 2(b) show area-specific random effects of the autoregressive and endemic components. Areas with a large number of cases have larger autoregressive random effects. However, no such structure was found in the endemic random effects. There was no evidence of correlation between endemic and autoregressive random effects (Appendix B Fig. 3) which supports our decision to use independent random effects in our model.

3.3. Endemic–epidemic models

We will first discuss the results of the simple model with DTaP vaccination coverage in 19–35 month old children in the epidemic component as the only covariate and no neighborhood component or other demographic covariates included. The estimate of the exponentiated autoregressive intercept from this model ($\exp \left({\alpha }_{AR}\right)=0.07$; 95% CI: 0.003, 1.43) is lower compared to the autoregressive intercept estimate from the ecological vaccine model ($\mathrm{e}\mathrm{x}\mathrm{p}\left({\alpha }_{AR}\right)=1.1$; 95% CI: 0.75, 1.52) and has very wide confidence intervals (Table 4). The endemic estimate from the endemic-epidemic model (2.23; 95% CI: 1.99, 2.47) is lower than that from the ecological vaccine model (3.16; 95%CI: 2.77, 3.50). There is no comparable estimate to vaccine effect $\varphi$ in this model. Here, the effect of DTaP vaccine coverage is estimated as $2^{\mathrm{e}\mathrm{x}\mathrm{p}\left({\beta }_v\right)}$. It is interpreted as for each doubling of pertussis undervaccination rate (or doubling of susceptible population), the epidemic risk of pertussis increases multiplicatively by 3.54 fold (95% CI: 1.65, 23.05). Thus, this model suggests that epidemic pertussis risk was statistically significantly associated with DTaP vaccination coverage in 19–35 month old children. There is no strong temporal or seasonal trend in the data (Hill et al., 2016). A statistically significant overdispersion parameter suggests that using the negative binomial distribution was a more appropriate choice for this model than Poisson distribution.

We fit a second endemic–epidemic model where we split the epidemic component into the neighborhood and autoregressive components and added sociodemographic factors in the endemic component to estimate the effect of sociodemographic factors on endemic pertussis risk. For this model, for each doubling of under-vaccination rate (or doubling of susceptible population), the epidemic risk of pertussis increased multiplicatively by 2.76 fold higher (95% CI: 1.44, 16.6) (Table 4). This result is consistent with the under-vaccination coverage estimate from the simple model. Household size and median income were statistically significantly associated with endemic pertussis risk. For every unit increase in proportion of households with more than 4 individuals, the endemic risk of pertussis increased by 15% (95% CI: 8.3%, 20.9% ), adjusting for other covariates. For every $10,000 increase in median income, endemic pertussis risk decreased by 45% (95% CI: 25%, 60%). These are ecological associations and may not apply at the individual level. No significant association was found between pertussis risk and proportion of residents who were White, foreign born, non-citizens, spoke a language other than English at home, and had education less than high school.

The endemic–epidemic model appears to fit the data well (Fig. 3). According to this model, 32.8% of the time-averaged mean pertussis risk is explained by the endemic component, 34% by the autoregressive component, and 33% by the neighborhood component. The model suggests that a large proportion of cases in Renton, Lake Washington, and Bellevue school districts come from neighboring areas. This might be because these school districts have many neighbors and are also highly populated. On the other hand, Vashon Island, being an island (Appendix B Fig 1), has most of its incidence explained by the autoregressive component.