Categorical Predictors
A categorical predictor is a column whose values are labels — "north", "variant_b", "opted_in" — rather than magnitudes. pathmc treats such columns as first-class terms in the DSL: it expands each one into indicator columns with treatment coding, gives every represented level a named coefficient, and freezes the level set at build time so that predict() and do() can re-encode labels safely. This page is the reference for that behaviour; the Categorical Predictors example is the worked tutorial.
When a column is treated as categorical
A plain term such as y ~ region is inspected against the data when model() or simulate() is called:
| Column dtype | Treated as | Notes |
|---|---|---|
pandas Categorical |
categorical | declared category order is respected; unobserved categories are dropped |
pandas StringDtype |
categorical | |
object holding strings or bytes |
categorical | values are inspected, not just the dtype |
object holding numbers |
continuous | e.g. integers stored as Python objects |
| integer or float | continuous | wrap in C() to declare labels |
| polars / other narwhals frames | as above | converted to pandas for inspection |
Wrapping a term in C() forces categorical treatment regardless of dtype. That is the only way to treat an integer-coded column (store_type as 1, 2, 3) as a set of labels.
Categorical resolution is part of building the model, not of parsing the spec. A data-free model (data=None) parses C(region) but has no levels; equations() shows levels=[]. Everything downstream — column names, coefficient shapes, valid intervention labels — is fixed by the frame you pass to model().
The C() syntax
C(variable, reference=<level>, prior='independent' | 'hierarchical')
| Argument | Default | Meaning |
|---|---|---|
| variable | required | a plain column name |
reference |
first level (see below) | the baseline level absorbed into the intercept; must be an observed level |
prior |
'independent' |
'hierarchical' places a shared population distribution over the level effects |
Level literals are parsed with Python semantics, so reference='south' gives a string and reference=2 gives an integer that must match the column’s values. Unknown arguments (for example contrast=) raise a ParseError; only treatment and cell-means coding are available.
C() cannot take a coefficient prefix (b*C(region) or 2*C(region)), because it expands to several coefficients rather than one. Refer to its coefficients by the vector name described below.
Coding schemes
Which indicator columns are built depends on whether the equation has an intercept.
Treatment coding (intercept present)
One level is the reference. Every other level receives an indicator column named variable[T.level], and its coefficient is the difference from the reference:
region value |
region[T.south] |
region[T.west] |
\mu (other terms omitted) |
|---|---|---|---|
| north (reference) | 0 | 0 | \beta_0 |
| south | 1 | 0 | \beta_0 + \beta_{\text{south}} |
| west | 0 | 1 | \beta_0 + \beta_{\text{west}} |
Cell-means coding (0 +)
With the intercept suppressed there is nothing to absorb a baseline, so every level gets a column named variable[level] and its coefficient is the expected outcome for that level (at zero for any continuous predictors):
region value |
region[north] |
region[south] |
region[west] |
|---|---|---|---|
| north | 1 | 0 | 0 |
| south | 0 | 1 | 0 |
| west | 0 | 0 | 1 |
Both schemes are full rank. The reference level is a reparametrization: pairwise contrasts such as \beta_{\text{west}} - \beta_{\text{south}} are identical whichever level is the reference, so choose whatever makes the coefficients easiest to read.
One baseline per equation
The design stays full rank because exactly one term in each equation supplies the baseline of the linear predictor. With an intercept, that is the intercept, and every independent-prior categorical is treatment coded. With 0 +, the first categorical is cell-means coded and any later categorical is treatment coded against it, so y ~ 0 + region + ch builds region[north], region[south], region[west], ch[T.b] rather than five collinear columns. A reference= argument on a cell-means term has nothing to select and raises a warning.
A hierarchical term also keeps one coefficient per level. The term that supplies the baseline is cell-means coded; a later hierarchical term keeps those coefficients but constrains them to sum to zero, so the design stays full rank. See Hierarchical.
Level order and the default reference
| Column type | Level order | Default reference |
|---|---|---|
strings / StringDtype / object |
sorted by string value | first sorted level |
pandas Categorical |
declared category order, restricted to observed categories | first observed category in that order |
integers via C() |
sorted by string rendering of the value | first in that order |
Distinct values that render to the same string (the integer 1 and the string "1") are rejected, because their coefficients would collide. Missing values are rejected: drop or impute before fitting.
Because C() on an integer column sorts by the string rendering, codes 1..12 are ordered 1, 10, 11, 12, 2, 3, ..., 9 and the reference is 1. This is harmless for fitting but matters whenever you write a per-level vector (priors, simulate() parameters): read the order from model.pymc_model.coords or declare the column as a pandas Categorical with the order you want.
Parameter names and coordinates
Categorical coefficients never enter the equation’s flat beta_{lhs} vector. They live in their own vector whose entries are indexed by level name:
| Object | Name | Dimension / shape |
|---|---|---|
| coefficient vector | beta_{lhs}_{variable} |
({lhs}_{variable}_levels,) |
| level coordinate | {lhs}_{variable}_levels |
non-reference levels (treatment) or all levels (cell means and hierarchical), as strings |
indicator matrix (internal pm.Data) |
_cat_{lhs}_{variable} |
(n_rows, n_columns) |
| hierarchical population mean | mu_beta_{lhs}_{variable} |
scalar; only when the term supplies the equation’s baseline (see Hierarchical) |
| hierarchical population scale | sigma_beta_{lhs}_{variable} |
scalar |
For y ~ x + region with reference north:
model.pymc_model.coords
# {'y_predictors': ['Intercept', 'x'], 'y_region_levels': ['south', 'west']}
idata.posterior["beta_y_region"].sel(y_region_levels="south")The same variable can appear in several equations; each equation owns its own vector (beta_m_region, beta_y_region) with its own coordinate.
Priors
Independent (default)
Every categorical coefficient defaults to Normal(0, 10), matching ordinary slopes. Override by name. A scalar prior is broadcast across levels; to give levels different priors pass vectors and name the level dimension, in coordinate order (here south, then west):
from pathmc import Prior
model.set_priors({"beta_y_region": Prior("Normal", mu=0, sigma=2)})
model.set_priors({
"beta_y_region": Prior(
"Normal", mu=[1.0, -0.5], sigma=[1.0, 1.0], dims=("y_region_levels",)
)
})Hierarchical
C(variable, prior='hierarchical') replaces the independent priors with a partially pooled structure over every level, with no reference:
\beta_{\ell} \sim \mathcal{N}(\mu_\beta, \sigma_\beta), \qquad \mu_\beta \sim \mathcal{N}(0, 10), \qquad \sigma_\beta \sim \text{HalfNormal}(1).
Levels with little data are shrunk towards \mu_\beta by an amount the data determines through \sigma_\beta. This is the recommended choice for many or sparsely observed levels; the example shows the effect on twelve stores with 3–80 observations each.
Pooling all levels is what keeps the result independent of the reference: if one level were absorbed into the intercept it would escape shrinkage, and changing the reference would change the model. Because the term already carries a population mean, it takes over the intercept’s role. sales ~ x + C(store, prior='hierarchical') compiles to x + C(store, ...) with mu_beta_sales_store as the grand intercept, and equations() and design() show no separate Intercept. Treatment-coded terms leave the baseline with the intercept, so a hierarchical term absorbs it regardless of where it appears: y ~ x + ch + C(region, prior='hierarchical') and y ~ x + C(region, prior='hierarchical') + ch build the same design. Only when an earlier categorical already supplies the baseline — a cell-means term under 0 +, or an earlier hierarchical term — do those coefficients sum to zero, \beta \sim \text{ZeroSumNormal}(0, \sigma_\beta), with no mu_beta_{lhs}_{variable}. The design then uses one sum-to-zero contrast column per non-reference level (ch[S.b]) rather than a second full dummy set, so the matrix stays full rank. A reference= argument has no effect and raises a warning. The same rule is applied to the formula alone, so priors(), equations(), and set_priors() on a data-free model describe the model that compiles once data are attached. The only exception is a 0 + equation where a string column recognised only from the data precedes an explicit C() term: the data-free view treats C() as the baseline supplier, so it omits a reference= that the fitted model will use, and a following hierarchical term shows a population mean the fitted model will not have.
Override the hyperpriors by name:
model = pathmc.model(
"sales ~ C(store, prior='hierarchical')",
data=df,
priors={
"mu_beta_sales_store": Prior("Normal", mu=0, sigma=5),
"sigma_beta_sales_store": Prior("HalfNormal", sigma=2),
},
)Introspection
Everything the compiler froze is visible before fitting:
| Method | Shows |
|---|---|
model.equations() |
C(region, reference='north', levels=['north', 'south', 'west']) in the linear predictor; the reference is omitted for cell-means and hierarchical terms, and prior='hierarchical' is shown when set |
model.design("y") |
the indicator columns with their exact names |
model.priors() |
beta_{lhs}_{variable} and any hierarchical hyperpriors |
model.graph() |
the categorical variable as an exogenous node |
model.summary() |
posterior rows labelled beta_y_region[south], beta_y_region[west] |
Interventions and queries by label
Treatment coding is internal. The query methods do(), ate(), prob(), cate(), counterfactual(), and refute_placebo() accept a level label for a categorical variable wherever they take an intervention (set=, values=, do=), a condition (condition=), or evidence (evidence=), and pathmc rebuilds the indicator rows against the frozen state. The interpret helpers do not accept labels yet; see Current limitations.
model.fit()
south = model.do(set={"region": "south"}) # all rows set to south
contrast = south - model.do(set={"region": "north"}) # DoResult arithmetic
model.ate("y", "region", values=("north", "south")) # same contrast, EstimandResult
model.ate("y", "region", values=("west", "south")) # any pair of observed labels
model.do(set={"region": "west", "x": 1.0}) # mix with numeric interventions
model.do(set={"region": "south"}, kind="predictive") # propagate residual noise
model.prob("y > 2", set={"region": "south"})
model.cate("y", "x", values=(0, 1), condition={"region": "south"})
model.counterfactual( # one unit's region flipped, holding its noise fixed
evidence={"x": 0.4, "region": "south", "y": 3.5}, do={"region": "west"}
)
model.refute_placebo("y", "region", values=("north", "south")) # permutes labelsSetting the reference level selects the all-zero indicator row; setting any other level activates its column. A label must be a single scalar — per-row arrays of labels are not accepted by do().
Because labels are not numbers, the extrapolation warning and scaling that apply to numeric interventions are skipped for categorical ones. counterfactual() always needs a label for every categorical predictor in evidence, even with allow_partial_evidence=True: a label has no population mean to fall back on.
Testing the DAG with categorical variables
test_implications() and falsify() test conditional independences with partial correlations. A categorical variable enters those tests as its indicator columns. In a conditioning set that adjusts for group membership exactly as a regression would. When the categorical is one of the two variables being tested, the test becomes the linear multivariate generalisation (Wilks’ lambda with Rao’s F approximation), which is the partial F-test of the indicator block when the other variable is continuous and reduces to the usual t-test when both are single columns. The partial_corr column of test_implications() then reports the largest canonical correlation, which is non-negative. Integer-coded labels declared with C() are encoded the same way. A single-level categorical contributes no indicator columns, so a conditioning set that contains one is rank-collapsed: the row is neither a pass nor a fail of the independence that was written down.
Prediction on new rows
predict(data=new_df) swaps the replacement predictor values into the model’s mutable data nodes and samples the posterior predictive. Categorical columns are encoded with the fitted levels and reference, so:
- rows may arrive in any order, and a level may be absent from the new frame;
- the frame must have the same number of rows as the fitting data (variable-length prediction is not supported yet) and contain every predictor column;
- a level that was not observed at fit time raises a
ValueErrornaming the predictor, the offending level, and the fitted levels.
The last point is deliberate. Recomputing the coding from new data — as pandas.get_dummies or a fresh patsy call would — silently changes what every coefficient means.
flowchart LR
D["fitting data"] --> F["resolve once<br/>levels · reference · columns"]
F --> C["compile<br/>pm.Data indicator matrix<br/>beta_{lhs}_{var}[levels]"]
C --> P["predict(data=…)<br/>re-encode → swap matrix"]
C --> I["do(set={var: label})<br/>encode label → swap matrix"]
P -. "unseen label" .-> E["ValueError"]
I -. "unseen label" .-> E
Scaling and numeric labels
When you pass scaling= to model() or simulate(), channel scaling is applied to numeric predictors only. A column declared categorical with C() — including integer-coded labels — is excluded from the fitted ScalingFactors, and a pre-fitted ScalingFactors object that happens to contain such a column has that entry dropped. Dividing a group code by its maximum would be meaningless.
Simulation
pathmc.simulate() uses the same resolution. Supply the coefficient vector under its usual name with one entry per represented level in frozen order; simulate_params_template() reports the required shapes:
pathmc.simulate(
"y ~ x + region",
data=exog,
params={"beta_y": [1.0, 0.5], "beta_y_region": [2.0, -1.0], "sigma_y": 0.5},
random_seed=0,
)Errors you may see
Each check runs at model() time (or at the offending predict()/do() call) and names both the problem and the fix.
| Situation | Exception | Message names |
|---|---|---|
C(region) but no region column |
ValueError |
the missing column |
reference= not among observed levels |
ValueError |
the available levels |
| fewer than two observed levels | ValueError |
the level count |
| missing values in the column | ValueError |
the column; drop or impute |
| distinct values rendering to the same label | ValueError |
both values and their types |
coefficient prefix on a categorical (b*region, 2*region) |
ParseError |
the prefix to remove |
unknown C() argument |
ParseError |
the valid arguments |
interaction with a categorical (region:x, C(region):x) |
NotImplementedError / ParseError |
the interaction term |
transform on a categorical (adstock(region, …)) |
NotImplementedError |
the term |
| same column used categorically and continuously (any equation) | NotImplementedError |
the column and equation |
| same categorical listed twice in one equation | ValueError |
the column and equation |
categorical column on the left of ~ |
NotImplementedError |
the outcome; suggests recoding + families= |
object column mixing numbers and labels |
ValueError |
the column; suggests C() or recoding |
| panel model with a categorical | NotImplementedError |
— |
| unseen level in predict(), do(), counterfactual(), or refute_placebo() | ValueError |
predictor, level, fitted levels |
array of labels in do(set=…) |
ValueError |
the predictor |
| labels passed for a numeric treatment in refute_placebo() | ValueError |
the treatment |
categorical predictor missing from counterfactual(evidence=…) |
ValueError |
the predictor |
reference= on a cell-means or hierarchical term |
UserWarning |
the predictor; the argument has no effect |
Current limitations
- Cross-sectional only. Panel models (
panel=) reject categorical terms; encode indicators manually if needed. - No user-selected contrasts. Independent-prior terms are treatment or cell-means coded. A later hierarchical term uses an internal sum-to-zero constraint so the design stays full rank. Helmert, polynomial, and ordered contrasts are not available.
- No categorical interactions.
region:xandC(region):xare rejected; fit per-level models or construct indicator × covariate columns yourself. - No categorical outcomes. A string column cannot be on the left of
~; recode and pick a family. - Interpret API. predictions(), comparisons(), slopes(), att(), atu(), and effect() do not yet accept labels for categorical variables. Use do(), ate(), prob(), cate(), counterfactual(), and refute_placebo() as shown above.
- Fixed row count for
predict(data=…).
See the Categorical Predictors example for a complete simulate → inspect → fit → intervene → predict walkthrough with figures.