Your First Hierarchical Model

Goal

Build, fit, and interpret a multi-market hierarchical model with partial pooling and optional CRE (Mundlak) correction using the DSAMbayes R API.

This page is a hands-on tutorial. For the broader methodological questions behind prior-setting, diagnostics, and interpretation, use the workflow pages:

Prerequisites

Dataset

This walkthrough uses the tracked data/geo_panel/demo_data_geo_panel.csv , a panel dataset with weekly observations across multiple geos. Key columns:

  • Response: revenue, weekly synthetic KPI per geo.
  • Group: geo, geo identifier.
  • Media: channel0_signal, channel1_signal, channel2_signal, channel3_signal.
  • Controls: t_scaled, sin52_1, cos52_1, events, newsletters.
  • Date: date, weekly date index.
library(DSAMbayes)
panel_df <- read.csv("data/geo_panel/demo_data_geo_panel.csv")
table(panel_df$geo)  # Check group counts

Step 1: Construct the hierarchical model

The (term | group) syntax tells DSAMbayes to fit random effects. Terms inside the parentheses get group-specific deviations from the population mean:

model <- blm(
  revenue ~
    t_scaled + sin52_1 + cos52_1 + events + newsletters +
    channel0_signal + channel1_signal + channel2_signal + channel3_signal +
    (1 + channel0_signal + channel1_signal + channel2_signal + channel3_signal | geo),
  data = panel_df
)

This specifies:

  • Population (fixed) effects for all terms, the average effect across markets.
  • Random intercepts and slopes for media terms by market, each market can deviate from the population average.

Step 2: Set boundaries

model <- model %>%
  set_boundary(
    channel0_signal > 0, channel1_signal > 0,
    channel2_signal > 0, channel3_signal > 0
  )

Boundaries apply to the population-level coefficients.

Step 3: (Optional) Add CRE / Mundlak correction

If you suspect that group-level spending patterns are correlated with unobserved market characteristics (e.g. high-spend markets also have higher baseline demand), CRE can model this correlation when the conditional mean of the group effect is adequately represented by the included group means:

model <- model %>%
  set_cre(vars = c(
    "channel0_signal", "channel1_signal",
    "channel2_signal", "channel3_signal"
  ))

This adds one cre_mean_<channel> term for each media variable as a fixed effect. Conditional on the CRE specification and other included controls, the original media coefficient represents the within-group temporal association. CRE does not solve omitted time-varying confounding or establish a causal effect.

See CRE / Mundlak for when and why to use this.

Step 4: Fit with MCMC

fitted_model <- model %>%
  fit(cores = 4, iter = 2000, warmup = 1000, seed = 42)

Hierarchical models are slower than BLM, expect 10–30 minutes depending on group count and data size. First-time Stan compilation of the hierarchical template adds 2–3 minutes.

Step 5: Check diagnostics

chain_diagnostics(fitted_model)

Pay special attention to Rhat and ESS for sd_* parameters (group-level standard deviations), which are often harder to estimate than population coefficients.

Step 6: Extract the posterior

post <- get_posterior(fitted_model)

For hierarchical models, coefficient draws from get_posterior() return vectors (one value per group) rather than scalars. The population-level (fixed-effect) estimates are averaged across groups.

Step 7: Group-level results

Fitted values and decomposition are returned per group:

# Fitted values, one row per observation, grouped by geo
fit_tbl <- fitted(fitted_model)
head(fit_tbl)

# Decomposition, per-group predictor contributions
tbl_map <- data.frame(
  actual.vars = c(
    "intercept", "t_scaled", "sin52_1", "cos52_1", "events", "newsletters",
    "cre_mean_channel0_signal", "cre_mean_channel1_signal",
    "cre_mean_channel2_signal", "cre_mean_channel3_signal",
    "channel0_signal", "channel1_signal", "channel2_signal", "channel3_signal"
  ),
  group = c(rep("BASE", 10), rep("media", 4)),
  ref_points = c("none", rep("mean", 13))
)
decomp_result <- decomp(fitted_model, tbl_map = tbl_map)
# One native decomposition result per retained group:
head(decomp_result$decomp[[1]]$DecompedData)

Step 8: Budget optimisation (population level)

Budget optimisation uses population-level (fixed-effect) beta draws, not group-specific totals:

# See Budget Optimisation docs for full scenario specification
result <- optimise_budget(fitted_model, scenario = my_scenario)

Key differences from BLM

Aspect BLM Hierarchical
Data structure Single market Panel (multiple groups)
Coefficient draws Scalars Vectors (one per group)
Fit time 2–5 min 10–30 min
Decomposition Direct Native result per retained geo; probabilistic media transforms are rejected
Forest/prior-posterior plots Direct Group-averaged population estimates
Stan template bayes_lm_updater_revised.stan general_hierarchical.stan (templated per group count)

Common pitfalls

Pitfall Symptom Fix
Too few groups Weak partial pooling; group SDs poorly estimated Need 4+ groups for meaningful hierarchical structure
Too few obs per group High Rhat on sd_* parameters Increase iterations; simplify random-effect structure
CRE with too many vars More CRE variables than groups Reduce CRE variable set; see identification warnings
CRE mean has zero variance scale=TRUE aborts with constant column error Use model.type: re (without CRE) or model.scale: false

Next steps