Introduction and Exploratory Data Analysis

Introduction

Pertussis, more commonly known as whooping cough, is a disease chacterized by a serious cough followed by an inhale that sounds like a “whoop”[1]. It is usually seen in unvaccinated people or children who have not recieved all doses of the vaccine. Whooping cough cases are typically pretty low with small fluctuations. However, in 2024, there was an outbreak of whooping cough. Experts believed the outbreak may have been caused by an increase in vaccination hesitation or disruptions to children’s normal vaccination schedule due to the COVID-19 pandemic[2].

Our goal is to analyze whooping cough data using ARMA and POMP modeling techniques. We focus specifically on whooping cough cases in the East North Central states, which include Michigan, Ohio, Indiana, Illinois, and Wisconsin. Our data was obtained from the CDC’s National Notifiable Diseases Surveillance System database (NNDSS)[3] and contains weekly data from 2017 to early 2025.

Exploratory Data Analysis

The graph below displays the weekly cases for all five east north central states from 2017-2025.

Prior to 2024, cases remained fairly low with some small peaks. We see a large outbreak in 2024, where cases were more than twice as high as the previously highest peak.

Unfortunately, our dataset is missing over two full years of data. All of 2022 is missing, and there is very little data for 2021 or the later half of 2020. In total, there are 127 missing weeks. Some of these were orignally included in the dataset with a missing value for the number of cases; while others did not have a corresponding row in the dataset at all. We believe this was due to a reporting error by the CDC, combined with pandemic-ero reporting issues. This issue is not limited to whooping cough, as other diseases in the NNDSS database are also missing for this time period.

We also collected data on births, deaths, population, and vaccination rates during this time period. The weekly death count was obtained from the provisional CDC tables[4]. Weekly data was available for 2018-2025. For 2017, we assumed deaths occured at a constant rate across the year, and divided the yearly number of deaths by 52. Vaccination data was obtained from the Michigan county immunization report cards[5]. We were unable to find sufficient vaccination data for the other four states, so we assume that rates were similar across all five states.

There is a clear decrease of the vaccine coverage rate in 2020, which persists through 2024. This is likely due to the COVID-19 pandemic, which upset typical medical appointments. The full sequence of whooping cough vaccines requires five doses and so it is unsuprising that the pandemic resulted in a decrease in the number of fully vaccinated people. Also, attitudes toward vaccines have also shifted away from vaccines. We will attempt to incorporate this data into our models.

We include exploratory analysis of births and deaths here as well, although we did not include them in our model due to time constraints. We do not expect this to have significantly changed our model. Weekly births were fairly consistent, with a similar cyclical trend in each year. There is a slight decrease over time. Deaths also remained fairly constant, again with a slight cyclical pattern, except for 2021 and 2022. This was also likely due to the COVID-19 pandemic which would have resulted in additional deaths. The drop seen at the end of the dataset is likely due to a delay in releasing provisional mortality counts for some states.

ARMA

ARMA model

We take the first difference to remove trends and center the series around a constant mean (≈0), which is required for valid ARMA modeling. The differenced series then behaves like a random walk with some conditional variance.

Below, we plot the ACF of the differenced pertussis time series.

After differencing, both the ACF and PACF show a single significant spike at lag-1, with most other lags lying within the confidence bounds.

ARMA Model selection

We calculated the Akaike Information Criterion (AIC) for each ARMA(p,q) models for p<5 and q<5[6].

MA0 MA1 MA2 MA3 MA4 MA5
AR0 2788.33 2758.98 2757.76 2753.54 2751.44 2750.57
AR1 2765.01 2756.23 2758.18 2753.41 2750.82 2752.32
AR2 2764.88 2758.12 2755.34 2739.09 2752.00 2742.69
AR3 2750.67 2751.74 2753.74 2745.95 2747.10 2750.01
AR4 2751.85 2753.74 2754.78 2746.99 2748.76 2750.75

Although ARMA(2,4) outperformed the other specifications on the differenced series, the estimation still experienced convergence problems—most likely a consequence of the series’ time‑varying variance.

While the ARMA residuals oscillate randomly around zero, they display non‑constant variance over time. Therefore, we’ll fit an ARCH model to explicitly capture and forecast this conditional volatility.

ARCH Model

An ARCH specification naturally captures time‑varying volatility by making the variance depend on past shocks and outbreak magnitude. By modeling conditional heteroskedasticity, we can more accurately forecast periods of heightened uncertainty and better inform resource allocation during outbreaks[7].

To capture the variability in pertussis case differences, we fit an ARCH(1) model in which the conditional variance depends both on the previous period’s shock and on last period’s case count as an external regressor. We assume that the time series of differenced pertussis cases, \(y_t\), can be decomposed as \[ y_t = \mu + \epsilon_t, \]

where \(\epsilon_t \sim N(0, \sigma_t^2)\) is an error term with time-varying variance. We model the conditional variance \(\sigma_t^2\) using an ARCH-X framework, in which the variance depends not only on past errors but also on an external predictor. Specifically, we assume \[ \sigma_t^2 = \omega + \alpha\,\epsilon_{t-1}^2 + \beta\,x_{t-1}, \]

where: - \(\omega > 0\) is a constant, - \(\alpha \ge 0\) is the ARCH parameter capturing the impact of the previous period’s squared error, - \(\beta\) quantifies the effect of the exogenous regressor \(x_{t-1}\), which is the number of pertussis cases at time t-1.

This specification allows the conditional variance of the model to adjust dynamically to both its own recent shock history and to external signals[7].

# x_exog generated as the lagged pertussis time series
x_exog <- c(NA, head(df$Cases, -1))
x_exog <- x_exog[-1]
pertussis_diff_adjusted <- pertussis_diff[-1]

spec <- ugarchspec(
  variance.model = list(
    model = "sGARCH",          
    garchOrder = c(1, 0),      
    external.regressors = matrix(x_exog, ncol = 1)  # using the lagged cases as exogenous
  ),
  mean.model = list(
    armaOrder = c(0, 0),       
    include.mean = TRUE         
  ),
  distribution.model = "norm" 
)

fit <- ugarchfit(spec = spec, data = pertussis_diff)
param_mat <- fit@fit$matcoef

param_df <- data.frame(
  Parameter    = rownames(param_mat),
  Estimate     = param_mat[, 1],
  `Std. Error` = param_mat[, 2],
  `t value`    = param_mat[, 3],
  `Pr(>|t|)`   = param_mat[, 4],
  row.names    = NULL,
  check.names  = FALSE
)

kable(
  param_df,
  caption = "ARCH Model Parameter Estimates",
  digits  = c(0, 4, 4, 2, 3),
  align   = c("l", "r", "r", "r", "r")
)
ARCH Model Parameter Estimates
Parameter Estimate Std. Error t value Pr(>|t|)
mu 0.2730 0.1488 1.83 0.067
omega 0.0000 0.0011 0.00 1.000
alpha1 0.5332 0.1102 4.84 0.000
vxreg1 2.1387 0.2257 9.48 0.000

Substituting our estimated parameters, we have \[ y_t = 0.273 + \varepsilon_t, \] with \[ \sigma_t^2 = 0.533\,\varepsilon_{t-1}^2 \;+\; 2.139\,x_t. \]

Since both the ARMA and ARCH models are fitted to the same differenced series, their log‑likelihoods can be compared directly.

Log‑Likelihood Comparison of Models
Model Log‑Likelihood
ARMA(2,4) -1368.0
ARCH(1) -1203.8

The ARCH model obtains a log-likelihood of -1203, marking a large improvement from our original ARMA models (-1368).

To evaluate our model fit, first look at the ACF of the residuals as well as the squared residuals.

We then applied the ARCH-LM test, which tests

\(H_0\): no remaining ARCH effects (homoskedastic residuals)

\(H_1\): presence of conditional heteroskedasticity[8]

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  z
## Chi-squared = 25.754, df = 12, p-value = 0.01163

We also ran the Ljung–Box test with

\(H_0\): residuals are independently distributed (no serial correlation)

\(H_1\): residuals are not independently distributed[9]

to verify that the mean model has captured all linear dependence.

## 
##  Box-Ljung test
## 
## data:  z
## X-squared = 55.22, df = 12, p-value = 1.653e-07

Despite relatively flat ACFs for both residuals and squared residuals, the ARCH-LM test (p=0.01163) and Ljung–Box test (p<0.05) both reject their null hypotheses, revealing persistent volatility clustering and serial correlation; this suggests the ARCH(1) model fails to capture the full mean–variance dynamics, likely due to unmodeled external factors.

Because pertussis counts are discrete and over‑dispersed, a continuous ARCH model can miss important epidemic dynamics, as our diagnostics indicate. Therefore, we implemented a compartmental SIR framework to better capture these features.

POMP

SIR Model

We first begin by building a simple SIR compartment model[6] for the outbreak of whooping cough that occurred during 2024 where:

  • \(S\): Susceptible
  • \(I\): Infectious
  • \(R\): Recovered
  • \(C\): Number of reported whooping cough cases
  • \(\mu_{SI}\): Rate at which individuals in S transition to I
  • \(\mu_{IR}\): Rate at which individuals in I transition to R
  • \(\rho\): Rate at which cases are reported
SIR Model Diagram

SIR Model Diagram

We assume that the population size, \(N\), remained constant at 48 million individuals. The number of people in each compartment is computed as follows: \[ \begin{aligned} S(t) &= S_0 - N_{SI}(t),\\[6pt] I(t) &= I_0 + N_{SI}(t) - N_{IR},\\[6pt] R(t) &= R_0 + N_{IR}(t) \end{aligned} \] We model the number of people transitioning from one compartment to the next as follows: \[ \begin{aligned} \frac{dN_{SI}}{dt} \sim Binomial(S(t), 1- exp(\frac{-\beta I(t)}{N}\Delta t))\\[6pt] \frac{dN_{IR}}{dt} \sim Binomial(I(t), 1- exp(\mu_{IR}\Delta t))\\[6pt] \end{aligned} \] Note that this is for the 2024 outbreak only—we are not attempting to model the smaller fluctuations that occured from 2017-2023. We chose to start with this smaller time period since it was more straightforward and allowed us to avoid the missing data issues present in our data.

We initially use the following parameters to simulate from our model based on general searches: \(\beta=4\), \(mu_{IR}\)=0.5, \(\rho=0.2\), \(\eta=0.135\)

The simulations show that we aren’t capturing the reported whooping cough cases very well. We then perform a local search for tuning our parameters to see if we can better estimate the outbreak.

SEIR Model

We then expand our model to include an exposure compartment and attempt to model all of the cases spanning from 2017 to 2025.

  • \(S\): Susceptible

  • \(E\): Exposed

  • \(I\): Infectious

  • \(R\): Recovered

  • \(C\): Number of reported whooping cough cases

  • \(\mu_{SI}\): Rate at which individuals in S transition to E

  • \(\mu_{EI}\): Rate at which individuals in E transition to I,

  • \(\mu_{IR}\): Rate at which individuals in I transition to R

  • \(\rho\): Rate at which cases are reported

We assume that the population size, \(N\), remained constant at 48 million individuals. The number of people in each compartment is computed as follows:

\[ \begin{aligned} S(t) &= S_0 - N_{SE}(t),\\[6pt] E(t) &= E_0 + N_{SE}(t) - N_{EI}(t),\\[6pt] I(t) &= I_0 + N_{EI}(t) - N_{IR},\\[6pt] R(t) &= R_0 + N_{IR}(t) \end{aligned} \] We model the number of people transitioning from one compartment to the next as follows:

\[ \begin{aligned} \frac{dN_{SE}}{dt} \sim Binomial(S(t), 1- exp(\frac{-\beta(t)I(t)}{N}\Delta t))\\[6pt] \frac{dN_{EI}}{dt} \sim Binomial(E(t), 1- exp(\mu_{EI}\Delta t))\\[6pt] \frac{dN_{IR}}{dt} \sim Binomial(I(t), 1- exp(\mu_{IR}\Delta t))\\[6pt] \end{aligned} \] We model the transmission rate \(\beta(t)\) as a function of time since we see a surge of cases in 2024.

We let \[ \beta(t) = \begin{cases} &\beta_{\text{no outbreak}}, & t < t_{\text{outbreak start}}, \\ &\beta_{\text{outbreak}}, & t \ge t_{\text{outbreak start}}, \end{cases} \] meaning \(\beta\) switches from \(\beta_{\text{no outbreak}}\) to \(\beta_{\text{ outbreak}}\) at time \(t_{\text{outbreak start}}\). We choose \(t_{\text{outbreak start}}\) to be April 7th 2024 (week 332) since this is before the number of cases begin to surge.

Using the following parameters,

  • \(\beta_{\text{no outbreak}} = 4.75\)
  • \(\beta_{\text{outbreak}} = 5.3\)
  • \(\mu_{\mathrm{IR}} = 0.5\)
  • \(\mu_{\mathrm{EI}} = 0.5\)
  • \(N = 48{,}000{,}000\)
  • \(\eta = 0.102\)
  • \(\rho = 0.35\)
  • \(k = 5\)

we obtain our simulation results. Some of our simulations show an outbreak that occurs at a similar time as the reported cases, though we are not estimating the number of cases correctly. We also see that we can simulate the number of cases for when we have missing data for 2022.

SEIR Model for ARCH Comparison

Since we are interested in benchmarking our POMP model against our ARCH model, we reran the SEIR model from before, but removed the first data point to account for the differencing data used in ARCH contains one less data point. This ensure that both models are using the same data which allows us to make fair comparisons.

Performing a global search across the same parameter space as before, we obtained a log likelihood of -1442. Simulating from the parameters we obtained for this log likelihood, we still see that the simulations fail to capture the surge in reported whooping cough cases.

## # A tibble: 1 × 10
##   base_beta outbreak_beta   rho   eta mu_EI mu_IR       N     k loglik loglik.se
##       <dbl>         <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl>  <dbl>     <dbl>
## 1      7.83          7.84 0.288 0.884  1.23  37.9  4.80e7  5.00 -1442.     0.181

Comparisons

ARMA vs. POMP

As our differenced series exhibited clear heteroskedasticity that our initial ARMA model could not capture, we next fitted an ARCH(1) model. This specification improved the log‑likelihood substantially (–1203 vs. –1366), but diagnostic tests still revealed residual dependencies and lingering volatility clustering—signs that important dynamics remained unmodeled. Moreover, pertussis case counts are discrete and often over‑dispersed, so a continuous ARCH model may overlook key epidemic behavior; thus, we explored implementing a compartmental SIR model to better accommodate those features.

For our POMP model, we started with a simple SIR model for a single outbreak (2024). We then modified this model by adding an exposure period, modeling the entire time period from 2017-2025, allowing \(\beta\) to vary with time, and adding vaccination data. All of these actually proved to be worse than the original SIR model. Unsurprisingly, it was much harder to model the entire time period than a single outbreak. We included an exposure period and allowed \(\beta\) to vary with time in a attempt to better model this time period, but it did not work well. Futhermore, the inclusion of vaccination data led the outbreak to die out quickly and did not help as we intended. The inclusion of vaccination data caused the outbreak to die out quickly.

Although the POMP model provides insights into disease transmission dynamics through its compartmental structure, our ARCH model demonstrated superior performance with higher log‑likelihood (–1203 vs. –1442). Since first differencing is a linear transformation with a constant Jacobian, we can directly compare these fits after accounting for the dropped initial term, confirming the ARCH model’s superior grounding. While the POMP framework attempts to model latent states and transmission processes explicitly, its poorer fit suggests potential misspecification of the underlying epidemic mechanics. The ARCH model’s ability to capture conditional variance patterns in the time series without imposing restrictive epidemiological assumptions proved more effective for forecasting pertussis case dynamics in our dataset. However, this advantage should be interpreted with caution, as the ARCH model’s statistical superiority comes at the expense of epidemiological realism and may not generalize beyond our dataset.

Discussion

To improve our ARCH model, future work could shift to a count‐based INGARCH (or negative‐binomial INGARCH) to enforce non‑negativity, and directly handle over dispersion[10]. Introducing nonlinear effects of lagged incidence and volatility dynamics would allow the model to distinguish outbreak from endemic phases more realistically. Finally, adopting asymmetric or heavy‑tailed error distributions could improve the fit for occasional large spikes in case counts.

Our POMP models were unable to model this data well. One idea for future work is to include birth and death data. We initially intended to do this, but ran out of time. We also only had vaccination coverage for children, which may have contributed to the outbreak dying out immediately when using that data. Future work could use vaccination rates for the whole population, which we were unable to find. The missing data from 2022 also made things difficult. We attempted to find alternative data sources to replace the missing values, but were unable to find any data sources that did not have the same issue.

Scholarship

There were two other final projects that looked at whooping cough. The first, from 2020, looked at pertussis cases in Michigan from 2006 to 2011[11]. There were no large outbreaks in this case as cases only ranged from 0 to 21 weekly cases. The second pertussis project, from 2016, looked at a similar population and time period[12]. We looked at a larger population (Michigan, Ohio, Wisconsin, Indiana, and Illinois) over a more interesting time period (2017-2024) that actually had an outbreak.

Both projects focused only on POMP models with slight modifications. The first started with an SIR model and added vaccination rates. The second started with an SEIR model and added birth, death, and vaccination data. We attempted something similar with our POMP models, but we also compared to ARMA models. Our ARMA models also went beyond what was covered in class, which was not commonly seen in other final projects. There were other projects that compared ARMA and POMP that inspired us but they did not explore additional ARMA modifications[13] [14][15] [16].

We also approached our POMP models differently than these previous projects. Since we had a longer time period with one outbreak, we experimented with modeling the entire time period versus just a portion of it. We also attempted to use different parameters for different time periods. Unsurprisingly, our models for the shorter time period that contained the outbreak worked much better than models for the entire time period. However, attempting to model a longer time period containing outbreaks and general low-level disease was not commonly seen in other final projects.

References

References.

1. Whooping cough. (n.d.). Mayo Clinic. https://www.mayoclinic.org/diseases-conditions/whooping-cough/symptoms-causes/syc-20378973
2. O’Neil, E. (2024). Worst u.s. Whooping cough outbreak in a decade has infected thousands. NPR. https://www.npr.org/2024/10/10/nx-s1-5127422/whooping-cough-outbreak-u-s-seattle
3. CDC national notifiable diseases surveillance system. (n.d.). https://data.cdc.gov/browse?category=NNDSS&sortBy=relevance&pageSize=20&page=1
4. Provisional mortality statistics. (n.d.). https://wonder.cdc.gov/mcd-icd10-provisional.html
5. County immunization report cards. (n.d.). Michigan Health; Human Services. https://www.michigan.gov/mdhhs/adult-child-serv/childrenfamilies/immunizations/data-statistics/countyimmsreportcards
6. Ionides, E. (2025). Notes for STATS/DATASCI 531, Modeling and Analysis of Time Series Data. https://ionides.github.io/531w25/
7. Ionides, E. L. (2025). Modeling and analysis of time series data DATASCI/STATS 531 chapter 17: A case study of financial volatility and a POMP model with observations driving latent dynamics. https://ionides.github.io/531w25/17/index.html
8. Engle’s ARCH test. (2024). https://www.mathworks.com/help/econ/engles-arch-test.html
9. Bobbitt, Z. (2020). Ljung-box test: Definition + example. https://www.statology.org/ljung-box-test/
10. Zhu, F. (2011). Modeling time series of counts with COM-poisson INGARCH models. https://www.sciencedirect.com/science/article/pii/S0895717711007540
11. Final project: Analysis of pertussis transmission by SIR model. (2020). https://ionides.github.io/531w20/final_project/
12. Final project: Exploratin of pertussis transmission dynamics and vaccination effect period in michigan. (2016). https://ionides.github.io/531w16/final_project/
13. Final project: Predator-prey model with POMP. (2024). https://ionides.github.io/531w20/final_project/
14. Final project: COVID-19 modeling comparison: ARIMA vs SEIR. (2024). https://ionides.github.io/531w24/final_project/project04/blinded.html
15. Final project: ARMA and POMP analysis on COVID-19 variants in the US. (2022). https://ionides.github.io/531w22/final_project/project21/blinded.html
16. Final project: Modeling COVID-19 cases in michigan: ARMA model v.s. SEIR POMP model. (2021). https://ionides.github.io/531w21/final_project/project11/blinded.html