run_hmm: Hidden Markov Model (HMM) for Path Dependence (Counts I and...

In bivarhr: Bivariate Hurdle Regression with Bayesian Model Averaging

Hidden Markov Model (HMM) for Path Dependence (Counts I and C)

Description

Fits a bivariate count Hidden Markov Model (HMM) to the series I and C with a self-contained, dependency-free base-R implementation of the Baum-Welch (EM) algorithm. Each latent state j emits I and C as conditionally independent counts given the state, with either Poisson or Negative Binomial emissions. With Poisson emissions the model coincides with an intercept-only Poisson depmix; the Negative Binomial family adds a per-state, per-series dispersion parameter and is the appropriate choice for overdispersed counts (the same motivation as the hurdle negative-binomial core of this package).

Usage

run_hmm(
  DT,
  nstates = 3,
  family = c("poisson", "nbinom"),
  seed = NULL,
  n_starts = 10,
  max_iter = 500,
  tol = 1e-08,
  dir_csv = NULL,
  dir_out = NULL
)

Arguments

DT	A data.frame or data.table containing at least the columns I and C, interpreted as non-negative integer count series observed over time.
nstates	Integer; number of latent Markov states (default 3).
family	Character; emission family, either "poisson" (default) or "nbinom" (Negative Binomial, mean–dispersion parametrization). Negative Binomial is recommended when the counts are overdispersed.
seed	Integer or NULL; seed controlling the random EM restarts. The first restart always uses a deterministic quantile-based initialization, so results are reproducible even when seed is NULL. The caller's random-number state is preserved (restored on exit).
n_starts	Integer; number of EM restarts (multi-start) used to mitigate convergence to a poor local optimum. The fit with the highest log-likelihood is returned (default 10).
max_iter	Integer; maximum EM iterations per restart (default 500).
tol	Numeric; relative convergence tolerance on the log-likelihood increase between EM iterations (default 1e-8).
dir_csv	Character scalar or NULL; directory where the state-sequence CSV ("hmm_states.csv") is written. If NULL (default), no CSV is written.
dir_out	Character scalar or NULL; directory where the fitted HMM object ("hmm_fit.rds") is saved. If NULL (default), the object is not saved to disk.

Details

Model. The parameters are an initial-state distribution \delta (length nstates), a row-stochastic transition matrix \Gamma, per-state means \mu^I_j, \mu^C_j and, for family = "nbinom", per-state dispersions \phi^I_j, \phi^C_j. The emission density at time t in state j factorizes as

b_j(t) = f(I_t; \mu^I_j, \phi^I_j)\, f(C_t; \mu^C_j, \phi^C_j),

with f the Poisson or Negative Binomial mass function. As \phi \to \infty the Negative Binomial converges to the Poisson.

Estimation. The Baum-Welch algorithm is used. A vectorized, log-space forward-backward pass computes the smoothed state and transition posteriors (E-step) using the log-sum-exp identity, so the recursion is numerically stable for long series and small emission probabilities. The M-step updates \delta and \Gamma in closed form; the means \mu are the posterior-weighted sample means; for Negative Binomial emissions the dispersions \phi are updated by a bounded one-dimensional maximization of the weighted emission log-likelihood. The EM iteration is monotone in the log-likelihood; convergence is declared on a relative increase below tol.

Robustness. Several safeguards are applied: multiple restarts (the first deterministic, the rest random) keep the best fit; empty or near-empty states retain their previous parameters instead of collapsing; the RNG state of the caller is preserved; and the estimated states are relabeled into a canonical order (increasing \mu^I) so that the output is invariant to the arbitrary labeling of hidden states. Model-selection criteria AIC and BIC are returned to help choose nstates.

When dir_csv is supplied, a CSV named "hmm_states.csv" is written with columns t and state; when dir_out is supplied, the fitted object is saved as "hmm_fit.rds".

Value

An object of class "bivarhr_hmm": a list with

• fit: a list with nstates, family, the estimated parameters (mu_I, mu_C, and, for "nbinom", phi_I, phi_C), delta, Gamma, the maximized logLik, npar, AIC, BIC, the smoothed state posteriors (posterior, a T \times nstates matrix), and convergence information (iterations, converged).

• states: integer vector of Viterbi-decoded latent states.

If estimation fails (e.g., degenerate data), the function returns NULL.

Examples

DT <- data.frame(
  I = rpois(40, lambda = 4),
  C = rpois(40, lambda = 3)
)

# 'n_starts' is kept small here for speed; the default (10) is more robust.
res_hmm <- run_hmm(DT, nstates = 2, seed = 1, n_starts = 2)
if (!is.null(res_hmm)) {
  print(res_hmm)
  table(res_hmm$states)
}


# Negative Binomial emissions for overdispersed counts:
res_nb <- run_hmm(DT, nstates = 2, family = "nbinom", seed = 1, n_starts = 2)

bivarhr documentation built on July 7, 2026, 1:06 a.m.