Integer-Valued Functional Data Analysis for Measles Forecasting

TL;DR

This paper introduces a novel integer-valued functional time series model for measles count forecasting, capturing complex seasonality and providing accurate uncertainty quantification, which improves multi-month ahead predictions and peak timing estimates.

Contribution

The paper develops a new Bayesian methodology modeling measles counts as an integer-valued functional time series with dynamic basis expansion and overdispersion handling, enhancing forecast accuracy and uncertainty quantification.

Findings

Improved multi-month ahead forecasts with tighter intervals

Accurate uncertainty quantification for peak timing

Enhanced modeling of dynamic seasonal patterns

Abstract

Measles presents a unique and imminent challenge for epidemiologists and public health officials: the disease is highly contagious, yet vaccination rates are declining precipitously in many localities. Consequently, the risk of a measles outbreak continues to rise. To improve preparedness, we study historical measles data both pre- and post-vaccine, and design new methodology to forecast measles counts with uncertainty quantification. We propose to model the disease counts as an integer-valued functional time series: measles counts are a function of time-of-year and time-ordered by year. The counts are modeled using a negative-binomial distribution conditional on a real-valued latent process, which accounts for the overdispersion observed in the data. The latent process is decomposed using an unknown basis expansion, which is learned from the data, with dynamic basis coefficients. The resulting framework provides enhanced capability to model complex seasonality, which varies dynamically from year-to-year, and offers improved multi-month ahead point forecasts and substantially tighter forecast intervals (with correct coverage) compared to existing forecasting models. Importantly, the fully Bayesian approach provides well-calibrated and precise uncertainty quantification for epi-relevent features, such as the future value and time of the peak measles count in a given year. An R package is available online.

Figures

Tables

Table 1. Table 1 : Median interval width and percent coverage for 95% forecasting intervals of peak time and peak value. The results are given for measles forecasting of weeks 10-52 and the simulated data with r = 1000 𝑟 1000 r=1000 and r = 10 𝑟 10 r=10 . The narrowest intervals achieving 95% coverage are in bold. The proposed integer-valued functional data models *NB-FTS and *Pois-FTS clearly outperform the Gaussian functional data model. *NB-FTS typically achieves the correct coverage, although in some cases produces wider intervals. Other methods are omitted, since they do not readily admit posterior predictive intervals for epi-relevant features.

Measles Peak Value (10-52)	*NB-FTS	*Pois-FTS	Gauss-FTS
Median Width	2535	1504	2717
Coverage Probability	0.97	0.84	0.68
Measles Peak Time (10-52)	*NB-FTS	*Pois-FTS	Gauss-FTS
Median Width	11	10	39
Coverage Probability	0.94	1.00	0.90
Simulated Peak Value ($r=1000$)	*NB-FTS	*Pois-FTS	Gauss-FTS
Median Width	50	41	42
Coverage Probability	0.95	0.94	0.88
Simulated Peak Time ($r=1000$)	*NB-FTS	*Pois-FTS	Gauss-FTS
Median Width	16	14	16
Coverage Probability	1.00	1.00	1.00
Simulated Peak Value ($r=10$)	*NB-FTS	*Pois-FTS	Gauss-FTS
Median Width	74	48	49
Coverage Probability	0.96	0.96	0.86
Simulated Peak Time ($r=10$)	*NB-FTS	*Pois-FTS	Gauss-FTS
Median Width	16	14	18
Coverage Probability	1.00	1.00	0.99

Table 2. Table 2 : Median interval width and percent coverage for 95% forecasting intervals. The results are separated by weeks forecasted (10-52 or 26-52) and pre-vaccine (before 1963) and post-vaccine (after 1964) for the measles data, and by r = 1000 𝑟 1000 r=1000 and r = 10 𝑟 10 r=10 for the simulated data. The narrowest intervals achieving 95% coverage are in bold. Intervals are substantially wider pre-vaccine. Notably, only the proposed method *NB-FTS is capable of simultaneously achieving the nominal 95% coverage while maintaining narrow intervals.

Pre-vaccine (10-52)	*NB-FTS	*Pois-FTS	Gauss-FTS	SARIMA	SARIMA-sqrt	NB-TS	Pois-TS
Median Width	543	444	1117	2947	1751	600	210
Coverage Probability	0.96	0.84	0.89	0.79	0.92	0.80	0.13
Post-vaccine (10-52)	*NB-FTS	*Pois-FTS	Gauss-FTS	SARIMA	SARIMA-sqrt	NB-TS	Pois-TS
Median Width	78	72	319	3759	1019	37	142
Coverage Probability	0.98	0.94	0.89	0.99	1.00	0.69	0.08
Pre-vaccine (26-52)	*NB-FTS	*Pois-FTS	Gauss-FTS	SARIMA	SARIMA-sqrt	NB-TS	Pois-TS
Median Width	244	202	353	2938	1189	729	154
Coverage Probability	0.96	0.87	0.97	1.00	0.99	0.94	0.09
Post-vaccine (26-52)	*NB-FTS	*Pois-FTS	Gauss-FTS	SARIMA	SARIMA-sqrt	NB-TS	Pois-TS
Median Width	32	33	100	3707	867	36	70
Coverage Probability	0.95	0.92	0.95	1.00	1.00	0.82	0.30
Simulated ($r=1000$)	*NB-FTS	*Pois-FTS	Gauss-FTS	SARIMA	SARIMA-sqrt	NB-TS	Pois-TS
Median Width	20	19	22	73	33	49	24
Coverage Probability	0.96	0.93	0.77	0.92	0.92	0.93	0.61
Simulated ($r=10$)	*NB-FTS	*Pois-FTS	Gauss-FTS	SARIMA	SARIMA-sqrt	NB-TS	Pois-TS
Median Width	24	20	26	97	42	42	21
Coverage Probability	0.96	0.91	0.77	0.92	0.90	0.94	0.56

Equations

$$[Z_i(\tau )|r,{\theta }_i(\cdot )]\stackrel{indep}{\sim }\text{mbox}NB\left\{r,\frac{\exp \left[{\theta }_i(\tau )\right]}{r+\exp \left[{\theta }_i(\tau )\right]}\right\},\tau \in \mathcal{T}$$

$$\mathbb{E}[Z_i(\tau )|r,{\theta }_i(\cdot )]$$

$$\text{mbox}Var[Z_i(\tau )|r,{\theta }_i(\cdot )]$$

$${\theta }_i(\tau )$$

$${\mu }_i(\tau )$$

$$\mathbb{E}[Z_i(\tau )|{\mu }_i(\cdot ),{\sigma }_{ϵ},r]=\exp \left[{\mu }_i(\tau )\right]\exp ({\sigma }_e^2/2)$$

$${\beta }_{k,i}={\mu }_k+{\varphi }_k({\beta }_{k,i-1}-{\mu }_k)+{\eta }_{k,i},{\eta }_{k,i}\stackrel{indep}{\sim }N(0,{\sigma }_{{\eta }_{k,i}}^2)$$

$$\text{mbox}Cov\left[{\theta }_i(\tau ),{\theta }_i(u)\right]=\text{mbox}Cov\left[{\vartheta }_{\tau +(i-1)m},{\vartheta }_{u+(i-1)m}\right]={\mathbf{f}}^{\prime }(\tau )\text{mbox}Cov\left({\mathbf{\beta }}_i\right)\mathbf{f}(u)+\mathbb{I}\{\tau =u\}{\sigma }_{ϵ}^2$$

$$\text{mbox}Cov\left[{\theta }_i(\tau ),{\theta }_{i-\ell }(u)\right]=\text{mbox}Cov\left[{\vartheta }_{\tau +(i-1)m},{\vartheta }_{u+(i-\ell -1)m}\right]={\mathbf{f}}^{\prime }(\tau )\text{mbox}Cov\left({\mathbf{\beta }}_i,{\mathbf{\beta }}_{i-\ell }\right)\mathbf{f}(u)$$

$$\text{mbox}Cov\left[{\theta }_i(\tau ),{\theta }_{i-\ell }(u)\right]=\sum _{k=1}^K\left[f_k(\tau )f_k(u){\varphi }_k^{\ell }{\sigma }_{{\eta }_k}^2/\left(1-{\varphi }_k^2\right)\right]$$

$${\beta }_{k,i}={\mu }_k+\sum _{j=1}^p{x}_{i,j}{\alpha }_{j,k}+{\gamma }_{k,i},{\gamma }_{k,i}={\varphi }_k{\gamma }_{k,i-1}+{\eta }_{k,i}$$

$$\mathbb{E}[Z_i(\tau )|{\mu }_i(\cdot ),{\sigma }_{ϵ},r]=\exp \left[\tilde{\mu }(\tau )\right]\left\{\prod _{j=1}^p\exp \left[x_{i,j}{\tilde{\alpha }}_j(\tau )\right]\right\}\exp \left[{\tilde{\gamma }}_i(\tau )\right]\exp ({\sigma }_e^2/2)$$

Full text

KEYWORDS: Bayesian modeling; disease; prediction; MCMC; public health; time series

1 Introduction

Forecasting the spread of infectious diseases is a fundamental goal of epidemiologists and healthcare suppliers. Accurate weeks- and months-ahead forecasts provide essential information for planning and allocation of medical resources and contribute to an informed and healthy community. Precise uncertainty quantification of future disease counts greatly aids medical and public health officials, and can save human lives. Other epi-relevant features, such as the peak number of cases and the time at which the peak occurs, may be equally important to forecast (Tabataba et al., , 2017). Indeed, despite the wide availability of vaccines for many infectious diseases in the United States, certain disease incidences are on the rise in many states, and often accompanied by fatalities (Van Panhuis et al., , 2013).

Among infectious diseases, measles presents an important and imminent challenge. Measles infections commonly occur among school-aged children, with incidence patterns that demonstrate complex seasonality (see Figure 1): in endemic localities, distinct intra-year seasonality persists both during and between epidemics, yet seasonality dissipates as the disease is eradicated (Durrheim et al., , 2014). Accounting for intra-year seasonality—which may evolve or vanish over time—is therefore essential for measles modeling and forecasting. Importantly, measles is the most infectious communicable disease known: an infected individual in a susceptible population generates 12-18 new cases on average (Durrheim et al., , 2014). As a result, vaccination rates must be exceptionally high to prevent measles outbreaks.

The state of Texas is notably problematic: it is second most populous state in the United States, with a rapidly declining measles vaccination rate. Texas state law allows for “reasons of conscience” (i.e., nonmedical) exemptions from immunization. The Texas Department of State Health Services reports that 1.07% of students—or around 50,000 children—were granted exemptions in 2018, far above the 0.45% who were granted exemptions in 2010 (Texas Department of State Health Services, 2018a, ). Lo and Hotez, (2017) carefully articulate the dangers of lapses in measles vaccination coverage in Texas, noting that consequences extend beyond school children and may adversely affect unvaccinated infants and adults, with implications for both public health and the statewide economy. The perilous combination of (i) a highly infectious disease and (ii) an expanding population of susceptible, nonvaccinated children is particularly concerning, and leads experts to warn of a serious and imminent risk for a measles outbreak in Texas in the coming years (Hotez, , 2016).

To improve preparedness for a measles outbreak, we develop new methodology that provides accurate weeks- and months-ahead forecasts with precise uncertainty quantification. Within a Bayesian framework, we model measles counts as integer-valued functional data: measles counts are a function of time-of-year with replicates (time-ordered) across years. This approach, illustrated in Figure 1, allows for flexible modeling of the intra-year seasonality, which varies dynamically from year-to-year—an important feature for modeling measles incidence. Here, the intra-year seasonality pattern is modeled as unknown and learned from the data, which provides the foundation for accurate long-term forecasting coupled with well-calibrated prediction intervals. The measles counts are modeled using a negative-binomial distribution conditional on a real-valued latent process, which offers several advantages: (i) it is well-defined for count data, which is important for modeling small counts that occur during low seasons and especially post-vaccine; (ii) it accounts for overdispersion observed in the data; (iii) it is compatible with an efficient Gibbs sampling algorithm for posterior inference; and (iv) by conditioning on a real-valued process, it provides a general framework for building upon the multitude of successful models for real-valued functional data.

Prediction of disease incidence has become an important statistical research topic in recent years (Unkel et al., , 2012; Nsoesie et al., , 2014). Successful models commonly assume dynamic and integer-valued distributions for disease counts. The susceptible-infected-recovered (SIR) model, which specifies a mechanistic model for the spread of infectious disease via differential equations, corresponds to a negative-binomial model with latent dynamics (Bjørnstad et al., , 2002; Dalziel et al., , 2016). For influenza forecasting, Osthus et al., (2018) improve upon the parametric limitations of the SIR model by incorporating additional local dynamics via a reverse-random walk model. Martínez-Bello et al., (2017) model dengue disease counts using a dynamic Poisson model, which includes covariates. These methods, including more classical seasonal time series models (Shumway and Stoffer, , 2000), focus primarily on local week-to-week time-dependence, with limited flexibility for intra-year seasonalities. Alternatively, Brooks et al., (2015) propose an empirical Bayesian approach for modeling influenza counts as smooth functions of week-of-year. The aggregate success of these methods serves as motivation for the proposed Bayesian integer-valued functional data model.

The paper is organized as follows. The integer-valued functional data model is defined in Section 2, with an accompanying MCMC algorithm for posterior inference in Section 3. The measles forecasting study is in Section 4, which also includes in-sample inference for the vaccine effect. A simulation study is in Section 5. We conclude in Section 6.

2 Modeling Integer-Valued Functional Data

2.1 The Model

Suppose we observe (non-negative) integer-valued functions $Z_{i}\!:\mathcal{T}\to\mathcal{N}=\{0,1,\ldots,\infty\}$ for $i=1,\ldots,n$ on some domain $\mathcal{T}\in\mathbb{R}^{D}$ for $D\in\mathbb{Z}^{+}$. For measles incidence data, $Z_{i}(\tau)$ corresponds to the measles counts in week $\tau$ of year $i$. The intra-year index $\tau\in\mathcal{T}$ is continuous (with $D=1$), so the model is well-defined for finer time scales, such as daily counts. The integer-valued functional data are modeled using a negative-binomial distribution conditional on a real-valued latent process, $\theta_{i}\!:\mathcal{T}\to\mathbb{R}$, and a scalar dispersion, $r>0$:

[TABLE]

The dispersion $r$ provides distributional flexibility, and in particular may incorporate overdispersion of the integer-valued data $Z_{i}(\cdot)$, while the real-valued process $\theta_{i}(\cdot)$ accounts for functional and other dependence. Explicitly, the negative-binomial density of $z\sim\mbox{NB}(r,\pi)$ for $r>0$ and $\pi\in[0,1]$ is $p(z|r,\pi)=\Gamma(r+z)/\left[\Gamma(r)z!\right]\pi^{z}(1-\pi)^{r}$ supported on $z\in\mathcal{N}$. The conditional expectation and variance of $Z_{i}(\tau)$ in (1) are related to $r$ and $\theta_{i}(\cdot)$ as follows:

[TABLE]

The process $\theta_{i}(\cdot)$ completely determines the conditional expectation of the integer-valued functional data $Z_{i}(\cdot)$, while the conditional variance exceeds the conditional expectation with greater overdispersion as $r>0$ decreases. Alternatively, model (1) may be derived via the Poisson-Gamma hierarchical model $[Z_{i}(\tau)|\lambda_{i}(\tau)]\stackrel{{\scriptstyle indep}}{{\sim}}\mbox{Poisson}\left[\lambda_{i}(\tau)\right]$ with $[\lambda_{i}(\tau)|r,\theta_{i}(\cdot)]\sim\mbox{Gamma}\left\{r,r\exp\left[\theta_{i}(\tau)\right]\right\}$ and marginalizing over $\lambda_{i}(\tau)$. In the limit $r\rightarrow\infty$, (1) converges to a Poisson distribution with mean parameter $\exp\left[\theta_{i}(\tau)\right]$.

$Z_{i}(\tau)$ is modeled as conditionally independent in (1), with all functional and other dependence, such as time dependence, relegated to the real-valued process $\theta_{i}(\tau)$. Importantly, the vast majority of Bayesian models for functional data are designed for real-valued processes, which suggests a wide range of options for $\theta_{i}(\cdot)$. We model the real-valued process using a functional basis expansion with Gaussian innovations:

[TABLE]

where $\mu_{i}(\tau)$ is the conditional expectation of $\theta_{i}(\tau)$, $\epsilon_{i}(\cdot)$ is a white noise process, $f_{k}\!:\mathcal{T}\to\mathbb{R}$ for $k=1,\ldots,K$ are (known or unknown) basis functions, and $\{\beta_{k,i}\}$ are the corresponding basis coefficients. Model (4)-(5) is designed to reproduce the classical setting for functional data analysis, where a real-valued function $\theta_{i}(\tau)$ is modeled as noisy observations of a smooth function $\mu_{i}(\tau)$, commonly via a basis expansion. For maximal flexibility, we model the basis functions $\{f_{k}\}$ as smooth yet unknown functions (subject to identifiability constraints), which produces a data-adaptive functional basis (Kowal, , 2018; Kowal and Bourgeois, , 2018). For modeling measles counts, learning $\{f_{k}\}$ corresponds to learning the intra-year seasonalities, with year-specific weights given by the coefficients $\{\beta_{k,i}\}$ (see Section 2.2).

The conditional expectation of $Z_{i}$ in (2) is non-smooth, since $\theta_{i}$ contains a white noise component $\epsilon_{i}$. A smoother representation is obtainable by marginalizing over $\epsilon_{i}$:

[TABLE]

where $\mu_{i}$ is smooth. Equation (4) may also include an additive offset or exposure, say $E_{i}(\tau)>0$, on the log-scale such that $\mathbb{E}[Z_{i}(\tau)|\mu_{i}(\cdot),\sigma_{\epsilon}]=\exp\left[\mu_{i}(\tau)+\log E_{i}(\tau)\right]\exp(\sigma_{e}^{2}/2)=E_{i}(\tau)\exp\left[\mu_{i}(\tau)\right]\exp(\sigma_{e}^{2}/2)$. In Section 4, we include an offset for the population of Texas.

The combination of a negative-binomial model (1) and a Gaussian model (4) has appeared previously in other contexts. In place of the basis expansion (5), Davis and Wu, (2009) propose a linear time series model, Zhou et al., (2012) use a regression model, and Klami, (2015) introduce a factor model, while each emphasizes distributional flexibility, in particular for overdispersion. The model framework of (1), (4), and (5) is compatible with a computationally efficient and convenient blocking structure for a Gibbs sampler, as in Zhou et al., (2012) and Klami, (2015): $\theta_{i}(\cdot)$ is sampled from a full conditional Gaussian distribution using a Pólya-Gamma data augmentation scheme (Polson et al., , 2013), and conditional on $\theta_{i}(\cdot)$, model (4)-(5) is simply a Gaussian functional data model with “data” $\theta_{i}(\cdot)$ and therefore may utilize existing samplers for $\sigma_{\epsilon}$ and the parameters comprising $\mu_{i}(\cdot)$.

Model (5) may be accompanied by various functional data models, such as function-on-scalars regression (Goldsmith et al., , 2015; Goldsmith and Kitago, , 2016), functional multi-level and mixed models (Morris and Carroll, , 2006; Di et al., , 2009; Zhu et al., , 2011), functional autoregressive models for time-ordered functional data (Kowal et al., 2017c, ), among others (e.g., Ramsay and Silverman, , 2005; Morris, , 2015). To incorporate time-dependence from year-to-year, we consider the following autoregressive model:

[TABLE]

Generalizations for multiple lags and vector autoregressions are available, yet less parsimonious. The model properties of (4), (5), and (7) are discussed in Section 2.3.

The priors for $\mu_{k}$ and $\eta_{k,i}$ control the shrinkage behavior of the model for $\beta_{k,i}$. Without adequate shrinkage on $\mu_{k}$ and $\eta_{k,i}$, the model may be overly sensitive to the number of basis functions $K$. To mitigate the impact of the choice of $K$, we introduce ordered shrinkage across coefficients $k=1,\ldots,K$ via a multiplicative gamma process (MGP) prior (Bhattacharya and Dunson, , 2011) for $\mu_{k}$ and $\eta_{k,i}$. The resulting variance for each of $\mu_{k}$ and $\eta_{k,i}$ is stochastically decreasing in $k$, so coefficients become a priori less important for larger $k$. Compared to alternative approaches that treat $K$ as unknown (Suarez and Ghosal, , 2017), the MGP approach is computationally scalable and does not require complex and intensive computing procedures such as reversible jump MCMC. Specifically, the priors are $\mu_{k}\stackrel{{\scriptstyle indep}}{{\sim}}N(0,\sigma_{\mu_{k}}^{2})$ with $\sigma_{\mu_{k}}^{-2}=\prod_{\ell\leq k}\delta_{\mu_{\ell}}$, $\delta_{\mu_{1}}\sim\mbox{Gamma}(a_{\mu_{1}},1)$, $\delta_{\mu_{\ell}}\sim\mbox{Gamma}(a_{\mu_{2}},1)$ for $\ell>1$, and $\sigma_{\eta_{k,i}}^{2}=\sigma_{\eta_{k}}^{2}/\xi_{\eta_{k,i}}$ with $\sigma_{\eta_{k}}^{-2}=\prod_{\ell\leq k}\delta_{\eta_{\ell}}$, $\delta_{\eta_{1}}\sim\mbox{Gamma}(a_{\eta_{1}},1)$, $\delta_{\eta_{\ell}}\sim\mbox{Gamma}(a_{\eta_{2}},1)$ for $\ell>1$, and $\xi_{\eta_{k,i}}\stackrel{{\scriptstyle iid}}{{\sim}}\mbox{Gamma}(\nu_{\eta}/2,\nu_{\eta}/2)$ with $\nu_{\eta}\sim\mbox{Uniform}(2,128)$ to induce heavier tails in the marginal distribution of $\eta_{k,i}$. The hyperpriors $a_{\mu_{1}},a_{\mu_{2}},a_{\eta_{1}},a_{\eta_{2}}\stackrel{{\scriptstyle iid}}{{\sim}}\mbox{Gamma}(2,1)$ allow the data to determine the rate of ordered shrinkage. For the autoregressive coefficient in (7), we assume the prior $\left(\phi_{k}+1\right)/2\stackrel{{\scriptstyle iid}}{{\sim}}\mbox{Beta}(a_{\phi},b_{\phi})$ to constrain $|\phi_{k}|<1$ for stationarity, and select $a_{\phi}=5$ and $b_{\phi}=2$ for measles and simulated data. These choices were successfully applied for Gaussian functional data in Kowal, (2018) and Kowal and Bourgeois, (2018).

2.2 Modeling the Basis Functions

While many options exist for the basis functions $f_{k}(\cdot)$ in (5), we follow the highly flexible approach of Kowal et al., 2017a for unknown $f_{k}$ with the computational improvements in Kowal, (2018) and Kowal and Bourgeois, (2018). Modeling each basis function $f_{k}$ as unknown (i) produces a data-adaptive functional basis and (ii) incorporates the uncertainty about $f_{k}$ into the posterior distribution. Specifically, we let $f_{k}(\tau)=\bm{b}^{\prime}(\tau)\bm{\psi}_{k}$ for $\bm{b}^{\prime}(\tau)=(b_{1}(\tau),\ldots,b_{L_{m}}(\tau))$ a vector of low-rank think plate splines and $\bm{\psi}_{k}$ a vector of unknown coefficients. Low-rank thin plate splines are flexible, smooth, and typically efficient within MCMC samplers (Crainiceanu et al., , 2005), and are well-defined for $\mathcal{T}\subset\mathbb{R}^{D}$ with $D\in\mathbb{Z}^{+}$. For smoothness, we assume the prior $\bm{\psi}_{k}\sim N(\bm{0},\lambda_{f_{k}}^{-1}\bm{\Omega}^{-1})$, where $\bm{\Omega}$ is a $L_{m}\times L_{m}$ known roughness penalty matrix, such as $\bm{\Omega}=\int\ddot{\bm{b}}(\tau)\ddot{\bm{b}}^{\prime}(\tau)d\tau$ for $\ddot{\bm{b}}$ the second derivative of $\bm{b}$, and $\lambda_{f_{k}}>0$ is an unknown smoothing parameter. Details on the construction of $\bm{b}(\cdot)$ and $\bm{\Omega}$, including useful reparametrizations, are given in Kowal, (2018). For identifiability, we enforce the matrix orthonormality constraint $\bm{F}^{\prime}\bm{F}=\bm{I}_{K}$, where $\bm{F}=(\bm{f}_{1},\ldots,\bm{f}_{K})$ is the $m\times K$ basis matrix and $\bm{f}_{k}=(f_{k}(\tau_{1}),\ldots,f_{k}(\tau_{m}))^{\prime}$ is the basis function $f_{k}$ evaluated at the observation points $\tau_{1},\ldots,\tau_{m}$. The matrix orthonormality constraint, coupled with the ordering implied by the MGP prior, is sufficient for identifiability, and may be leveraged to improve computationally efficiency in sampling the coefficients $\{\beta_{k,i}\}$ as in Kowal, (2018) and Kowal and Bourgeois, (2018).

2.3 Properties of the Model

To understand the implications of a functional time series model for seasonal time series data, we consider the covariance properties of the real-valued process $\theta_{i}(\tau)$ implied by (4)-(5). Naturally, the properties of $\theta_{i}(\tau)$ are important for the integer-valued functional data $Z_{i}(\tau)$, and for real-valued functional data we may simply replace $\theta_{i}(\tau)$ with the functional observations. Let $\{\vartheta_{t}\}_{t=1}^{T}$ denote a seasonal time series with seasonality $m$ and $n$ seasons, so $T=nm$. To model the seasonal time series as a functional time series, we map $\theta_{i}(\tau)=\vartheta_{\tau+(i-1)m}$ for each $\tau=1,\ldots,m$. In our application, the seasonality is $m=52$ weeks and the seasons are years $i=1,\ldots,n$. Consequently, the covariance of the process $\theta_{i}(\tau)$ directly determines the seasonal covariance of $\vartheta_{t}$.

Let $\bm{f}^{\prime}(\tau)=(f_{1}(\tau),\ldots,f_{K}(\tau))$ and $\bm{\beta}_{i}^{\prime}=(\beta_{1,i},\ldots,\beta_{K,i})$.

Proposition 1**.**

The contemporaneous covariance function of $\theta_{i}(\tau)$ in (4), conditional on $\{f_{k}\}$, is full rank with

[TABLE]

The covariance function in (8) is the covariance between weeks $\tau$ and $u$ of the same year $i$. If the coefficients $\{\bm{\beta}_{i}\}$ are stationary, then the covariance of $\bm{\beta}_{i}$ is constant and the within-year covariance of $\vartheta_{t}$ does not vary from year-to-year. The week-to-week dependence is primarily governed by the basis functions $\{f_{k}\}$, which are smooth, implying that the seasonal covariance in (8) is also smooth. For weeks in different years, we have the following:

Proposition 2**.**

The lag-$\ell$ autocovariance function of $\theta_{i}(\tau)$ in (4), conditional on $\{f_{k}\}$, is rank $K$ with

[TABLE]

For a seasonal time series, (9) determines the covariance of $\vartheta_{t}$ between different seasons of different years. By modeling the dynamics of $\bm{\beta}_{i}$ (with respect to $i$), we may incorporate a flexible model for time-varying seasonality. In the special case of (7), the covariance in (9) simplifies to

[TABLE]

Seasonality is determined by the basis functions $f_{k}$, which are learned from the data, while year-to-year dependence is controlled by the autoregressive coefficients $\{\phi_{k}\}$. For the same week $\tau$ in different years $i$ and $i-\ell$, (10) simplifies to $\mbox{Cov}\left[\theta_{i}(\tau),\theta_{i-\ell}(\tau)\right]=\sum_{k=1}^{K}\left[f_{k}^{2}(\tau)\phi_{k}^{\ell}\sigma_{\eta_{k}}^{2}/\left(1-\phi_{k}^{2}\right)\right]$. Naturally, the year-over-year covariance depends on the time-of-year via $f_{k}^{2}(\tau)$: for example, the seasonal behavior during the peak months in the late spring may differ from the seasonal behavior during low seasons in the fall (see Figure 1).

3 MCMC Sampling Algorithm

We develop a computationally efficient MCMC sampling algorithm for model (1), (4), and (5). A fundamental observation is that, conditional on $\theta_{i}(\tau)$, the remaining parameters in $\mu_{i}(\tau)$ and $\sigma_{\epsilon}$ may be sampled using MCMC methods for Gaussian functional data, such as Morris and Carroll, (2006), Di et al., (2009), Zhu et al., (2011), Kowal, (2018), and Kowal and Bourgeois, (2018). Therefore, the primary consideration is obtaining a sampler for the negative-binomial parameters $r$ and $\theta_{i}(\tau)$ in (1). For the dispersion parameter $r>0$, we assume the half-Cauchy prior distribution $r\sim C^{+}(0,10)$ and sample from the full conditional distribution using the univariate slice sampler (Neal, , 2003).

To draw from the full conditional distribution of $\{\theta_{i}(\cdot)\}$, we employ a Pólya–Gamma data augmentation scheme (Polson et al., , 2013), which produces a computationally efficient sampling algorithm that does not require any tuning. The likelihood for $\theta_{i}(\tau)$ is proportional to $p(Z_{i}(\tau)|r,\theta_{i}(\cdot))\propto\left\{\exp\left[\theta_{i}(\tau)-\log r\right]\right\}^{Z_{i}(\tau)}/\left\{1+\exp\left[\theta_{i}(\tau)-\log r\right]\right\}^{r+Z_{i}(\tau)}$, which is proportional to a $Z$-distribution for $\left[\theta_{i}(\tau)-r\right]$. The likelihood for $\theta_{i}(\tau)$ is therefore proportional to the likelihood implied by $Z_{i}^{PG}(\tau)=\theta_{i}(\tau)+\nu_{i}(\tau)$ with $Z_{i}^{PG}(\tau)\equiv\left[Z_{i}(\tau)-r\right]/\left[2\xi_{i}(\tau)\right]+\log r$ and $[\nu_{i}(\tau)|\xi_{i}]\stackrel{{\scriptstyle indep}}{{\sim}}N(0,\xi_{i}^{-1}(\tau))$ for $\xi_{i}(\tau)\stackrel{{\scriptstyle iid}}{{\sim}}PG(Z_{i}(\tau)+r,0)$ a Pólya-Gamma random variable (Kowal et al., 2017b, ). Sampling proceeds by drawing $\theta_{i}(\cdot)$ from its Gaussian full conditional distribution and the auxiliary Pólya-Gamma random variables from the full conditional distribution $\xi_{i}(\tau)\stackrel{{\scriptstyle iid}}{{\sim}}PG\left[Z_{i}(\tau)+r,\theta_{i}(\tau)-\log r\right]$. Additional background on Pólya-Gamma augmentation is available in Zhou et al., (2012), Polson et al., (2013), and Klami, (2015).

Let $\tau_{1},\ldots,\tau_{m}$ denote the observation points of $Z_{i}(\cdot)$. While these points are assumed to be common for all $i=1,\ldots,n$, sparsely- or irregularly-sampled functional data may be accommodated via an imputation step. An outline of the Gibbs sampling algorithm is below:

Imputation: Sample $\left[Z_{i}(\tau^{*})|\cdots\right]\stackrel{{\scriptstyle indep}}{{\sim}}\mbox{NB}\left(r,\exp\left[\theta_{i}(\tau^{*})\right]/\left\{r+\exp\left[\theta_{i}(\tau^{*})\right]\right\}\right)$ for all unobserved $Z_{i}(\tau^{*})$; 2. 2.

Dispersion: Sample from $[r|\{Z_{i}(\cdot),\theta_{i}(\cdot)\}]$ using a slice sampler (Neal, , 2003). 3. 3.

Parameter Expansion: Sample $[\xi_{i}(\tau_{j})|\cdots]\stackrel{{\scriptstyle indep}}{{\sim}}PG\left[Z_{i}(\tau_{j})+r,\theta_{i}(\tau_{j})-\log r\right]$ using Polson et al., (2013); 4. 4.

Process: Sample $[\theta_{i}(\tau_{j})|\cdots]\stackrel{{\scriptstyle indep}}{{\sim}}N\left(Q_{\theta_{j,i}}^{-1}\ell_{\theta_{j,i}},Q_{\theta_{j,i}}^{-1}\right)$ where $Q_{\theta_{j,i}}=\xi_{i}(\tau_{j})+\sigma_{\epsilon}^{-2}$ and $\ell_{\theta_{j,i}}=\xi_{i}(\tau_{j})Z_{i}^{PG}(\tau_{j})+\sigma_{\epsilon}^{-2}\mu_{i}(\tau_{j})$ with $Z_{i}^{PG}(\tau_{j})\equiv\left[Z_{i}(\tau)-r\right]/\left[2\xi_{i}(\tau)\right]+\log r$; 5. 5.

Rest: given $\{\theta_{i}(\tau_{j})\}$, sample the parameters of the Gaussian functional data model implied by (5), including $\sigma_{\epsilon}$, $\{\beta_{k,i}\}$, and $\{f_{k}\}$.

For model (7), the final Gibbs block above uses the efficient sampler in Kowal, (2018). While we omit details for brevity, the Kowal, (2018) sampler iteratively draws (i) the basis functions $\{f_{k}\}$ using a Bayesian backfitting sampler for Bayesian splines (Hastie and Tibshirani, , 2000), (ii) the (dynamic) basis coefficients $\{\beta_{k,i}\}$ using a state space sampler (Durbin and Koopman, , 2002), and (iii) the variance components $\sigma_{\epsilon}$ and the MGP shrinkage parameters with known full conditional distributions. An appealing feature of the proposed Gibbs sampler is its modularity: substituting other models in (5) only impacts the final Gibbs block above, which may directly incorporate existing samplers for Gaussian functional data models.

Draws from the posterior predictive distribution, say $\tilde{Z}_{i}(\tau^{*})$, are generated using step 1. Inference on epi-relevant features proceeds by evaluating $h(\tilde{Z}_{i})$ for each draw from the posterior predictive distribution, where $h$ computes the epi-relevant feature, such as the maximum count over $\tau^{*}$ in year $i$.

4 Measles Counts in Texas

We model and forecast weekly measles counts in Texas using publicly-available data from Project Tycho (Van Panhuis et al., , 2013). Data are available from 1928-2002 with approximately 10% of counts missing, and exhibit substantial overdispersion: the sample mean and sample variance of the observed counts are 355 and 628,746, respectively. Figure 1 plots the observed counts, and demonstrates how we construct a functional time series using seasonal time series: the counts are re-organized as functions of week-of-year and time-ordered by year. The annual seasonality is evident, with maximal variability between March and June. The amplitude varies from year-to-year, and vanishes altogether shortly following the introduction of the vaccine in 1963 (Plotkin, , 2014). In more recent years, measles counts are comparatively low: from 2006-2016 there were 53 total cases, 27 of which occurred in 2013 alone (Texas Department of State Health Services, 2018b, ). However, as noted previously, declining vaccination rates in Texas imply a rapidly growing susceptible population, which continuously increases the risk of an outbreak.

4.1 Measles Counts: Forecasting

For optimal preparedness and resource allocation, multi-weeks–ahead forecasts with uncertainty quantification are essential. However, adequate forecasting performance on post-vaccine data alone may be insufficient: in the event of a measles outbreak, counts may reflect pre-vaccine seasonality patterns and volume. Therefore, it is imperative that forecasting technologies be capable of performing adequately both pre- and post-vaccine.

Let $Z_{i}(\tau)$ denote the measles incident count for week $\tau=1,\ldots,m=52$ within year $i=1,\ldots,n$. Equivalently, we may write the observed measles counts more traditionally as a time series of counts $z_{t},t=1,\ldots,T=nm$ with $Z_{i}(\tau)=z_{\tau+(i-1)m}$. To accompany the measles counts, we include annual population data for Texas (U.S. Bureau of the Census, Federal Reserve Bank of St. Louis, , 2018) as an offset $E_{i}(\tau)=E_{i}$ in (4)-(5), so (4) becomes $\theta_{i}(\tau)=\log E_{i}(\tau)+\mu_{i}(\tau)+\epsilon_{i}(\tau)$, which implies the conditional expectation (6) may be expressed as a rate, $\mathbb{E}[Z_{i}(\tau)/E_{i}(\tau)|\mu_{i}(\cdot),\sigma_{\epsilon},r]=\exp\left[\mu_{i}(\tau)\right]\exp(\sigma_{e}^{2}/2)$.

4.1.1 Forecasting Design

For each year from 1950-1980, we compute multi-weeks-ahead forecasts given (i) historical measles counts and population data from 1928 up to the forecasting year and (ii) measles counts and population data from the first $m_{0}$ weeks of the forecasting year. By observing a small number $m_{0}$ of weeks at the beginning of the year, the forecasting models may adapt to each year’s distinct seasonality pattern in order to forecast the remaining $m-m_{0}$ weeks of the year. We consider $m_{0}=9$ weeks and $m_{0}=25$ weeks: the former case $m_{0}=9$ is more challenging, since it constrains the partial observations almost entirely to January-February and therefore requires forecasts during the peak months of March-June as well as the subsequent decline. The forecasting period 1950-1980 includes pre- and post-vaccine years, and terminates in 1980 due to substantial missingness thereafter (see Web Figure 12).

An illustration of the forecasting design is given in Figure 2 for 1961, which immediately precedes the introduction of the measles vaccine. Given observations from the first $m_{0}=9$ weeks (yellow region), we use the posterior predictive distribution under model (1), (4), (5), and (7) to compute (i) the posterior predictive expected value (6) for the future counts in weeks 10-52 (blue line), (ii) 95% posterior predictive intervals for the future counts in weeks 10-52 (gray region), and (iii) the posterior predictive distribution for the future peak measles count and the time at which it occurs. This figure appears in color in the electronic version of this article, and color refers to that version. The measles forecasts are accurate, with prediction intervals that widen during peak times from March-July, narrow for July-November, and widen again in December, which matches the observed pattern of variability in Figure 1. For both peak value and peak time, the corresponding posterior predictive distribution is centered around the future observed value. An analogous plot given observations from the first $m_{0}=25$ weeks in 1961 is in Web Figure 9: the forecasting intervals are far more narrow, yet the posterior distribution of peak time has an interesting bimodality, which assigns large probability to the peak time occurring in either July or December.

4.1.2 Forecasting Methods

We consider a variety of functional, seasonal, and time series forecasting methods with different distributional assumptions. For each forecasting year $t$, we are interested in computing point and interval forecasts for an unobserved (future) count $Z_{t}(\tau^{*})$ for weeks $\tau^{*}>m_{0}$ with known offset $E_{t}(\tau^{*})$ based on (i) historical data from previous years $\{Z_{i}(\cdot),E_{i}(\cdot):i<t\}$ and (ii) partial observations from the current year $\{Z_{t}(\tau):\tau\leq m_{0}\}$. While accurate point forecasts are undoubtedly important, well-calibrated and precise interval forecasts provide useful information for planning and allocation of resources. In addition to forecasting intervals for future measles counts, we consider 95% prediction intervals for the peak count, $\max\{Z_{t}(\tau^{*}):\tau^{*}>m_{0}\}$, and peak time, $\arg\max_{\tau^{*}}\{Z_{t}(\tau^{*}):\tau^{*}>m_{0}\}$, for each year $t$ from 1950-1980. Note that these quantities are not readily available for all methods.

We consider two variations of the proposed integer-valued functional time series model: a negative binomial distribution with unknown $r$ (NB-FTS) and an approximate Poisson distribution (Pois-FTS) with $r=1000$. Forecasts are computed by treating $Z_{t}(\tau^{*})$ as missing data for $\tau^{*}>m_{0}$ and imputing from the posterior predictive distribution (see Section 3). We select $K=6$ and run the MCMC for 30,000 iterations, discarding the first 5,000 iterations as a burn-in and retaining every 5th simulation to reduce autocorrelation.

As a direct competitor, we include the Gaussian functional time series model (Gauss-FTS) of Kowal, (2018), which uses the same model for (5) and (7), but replaces $\theta_{i}(\tau)$ in (4) with the observed functional data. Although not well-defined for integer-valued data, Gauss-FTS uses the variance-stabilizing transformation for the Poisson distribution: the input data are $Y_{i}(\tau)=\sqrt{Z_{i}(\tau)/E_{i}(\tau)}$ and the posterior predictive distribution of $E_{t}(\tau^{*})\left[Y_{t}(\tau^{*})\right]^{2}$ is computed for forecasting and inference. We use the same choice of $K$ and MCMC specifications as for NB-FTS and Pois-FTS. Among other functional data methods, we compute forecasts based on (i) the pointwise sample means of the functional observations rescaled by the offsets (Mean-FDA), $\hat{Z}_{t}(\tau^{*})=E_{t}(\tau^{*})\left[\frac{1}{t-1}\sum_{i<t}Z_{i}(\tau^{*})/E_{i}(\tau^{*})\right]$, and (ii) the previous functional observation rescaled by the offsets (RW-FDA), $\hat{Z}_{t}(\tau^{*})=E_{t}(\tau^{*})\left[Z_{t-1}(\tau^{*})/E_{t-1}(\tau^{*})\right]$, i.e., a random walk forecast. While these are nominally functional data estimators, they are also seasonal time series estimators, and therefore provide an important baseline.

We also include more classical seasonal and integer-valued time series models, which are broadly popular for disease forecasting (Shumway and Stoffer, , 2000; Martinez et al., , 2011). Using the time series data $z_{1},\ldots,z_{m_{0}+(t-1)m}$ up to week $m_{0}$ of the current year $t$, we fit a seasonal autoregressive moving average (SARIMA) model using the auto.arima() function in the forecast package in R (Hyndman et al., , 2018; Hyndman and Khandakar, , 2008) with the order of the autoregressive, moving average, seasonal autoregression, and seasonal moving average components (based on $m=52$ week seasonality) selected by AIC. Similarly, we fit a SARIMA model to the square-root transformed count data (sqrt-SARIMA) to stabilize the variance. In both cases, the counts are scaled by the offset before fitting and rescaled for forecasting and inference. Among time series models for count data, we include negative-binomial (NB-TS) and Poisson (Pois-TS) models implemented using the tsglm() function in the tscount package in R (Liboschik et al., 2017b, ; Liboschik et al., 2017a, ). We select a log-link function and include past observations and past means from lags 1,2, 52, and 53, which incorporates local time-dependence and annual seasonality.

4.1.3 Forecasting Evaluation Metrics

Forecasting performance is evaluated based on (i) point forecasts, (ii) interval forecasts, and (iii) interval forecasts of epi-relevant features, namely the peak count and the peak time. In each case, results are divided into pre- and post-vaccine eras, with 14 years pre-vaccine and 16 years post-vaccine (we omit 1964 due to substantial missingness). For each year, we forecast $m-m_{0}$ weeks of measles counts, totaling 602 weeks pre-vaccine and 688 weeks post-vaccine for $m_{0}=9$, and 378 weeks pre-vaccine and 432 weeks post-vaccine for $m_{0}=25$.

For a point forecast, say $\hat{Z}_{i}(\tau^{*})$, we consider two definitions of mean absolute error (MAE). We compute the MAE averaged across all weeks $\tau^{*}>m_{0}$ in each year $i$, $\mbox{MAE}_{i}=(m-m_{0})^{-1}\sum_{\tau^{*}=m_{0}+1}^{m}|Z_{i}(\tau^{*})-\hat{Z}_{i}(\tau^{*})|$, which illustrates how forecasting performance varies from year-to-year. To compare forecasting performance for different forecasting horizons, we also compute the MAE averaged across all years in each week $\tau^{*}$, $\mbox{MAE}(\tau^{*})=(n-n_{0})^{-1}\sum_{i=n_{0}+1}^{n}|Z_{i}(\tau^{*})-\hat{Z}_{i}(\tau^{*})|$ where $n_{0}$ is the number of years used for estimation and $n$ is the total number of years.

Interval forecasts are evaluated using empirical coverage probability (ECP) and median interval width (MIW): for a prediction interval $(\hat{Z}_{i}^{lower}(\tau^{*}),\hat{Z}_{i}^{upper}(\tau^{*}))$, \mbox{ECP}=(n-n_{0})^{-1}(m-m_{0})^{-1}\sum_{i=n_{0}+1}^{n}\sum_{\tau^{*}=m_{0}+1}^{m}\mathbb{I}\Big{\{}\hat{Z}_{i}^{lower}(\tau^{*})\leq Z_{i}(\tau^{*})\leq\hat{Z}_{i}^{upper}(\tau^{*})\Big{\}} and $\mbox{MIW}=\mbox{median}\{\hat{Z}_{i}^{upper}(\tau^{*})-\hat{Z}_{i}^{lower}(\tau^{*}):i=n_{0}+1,\ldots,n,\tau^{*}=m_{0}+1,\ldots,m\}$. Both metrics are important: ECP assesses calibration of the interval, while MIW describes the precision of the uncertainty quantification. The goal is to obtain the most narrow forecasting intervals that achieve the 95% nominal coverage.

4.1.4 Forecasting Results

A comparison of the point forecasts is given in Figure 3, which includes both MAE metrics for forecasting weeks 10-52. The proposed methods NB-FTS and Pois-FTS perform the best across all years, with notable improvements over competing methods in the pre-vaccine era. The specific time-of-year improvements are illustrated in the $\mbox{MAE}(\tau^{*})$ plots in the bottom panels: the NB-FTS and Pois-FTS forecasts are substantially more accurate during weeks 10-20, which corresponds to peak measles counts, especially pre-vaccine. Unsurprisingly, Mean-FDA performs worst post-vaccine, since the forecasts are based on pre- and post-vaccine averages. The results for weeks 26-52 are in Web Figure 10: the proposed methods NB-FTS and NB-Poisson remain superior, albeit by smaller margins.

The most striking results are the forecasting interval comparisons for weeks 10-52 (Figure 4), weeks 26-52 (Web Figure 11), as well as simulated data from Section 5 (Web Figure 8), with exact ECP and MIW values given in Web Table 2. Notably, only the proposed NB-FTS model provides forecasting intervals that consistently achieve the 95% nominal coverage. Naturally, correct coverage is essential for utility and interpretability of the intervals. In addition to correct coverage, NB-FTS produces substantially more narrow forecasting intervals, especially among methods that achieve nearly the nominal 95% coverage, which implies greater precision in the forecasting intervals. For example, in the pre-vaccine era, NB-FTS forecast intervals achieve 96% coverage with a median width of 543 counts, while the most competitive alternative, sqrt-SARIMA, has 92% coverage with a median width of 1751 counts—which is more than three times larger. Pois-FTS is competitive with NB-FTS, but typically suffers from undercoverage. Interestingly, the integer-valued distribution appears to be important: Gauss-FTS, despite including a variance-stabilizing transformation, either fails to provide adequate coverage or produces much wider forecast intervals than NB-FTS. Clearly, NB-FTS provides superior forecast intervals, which are both correctly calibrated at 95% and more precise than competing methods.

Similarly, Table 1 compares the 95% forecasting intervals for the epi-relevant features corresponding to the peak time and peak values of the measles counts over the weeks 10-52, as well as for the simulated data in Section 5. For resource allocation and planning, it is important to know when the measles season will peak, and how many people will be affected at the peak. For the Bayesian (integer-valued and Gaussian) functional data models, inference for these quantities are readily available via the posterior predictive distribution (as in Figure 2). Impressively, the proposed integer-valued functional data models provide the correct nominal coverage with much narrower intervals than the Gaussian model.

These results cumulatively demonstrate the clear advantages of using the proposed integer-valued functional data model for point, interval, and epi-relevant forecasts of measles counts. The integer-valued distribution offers substantial improvements for NB-FTS and Pois-FTS relative to Gauss-FTS, including both point estimation and interval coverage, while the functional data approach outperforms seasonal and integer-valued time series models, especially among forecast intervals and pre-vaccine point forecasts.

4.2 Measles Counts: Inference

The proposed modeling framework may be used for in-sample estimation, inference, and imputation of missing values. Consider the introduction of the measles vaccine in 1963 in Figure 1: while there is a notable change in the measles incidence pattern, the decline in measles counts is non-monotone due to the presence of seasonality. A functional time series approach offers an intuitive decomposition: the curves $f_{k}(\tau)$ model intra-year seasonality, while the basis coefficients $\beta_{k,i}$ capture how the level and seasonality vary from year-to-year. We propose the following extension of (7) to include predictors $\{x_{i,j}\}_{j=1}^{p}$ in year $i$, such as yearly time trends and a vaccination effect:

[TABLE]

with $\eta_{k,i}$ distributed as in (7), which is a multiple linear regression with autoregressive errors for each basis coefficient $k=1,\ldots,K$. Model (11) is utilized in Kowal, (2018) for Gaussian function-on-scalars regression, but is easily adaptable to the integer-valued setting via the MCMC algorithm in Section 3. For interpretability, we may rewrite the conditional mean of the count observations (6) as follows:

[TABLE]

where $\tilde{\mu}(\tau)=\sum_{k}f_{k}(\tau)\mu_{k}$ is an intercept function, $\tilde{\alpha}_{j}(\tau)=\sum_{k}f_{k}(\tau)\alpha_{j,k}$ is the regression function for predictor $j$, and $\tilde{\gamma}_{i}(\tau)=\sum_{k}f_{k}(\tau)\gamma_{k,i}$ is the year-specific autoregressive random effect. The representation in (12) is analogous to standard log-linear models, in particular for Poisson and negative-binomial distributions with log-link functions. The seasonal effect of predictor $j$ is capture by $\tilde{\alpha}_{j}(\tau)$; when $\tilde{\alpha}_{j}(\tau)=0$ for all weeks $\tau$, the effect is null. Following Kowal, (2018), we assume nested horseshoe priors on the regression coefficients: $\alpha_{j,k}\stackrel{{\scriptstyle indep}}{{\sim}}N(0,\sigma_{\alpha_{j,k}}^{2})$ with $\sigma_{\alpha_{j,k}}\stackrel{{\scriptstyle indep}}{{\sim}}C^{+}(0,\lambda_{j})$, $\lambda_{j}\stackrel{{\scriptstyle iid}}{{\sim}}C^{+}(0,\lambda_{0})$, and $\lambda_{0}\sim C^{+}(0,1/\sqrt{p})$, which provides a hierarchy of local $(j,k)$ shrinkage, predictor-specific $j$ shrinkage, and global shrinkage.

We implement model (1), (4), (5), and (11) for the entire measles counts dataset from 1928-2002. For predictors, we include a linear time trend $i$, an indicator of post-vaccine years $\mathbb{I}\{i\geq 1963\}$, and a vaccine-year interaction $i\mathbb{I}\{i\geq 1963\}$. We select $K=6$ and run the MCMC for 30,000 iterations, discarding the first 5,000 iterations as a burn-in and retaining every 5th simulation. Traceplots demonstrate good mixing and suggest convergence, and effective samples sizes are sufficiently large. Missing values are automatically imputed within the Gibbs sampler (see Web Figure 12). The model clearly recognizes overdispersion: the 95% highest posterior density (HPD) interval for the dispersion parameter $r$ is $[4.72,6.92]$. In addition, the year-to-year dynamics—even after adjusting for predictors—are substantial: posterior credible intervals for autoregressive coefficients $\phi_{k}$ in (11) exclude zero for $k=1$ (95% HPD interval $[0.45,0.89]$) and $k=2$ (95% HPD interval $[-0.527,-0.003]$).

Inference for the regression coefficients $\tilde{\alpha}_{j}(\tau)$ is in Figure 5, including posterior means and 95% posterior pointwise intervals and simultaneous bands. There is a slightly positive pre-vaccine linear year-to-year trend, particularly in the fall. More importantly, there is a substantial linear year-to-year decline in measles counts post-vaccine, with the largest effect during the peak months of March-June. These effects are visually confirmed in Figure 6, which displays the learned seasonalities via $f_{k}$ and the year-to-year weights $\beta_{k,i}$. The dynamic coefficients $\beta_{k,i}$ show a slight linear trend pre-vaccine and a substantial trend in the oppose direction post-vaccine, especially for $k=1$. Interestingly, the learned basis functions $f_{1}$ and $f_{2}$ are simple: $f_{1}$ is sinusoidal and $f_{2}$ is linear. However, $f_{k}$ for $k>2$ are more complex, with $f_{3}$ partly capturing the peak months. Note that the coefficients $\{\beta_{k,i}\}$ for $k\leq 4$ are visibly distinct from zero, suggesting that the more complex $f_{k}$ curves are important.

5 Simulations

We conducted a simulation study to validate the results from the forecasting comparison and evaluate the model performance under different distributions. We simulate data according to equations (1), (4), (5), and (7) for $n=50$ integer-valued functions observed at $m=50$ points, and introduce missingness for 10% of the observations. For the dispersion parameter in (1), we consider $r=1000$, which approximates the Poisson distribution, and $r=10$, which provides substantial overdispersion similar to the measles data.

For equally-spaced observation points $\tau\in[0,1]$, we introduce functional (or seasonal) dependence using the basis $f_{1}^{*}(\tau)=1/\sqrt{m}$ and $f_{k}^{*}$ an orthogonal polynomial of degree $k$ for $k=2,3,4$. For time-dependence, let $\beta_{k,i}^{*}=\sqrt{m}+\frac{\sqrt{m}}{k}\gamma_{k,i}^{*}$ with $\gamma_{k,i}^{*}=0.8\gamma_{k,i-1}^{*}+\eta_{k,i}^{*}$ and $\eta_{k,i}^{*}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,\sqrt{1-0.8^{2}})$ for $i=1,\ldots,n=50$. The real-valued process is simulated as $\theta_{i}^{*}(\tau)=\mu_{i}^{*}(\tau)+\sigma^{*}\epsilon_{i}^{*}(\tau)$, where $\mu_{i}^{*}(\tau)=\sum_{k=1}^{4}f_{k}^{*}(\tau)\beta_{k,i}^{*}$, $\sigma^{*}=\mbox{sd}[\mu_{i}^{*}(\tau)]/\mbox{RSNR}$ for root signal-to-noise ratio RSNR = 10, and $\epsilon_{i}^{*}(\tau)\stackrel{{\scriptstyle iid}}{{\sim}}N(0,1)$. The observed integer-valued functional data are simulated from $Z_{i}(\tau)\stackrel{{\scriptstyle indep}}{{\sim}}\mbox{NB}\left(r,\exp\left[\theta_{i}^{*}(\tau)\right]/\left\{r+\exp\left[\theta_{i}^{*}(\tau)\right]\right\}\right)$, and 10% of the observations $Z_{i}(\tau)$ are deleted at random to induce missingness. We define the true curves to be the conditional expectations $Z_{i}^{*}(\tau)=\exp\left[\mu_{i}^{*}(\tau)\right]\exp\left[(\sigma^{*})^{2}/2\right]$ akin to (6). Given all observed counts $\{Z_{i}(\tau\}_{i=1}^{n-1}$ prior to time $n$ and the first 30 observation points at time $n$ $\{Z_{n}(\tau_{j})\}_{j=1}^{30}$, the goal is to forecast the true curves $\{Z_{n}^{*}(\tau_{j})\}_{j=31}^{50}$ for the remaining 20 observation points. The same competing methods are considered as in Section 4.1 and the simulations are repeated 100 times. Note that empirical coverage probability (ECP) is evaluated based on future counts $\{Z_{n}(\tau_{j})\}_{j=31}^{50}$ rather than the true curves $\{Z_{n}^{*}(\tau_{j})\}_{j=31}^{50}$.

An illustration of the forecasting design, analogous to Figure 2, is provided in Web Figure 13 for $r=1000$ and Web Figure 14 for $r=10$. The point and interval forecasts perform well in both cases, although the $r=10$ case is far more challenging due to the increased variability. Aggregating across all simulations, the MAEs defined in Section 4.1 are displayed in Figure 7. The functional models NB-FTS, Pois-FTS, and Gauss-FTS perform similarly for $\mbox{MAE}_{i}$, while the proposed integer-valued methods NB-FTS and Pois-FTS offer clear improvements in $\mbox{MAE}(\tau^{*})$ across all forecast horizons. Figure 8 provides ECPs and MIWs for comparing forecast intervals. As for the measles data, NB-FTS achieves the correct nominal coverage with narrower intervals than competing methods. Notably, Gauss-FTS produces similar interval widths, yet suffers from substantial undercoverage. Lastly, coverage for the peak count and peak time are in Table 1, with NB-FTS and Pois-FTS again outperforming Gauss-FTS. Clearly, the integer-valued distribution is important for well-calibrated and precise uncertainty quantification.

6 Discussion

The proposed framework for integer-valued functional data offers clear improvements in modeling and forecasting measles counts relative to existing methods. The methodology is fully Bayesian with an integer-valued distribution for the observed data, which produces correctly calibrated and accurate uncertainty quantification for measles counts and epi-relevant features, including the peak number of counts and the time at which the peak occurs. The intra-year seasonality is learned from the data, and includes year-to-year dynamic dependence. Importantly, the model is sufficiently flexible and accurate in both pre- and post-vaccine eras, and therefore offers utility in the modern era in which vaccines are available yet vaccination rates are declining rapidly. More broadly, the model and accompanying Gibbs sampler are designed to build upon methods and algorithms for Gaussian functional data, such as functional autoregressions and function-on-scalar regressions, which suggests a widely generalizable modeling framework.

Given the success of the proposed methodology, there are a number of promising extensions. First, we may consider disease counts for multiple related diseases or multiple states concurrently, which requires a multivariate functional data model (Kowal et al., 2017a, ). Similarly, external data sources such as Google search trends, Twitter, or Wikipedia have proven useful for influenza forecasting (Paul et al., , 2014; Hickmann et al., , 2015), and may be incorporated for measles forecasting via the regression model (11) or suitable modifications. Lastly, methodological extensions for binomial functional data are readily available via a similar Pólya-Gamma data augmentation scheme.

Supporting Information

Additional supporting information may be found online in the Supporting Information section at the end of the article, including Web Figures, Tables, and R code.

Web Appendix A