Your First BLM Model

Goal

Build, fit, and interpret a single-market Bayesian linear model (BLM) using the DSAMbayes R API.

Prerequisites

  • DSAMbayes installed locally (see Install and Setup).
  • Familiarity with R and lm()-style formulas.

Dataset

This walkthrough uses the synthetic dataset shipped at data/synthetic_dsam_example_wide_data.csv. It contains weekly observations for a single market with columns for:

  • Response: kpi_os_hfb01_value — weekly KPI (e.g. revenue).
  • Media: m_tv, m_search, m_social, m_display, m_ooh, m_email, m_affiliate.
  • Controls: t_scaled (trend), sin52_1/cos52_1/sin52_2/cos52_2 (Fourier seasonality), price_index, distribution, bm_ikea_trust_12r (brand metric).
  • Holidays: h_black_friday, h_christmas, h_new_year, h_easter, h_summer_sale.
library(DSAMbayes)
df <- read.csv("data/synthetic_dsam_example_wide_data.csv")
str(df)

Step 1: Construct the model

blm() creates an unfitted model object. No fitting happens yet.

model <- blm(
  kpi_os_hfb01_value ~
    t_scaled + sin52_1 + cos52_1 + sin52_2 + cos52_2 +
    h_black_friday + h_christmas + h_new_year + h_easter + h_summer_sale +
    bm_ikea_trust_12r + price_index + distribution +
    m_tv + m_search + m_social + m_display + m_ooh + m_email + m_affiliate,
  data = df
)

Inspect defaults:

peek_prior(model)      # normal(0, 5) for all terms
peek_boundary(model)   # (-Inf, Inf) — unconstrained

Step 2: Set boundaries

Media channels should have non-negative effects. Use inequality notation:

model <- model %>%
  set_boundary(
    m_tv > 0, m_search > 0, m_social > 0,
    m_display > 0, m_ooh > 0, m_email > 0, m_affiliate > 0
  )

Step 3: (Optional) Override priors

Default priors are weakly informative. Override only with domain knowledge:

model <- model %>%
  set_prior(price_index ~ normal(-0.2, 0.1))

See Minimal-Prior Policy for guidance.

Step 4: Fit with MCMC

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

First-time Stan compilation takes 1–3 minutes. Subsequent runs use a cached binary. With 2 chains on synthetic data, sampling typically completes in under 2 minutes.

Step 5: Sampler diagnostics

chain_diagnostics(fitted_model)
Metric Good Concern
Max Rhat < 1.01 > 1.05 means chains have not converged
Min ESS (bulk) > 400 < 200 means too few effective samples
Divergences 0 Any non-zero count warrants investigation

Step 6: Extract the posterior

post <- get_posterior(fitted_model)

post is a tibble with one row per draw containing coef (named coefficient list), yhat (fitted values), noise_sd, r2, rmse, and smape.

Summarise coefficients:

library(dplyr); library(tidyr)

coef_summary <- post %>%
  select(coef) %>%
  unnest_wider(coef) %>%
  pivot_longer(everything(), names_to = "term") %>%
  group_by(term) %>%
  summarise(
    mean = mean(value), median = median(value), sd = sd(value),
    ci_low = quantile(value, 0.025), ci_high = quantile(value, 0.975),
    .groups = "drop"
  )
print(coef_summary, n = 30)

What to look for:

  • Media coefficients should be positive (boundaries enforce this).
  • Wide credible intervals mean the prior dominates — the data cannot identify the effect precisely.

Step 7: Assess model fit

fit_tbl <- fitted(fitted_model)
cat("Median R²:", median(r2(fitted_model)), "\n")
cat("Median RMSE:", median(rmse(fitted_model)), "\n")

For a well-specified MMM on weekly data, in-sample R² above 0.85 is typical.

Step 8: Response decomposition

decomp_tbl <- decomp(fitted_model)
head(decomp_tbl)

Shows each term’s contribution (coefficient × design-matrix column) to the predicted KPI at each time point — the foundation for media contribution and ROI reporting.

Step 9: MAP for rapid iteration

During development, use MAP for fast point estimates:

map_model <- model %>% optimise(n_runs = 10)
get_posterior(map_model)

Use MCMC for final reporting; MAP for formula iteration.

Common pitfalls

Pitfall Symptom Fix
Forgetting to set boundaries Media coefficients go negative Add set_boundary(m_x > 0)
Too few iterations High Rhat, low ESS Increase iter and warmup
Missing controls High residual autocorrelation Add trend, seasonality, or holiday terms
Scaling confusion Coefficients look wrong model.scale: true is default; get_posterior() back-transforms automatically

Next steps