Model Classes

Purpose

dsambayesplus provides four model classes for Bayesian marketing mix modelling. Each class targets a different data structure and pooling strategy. This page describes the constructor pathways, fit support, and practical limitations of each class so that an operator can select the appropriate model for a given dataset.

The fixed-effects (FE) class is a separate within-unit estimator. A hierarchical model with unit random intercepts remains a random-effects model, not an FE estimator.

Use this page to choose the right modelling surface for your data structure and decision problem. Do not use it as a substitute for the broader workflow: class selection does not settle prior design, computational trustworthiness, model adequacy, or causal interpretation. For that framing, start with What Principled Means, Stage 2: Model and Priors, and Stage 5: Model Adequacy.

Selection discipline

• Choose the simplest class that matches the real data structure.
• Do not move to hierarchical or pooled models only because they look more advanced.
• Treat model class as a structural choice, not proof that the resulting model is decision-ready.
• After choosing a class, return to the workflow pages for prior-setting, diagnostics, and interpretation discipline.

Class summary

Fixed effects (fixed_effects)

Construction

model <- fixed_effects(kpi ~ m_tv + m_search + trend, panel_df, market) |>
  set_date(week)

FE v2 removes unrestricted unit intercepts with deterministic orthonormal within-unit contrasts. It estimates common slopes and residual noise through MCMC. Formula terms must be bare numeric columns, and retained unit-date keys must be unique. Unequal unit sizes and irregular date spacing are allowed.

Output boundary

Use get_posterior() for reporting-scale slope and residual-noise draws and within-contrast generated quantities. After MCMC, use fixed_effects_intercepts() for conditional retained-unit intercept draws, fitted() and observed() for retained level data, and predict() for new rows from retained units. Its runner also writes an original-row conditional level-diagnostic table for retained fitting rows. FE v2 rejects MAP, unseen-unit prediction, generic residual diagnostics, PSIS model selection, decomposition, counterfactual analysis, optimisation, forecasting, and deployment. The YAML runner supports validation, dry-run, and bounded MCMC fitting with dedicated coefficient and contrast-space artefacts. It also supports opt-in seen-unit expanding-window validation: each holdout unit must be present in that fold’s training panel, and each training unit must have at least two observations. It does not write generic level-scale output files.

P4-J’s bounded live runner, explicit-prior prior-only, one-replication recovery, and sampler gates passed. These checks do not estimate repeated-sampling coverage or production reliability. The implementation is not production-qualified, so FE remains implemented_not_qualified. See Estimator Capabilities for the full contract.

BLM (blm)

Construction

model <- blm(kpi ~ m_tv + m_search + trend + seasonality, data = df)

blm() dispatches on the first argument. When passed a formula, it creates a blm object with default priors and boundaries. When passed an lm object, it creates a bayes_lm_updater whose priors are initialised from the OLS coefficient estimates and standard errors.

Fit support

Post-fit accessors

• get_posterior(), coefficient draws, fitted values, metrics
• fitted(), predicted response on the original scale
• decomp(), native predictor-level decomposition. It returns a dsambayesplus_decomposition S3 result for BLM objects; use $DecompedData for the long-form contribution table. Hierarchical models return a keyed tibble with one such result in each row’s decomp list-column.
• optimise_budget(), decision-layer budget allocation

Decomposition migration

decomp() no longer returns the external DSAMdecomp S4 calculators. This is a breaking development-version change. Replace S4 slot access such as result@DecompedData with result$DecompedData; use result$DecompGroupProps for grouped totals. The DSAMdecomp-only plotting and post-processing helpers are not part of the native result. Build any required downstream summaries from these two tables.

Limitations

• No group structure. For multi-market data, use the hierarchical class.

Hierarchical (hierarchical)

Construction

The hierarchical class is created automatically when blm() detects random-effects syntax (|) in the formula:

model <- blm(
  kpi ~ m_tv + m_search + trend + (1 + m_tv + m_search | market),
  data = panel_df
)

Terms to the left of | become random slopes; the variable to the right defines the grouping factor. Multiple grouping terms are supported.

The ordinary random-effects specification partially pools group effects and uses both within-group and between-group information for population coefficients. Its interpretation requires the group effects to be conditionally mean-independent of the regressors. If this assumption is not credible, use a justified CRE specification or another design. Neither option solves omitted time-varying confounding by itself.

CRE / Mundlak extension

For correlated random effects, call set_cre() after construction:

cre_model <- blm(
  kpi ~ m_tv_signal + m_search_signal + trend + (1 | market),
  data = panel_df
) |>
  set_cre(vars = c("m_tv_signal", "m_search_signal"))

CRE v1 is restricted to a selected random-intercept group with no random slope in that group. Each CRE variable must be one simple numeric main effect and must not overlap an internal draw-dependent media transform. Fit preparation calculates the Mundlak means from the exact retained estimation frame. See CRE / Mundlak for details.

Fit support

Post-fit accessors

Same as BLM. Coefficient draws from get_posterior() return vectors (one value per group) rather than scalars. Budget optimisation uses the population-level (fixed-effect) coefficient draws from the beta parameter.

Limitations

• Stan template compilation uses a templated source (general_hierarchical.stan) rendered per number of groups and parameterisation mode. First compilation is slow; subsequent runs use a cached binary.
• Response decomposition is native and returns one result per retained group. Probabilistic media transforms remain unsupported for decomposition and are rejected explicitly.
• Formula offsets enter native decomposition as coefficient-one baseline components and must be included in the decomposition mapping table.
• Posterior forest and prior-vs-posterior plots average group-specific draws to produce a single population-level estimate.
• Offset support in the hierarchical Stan template is handled via stats::model.offset() within build_hierarchical_frame_data().

Pooled (pooled)

Construction

The pooled class is created by converting an existing BLM object with pool():

base <- blm(kpi ~ m_tv + m_search + trend + seasonality, data = df)
model <- pool(base, grouping_vars = c("channel"), map = pooling_map)

The map is a data frame with a variable column mapping formula terms to pooling dimension labels. Exact formula-term labels are preferred; raw variable names are accepted only when they resolve unambiguously to a single non-offset formula term. Priors and boundaries are reset to defaults when pool() is called.

Fit support

MAP fitting (fit_map) is not currently implemented for pooled models.

Post-fit accessors

Same as BLM. The design matrix is split into base terms (intercept + non-pooled) and media terms (pooled). The Stan template uses a per-dimension coefficient structure.

Limitations

• MAP fitting is not available.
• extract_stan_design_matrix() may return a zero-row matrix, which causes VIF computation to be skipped.
• The pooled Stan cache key includes sorted grouping variable names to avoid collisions between different pooling configurations.
• Time-series cross-validation is not supported for pooled models.

Class selection guide

In practice, the class decision should usually be driven by three questions:

1. Is the dataset a single time series or a grouped panel?
2. Do you need partial pooling across real groups, or pooling across labelled coefficient dimensions?
3. Is the added structure necessary for the business question, or are you adding complexity without a clear identifiability benefit?

Fit method selection

For inferential runs where diagnostics and uncertainty quantification matter, MCMC is the recommended fit method. MAP is useful for rapid iteration during model development. Fit method alone does not qualify an estimator for production use.

MAP returns one selected optimum, not posterior draws. Do not derive credible intervals or MCMC diagnostics from it; its point estimate can understate uncertainty, especially for hierarchical variance components. fit_map() retains the restart results for inspection (and runner outputs include optimisation_runs.csv), so materially different restart objectives should be treated as optimisation instability or competing local optima, not as a substitute for posterior uncertainty.

Cross-references

• Model Object Lifecycle, state transitions from construction to post-fit
• Priors and Boundaries, prior and constraint specifications
• CRE / Mundlak, correlated random effects for hierarchical models
• Diagnostics Gates, post-fit quality checks
• Stage 2: Model and Priors, workflow framing for structural choices and prior design
• Stage 5: Model Adequacy, what must still be checked after selecting a class

Model Object Lifecycle

dsambayesplus model objects (blm, hierarchical, pooled) are mutable S3 lists that progress through a well-defined sequence of states. Understanding these states helps avoid calling post-fit accessors on an unfitted object or forgetting to compile before fitting.

This page is an API/runtime reference. It explains how dsambayesplus model objects move through construction, compilation, fitting, and post-fit access. It is not the main guide for prior-setting, diagnostics meaning, or model adequacy. For that, use the Principled Bayesian Workflow.

State-machine diagram

                        ┌─────────────────────────────────────────┐
                        │           CREATED                       │
                        │  blm(), blm.formula(), blm.lm(),        │
                        │  as_bayes_lm_updater()                  │
                        │  Fields set: .formula, .original_data,  │
                        │    .prior, .boundaries                  │
                        └──────────────┬──────────────────────────┘
                                       │
            ┌──────────────────────────┼──────────────────────────┐
            ▼                          ▼                          ▼
    set_prior(obj, …)         set_boundary(obj, …)        set_date(obj, …)
    Mutates .prior            Mutates .boundaries         Sets .date_var
            │                          │                          │
            └──────────────────────────┼──────────────────────────┘
                                       │
                    (optional: pool() transitions blm → pooled,
                     resets .prior/.boundaries, adds .pooling_vars/.pooling_map)
                                       │
                                       ▼
                        ┌─────────────────────────────────────────┐
                        │           CONFIGURED                    │
                        │  Priors, boundaries, date variable are  │
                        │  set (may still use defaults).          │
                        └──────────────┬──────────────────────────┘
                                       │
                                       ▼
                        ┌─────────────────────────────────────────┐
                        │           COMPILED                      │
                        │  compile_model(obj)                     │
                        │  Sets .stan_model                       │
                        │  (pre_flight_checks auto-compiles if    │
                        │   .stan_model is NULL)                  │
                        └──────────────┬──────────────────────────┘
                                       │
                                       ▼
                        ┌─────────────────────────────────────────┐
                        │           PRE-FLIGHTED                  │
                        │  pre_flight_checks(obj, data)           │
                        │  Validates formula/data compatibility,  │
                        │  auto-compiles and auto-sets date_var   │
                        │  if missing. Sets .response_transform,  │
                        │  .response_scale.                       │
                        └──────────────┬──────────────────────────┘
                                       │
                                       ▼
                        ┌─────────────────────────────────────────┐
                        │           FITTED                        │
                        │  fit(obj) / fit_map(obj)                │
                        │  Calls pre_flight_checks internally,    │
                        │  then prep_data_for_fit → rstan.        │
                        │  Sets .stan_data, .date_val, .posterior │
                        └──────────────┬──────────────────────────┘
                                       │
                   ┌───────────────────┼───────────────────┐
                   ▼                   ▼                   ▼
           get_posterior(obj)    fitted(obj)         decomp(obj)
           Returns tibble of    Predicted values    Decomposition
           posterior draws      (yhat)              native dsambayesplus result
                   │                                       │
                   ▼                                       ▼
           optimise_budget(obj, …)                         Further analysis
           Budget allocation
           (requires fitted model)

States and key fields

Post-fit accessors

These functions require a fitted model (.posterior is not NULL):

Guards and auto-transitions

• pre_flight_checks() auto-compiles via compile_model() if .stan_model is NULL, and auto-sets .date_var to "date" if not already set.
• fit() and fit_map() call pre_flight_checks() internally, so explicit compilation is optional.
• get_posterior() aborts with a clear error if .posterior is NULL.

Object field reference

All fields are initialised by model_object_schema_defaults() in R/model_schema.R. The canonical field list:

For hierarchical MCMC, fit.mcmc.parameterization.positive_priors sets this field to centered (the default) or noncentered. Compilation uses it to select the positive-parameter prior parameterisation in the Stan template. The non-centred option requires lognormal_ms priors for positive hierarchical parameters.

Runner-injected fields

These are set by the YAML/CLI runner (run_from_yaml()) for artifact writing and are not part of the core modelling API:

• .runner_config, .runner_kpi_type, .runner_identifiability
• .runner_time_components, .runner_budget_optimisation
• .runner_model_selection, .runner_model_type

Priors and Boundaries

For the workflow guidance behind these controls, start with Stage 2: Model and Priors. This page is the technical contract for dsambayesplus prior and boundary behaviour.

Use this page when you need exact dsambayesplus semantics: supported prior families, override syntax, default generation rules, and scaling behaviour. Do not use it as the main argument for why a prior is reasonable. That reasoning belongs in the workflow pages and in your modelling rationale.

Purpose

This page defines how dsambayesplus specifies, defaults, overrides, and scales coefficient priors and parameter boundaries for all model classes. It covers the prior schema, supported families, default-generation logic, YAML override contract, and the interaction between priors, boundaries, and the scale=TRUE pathway.

How to read this page

• Use Minimal-Prior Policy if you want the short recommended operating rule.
• Use this page when you need to know exactly how dsambayesplus will interpret a prior or boundary specification.
• Return to Stage 2: Model and Priors if the question is whether a custom prior should be added at all.

Prior schema

Each model object carries a .prior tibble with one row per parameter. The columns are:

Supported prior families

All coefficient priors use normal(). The lognormal_ms family is parameterised by the mean and standard deviation on the original (non-log) scale; dsambayesplus converts these internally to log-space parameters.

Default prior generation

BLM and hierarchical (population terms)

default_prior.blm() calls standard_prior_terms(). On the reporting scale, the intercept prior is normal(ybar, sy), each slope prior is normal(0, sy / sx), and the residual-SD prior is normal(0, sy), where sy and sx are the response and term standard deviations. These defaults therefore map to unit-scale priors after internal standardisation and respond coherently when measurement units change.

Hierarchical (group-level standard deviations)

default_prior.hierarchical() additionally generates sd_<idx>[<term>] rows for each group factor. Each constrained group SD receives a zero-centred normal prior whose scale matches the corresponding coefficient scale: sy for a group intercept and sy / sx for a random slope. Because the SD parameter is constrained positive in Stan, this acts as a half-normal prior. It remains positive when observed group outcome means happen to be equal, and random-slope defaults change coherently when predictor measurement units change.

Specify hierarchical standard-deviation priors in reporting-scale units. With scale = TRUE, dsambayesplus transforms them to the Stan scale: intercept heterogeneity is divided by the response standard deviation, while slope heterogeneity is multiplied by the predictor-to-response standard-deviation ratio. The group intercept is defined at the internally centred predictor reference point. It is not a group intercept at raw predictor values of zero. With the defaults above, these transformations produce unit-scale group-SD priors in Stan space. Explicit user overrides retain their supplied reporting-scale meaning and are transformed by the same rules.

Default population and group-SD scales use the same population-model complete-case frame. This keeps the reporting-to-Stan transformation exact when a population term contains missing values. Random-slope terms must be numeric. Group-only random-slope terms are evaluated on the retained population rows. dsambayesplus aborts rather than deriving a group-SD default when a group term would drop additional rows or when a factor random slope is unsupported.

BLM from lm (Bayesian updating)

default_prior.bayes_lm_updater() initialises coefficient priors from the OLS point estimates (mean) and standard errors (sd), enabling informative Bayesian updating.

Pooled

default_prior.pooled() uses the BLM defaults for non-pooled terms (intercept, base regressors, noise_sd) and normal(0, 5) for each dimension-level pooled coefficient. Default pooled boundaries remain unconstrained; add explicit boundaries if a pooled dimension should be sign-restricted.

Boundary schema

Each model object carries a .boundaries tibble with one row per parameter:

Default boundaries are lower = -Inf, upper = Inf for all terms. No sign constraints are imposed by default.

YAML override contract

Prior overrides

priors:
  use_defaults: true
  overrides:
    - { parameter: m_tv, mean: 0.5, sd: 0.2 }
    - { parameter: price_index, mean: -0.2, sd: 0.1 }

Each override replaces the distribution call for the named parameter with normal(mean, sd). Overrides are sparse: only the listed parameters are changed; all other parameters keep their defaults.

In M1, use_defaults must remain true. The v2 runner is default-first: it always starts from the generated prior table, then applies sparse grouped aliases and explicit overrides.

The friendly YAML surface also accepts an explicit alias style:

priors:
  media_beta:
    distribution: HalfNormal
    sigma: 1

or:

priors:
  intercept:
    distribution: Normal
    mu: 120000
    sigma: 15000

HalfNormal is implemented by compiling to normal(0, sigma) plus an implied lower bound of 0 on targeted parameters that are otherwise unconstrained. The priors.likelihood.sigma alias compiles to the dsambayesplus noise_sd prior family.

Boundary overrides

boundaries:
  overrides:
    - { parameter: m_tv, lower: 0.0, upper: .Inf }
    - { parameter: competitor_discount, lower: -.Inf, upper: 0.0 }

Each override replaces the boundary entry for the named parameter. YAML infinity tokens (.Inf, -.Inf) are coerced during config resolution.

Scale semantics (scale = TRUE)

When model.scale: true (the default), the response and predictors are standardised before Stan fitting. Centring applies only to population formulae with an intercept. Through-origin formulae divide by the same standard deviations without centring. This affects both priors and boundaries.

Coefficient prior scaling

For BLM and hierarchical population terms, the fitted intercept refers to centred predictors on the standardised model-response scale. The reported population intercept is back-transformed to original predictor zero. For a log-response model, the unstandardised model response is log-KPI, not KPI.

Let cy and sy be the centre and SD of the model response, and cx[j] and sx[j] the corresponding predictor quantities. For authored normal locations m_alpha, m_beta[j] and prior SDs sd_alpha, sd_beta[j], scale_.prior() implements:

intercept location in Stan = (m_alpha - cy + sum(cx * m_beta)) / sy
intercept prior SD in Stan = sd_alpha / sy
slope location in Stan[j]  = m_beta[j] * sx[j] / sy
slope prior SD in Stan[j]  = sd_beta[j] * sx[j] / sy

These locations are the parameters of normal kernels; with finite bounds they are not generally the means of the truncated distributions. Stan assigns independent normal priors to the population coefficients, restricted by their bounds. For noise_sd, the normal location and SD, or the lognormal_ms mean and SD, are divided by sy.

The reported coefficients use the inverse map:

beta[j] = sy * beta_stan[j] / sx[j]
alpha   = sy * alpha_stan + cy - sum(cx * beta)

The current conversion does not generally preserve an independently authored raw-intercept prior jointly with the slopes. The inverse intercept depends on the slope draws, but prior scaling divides its SD only by sy and introduces no intercept-slope covariance. For example, normal(10, 2) for the intercept and normal(3, 4) for one slope, with cy=10, sy=2, cx=5, sx=1, give independent internal priors normal(7.5, 1) and normal(1.5, 2). Without bounds, the reported intercept then has variance 404 rather than the authored 4, and covariance -80 with the reported slope. This is an algebraic consequence of the implementation, not a fitted result.

See R/scale.R (scale_.prior() and unscale_coef_vector()), the preparation paths in R/blm.R and R/hierarchy.R, and the Stan population-prior blocks.

Boundary scaling

• Zero slope boundaries (0) are preserved by positive multiplicative scaling.
• Infinite boundaries (±Inf) are invariant under scaling.
• Finite non-zero boundaries for slope terms are scaled using scale_boundary_for_parameter(), which applies the same sx / sy ratio used for slope priors.
• Population-intercept boundaries in BLM and hierarchical models: with scale = TRUE, finite non-zero bounds are rejected. Zero and infinite bounds pass through unchanged. A zero lower bound constrains the internal centred intercept, not necessarily the reported raw-zero intercept. In the example above, alpha_stan=0.1 and beta_stan=2 satisfy that internal bound but back-transform to alpha=-9.8.
• If a finite non-zero boundary is specified for a parameter without a matching scale factor in the design matrix, dsambayesplus aborts with a validation error.

Practical implication

State the model-response units, scaling flag and intercept reference explicitly. Unstandardised coefficient reporting does not imply that centring preserved the authored joint prior or an intercept bound at raw predictor zero. With scale = FALSE, BLM and hierarchical population priors and bounds are passed to the unstandardised fit without these conversions. Through-origin models have no population intercept and retain zero centres when scaling.

Hierarchical group-intercept SDs have the separate centred-reference convention described above. Do not apply the population-intercept interpretation to group heterogeneity or generalise these boundary rules to pooled dimension terms.

Interaction with model classes

Programmatic API

Inspect priors and boundaries

peek_prior(model)
peek_boundary(model)

Override priors

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

Override boundaries

model <- model %>%
  set_boundary(
    m_tv > 0,
    competitor_discount < 0
  )

Minimal-prior policy

The recommended operating profile for MMM is documented in Minimal-Prior Policy. The policy keeps priors weak by default and uses hard constraints only when there is structural business knowledge.

Cross-references

• Model Classes, constructor and fit support per class
• Minimal-Prior Policy, governance guidance for prior specification
• Response Scale Semantics, link vs KPI scale behaviour
• Config Schema, YAML prior and boundary keys

Minimal-Prior Policy

This page is the short operating rule for prior-setting in dsambayesplus. Use it when you want a compact default policy. For the full workflow logic, see Stage 2: Model and Priors. For the mechanics of YAML and API prior specification, see Priors and Boundaries.

Purpose

Use a principled but low-friction prior setup that avoids specification-hunting while preserving identifiability in short, collinear MMM datasets.

Policy

1. Default-first: keep priors.use_defaults: true.
2. Sparse overrides: only add priors.overrides for high-conviction terms.
3. Selective bounds: add boundaries.overrides only for structural signs.
4. No blanket constraints: do not force all controls/media to one sign by default.
5. Diagnose before tightening: use pre-flight and diagnostics gates first, then add priors/bounds if uncertainty is still unstable.

YAML mapping

priors:
  use_defaults: true
  overrides:
    # Optional high-conviction override examples:
    # - { parameter: price_index, mean: -0.2, sd: 0.1 }
    # - { parameter: distribution, mean: 0.15, sd: 0.1 }

boundaries:
  overrides:
    # Optional structural sign constraints:
    # - { parameter: m_tv, lower: 0.0, upper: .Inf }
    # - { parameter: competitor_discount, lower: -Inf, upper: 0.0 }

When to override defaults

• Do override when domain mechanism is stable and defensible.
• Do not override only to improve one run’s fit metrics.
• Do not add bounds if sign can plausibly flip under promotion, pricing, or substitution effects.

Review checklist

• Are overrides fewer than the number of major business assumptions?
• Is each bound tied to a concrete causal rationale?
• Did diagnostics indicate a real identifiability problem before tightening?

When to use this page

• Use this page when you need a concise prior-setting policy for routine MMM work.
• Use Stage 2: Model and Priors when you need the reasoning behind that policy.
• Use Priors and Boundaries when you need exact dsambayesplus syntax, scaling rules, or boundary mechanics.

Response Scale Semantics

Purpose

dsambayesplus models can operate on an identity (level) or log response scale. This page defines how response scale is detected, stored, and used for post-fit reporting, so that operators understand which scale their outputs are on and how KPI-scale conversions work.

Response scale detection

Response scale is determined at construction time by detect_response_scale(), which inspects the left-hand side of the formula:

The detected value is stored in two model-object fields:

• .response_transform: "identity" or "log". Describes the mathematical transform applied to the response before modelling.
• .response_scale: "identity" or "log". Used as a label when reporting whether outputs are on the model scale or the KPI scale.

Both fields are set by the constructor and confirmed by pre_flight_checks().

Model scale vs KPI scale

For identity-response models, model scale and KPI scale are identical. For log-response models, fitted values and residuals on the model scale are in log units and must be exponentiated to obtain KPI-scale values.

This is a semilog model, not a log-log model. In dsambayesplus, a coefficient from log(kpi) ~ x means:

So for a one-unit increase in x, the exact KPI-scale percentage change is:

The common shortcut 100 * beta is only a small-coefficient approximation.

Interpreting log-response models

This is the section to use when an analyst asks, “what does the coefficient actually mean on the KPI scale?”

1. Coefficients stay on the model scale

For a model written as:

the coefficient beta returned by get_posterior() and summarised in posterior_summary.csv is a log-KPI coefficient. dsambayesplus does not silently convert coefficient tables into KPI-scale percentage effects.

2. The exact KPI-scale effect depends on the predictor change

For a change of \Delta x in a predictor, the model implies:

Special cases:

• If \Delta x = 1, the exact percent change is 100 * (exp(beta) - 1).
• If x is a binary indicator changing from 0 to 1, use the same exact formula.
• The shortcut 100 * beta is only acceptable when beta * \Delta x is small enough that the approximation error is negligible for the use case.

3. This is not automatically an elasticity

log(kpi) ~ x is a semilog model. The coefficient is an elasticity only if the predictor is also logged, for example log(kpi) ~ log(x).

So in dsambayesplus:

• log(kpi) ~ x gives a semilog coefficient.
• log(kpi) ~ log(x) would be interpreted as an elasticity.

4. Coefficients attach to the modeled column, not necessarily raw spend

dsambayesplus coefficients describe the predictor that actually enters the model matrix.

That matters because in MMM workflows the modeled term is often:

• an adstocked media signal,
• a saturated transform,
• a normalised exposure metric,
• or another user-authored transformed column.

So if your YAML media block points to transformed signal columns, the coefficient is per unit of that transformed signal, not per unit of raw spend. The same caution applies to interactive formula workflows.

5. Use the right output for the question

Use these surfaces consistently:

• posterior_summary.csv and get_posterior() for coefficient summaries on the model scale.
• fitted.csv and observed.csv for fitted and observed values on the model scale.
• fitted_kpi.csv, observed_kpi.csv, and fitted_kpi() for business-facing values on the KPI scale.

For log-response models, posterior_summary.csv is therefore the wrong place to read off a KPI-scale uplift directly. It is the right place to get beta, which you then interpret with 100 * (exp(beta * \Delta x) - 1).

6. dsambayesplus labels KPI-scale conversions explicitly

When dsambayesplus writes KPI-scale outputs for log-response models, it records:

• source_response_scale = "log"
• response_scale = "kpi"
• conversion_method

This is intended to make it obvious that the values have been back-transformed and to distinguish the default lognormal-mean conversion from the simpler pointwise exp() median-style conversion.

Post-fit accessors and scale behaviour

fitted(), model scale

fitted() returns predicted values on the model scale. For identity-response models this is the KPI scale. For log-response models this is the log scale.

fit_tbl <- fitted(model)
# fit_tbl$fitted is on model scale

fitted_kpi(), KPI scale

fitted_kpi() applies the inverse transform draw-wise before summarising. For log-response models the default conversion uses the conditional-mean estimator:

This is the bias-corrected back-transform that accounts for the log-normal variance term. The simpler exp(mu) estimator corresponds to the conditional median on the KPI scale. To request that estimator, pass log_response = "median":

# Default: conditional mean, bias-corrected
kpi_tbl <- fitted_kpi(model)

# Explicit conditional median
kpi_tbl <- fitted_kpi(model, log_response = "median")

The output includes source_response_scale (the model’s response scale), response_scale = "kpi", and conversion_method ("lognormal_mean" or "point_exp") to label the result.

observed(), model scale

observed() returns the observed response on the model scale after unscaling (if scale=TRUE).

observed_kpi(), KPI scale

observed_kpi() returns the observed response on the KPI scale. For log-response models, this applies exp() to the model-scale observed values.

to_kpi_scale() helper

The internal function to_kpi_scale(x, response_scale) implements the conversion:

• If response_scale == "log": returns exp(x).
• Otherwise: returns x unchanged.

This function is used consistently by fitted_kpi(), observed_kpi(), and runner artefact writers.

Runner artefact scale conventions

Runner artefact writers use the response scale metadata to determine which scale to report:

Interaction with scale = TRUE

The scale flag and response scale are orthogonal:

• scale = TRUE divides predictors and response by their standard deviations before Stan fitting. It centres them only when the population formula includes an intercept. No-intercept formulae (0 + or - 1) use zero centres, preserving their through-origin specification, including explicit offsets and group effects. Coefficients and fitted values are back-transformed to the original scale by get_posterior().
• Response scale determines whether the original scale is levels (identity) or logs (log).

Both transformations compose: a log-response model with scale=TRUE first takes the log of the response (via the formula), then standardises the logged values. Post-fit, draws are first unscaled, then (for KPI-scale outputs) exponentiated.

Likelihood scores and intervals

LOO and WAIC use density on the unstandardised model-response scale. For scaled fits the shared extractor subtracts log(s_y) from each Stan pointwise log likelihood. A log-response model is scored per unit of log KPI; there is no further Jacobian to the KPI scale. This follows the change-of-variables rule in the Stan User’s Guide. Compare only the same response transformation and observation set. Correcting standardisation does not make a levels score comparable to a log-response score.

predict() and fitted() summarise conditional-mean draws. They do not add observation noise. In fitted.csv, p2_5, p50 and p97_5 describe that mean; predictive_p2_5, predictive_p50 and predictive_p97_5 describe replicate observations from y_rep, including observation noise. Both use the model response scale. Predictive intervals are unavailable when those draws are absent. MAP reports point estimates with NA interval quantiles.

For a log-response model, the fitted_kpi.csv credible interval summarises exp(mu + sigma^2 / 2) draw by draw. It describes uncertainty in the KPI conditional mean. It is not an interval for a future KPI observation.

Earlier scaled no-intercept fits were centred incorrectly and require refitting. Loading an old object or regenerating its predictions cannot correct its fitted posterior. Other retained full fit objects can have LOO/WAIC scores regenerated with the corrected extractor; old score CSVs are not rewritten automatically.

Jensen’s inequality and draw-wise conversion

When converting log-scale posterior draws to KPI scale, dsambayesplus applies exp() to each draw individually before computing summaries (mean, median, credible intervals). This is the correct Bayesian approach because:

• E[exp(X)] ≠ exp(E[X]) when X has non-zero variance (Jensen’s inequality).
• Draw-wise conversion preserves the full posterior distribution on the KPI scale.
• Summary statistics (mean, quantiles) computed after conversion correctly reflect KPI-scale uncertainty.

Practical guidance

• Use identity-response models when the KPI is naturally additive and coefficients should represent unit changes.
• Use log-response models when the KPI is naturally multiplicative, when variance scales with level, or when the response must remain positive.
• Always check response_scale_label(model) before interpreting coefficient magnitudes.
• Do not call log-response coefficients elasticities unless the predictor is also logged. In log(kpi) ~ x, they are semilog coefficients.
• For KPI-scale percentage interpretation, use 100 * (exp(beta) - 1): not 100 * beta, unless the coefficient is small enough that the approximation is acceptable.
• Use fitted_kpi() for business reporting; use fitted() for diagnostics.
• Do not manually exponentiate posterior means from log-response models. Use fitted_kpi() or to_kpi_scale() on individual draws.

Cross-references

• Model Classes, constructor and formula conventions
• Priors and Boundaries, scale interaction with priors
• Diagnostics Gates, diagnostics computed on model scale
• Config Schema, target.*, media, controls, and model.scale keys

Run and public LOO comparisons require matching response and observation provenance, including density units. The runner leaves ranks and deltas unavailable when provenance is absent or incompatible; loo_compare_models() rejects such comparisons. This guard does not establish the reliability of the score or the validity of the estimator. See Compare Runs.

CRE / Mundlak

Purpose

The correlated random effects (CRE) pathway, implemented as a Mundlak device, augments hierarchical dsambayesplus models with group-mean terms. This separates within-group variation from between-group variation for selected regressors. It addresses correlation between group effects and regressors when the conditional mean of the selected random intercept is adequately represented by the included group means. It does not correct correlated random slopes, time-varying endogeneity, or omitted time-varying confounders.

When to use CRE

Use CRE when:

• The model is hierarchical (panel data with (term | group) syntax).
• Time-varying regressors (e.g. externally transformed media signals) have group-level means that may be correlated with the selected group intercept.
• You want to decompose effects into within-group (temporal) and between-group (cross-sectional) components.

Do not use CRE when:

• The model is BLM or pooled (CRE requires hierarchical class).
• The panel has only one group (no between-group variation exists).
• All regressors of interest are time-invariant (CRE mean terms would be constant).
• The selected group needs random slopes, or a CRE variable uses an internal draw-dependent adstock or Hill transform.

Construction

CRE is applied after model construction via set_cre():

model <- blm(
  kpi ~ m_tv_signal + m_search_signal + trend + (1 | market),
  data = panel_df
)
model <- set_cre(model, vars = c("m_tv_signal", "m_search_signal"))

What set_cre() does

1. Resolves the grouping variable. If the formula has one group factor, it is used automatically. If multiple group factors exist, the group argument must be specified explicitly.

2. Generates group-mean column names. For each variable in vars, a mean-term column is named cre_mean_<variable> (configurable via prefix).

3. Preserves the source data. set_cre() does not add generated columns to .original_data. It creates a provisional prior view for the generated terms; fit preparation resolves the final priors from the retained frame.

4. Updates the formula. The generated mean terms are appended to the population formula as fixed effects.

5. Builds authoritative means at fit time. dsambayesplus reconciles all population, random-effect, group, date, offset, and transform dependencies; excludes incomplete or non-finite rows; then calculates group means over exactly the rows and fitted covariate basis passed to Stan.

6. Extends priors and boundaries. Default prior and boundary entries are added for each new mean term. Explicit prior overrides remain unchanged.

YAML runner configuration

When using the runner, CRE is configured via:

model:
  type: cre

hierarchy:
  group: market
  random_intercept: true
  random_slopes: []
  cre_variables: [m_tv, m_search, m_social]
  cre_prefix: cre_mean_

The runner calls set_cre() during model construction for model.type: cre.

Mundlak decomposition

For a regressor $x_{gt}$ (group $g$, time $t$), dsambayesplus fits the equivalent parameterisation

where:

• Within-group effect: $\beta$, the coefficient on $x_{gt}$ after conditioning on the group mean.
• Contextual difference: $\theta$, the coefficient on $\bar{x}_g$. It is the difference between the between-group and within-group slopes.
• Total between-group slope: $\beta + \theta$.

The original coefficient on $x_{gt}$ in a standard random-effects model conflates both sources. Adding $\bar{x}_g$ as a fixed effect separates them.

Validation and identification warnings

Input validation

set_cre() validates:

• The model is hierarchical (aborts for BLM or pooled).
• The selected group exists and contains a random intercept but no random slope.
• Each CRE variable is one simple numeric population main effect. Factors, interactions, inline transforms, polynomials, and splines involving a CRE variable are rejected.
• No CRE variable overlaps an enabled internal draw-dependent media transform.
• No generated CRE mean term appears in a random-effects block.
• Generated column names do not collide with user-authored source columns.

The same typed structural checks run in direct data preparation, pre-flight, fit(), fit_map(), authored v2 config validation, compiled v1 validation, and runner preparation. pre_flight_strict = FALSE cannot downgrade them.

Identification warnings

warn_cre_identification() checks two conditions after CRE setup:

1. More CRE variables than groups. If length(vars) > n_groups, between-effect estimates may be weakly identified. The function emits a warning.

2. Near-zero within-group variation. For each CRE variable, the within-group residual ($x_{gt} - \bar{x}_g$) standard deviation is checked. If it is effectively zero, within-effect identification is weak. The function emits a per-variable warning.

Zero-variance CRE mean terms

If the combined between-group design is rank-deficient, pre-flight reports the problem and strict runner fitting aborts. Reduce hierarchy.cre_variables, add credible group variation, or use RE only when its independence assumption is defensible. Disabling scaling does not repair a rank-deficient CRE design.

Panel assumptions

• Balanced panels are not required. Each retained observation contributes once to its group mean, so unequal group sizes are supported.
• Missing and non-finite dependencies are reconciled jointly. A row omitted from the likelihood cannot influence a fitted CRE mean. Mean calculation does not use na.rm = TRUE after the retained frame is frozen.
• Generated columns are owned by dsambayesplus. Fit preparation rebuilds them from source variables. A user-authored column with a reserved generated name aborts instead of being overwritten.
• Prediction uses fitted means. Holdout, cross-validation, and deployment paths use training-only lookups derived from the fitted estimation frame.

Estimator evidence smoke

The default suite runs deterministic checks of the repository-owned RE/CRE DGP and evidence aggregation contracts. An explicit Stan smoke checks that one fixed RI-02 replication can be fitted through the public RE and CRE APIs:

Rscript scripts/check.R --estimator-evidence-smoke

This smoke asserts only that usable posterior coefficient draws exist. It does not measure recovery, calibration, bias, relative estimator performance, or production readiness. The immutable manifest, terminal-record, and aggregation workflow under scripts/estimator-evidence/ owns pilot and confirmatory Monte Carlo evidence. Do not treat a smoke result as estimator qualification or causal validation.

Retained diagnostics evidence

Runner fits retain two CRE-specific diagnostic artefacts under 40_diagnostics/:

• estimator_checks.csv records the typed structural-contract result and the pre-flight CRE configuration result;
• row_reconciliation.csv records input, retained, and excluded counts plus excluded source-row IDs and reason codes. It does not contain response values.

These artefacts establish which checks and rows the fit used. They do not establish causal identification or Monte Carlo calibration.

Fixed-effects status

dsambayesplus implements a coefficient-only fixed-effects (FE) estimator through model.type: fe with fit.method: mcmc. It uses an orthonormal Helmert contrast basis and a dedicated artefact contract. Its status is implemented_not_qualified: it passes structural checks and fits, but it has not been qualified through repeated-dataset Monte Carlo evidence.

Decomposition and reporting

CRE mean terms appear as ordinary fixed-effect terms in the population formula. This means:

• Posterior summary includes CRE mean-term coefficients alongside other population coefficients.
• Response decomposition via decomp() attributes fitted-value contributions to CRE mean terms separately from their within-group counterparts.
• Plots (posterior forest, prior-vs-posterior) include CRE mean terms.

Interpretation note: the CRE mean-term coefficient is a contextual difference, not the total between-group slope. For a model containing $x_{gt}$ and $\bar{x}_g$, calculate the total between-group slope as the sum of their two coefficients. Treat all three quantities as associational unless the research design supports a causal interpretation.

Cross-references

• Model Classes, hierarchical class construction
• Estimator Capabilities, current FE, RE, and CRE support boundaries
• Priors and Boundaries, prior schema for CRE terms
• Config Schema, model.type: cre and hierarchy.* YAML keys
• Diagnostics Gates, within-group variation check

Time Components

Purpose

dsambayesplus provides managed time-component generation through the effects.holidays config section. When enabled, the runner deterministically generates holiday feature columns from a calendar file and appends them to the compiled model formula. This page defines the configuration contract, generation logic, naming conventions, and audit properties.

Overview

Time components in dsambayesplus cover:

• Holidays: deterministic weekly indicator features derived from an external calendar file.
• Trend and seasonality: specified directly in the model formula (e.g. t_scaled, sin52_1, cos52_1). These are not generated by the time-components system; they are user-supplied columns in the data.

The managed-effects system is responsible only for holiday feature generation.

Time-index and panel assumptions

dsambayesplus models time using an ordered modelling-wave index, not elapsed calendar time. The standard MMM workflow assumes:

• observations are conditionally independent given the modelled mean, with a common residual standard deviation;
• the standard panel workflow assumes balanced data, so every group has the same ordered modelling waves; low-level hierarchical and CRE constructors do not enforce this, so upstream preparation must verify it before a production panel run;
• trend, seasonality and media carry-over are defined per modelling wave; and
• the date column identifies wave order and supports reporting, blocked cross-validation and holiday alignment.

Under this contract, one step in t_scaled and one adstock recursion step both mean one modelling wave. A longer calendar gap between two observed rows does not rescale either component. Elapsed calendar duration is outside the model contract.

For a panel, construct the wave variables once from the common date index and map them back to every group. Do not create a separate sequential index within each stacked group:

waves <- sort(unique(panel_df$date))
wave_index <- seq_along(waves)
wave_lookup <- match(panel_df$date, waves)

panel_df$t_scaled <- as.numeric(scale(wave_index))[wave_lookup]
panel_df$sin52_1 <- sin(2 * pi * wave_index[wave_lookup] / 52)
panel_df$cos52_1 <- cos(2 * pi * wave_index[wave_lookup] / 52)

YAML configuration

effects:
  holidays:
    enabled: true
    path: ../data/holidays.csv
    date_col: null
    label_col: holiday
    country: gb
    country_col: country
    date_format: null
    week_start: monday
    timezone: UTC
    prefix: holiday_
    window_before: 0
    window_after: 0
    aggregation_rule: count
    overlap_policy: count_all
    overwrite_existing: false

Key definitions

Calendar file contract

The holiday calendar is a CSV (or data frame) with at minimum:

Date column detection

If date_col is null, the system tries column names in order: date, ds, event_date. If none is found, validation aborts.

Label normalisation

Holiday labels are normalised to lowercase, alphanumeric-plus-underscore form via normalise_holiday_label(). For example:

• Black Friday → black_friday
• New Year's Day → new_year_s_day
• Empty labels → unnamed

The generated feature column name is {prefix}{normalised_label}, e.g. holiday_black_friday.

Generation pipeline

The runner calls build_weekly_holiday_features() with the following steps:

1. Parse and validate the calendar. validate_holiday_calendar() checks column presence, date parsing, and label completeness.

2. Expand holiday windows. expand_holiday_windows() replicates each event row across the [event_date - window_before, event_date + window_after] range.

3. Align to weekly index. Each expanded event-day is mapped to its containing week using week_floor_date() with the configured week_start.

4. Aggregate per week. Events are counted per week per feature. Under aggregation_rule: any, counts are collapsed to binary (0/1). Under overlap_policy: dedupe_label_date, duplicate label-date pairs within a week are removed before counting.

5. Join to model data. The generated feature matrix is left-joined to the model data by the date column. Weeks with no events receive zero.

6. Append to formula. Generated feature columns are appended as additive terms to the compiled population formula.

Weekly anchoring

All weekly alignment uses week_floor_date(), which computes the most recent occurrence of week_start on or before each date. The model data’s date column must contain week-start-aligned dates; normalise_weekly_index() validates this and aborts if dates are not aligned.

Supported week-start values

monday, tuesday, wednesday, thursday, friday, saturday, sunday.

Timezone handling

• Calendar dates are parsed using the configured timezone (default UTC).
• If the calendar contains POSIXt values, they are coerced to Date in the configured timezone.
• Character dates are parsed as ISO 8601 by default, or using date_format if specified.

Generated-term audit contract

Generated holiday terms are tracked for downstream diagnostics and reporting:

• The list of generated term names is stored in model$.runner_time_components$generated_terms.
• The identifiability gate in R/diagnostics_report.R uses this list to auto-detect baseline terms (via detect_baseline_terms()), so generated holiday terms are included in baseline-media correlation checks without requiring explicit configuration.

Feature naming collision

If two different holiday labels normalise to the same feature name, build_weekly_holiday_features() aborts with a collision error. Ensure calendar labels are distinct after normalisation.

Interaction with existing data columns

• If overwrite_existing: false (default), the runner aborts if any generated column name already exists in the data.
• If overwrite_existing: true, existing columns with matching names are replaced by the generated features.

Practical guidance

• Start with aggregation_rule: count to capture multi-day holiday effects (e.g. a holiday spanning two days in one week produces a count of 2).
• Use window_before and window_after for events with known anticipation or lingering effects (e.g. window_before: 7 for pre-Christmas shopping).
• Use aggregation_rule: any when you want binary holiday indicators regardless of how many event-days fall in a week.
• Check generated terms in the resolved config (config.resolved.yaml) and posterior summary to confirm which holidays entered the model.

Cross-references

• Config Schema, effects.holidays.* YAML keys
• Diagnostics Gates, identifiability gate uses generated terms
• Priors and Boundaries, priors for generated holiday terms follow default schema
• Output Artefacts, resolved config records generated terms

Diagnostics Gates

For the workflow interpretation of these checks, start with Stage 4: Computation and Sampler and Stage 5: Model Adequacy. This page is the threshold and policy reference.

Use this page when you need exact gate thresholds, status aggregation, or YAML policy semantics. It does not replace substantive model review: a model can clear threshold tables and still be a poor basis for interpretation.

Geographic FE/RE/CRE models also have separate pre-fit panel design diagnostics, including same-date geographic contrasts and conditional design support.

Editable gate profiles

Edit config/diagnostic_gates.yml to define team reporting thresholds. The shipped BLM, RE and CRE examples reference this complete profile:

diagnostics:
  enabled: true
  gates_file: diagnostic_gates.yml

Relative paths resolve from the model YAML directory, including when runme.R applies CLI data or holiday overrides. The file is validated before compilation, including --mode=validate. Each profile requires schema_version: 1, a non-empty profile name and description, and the complete thresholds mapping. For example, one entry in that mapping is:

thresholds:
  design_max_vif: {warn: 10, fail: 30}

The profile controls reporting, not the likelihood, priors, numerical design checks or estimator support. Strict > and < comparisons are used; equality does not breach a threshold, and failure takes precedence. Equal warn/fail thresholds have no warning band. fail: null makes a finite result warning-only. Non-finite evidence for an evaluated check cannot pass; checks that do not apply remain skipped with reasons. MCMC gates are skipped for MAP. Empty or entirely skipped evidence is not PASS and has unknown publishability.

Unknown or duplicate keys, incomplete mappings, non-finite thresholds, invalid domains and reversed warn/fail ordering are errors. warn is required and numeric; fail must be explicitly numeric or null. Do not combine a gate file with inline pareto_k_warn/pareto_k_fail or abs_corr_warn/abs_corr_fail.

An explicit file controls all configurable gate classifications, including failures in publish mode. Without a file, the package resolves the existing mode and inline thresholds into a named default profile. policy_mode still controls existing structural/preflight policies that are outside this file. A gate file cannot relax estimator-specific rank or support requirements.

The runner freezes the resolved profile before fitting and records its source, content hash and hash algorithm in 40_diagnostics/diagnostic_gates.resolved.yaml. The hash uses rlang::hash; it is labelled rlang_hash, not SHA-256. Editing a profile during sampling does not change the active run. Version-control the profile with the model configuration; change its name when issuing a new team policy. The resolved snapshot contains provenance fields and is an output artifact; use the four authored fields above when creating a new input profile.

VIF interpretation

VIF is a collinearity diagnostic, not a hypothesis test. The default warning above 10 and failure above 30 reuse the existing preflight thresholds. The report records the maximum and limiting term; diagnostic_vif.csv retains per-term values, statuses, limits and design scope. vif_bar.png uses the same resolved rule. Fewer than two population predictors is not applicable; singular or constant predictors cannot produce a passing VIF gate.

This population-design VIF is distinct from CRE between-design VIF in cre_estimability_vif.csv. Those estimator-specific diagnostics and safeguards remain unchanged. The standard diagnostic-report API remains unsupported for FE; FE runs retain their existing dedicated estimator artifacts. FE configs reject diagnostics.gates_file and do not produce a standard gate snapshot. No generic FE gate packet or new estimator qualification is implied.

Terminal output

Rscript runme.R and the package CLI print grouped diagnostic results, observed values, threshold rules, limiting terms/parameters/chains where available, skipped reasons and recovery guidance. The report is derived from the same structured checks saved to CSV. Colours supplement explicit PASS/WARN/FAIL/ SKIPPED words; plain and redirected output remain readable. A diagnostic FAIL is distinct from a failed fit: retain and inspect the fitted model and reports. Existing CLI codes remain 0 for completion without diagnostic failure, 1 for completed runs needing failure review, and 2 for execution/configuration errors.

Model selection for time-ordered data

PSIS-LOO and WAIC treat pointwise observations as conditionally exchangeable. That assumption is not generally appropriate for time-ordered MMM data, where nearby weeks can remain dependent after conditioning on the fitted model.

Use the runner’s expanding-window blocked CV or leave-future-out CV as the primary evidence when selecting among time-series MMM specifications. Treat PSIS-LOO, WAIC, and Pareto-k outputs as supplementary fit and influence diagnostics. They do not establish future-period predictive performance, causal validity, or a publish-gate pass on their own.

Purpose

dsambayesplus runs a deterministic diagnostics framework after model fitting. Each diagnostic check produces a pass, warn, or fail status. The resolved gate profile controls the configurable thresholds. This page defines the check taxonomy, threshold tables, policy modes, identifiability gate, and the overall status aggregation rule.

How to use this page

• Use Stage 4: Computation and Sampler to understand which checks are non-negotiable before trusting the posterior.
• Use this page to see the exact dsambayesplus thresholds and artifact semantics.
• Use Stage 5: Model Adequacy before treating a passing diagnostics table as permission for decomposition, model comparison, or optimisation.

Policy modes

Without an explicit gate file, the diagnostics framework resolves these preset modes from diagnostics.policy_mode:

The mode is resolved by diagnostics_policy_thresholds(mode) in R/diagnostics_thresholds.R.

Check taxonomy

Checks are organised into phases:

Each check row includes:

Design checks

Condition number thresholds by mode

P0 sampler checks (MCMC only)

Mode adjustments for sampler checks

dsambayesplus treats any Rhat above 1.01 as a failure in publish and strict modes, following the rank-normalised, folded Rhat guidance in Vehtari et al. (2021) and the Stan warnings guide. In explore mode, the fail threshold is deliberately relaxed to 1.10, while the warning threshold remains 1.01.

P1 residual checks

Mode adjustments for residual checks

In publish mode, values beyond the raw residual fail boundary remain warn. The persisted threshold field states that failure is disabled and retains the raw boundary for triage. Non-residual P0 failures can still make the overall run fail.

P1 boundary hit check

In explore mode, boundary hits cannot fail. Publish mode retains the raw fail > 0.20 boundary for triage but records values beyond it as warn. Failure is disabled for this check in publish mode. In strict mode, thresholds tighten to warn > 0.02, fail > 0.10.

P1 within-group variation check

This check applies to hierarchical models and flags groups where within-group variation is extremely low relative to between-group variation. In explore mode, the fail threshold is zero (cannot fail). Publish mode retains the raw fail < 0.05 boundary for triage but records values below it as warn; failure is disabled for this check. Strict mode applies warn < 0.15 and fail < 0.10.

Identifiability gate

The P1 pre_identifiability_baseline_media_corr check measures the maximum absolute correlation between baseline terms and media terms in the design matrix. It is configured via diagnostics.identifiability in YAML:

diagnostics:
  identifiability:
    enabled: true
    media_terms: [m_tv, m_search, m_social]
    baseline_terms: [trend, seasonality]
    baseline_regex: ["^h_", "^sin", "^cos"]
    abs_corr_warn: 0.80
    abs_corr_fail: 0.95

Term detection

• Media terms: explicitly listed in media_terms.
• Baseline terms: union of baseline_terms, generated time-component terms, and matches from baseline_regex patterns.
• Both sets are intersected with actual design-matrix columns and filtered to remove constant columns.

Thresholds by mode

Skip conditions

The identifiability gate reports skipped when:

• identifiability.enabled: false
• No configured media terms found in the design matrix
• No baseline terms detected from configured terms/regex
• All resolved baseline or media terms are constant

Overall status aggregation

The overall diagnostics status is determined by diagnostics_overall_status():

1. If any check has status == "fail" → overall status is fail.
2. If any check has status == "warn" (and none fail) → overall status is warn.
3. Otherwise → overall status is pass.

Checks with status == "skipped" do not affect the overall status.

P2 rows use the emitted post_model_selection_method, post_psis_loo, post_psis_loo_options_ignored, post_psis_loo_time_dependence, post_psis_loo_elpd, post_psis_loo_looic and post_psis_loo_pareto_k identifiers. They record supplementary predictive evidence and do not override P0 or P1 failures.

The runner also emits post_sampler_params_available when sampler parameters can be inspected and post_sampler_unavailable_for_optimise for the explicit MAP not-applicable path. These identifiers are part of the persisted report; they are not substitute MCMC diagnostics for a MAP fit.

Runner artefact output

The diagnostics framework produces:

Interpretation guidance

• skipped: no evaluated evidence; this is not a gate pass.
• pass: no enabled check breached its configured threshold. This does not establish substantive adequacy, causal validity or suitability for the intended decision.
• warn: review recommended; the model may have quality concerns but does not block the configured policy.
• fail: at least one enabled check breached a fail threshold; resolve it before using the result for decisions.

Common remediation actions

Cross-references

• Model Classes, which diagnostics apply to each class
• Response Scale Semantics, residuals are computed on model scale
• Config Schema, diagnostics.* YAML keys
• Output Artefacts, diagnostics artefact paths
• Diagnostics Plots, visual diagnostic outputs

Estimator Capabilities

Current implementation status

This document is the approved Gate A contract for hardening the random-effects (RE), correlated-random-effects (CRE), and fixed-effects (FE) estimators. Approval fixes the intended scope. It does not mean that every target capability is implemented or production qualified.

The current implementation status is authoritative until the relevant implementation and evidence gates have passed:

No estimator is production qualified by this document alone. Production status requires the implementation, package validation, Stan review where applicable, and predeclared Monte Carlo evidence.

Shared interpretation boundary

Estimator choice changes how time-invariant group heterogeneity is handled. It does not establish causal identification. RE, CRE, and FE can all be biased by omitted time-varying confounding, reverse causality, measurement error, or a misspecified response function. Applied work must assess those risks separately.

The package will fail closed for a structurally unsupported estimator and design combination. Advisory diagnostics will remain distinct from structural errors and will report their metric, threshold, and recovery action.

Estimator design-diagnostic policy

diagnostics_policy_thresholds() owns one nested estimator_design policy. This keeps estimator-specific design checks separate from the existing general design-matrix and fitted-model diagnostics. In particular, the estimator condition-number rule does not replace the existing kappa_* policy, and the within-share rule does not replace the existing within_var_* policy.

The version 1 estimator-design policy fixes these boundaries:

The group-level boundary is an operational small-support warning, not a claim that four groups guarantee reliable variance-component estimation. The VIF, condition-number, within-share, residual-degrees-of-freedom, and group-support boundaries are conservative screening rules rather than universal statistical constants.

For YAML runner fits, diagnostics.policy_mode is retained before fitting and selects this estimator-design policy during CRE pre-flight and preparation. Consequently, strict mode rejects a full-rank combined between-design with fewer than two residual degrees of freedom before Stan compilation. Direct API preparation uses the publish policy.

P2-A defines and tests this policy. P2-B through P2-E implement the numerical reports, apply the policy to the exact retained frame, and retain the resulting evidence. The Stan model and estimation algebra are unchanged. Structurally invalid retained designs now fail before fitting; finite support limitations remain advisory. These diagnostics do not qualify an estimator for production use.

Capability matrix

The following matrix records the implemented FE v2 boundary and the approved RE and CRE scope. Implementation does not imply production qualification.

Random effects

Target estimator

RE remains the existing Gaussian hierarchical estimator. Group coefficients are partially pooled through the declared random-effects blocks. This can be a useful variance model when its conditional independence assumptions are credible.

Retained support diagnostics

dsambayesplus retains diagnostics for each random-effects block:

• group count and group-size distribution;
• random-effects design dimensions and rank by group;
• all-zero and non-finite random-effect columns;
• within-group variation for random slopes; and
• covariance dimension relative to available groups.

All-zero columns, non-finite matrices, fewer than two group levels, and invalid dimensions are structural errors. The response is also rejected as a random-effect design term or grouping key. Weak within variation, per-group rank loss, and small group counts are advisory. Covariance dimension and the groups-per-covariance-parameter ratio are reported without an identification threshold. A group-constant random slope is not automatically described as mathematically unidentified when between-group information may still inform its covariance.

Assumption boundary

Plain RE does not protect population slopes when the random intercept is correlated with included regressors. In that case, use CRE only if its narrower conditional-mean specification is credible, or use another design.

Correlated random effects

Target estimator

CRE v1 is a random-intercept Mundlak specification. For each selected numeric regressor, dsambayesplus adds its mean over the exact retained estimation rows for the selected group. The coefficient on the original regressor is the within slope. The coefficient on its group mean is the contextual difference. Their sum is the between slope.

Supported scope

CRE v1 requires:

• one explicitly selected hierarchy group;
• a random intercept and no random slope in that selected group;
• a CRE variable that maps to one simple numeric population-design column;
• no overlap between that variable and an enabled internal media transform;
• group means computed from the same retained rows and fitted covariate basis used by Stan; and
• holdout, cross-validation, and deployment lookups derived from the fitted training frame only.

Other hierarchy blocks may retain their ordinary RE semantics. Repeated selected-group/date pairs are allowed when another hierarchy level provides multiple observations in the same period.

Rejected specifications and recovery

Retained-sample rule

set_cre() will become configuration-only. At fit time, dsambayesplus will build one explicit completeness and finiteness mask across the response, population terms, random-effects terms, grouping keys, date, offset, and transformation inputs. CRE means, default priors, design diagnostics, Stan data, and retained lookups will all derive from that frozen frame. na.rm = TRUE will not define the fitted group mean.

CRE does not correct omitted time-varying confounding or guarantee that the conditional mean of the group effect is linear in the included means.

Fixed effects

Implemented estimator

FE v2 is a separate fixed_effects model class constructed with fixed_effects(formula, data, unit). Before fitting, the caller must supply exactly one date mapping through the existing set_date() lifecycle.

The non-dry YAML path uses the same estimator, prior, boundary, scaling, and sampling contracts as the direct R API. It bypasses generic level-scale diagnostics, PSIS model selection, scenario analysis, optimisation, forecasting, and deployment. It can opt into seen-unit expanding-window validation through diagnostics.time_series_selection; each fold must retain complete rows, at least two training observations per unit, and only holdout units already seen in that fold’s training window. Its optional ELPD-by-fold PNG requires both diagnostics.time_series_selection.save_png and outputs.save_diagnostics_png; the latter does not enable generic FE diagnostic plots. A completed runner result has no publishability claim, and its factual FE diagnostics summary records qualification_status: not_assessed.

For each retained unit, FE builds a deterministic orthonormal Helmert contrast matrix Q_i and fits:

Q_i y_i ~ Normal(Q_i X_i beta, sigma)

This removes unrestricted unit intercepts without estimating shrinkage priors for them. Unequal unit sizes are allowed. Singleton units, duplicate retained unit-date keys, non-finite values, time-invariant regressors, and rank-deficient within designs will abort before sampling.

Version 2 output contract

FE v2 supports slope and noise inference, within-design diagnostics, contrast-space posterior predictive checks, and conditional retained-unit level output. For each retained unit and posterior draw it reconstructs:

alpha_i | beta, sigma, y_i ~ Normal(mean(y_i - X_i beta), sigma / sqrt(T_i))
mu_i,t = alpha_i + x_i,t beta

fixed_effects_intercepts() returns the conditional intercept draws. fitted() and observed() return retained level data keyed by unit and date. predict(newdata) returns level means only for units retained at fit time and preserves newdata order. The optional seed controls the conditional intercept reconstruction. The strongest possible release label remains production for coefficient inference only.

FE v2 rejects:

• unseen-unit prediction;
• observation-noise simulation and generic residual diagnostics;
• decomposition and counterfactual response;
• budget optimisation;
• generic time-series cross-validation; the runner supports only the documented seen-unit expanding-window validation path;
• pointwise LOO, WAIC, and run comparison;
• MAP fitting; and
• internal media transformations.

Contrast-space fitted values remain diagnostics. They are distinct from the conditional level means returned by the direct API. The YAML runner does not yet write generic level-scale output files. Repeated-sampling coverage and a counterfactual specification remain separate work.

The P4-G bounded synthetic integration fit proves only that the direct API can compile, sample, extract its declared output, and survive a same-environment RDS round-trip. P4-I proves the bounded runner orchestration and artefact contract with a mocked Stan boundary. P4-J adds bounded live runner and explicit-prior prior-only checks plus balanced and unbalanced coefficient, noise, and sampler regression gates. All passed their predeclared technical checks.

The two P4-J recovery fits are deterministic one-replication regression evidence. They do not estimate repeated-sampling coverage or production reliability. They also do not estimate bias, root mean squared error, failure rates, or robustness. FE therefore remains implemented_not_qualified. Repeated Monte Carlo evidence and its separate review are post-release work.

Qualification requirements

Each estimator receives its own release decision. A package-wide pass cannot hide a failed estimator cell.

Qualification requires:

1. deterministic design and algebra oracles;
2. typed structural failures before Stan compilation;
3. four-chain sampler diagnostics on every declared estimand and likelihood scale parameter;
4. repeated paired Monte Carlo evidence under RE-valid, RE-invalid/CRE-valid, unbalanced, weak-design, and misspecified scenarios;
5. conditional coverage and unconditional success-and-cover rates, so failed fits remain visible;
6. Monte Carlo uncertainty for bias, coverage, and paired loss comparisons;
7. bounded prior-sensitivity evidence; and
8. human statistical and Stan review where model code changes.

Possible release labels are production, experimental, and unsupported. Any failed core gate blocks production. An inconclusive stress cell narrows the supported capability and remains visible in the evidence report.

Related pages

• Model Classes
• CRE / Mundlak
• Diagnostics Gates
• Model Object Lifecycle

Counterfactual Response

Purpose

counterfactual_response() evaluates one fitted MCMC model under two aligned predictor paths. It returns posterior draws of the fitted response location under the proposed scenario, the reference response, and their difference.

The result is a model-implied scenario contrast. It is not, by itself, a causal effect or an attribution estimate. Causal interpretation requires a credible identification argument for the fitted model and the intervention being represented.

Usage skeleton

The following code assumes that fitted_model is a fitted MCMC object and future_data contains every required model variable and modelling wave.

scenario <- future_data
scenario$paid_search <- scenario$paid_search * 1.10

reference <- future_data

contrast <- counterfactual_response(
  fitted_model,
  scenario = scenario,
  reference = reference,
  scale = "kpi",
  interval = 0.9
)

contrast$summary
contrast$draws

Both data frames must resolve to the same retained rows. Dates, hierarchical group keys and row order must match. Predictor values may differ. The reported difference is scenario - reference. For balanced panels, supply the common ordered modelling waves used by the fitted model.

Returned quantities

The object has three components:

These are posterior draws of the fitted response location. They do not include new observation noise. On KPI scale, a log-response model can return a conditional median or a lognormal conditional mean. Use the draws to describe uncertainty in the fitted response surface, not the full range of future observations.

Use draw_ids to select existing .sample identifiers when the complete draw-by-row result would be too large:

contrast <- counterfactual_response(
  fitted_model,
  scenario,
  reference,
  draw_ids = seq(1, 2000, by = 4)
)

draw_ids reduces the returned draw-by-row object. The current implementation still materialises the fitted posterior before selecting those identifiers, so it does not remove the peak memory cost of posterior extraction.

Aggregate a decision horizon

Use aggregate_counterfactual_response() to sum responses within each posterior draw before calculating intervals:

horizon <- aggregate_counterfactual_response(contrast)

horizon$summary
horizon$draws

The result includes draw-wise scenario totals, reference totals, differences and percentage lift. Its summary includes posterior means, medians, central intervals and probability_difference_positive. That probability is the proportion of retained posterior draws where the aggregated scenario - reference difference is strictly above zero. It is not the probability that an intervention has a causal effect.

Group by a retained panel or date label when separate aggregates are needed:

by_market <- aggregate_counterfactual_response(contrast, by = "market")
by_wave <- aggregate_counterfactual_response(contrast, by = "date")

Aggregation is draw-wise. dsambayesplus first sums each draw over the relevant rows and then calculates posterior intervals. It does not sum row-level interval bounds.

Percentage lift is defined draw by draw as:

100 * difference_total / reference_total

It requires totals on the KPI or identity-response scale and a strictly positive reference total in every retained draw and group. For log-response totals, first call counterfactual_response(..., scale = "kpi"). Set include_percent_lift = FALSE when only absolute totals and differences are required.

YAML runner integration

Set scenario_analysis.enabled: true to evaluate two file-backed paths after a successful MCMC fit. The runner writes row summaries, aggregate summaries and metadata under 60_scenario_analysis/. Draw-level aggregate output is written only when scenario_analysis.save_draws: true.

scenario_analysis:
  enabled: true
  scenario_path: ../data/scenarios/planned.csv
  reference_path: ../data/scenarios/reference.csv
  scale: kpi
  log_response: mean
  interval: 0.9
  aggregate_by: []
  include_percent_lift: true
  save_draws: false

The scenario and reference files must follow the same alignment and input contracts described above. A failed requested scenario analysis does not erase the successful model fit. The returned runner result instead records outcome = "completed_with_scenario_analysis_fail" and the failure reason. See Config Schema and Output Artefacts.

Response and KPI scales

scale = "response" returns the model response scale. For log(y) models, this is the log scale.

scale = "kpi" converts both absolute paths before calculating their difference. For a log-response model:

• log_response = "median" uses exp(eta);
• log_response = "mean" uses the draw-specific lognormal adjustment exp(eta + noise_sd^2 / 2).

The second form estimates a conditional expected KPI level under the fitted Gaussian log-response model. It does not correct model misspecification or confer a causal interpretation.

Media transforms and carry-over

For a model fitted with probabilistic adstock and Hill saturation, the engine uses each posterior draw’s fitted decay, half-saturation and media coefficient. It uses the fixed Hill shape and scaling metadata stored with the fit. The full Stan posterior must therefore remain available. A materialised posterior table or compact deployment artefact does not contain enough information and is rejected.

Carry-over starts at the first supplied row. dsambayesplus does not invent media history before the scenario window. Hierarchical media transforms reset at the fitted panel boundary. The current pooled Stan model uses one media panel, so its carry-over continues across the complete supplied row sequence. Arrange pooled inputs to match the fitted row contract and review this limitation before using pooled transformed-media contrasts.

If the scenario continues from a known observed media history, prepend that common history to both scenario and reference. Keep the historical values identical in both paths. This reconstructs carry-over at the decision-horizon boundary. After evaluation, use .input_row or the date and panel keys to keep only the decision-horizon rows. Omitting known pre-horizon media resets carry-over and answers a different scenario question.

Changing spend early in the supplied window can affect later rows through carry-over. Scenario and reference paths should therefore cover the complete decision horizon rather than isolated rows.

Hierarchical paths are evaluated in the model’s group-major order. The summary$.input_row column maps every processed row back to its supplied data row. Draw-vector positions follow summary$.row; do not bind them back to the input data by position alone.

Supported and rejected fits

The first implementation supports fitted MCMC BLM, hierarchical and pooled models. It rejects:

• unfitted models;
• MAP fits, because one optimum is not a posterior sample;
• deployment artefacts, because they contain point estimates;
• transformed-media objects that no longer retain Stan alpha and k draws;
• unseen hierarchical levels;
• duplicate date/panel row keys;
• transformed paths that are not strictly ordered by modelling wave;
• transformed hierarchical paths whose panels do not contain identical wave sets;
• scenario and reference paths with different retained dates, panels or row counts.

optimise_budget() remains a separate scenario-authored response-surface workflow. It does not yet call this engine and does not automatically inherit fitted adstock or Hill parameters.

Budget Optimisation

Purpose

dsambayesplus provides a decision-layer budget optimisation engine that operates on fitted model posteriors. Given a channel scenario with spend bounds, response-transform specifications, and an objective function, the engine searches for the allocation that maximises the chosen objective while respecting channel-level constraints. This page defines the inputs, objectives, risk scoring, response-scale handling, and output structure.

Overview

Budget optimisation is separate from parameter estimation. It takes a fitted model and a scenario specification, then:

1. Extracts posterior coefficient draws for the scenario’s channel terms.
2. Generates feasible candidate allocations within channel bounds that sum to the total budget.
3. Evaluates each candidate across all posterior draws to obtain a distribution of KPI outcomes.
4. Ranks candidates by the configured objective and risk scoring function.
5. Returns the best allocation, channel-level summaries, response curves, and impact breakdowns.

Response-surface contract

The allocator’s response curves are scenario-authored: each channel’s response specification supplies the identity, atan, log1p, or Hill curve used to score allocations. optimise_budget() records this as response_surface = "scenario_authored" in the returned object and its summary artefact.

For models fitted with probabilistic adstock/Hill media transforms, these decision-layer curves are not the fitted Stan response. They do not reuse posterior adstock decay, posterior Hill half-saturation, historical pacing, or carry-over state. Treat them as an explicit scenario model, not as a model-sourced marginal-response estimate.

Entry point

result <- optimise_budget(model, scenario, n_candidates = 2000L, seed = 123L)

The optimize_budget() alias is also available for American English convention.

Scenario specification

The scenario is a structured list with the following top-level keys:

channels

A list of channel definitions, each containing:

Channel names and terms must be unique across the scenario.

budget_total

Total budget to allocate across all channels. All feasible allocations sum to this value.

reference_spend

Optional named list of per-channel reference spend values. If not provided, reference spend is estimated from the mean of the spend_col in the model’s original data.

objective

Defines the optimisation target and risk scoring:

Response transforms

Each channel can specify a response transform that maps raw spend to the transformed value used in the linear predictor. Supported types:

The response transform is applied within response_transform_value() and determines the shape of the channel’s response curve.

Objective functions

kpi_uplift

Maximises the expected change in KPI relative to the reference allocation. The metric for each candidate is:

evaluated across posterior draws $d$.

profit

Maximises expected profit, defined as:

where $\Delta\text{spend} = \text{candidate total} - \text{reference total}$.

Risk-aware scoring

The risk scoring function determines how the distribution of objective draws is summarised into a single score for ranking candidates:

Coefficient extraction

BLM and pooled models

Coefficient draws are extracted via get_posterior() and indexed by the scenario’s channel terms.

Hierarchical models

For hierarchical MCMC models, the population-level (fixed-effect) beta draws are extracted directly from the Stan posterior. If the model was fitted with scale=TRUE, draws are back-transformed to the original scale before optimisation. This ensures that optimisation operates on the population effect rather than group-specific random-effect totals.

Draw thinning

If max_draws is specified, a random subsample of posterior draws is used for computational efficiency. The subsampling uses the configured seed for reproducibility.

Work-size controls and practical limits

n_candidates must be a finite integer from 10 through .Machine$integer.max. scenario$posterior$draws must be a finite integer from 1 through .Machine$integer.max; it defaults to 500. Integer-valued doubles such as 2000 are accepted and normalised. Fractional, missing, non-finite and overflowing values are rejected before candidate or draw allocation begins.

Runtime grows approximately with candidate allocations multiplied by retained posterior draws. An efficient frontier repeats that work for every feasible budget multiplier. The optimiser keeps scalar search results for each candidate and retains full draw vectors only for the winning candidate, so candidate count no longer multiplies retained draw storage.

The defaults of 2,000 candidates and 500 posterior draws are the routine starting point. Benchmark the deployment host before exceeding any of the following review thresholds:

• 10,000 candidates;
• 2,000 retained posterior draws; or
• seven efficient-frontier levels.

These are measurement thresholds, not hard caps. Channel count, response transforms and host hardware also affect runtime. Run the bounded benchmark from the repository root:

source scripts/r-library-path.sh
dsambayesplus_set_r_library host
Rscript scripts/benchmark_budget_optimiser.R \
  --profile=review \
  --repetitions=3

The benchmark runs each case in a fresh process and writes results under the ignored results/ tree. On Linux it reports whole-process peak resident set size (RSS) from /usr/bin/time -v. Peak RSS includes package-load overhead, so compare runs on the same host and software environment; do not interpret it as an optimiser-only heap estimate.

Response-scale handling

Budget optimisation handles both identity and log response scales:

• Identity response: $\Delta\text{KPI}$ is the difference in linear-predictor draws between candidate and reference allocations.
• Log response: $\Delta\text{KPI}$ is computed via kpi_delta_from_link_levels(), which correctly accounts for the exponential back-transformation. If kpi_baseline is available, the delta is expressed in absolute KPI units; otherwise, it is expressed as a relative change.

The delta_kpi_from_link() and kpi_delta_from_link_levels() functions ensure Jensen-safe conversions by operating draw-wise.

Feasible allocation generation

sample_feasible_allocation() generates random allocations that:

1. Respect per-channel lower bounds.
2. Respect per-channel upper bounds.
3. Sum exactly to budget_total.

Allocation is performed by distributing remaining budget (after lower bounds) using exponential random weights, iteratively filling channels until the budget is exhausted. project_to_budget() ensures exact budget equality via proportional adjustment.

Output structure

optimise_budget() returns a budget_optimisation object containing:

Runner integration

When allocation.enabled: true in YAML, the runner calls optimise_budget() after fitting and writes artefacts under 80_optimisation/:

Constraints and guardrails

• Budget feasibility: if channel lower bounds sum to more than budget_total, the engine aborts.
• Upper bound capacity: if channel upper bounds cannot accommodate the full budget, the engine aborts.
• Missing terms: if a scenario term is not found in the posterior coefficients, the engine aborts with a descriptive error.

Cross-references

• Model Classes, fit support and posterior extraction per class
• Response Scale Semantics, KPI-scale conversion for log models
• Priors and Boundaries, boundary interaction with optimisation constraints
• Config Schema, allocation.* YAML keys
• Optimisation Plots, visual outputs from budget optimisation