Markov-Switching Mixed-Effects Hurdle Models

• Markov-switching mixed-effects hurdle models are a unified framework that separate latent disease-on/off regimes from structural zeros in overdispersed count data.
• The model integrates logistic regression for state transitions with autoregressive and endemic seasonal components to account for spatial and temporal heterogeneity.
• Bayesian MCMC and EM algorithm implementations provide robust estimation and superior predictive performance compared to traditional zero-inflated models.

Markov-switching mixed-effects hurdle models constitute a unified statistical framework designed to jointly model state-switching dynamics, overdispersed spatio-temporal count data, and the structural zero-generating process. These models are especially suited to epidemiological contexts, where disease presence and absence, subject to spatial and temporal covariates and overdispersion, must be distinctly characterized. The Markov-switching negative binomial hurdle model with mixed effects (“ZS-MSNBH model”) provides explicit representation of (i) latent disease-on/-off regimes, (ii) perfect-detection hurdle at zero, (iii) areal random effects, and (iv) autoregressive and endemic (seasonal and spatial) structure (Xu et al., 2023).

1. Model Structure and Key Components

The model posits a latent two-state Markov chain $S_t\in\{0,1\}$ for each spatial unit (or time series) $i$, with transitions governed by covariate-dependent transition probabilities. At each time point $t$, the observed count $Y_t$ is:

$$Y_t|S_t \sim \begin{cases} 0, & S_t=0\ \text{(absence)} \ \mathrm{ZTNB}(\mu_t,\phi), & S_t=1\ \text{(presence)} \end{cases}$$

where $\mathrm{ZTNB}$ denotes the zero-truncated negative binomial distribution. Thus, zeros arise strictly from the latent absence state ($S_t=0$), and positive counts (with all zeros strictly excluded) emanate from the presence state ($S_t=1$), enforcing a “hurdle” at zero.

Transition probabilities are specified through logistic regressions incorporating both covariates and spatial lag terms:

$$\begin{aligned} \text{logit}(p_{01,t}) &= \alpha_0^{(0)} + \mathbf{g}_t^\top\boldsymbol{\alpha}^{(0)} + \gamma_1 \sum_{j\in\text{Nei}(i)} S_{j,t-1} \ \text{logit}(p_{11,t}) &= \alpha_0^{(1)} + \delta^{(1)}\log(y_{t-1}+1) + \mathbf{g}_t^\top\boldsymbol{\alpha}^{(1)} + \gamma_2 \sum_{j\in\text{Nei}(i)} S_{j,t-1} \end{aligned}$$

with $p_{01,t}$ encoding the probability of reemergence and $p_{11,t}$ persistence. Covariates $\mathbf{g}_t$ may include population statistics, environmental measures, or socioeconomic indicators.

2. Hierarchical Regression: Autoregressive, Endemic, and Random Effects

The presence-state ($S_t=1$) count distribution relies on a structured mean formulation reflecting both short-term memory and long-term (seasonal/spatial) patterns:

$\mu_t = \mu_t^{AR} y_{t-1} + \mu_t^{EN}$

• Autoregressive term: $$\mu_t^{AR} = \exp(b_{0,i} + \mathbf{g}_t^\top\boldsymbol{\beta}^{AR})$$, with $b_{0,i} \sim N(\beta_0^{AR},\sigma_{b0}^2)$, enables area-specific short-term dynamics and covariate effects.
• Endemic term: $$\mu_t^{EN} = \exp(b_i + \beta_2^{EN}\sin(2\pi t/52) + \beta_3^{EN}\cos(2\pi t/52))$$, with $$b_i\sim N(\beta_0^{EN}+\beta_1^{EN}\log(\mathrm{pop}_i),\sigma_b^2)$$, accounts for baseline differences (linked to population) and cyclical seasonality.

Random effects $b_{0,i}, b_i$ encode unobserved heterogeneity across spatial units.

Furthermore, overdispersion parameter $\phi_t$ may itself be modeled as a log-linear function of covariates and lagged counts:

$$\log(\phi_t) = \alpha_0^{(2)} + \mathbf{g}_t^\top\boldsymbol{\alpha}^{(2)} + \delta^{(2)}\log(y_{t-1}+1)$$

This flexibility enables robust modeling of heteroskedasticity in the data.

3. Likelihood, Bayesian Specification, and Priors

The joint model likelihood combines (a) the Markov-switching latent process and (b) the hurdle count component:

$$p(S_{1:T},Y_{1:T}|\Theta) = P(S_1) \prod_{t=2}^{T}P(S_t|S_{t-1},\mathbf{g}_t) \prod_{t=1}^{T} p(Y_t|S_t,\mu_t,\phi_t)$$

• $P(S_t|S_{t-1},\mathbf{g}_t)$ as above.
• $p(Y_t|S_t)$ is degenerate at zero for $S_t=0$, and $\mathrm{ZTNB}(\mu_t,\phi_t)$ for $S_t=1$.

Priors are typically weakly-informative Gaussian or uniform, e.g.

$$\boldsymbol{\beta}^{AR},\;\boldsymbol{\alpha}^{(\cdot)},\;\beta_{k}^{EN} \sim N(0,100^2), \quad \sigma_{b0}^2\sim\text{Inv-Gamma}(0.1,0.1),\quad \sigma_b\sim\text{Uniform}(0,10)$$

An uninformative $\text{Bernoulli}(0.5)$ prior initializes the latent state $S_1$.

4. Estimation: Bayesian MCMC and Maximum Likelihood

The model admits both Bayesian and frequentist implementations.

• Bayesian inference: Posterior samples are obtained via Gibbs sampling, with Metropolis–Hastings updates for logistic and over-dispersion parameters, and sequential updates for (1) hidden states $S_{1:T}$ (forward-backward algorithms), (2) random effects, and (3) regression parameters. Convergence is assessed using the Gelman–Rubin statistic ($\hat R<1.05$) and effective sample sizes.
• Maximum likelihood (EM algorithm): The latent switching sequence $S_{1:T}$ is treated as missing data. The E-step uses forward–backward recursion to compute $Q_t(i)=P(S_t=i|Y_{1:T},\Theta^{(m)})$, and the M-step maximizes the expected complete-data log-likelihood. The negative binomial GLM routines and logistic transition regression updates are executed iteratively (Xu et al., 2023).

Implementations are straightforward in probabilistic programming frameworks such as NIMBLE or JAGS.

5. Practical Application and Implementation Guidelines

Key recommendations for application include:

• Centering and scaling continuous covariates to improve MCMC mixing.
• Inclusion/testing of spatial lag terms $\sum_{j\in\text{Nei}(i)}S_{j,t-1}$ to accommodate spatial autocorrelation.
• Monitoring persistence probabilities (posterior mean of $p_{11,t}$) to guard against excessive regime switching.
• Evaluating predictive adequacy via out-of-sample prediction, WAIC, or ranked probability scoring.

A hurdle at zero (in contrast to zero-inflated models) is appropriate when perfect-detection can be assumed, as in diseases with high reporting rates or analyses focusing on true extinction/persistence.

6. Comparison and Theoretical Motivation

The Markov-switching negative binomial hurdle model distinctly separates the mechanisms generating observed zeros (structural/no-disease) from those governing positive counts, further differentiating between persistence and reemergence transitions—each modulated by its own set of covariates and lag structure. This structure allows for (i) explicit modeling of “true absence” periods versus “disease-on” periods, and (ii) the accommodation of spatial, temporal, and areal heterogeneity in both disease presence and case counts.

Compared to zero-inflated alternatives, hurdle models enforce perfect detection in the “on” state, while zero-inflated specifications allow for under-reporting. Analysis of fit on real-world data indicates that both classes of Markov-switching models provide substantial improvements in out-of-sample prediction over non-switching hurdle and zero-inflated benchmarks; zero-inflated Markov-switching models may exhibit superior predictive performance in some cases (Xu et al., 2023).

7. Synthesis and Domain Relevance

Markov-switching mixed-effects hurdle models offer a comprehensive stochastic framework unifying regime-switching, perfect-detection structural zeros, multi-level random effects, autoregressive memory, and endemic-seasonal patterns for count data. Their utility in infectious disease epidemiology derives from their capacity to (a) model extinction, (b) capture heterogeneous reemergence/persistence dynamics, and (c) flexibly incorporate spatial and temporal covariate effects. This approach is particularly advantageous for high-resolution spatio-temporal surveillance and intervention evaluation in settings with strongly patterned incidence and extinction dynamics (Xu et al., 2023).