from _api_doc_utils import *ABCOLS
Abundance-based constrained OLS for categorical modifiers
1 Where it fits
Group: Regression
ABCOLS is an OLS reparameterization for categorical main effects and categorical modifiers using abundance-based constraints / weighted effect coding.
Instead of treating one level as the omitted baseline, it estimates an overcomplete dummy/interactions design under linear constraints that force categorical effects to average to zero under empirical level frequencies. This makes the intercept and continuous main slopes sample-abundance-weighted averages rather than reference-category coefficients.
2 Criterion and identification
Let \(X_f\) be the full, overcomplete design containing the intercept, continuous regressors, every observed indicator level, and all requested continuous-by-category and category-by-category terms. The class constructs a constraint matrix \(A\) from empirical category abundances and solves
\[ \min_\theta \|y-X_f\theta\|_2^2 \quad\text{subject to}\quad A\theta=0. \]
For a categorical main effect, the constraint is the abundance-weighted zero mean \(\sum_\ell \hat p_\ell\theta_\ell=0\). Continuous-by-category coefficients use the same weighted-zero constraint over category levels. A category-by-category table receives weighted zero-margin restrictions along both dimensions, with one redundant margin removed.
The implementation eigendecomposes \(A'A\) and collects an orthonormal null-space basis \(Q\) satisfying \(AQ=0\). It then parameterizes \(\theta=Q\phi\) and solves
\[ \hat\phi=\arg\min_\phi\|y-X_fQ\phi\|_2^2, \qquad \hat\theta=Q\hat\phi. \]
This identifies one empirically centered representative from the otherwise non-unique full-indicator coefficient system. The constraints depend on the estimation sample, so coefficients need not be invariant to changing the sample’s category proportions.
3 Implementation walkthrough
The abundance-based parameterization is constructed explicitly; it is not implemented by choosing omitted levels and recoding coefficients afterward.
Continuous columns are optionally centered by their estimation-sample means, which are retained for prediction. Each categorical column must use contiguous observed codes \(0,\ldots,L_k-1\); a gap is rejected rather than silently treated as an empty level.
The code computes empirical marginal level shares and, for every requested category-by-category interaction, empirical cell shares. Every requested interaction cell must occur. These sample shares become coefficients in the restrictions, so the constraints and the coefficient interpretation are data-dependent.
The full design is materialized in a deterministic order: intercept, continuous main effects, all indicator columns for each category, requested continuous-by-category products, and requested category-by-category cell indicators. Metadata for every column is stored alongside human-readable names.
Constraint rows are assembled against those column metadata. Main categorical effects and continuous-by-category blocks receive abundance-weighted zero-sum rows. A two-category table receives weighted margin restrictions in both directions, with one redundant margin omitted before numerical null-space construction.
The code forms \(A'A\), performs a symmetric eigendecomposition, and defines a numerical null direction when its eigenvalue is no larger than
\[ 100\max(m,p)\epsilon\max\{\lambda_{\max}(A'A),1\}. \]
The selected eigenvectors are the orthonormal columns of \(Q\). This avoids choosing an arbitrary reference cell but squares the condition number of \(A\).
It materializes \(Z=X_fQ\), checks that at least one free coordinate and positive residual degrees of freedom remain, solves \(Z\phi\approx y\) by least squares, and maps back with \(\theta=Q\phi\). Fitted values use the original full design and constrained coefficients.
Prediction centers new continuous values with the stored training means, rejects unseen categorical codes, rebuilds the same full column layout, and multiplies by \(\hat\theta\). The reported maximum violation is \(\|A\hat\theta\|_\infty\) and is a direct numerical audit of the reparameterization.
The null-space route makes the estimand transparent and keeps the least-squares solve unconstrained. Its cost is dense expansion plus a tolerance-sensitive eigendecomposition; a sparse QR of \(A'\) would be a more scalable construction for high-cardinality systems.
4 Inference
Inference is homoskedastic and conditional on the constructed design and empirical constraints. With \(Z=X_fQ\), \(r=\operatorname{rank}(Z)\), and \(e=y-Z\hat\phi\),
\[ \hat\sigma^2=\frac{e'e}{n-r}, \qquad \widehat{\operatorname{Var}}(\hat\phi) = \hat\sigma^2(Z'Z)^{-1}, \qquad \widehat{\operatorname{Var}}(\hat\theta) = Q\widehat{\operatorname{Var}}(\hat\phi)Q'. \]
The reported standard errors are the square roots of the final diagonal. The class does not implement heteroskedastic, clustered, serial-correlation-robust, or bootstrap covariance. It also treats the abundance constraints as fixed rather than accounting for their sampling variation.
5 Performance and numerical behavior
The main risk is feature expansion. If \(p_f\) is the width of the full design, the constraint eigendecomposition and covariance are dense \(p_f\times p_f\) operations with roughly cubic factorization cost and quadratic storage. Category-by-category interactions can make \(p_f\) grow as the product of level counts, and the fit requires every requested interaction cell to be represented. Prediction rebuilds the same dense full design and rejects unseen or inconsistently ordered category levels. This estimator is appropriate for moderate categorical systems where the centered coefficient parameterization is substantively useful, not high-cardinality sparse designs.
6 Python API
Constructor: cm.ABCOLS
Call fit(y, x, categories, cont_cat_interactions=None, cat_cat_interactions=None, center_continuous=True). Categorical inputs are zero-based dense uint32 codes. predict(x, categories) returns fitted means for new rows under the same coding scheme. summary() reports constrained coefficients, standard errors, column names, constraint names, residual variance, residual degrees of freedom, rank, and a maximum-constraint-violation diagnostic.
print(inspect.signature(cm.ABCOLS))()cls = cm.ABCOLS
display(HTML(html_table(["Public method"], public_methods(cls))))| Public method |
|---|
column_names(self, /) |
constraint_matrix(self, /) |
design_matrix(self, /) |
fit(self, /, y, x, categories, cont_cat_interactions=None, cat_cat_interactions=None, center_continuous=True) |
fitted_values(self, /) |
predict(self, /, x, categories) |
residuals(self, /) |
summary(self, /) |
7 Minimal example
rng = np.random.default_rng(2026)
group = np.repeat(np.array([0, 1, 2], dtype=np.uint32), [36, 54, 30])
sex = np.tile(np.array([0, 1], dtype=np.uint32), len(group) // 2)
categories = np.column_stack([group, sex]).astype(np.uint32)
x_raw = rng.normal(size=len(group))
x = x_raw[:, None]
x_centered = x_raw - x_raw.mean()
y = 1.25 + 1.1 * x_centered + np.array([-0.75, 0.15, 0.95])[group] + np.array([0.45, -0.2, 0.1])[group] * x_centered + 0.35 * sex + rng.normal(scale=0.08, size=len(group))
model = cm.ABCOLS()
model.fit(y, x, categories, cont_cat_interactions=[(0, 0)], cat_cat_interactions=[(0, 1)])
print(model.summary()['column_names'])
print(model.summary()['coef'][:6])
print(model.predict(x[:3], categories[:3]))['Intercept', 'x0', 'c0[0]', 'c0[1]', 'c0[2]', 'c1[0]', 'c1[1]', 'x0:c0[0]', 'x0:c0[1]', 'x0:c0[2]', 'c0[0]:c1[0]', 'c0[0]:c1[1]', 'c0[1]:c1[0]', 'c0[1]:c1[1]', 'c0[2]:c1[0]', 'c0[2]:c1[1]']
[1.5119991962696528, 1.175173399782689, -0.8346135162279964, 0.07818245705941185, 0.8608077967666543, -0.1750127878211632]
[-0.8050577138355097, 1.1447560278107183, -2.524068516552255]8 summary() contract
The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.
rng = np.random.default_rng(2126)
group = np.repeat(np.array([0, 1, 2], dtype=np.uint32), [24, 30, 18])
sex = np.tile(np.array([0, 1], dtype=np.uint32), len(group) // 2)
categories = np.column_stack([group, sex]).astype(np.uint32)
x_raw = rng.normal(size=len(group))
x = x_raw[:, None]
x_centered = x_raw - x_raw.mean()
y = 0.8 + 0.9 * x_centered + np.array([-0.4, 0.1, 0.6])[group] + np.array([0.2, -0.1, 0.05])[group] * x_centered + 0.25 * sex + rng.normal(scale=0.1, size=len(group))
model = cm.ABCOLS()
model.fit(y, x, categories, cont_cat_interactions=[(0, 0)], cat_cat_interactions=[(0, 1)])
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))| summary() key | shape |
|---|---|
coef |
(16,) |
se |
(16,) |
column_names |
(16,) |
constraint_names |
(7,) |
sigma2 |
() |
df_resid |
() |
rank |
() |
max_constraint_violation |
() |
continuous_means |
(1,) |
n_levels |
(2,) |