bivarhr: When Two Count Series Need to Talk — A Bayesian Framework for Bivariate Hurdle Models

You can also find this library at CRAN and download it directly from R and RStudio.

LISTEN TO THIS POST AS A PODCAST:

An introduction to the R package that brings together bivariate hurdle regression, horseshoe regularization, and multi-method causal inference under one roof.

The Problem Nobody Teaches You in Grad School

You have two count variables measured over time. Maybe they are insurgent attacks and counterinsurgent operations in a conflict zone. Maybe they are weekly disease case counts and mortality figures. Maybe they are criminal incidents and arrests in a city.

You know the standard toolkit: Poisson regression, negative binomial if there is overdispersion, maybe a zero-inflated model if zeros pile up. But then the real world intervenes:

• Both series have far too many zeros. Not the kind of “a few extra zeros” that zero-inflated models handle — the kind where 60–80% of observations are zero, and the rest are counts. This is the hallmark of hurdle data: a process that first decides whether anything happens at all, and then, if it does, decides how much.
• The two series influence each other over time. Do counterinsurgent operations today predict insurgent attacks next month? Does disease incidence Granger-cause mortality with a lag? You need cross-lagged dynamics, not two separate regressions.
• You have many potential predictors — economic indicators, population density, climate variables, regime dummies — and your sample size is modest. Overfitting is a real danger.
• You want causal claims, not just correlations, and a single method is not enough. You need converging evidence.

If this sounds familiar, the R package bivarhr was built for exactly your situation.

What bivarhr Does

bivarhr is an open-source R package (MIT-licensed, available on GitHub) that provides a unified workflow for:

1. Bivariate hurdle negative binomial regression — jointly modeling two zero-heavy count series with separate zero-generating and count-generating processes.
2. Horseshoe priors — automatic Bayesian regularization that shrinks irrelevant coefficients toward zero while letting strong signals survive.
3. Bayesian Model Averaging (BMA) via stacking — combining predictions across many candidate models (different lag orders, different regularization strengths) instead of betting everything on one.
4. Multi-method causal inference — transfer entropy, Dynamic Bayesian Networks, Hidden Markov Models, VARX models, synthetic control, and sensitivity analysis.
5. Rigorous validation — temporal placebo tests, rolling out-of-sample evaluation, extreme bounds analysis, and counterfactual average treatment effects.

All of this runs on top of Stan via the cmdstanr interface, which means full Hamiltonian Monte Carlo sampling with the No-U-Turn Sampler (NUTS) — the gold standard for Bayesian computation.

The Core Model: Bivariate Hurdle Negative Binomial

Let us unpack what the model actually does, without unnecessary jargon.

The Hurdle Idea

A hurdle model splits the data-generating process into two questions:

1. Did anything happen? (Zero vs. non-zero.) This is modeled with a Bernoulli distribution and a logit link — essentially a logistic regression.
2. If something happened, how much? (Positive counts only.) This is modeled with a truncated negative binomial distribution and a log link.

Mathematically, for each time point $t$ and series $I$:

$P(Y_{I,t} = 0) = 1 – \pi_{I,t}$

$P(Y_{I,t} = y \mid y > 0) = \pi_{I,t} \cdot \frac{f_{\text{NB}}(y \mid \mu_{I,t}, \phi_I)}{1 – f_{\text{NB}}(0 \mid \mu_{I,t}, \phi_I)}$

where $\pi_{I,t}$is the probability of crossing the hurdle (a non-zero count), $\mu_{I,t}$is the conditional mean of the negative binomial, and $\phi_I$ is the dispersion parameter.

The key distinction from zero-inflated models: a hurdle model does not distinguish between “structural zeros” and “sampling zeros.” It simply asks whether the count is zero or not, and if not, how large. This is conceptually cleaner and often more appropriate for event count data.

Bivariate Means Cross-Lags

The “bivariate” part means that the model is estimated jointly for two series, $I$ and $C$. What makes this powerful is that the design matrices for $I$ can include lagged values of $C$, and vice versa. Four specifications are available:

By comparing the predictive performance of these specifications, you can formally test whether one series has leading information about the other — a Granger-causality test conducted within a fully Bayesian framework.

Horseshoe Priors: Letting the Data Decide

When you have many candidate predictors and a modest sample size, ordinary maximum likelihood will happily overfit. The horseshoe prior solves this elegantly.

Each coefficient $\beta_j$ gets a prior that is a mixture of a tight “spike” near zero and a broad “slab” away from zero. The balance between spike and slab is governed by:

• A global shrinkage parameter $\tau$ that controls overall sparsity — how many coefficients the model expects to be non-zero.
• A local shrinkage parameter $\lambda_j$ for each coefficient — allowing individual signals to escape shrinkage if the data support it.

The regularized version used by bivarhr (following Piironen & Vehtari, 2017) adds a slab regularization term that prevents the unbounded coefficient estimates that can afflict the original horseshoe. The result: noise is suppressed, genuine signals are preserved, and you do not need to manually select variables.

In practice, you set a prior guess $\tau_0$for the fraction of coefficients you expect to be non-zero (e.g., 0.1 for strong sparsity, 0.5 for moderate). The model does the rest.

Bayesian Model Averaging: Stop Picking One Model

One of the most consequential analytical choices is the lag order: how many past time steps of $C$ should enter the equation for $I$? Different lag orders can tell different stories.

Rather than picking one, bivarhr uses stacking (Yao et al., 2018) to combine predictions across many models. The algorithm finds optimal weights $w_1, \ldots, w_M$such that the combined predictive distribution:

$p_{\text{stack}}(y_t) = \sum_{m=1}^{M} w_m \cdot p(y_t \mid y_{1:t-1}, \mathcal{M}_m)$

maximizes the expected log predictive density (ELPD), estimated via Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO).

This means the final inference is not tied to a single model. Models with good out-of-sample predictive performance get high weight; poor models are down-weighted automatically. The stacking weights themselves are informative: if specification A ($C → I$) consistently gets higher weight than specification D (no cross-lags), that is evidence of a genuine cross-series effect.

The Causal Inference Toolkit

A bivariate regression with cross-lags is a necessary but not sufficient condition for causal claims. bivarhr therefore provides six complementary causal inference methods:

Transfer Entropy

Transfer entropy measures how much knowing the past of series $I$ reduces your uncertainty about the current value of series $C$, beyond what $C$‘s own past already tells you. It is an information-theoretic generalization of Granger causality that makes no linearity assumptions.

bivarhr computes transfer entropy for three data transformations — raw counts, rates (counts per unit exposure), and binary presence/absence — each pre-whitened via an appropriate GLM to remove confounding trends. Permutation-based significance tests (with Benjamini-Hochberg correction) control the false discovery rate.

Dynamic Bayesian Networks (DBN)

A DBN learns a directed acyclic graph over discretized versions of $I$, $C$, and any regime variable, with edges restricted to flow from $t-1$ to $t$. The learned structure reveals which variables are direct parents of which — does yesterday’s regime predict today’s attacks, or only through yesterday’s operations?

Hidden Markov Models (HMM)

An HMM with multivariate Poisson emissions assumes that both $I$ and $C$ are driven by a latent (hidden) state that evolves over time. The inferred state sequence can reveal regimes — a “calm” state with low counts for both series, an “escalation” state with high counts, and so on — without you having to define them in advance.

VARX Models

A Vector Autoregression with Exogenous variables provides a classical time-series check. bivarhr fits bivariate VAR models and reports stability, serial correlation, normality, and ARCH diagnostics — a complementary perspective to the Bayesian hurdle model.

Sensitivity Analysis

Using the framework of Cinelli and Hazlett (2020), bivarhr quantifies how robust an OLS regression result is to an unobserved confounder. The output is a “robustness value”: how strong would an omitted confounder need to be (in terms of partial $R^2$) to explain away the estimated effect? This gives a concrete, interpretable measure of fragility.

Synthetic Control

A Bayesian Structural Time Series (BSTS) model constructs a synthetic counterfactual: what would series $I$ have looked like in the absence of a treatment (e.g., a policy intervention, a conflict escalation)? The difference between observed and counterfactual is the estimated causal effect, with full posterior uncertainty intervals.

Validation: Proving Your Results Are Not Artifacts

Temporal Placebo Test

Shuffle the time indices of your data, refit the model, and compare the ELPD to the original. If the original model’s ELPD is substantially higher, the temporal structure is genuine — not an artifact of overfitting.

Rolling Out-of-Sample Evaluation

Split the data at multiple cut points (60%, 70%, 80%, 90% of the sample), fit the model on the training portion, and forecast the remaining observations. The RMSE of these forecasts gives an honest assessment of predictive accuracy.

Extreme Bounds Analysis (EBA)

Re-fit the model with every combination of control variables and check whether the key coefficients (e.g., cross-lag effects) remain stable. If a coefficient flips sign or loses significance when you add or remove a single control, it is fragile.

G-Computation for Average Treatment Effects

Counterfactual ATEs answer the question: “What would happen to $I$ if we set the cross-lags and transition variables to zero?” This is the causal estimand that many applied researchers actually care about, computed via posterior predictive simulation.

Getting Started

Installation

# Install the package from GitHub
devtools::install_github("isadorenabi/bivarhr")
# Install cmdstanr (required for model fitting)
install.packages("cmdstanr",
                 repos = c("https://stan-dev.r-universe.dev",
                           getOption("repos")))
cmdstanr::install_cmdstan()
# Optional: install packages for full causal inference functionality
install.packages(c("RTransferEntropy", "bnlearn", "sensemakr",
                   "CausalImpact", "vars", "openxlsx"))

Minimal Example

library(bivarhr)
library(data.table)
# Simulate data
n <- 100
DT <- data.table(
  I = rpois(n, 5),
  C = rpois(n, 3),
  zI = as.integer(rpois(n, 5) > 0),
  zC = as.integer(rpois(n, 3) > 0),
  t_norm = seq(-1, 1, length.out = n),
  t_poly2 = seq(-1, 1, length.out = n)^2,
  Regime = factor(sample(c("A", "B"), n, TRUE)),
  trans_PS = sample(0:1, n, TRUE),
  trans_SF = sample(0:1, n, TRUE),
  trans_FC = sample(0:1, n, TRUE),
  log_exposure50 = rep(0, n)
)
# Fit a single model with bidirectional cross-lags and 2 lags
fit <- fit_one(DT, k = 2, spec = "C",
               iter_warmup = 500, iter_sampling = 500, chains = 2)
# Inspect
print(fit$fit$summary())

Bayesian Model Averaging Across Lags and Hyperparameters

# Define a grid of horseshoe hyperparameters
hs_grid <- expand.grid(
  hs_tau0 = c(0.1, 0.5, 1.0),
  hs_slab_scale = c(1, 5),
  hs_slab_df = 4
)
# Run BMA over lag orders 0 through 3
bma_results <- select_by_bma(
  DT = DT, spec = "C", k_grid = 0:3, hs_grid = hs_grid,
  iter_warmup = 1000, iter_sampling = 1200, chains = 4,
  use_parallel = TRUE
)
# View model ranking by ELPD
print(bma_results$table)

Transfer Entropy

te_results <- run_transfer_entropy(DT, lags = 1:3, shuffles = 200, seed = 123)
print(te_results)

Where This Matters

The applications extend well beyond conflict analysis. Any domain with paired zero-heavy count time series is a candidate:

• Public health: Disease incidence and mortality; hospital admissions and readmissions; vaccination rates and outbreak counts.
• Criminal justice: Criminal incidents and law enforcement responses; drug seizures and overdose deaths.
• Economics: Firm entry and exit counts; patent applications and citations; trade flow counts between countries.
• Ecology: Predator and prey counts; species occurrence and extinction events; colonization and local disappearance.
• Social science: Protest events and government responses; legislative proposals and vetoes; online misinformation posts and fact-checking activity.

Under the Hood: Technical Details

For those who want to look inside:

The Stan model code is fully transparent — you can inspect it via get_hurdle_model() — and the generated quantities block produces posterior predictive checks, pointwise log-likelihoods, and fitted values out of the box.

Diagnostic Checklist

Before trusting your results, bivarhr encourages you to verify:

The Bigger Picture

What makes bivarhr distinctive is not any single component — hurdle models exist, horseshoe priors exist, transfer entropy exists — but the integration. The package recognizes that in applied work, the bottleneck is not any one method’s sophistication. It is the workflow: getting from raw data to a validated causal claim, with uncertainty quantification, in a reproducible pipeline.

By combining a carefully specified bivariate hurdle model with automatic regularization, principled model averaging, and six complementary causal inference methods — all validated through placebo tests and out-of-sample evaluation — bivarhr gives applied researchers a framework that is both rigorous and practical.

If your work involves paired count time series with excess zeros and you care about causal inference, this package deserves a serious look.

bivarhr is open-source (MIT License) and available at github.com/isadorenabi/bivarhr.

Author: José Mauricio Gómez Julián. ORCID.

Citation:

@software{gomez2025bivarhr,
  author = {Gómez Julián, José Mauricio},
  title = {bivarhr: Bivariate Hurdle Regression with Bayesian Model Averaging},
  year = {2025},
  url = {https://github.com/isadorenabi/bivarhr}
}

Built on Stan, cmdstanr, and the broader Bayesian ecosystem. The author thanks the Stan development team for their foundational work.