ABCOLS

Abundance-based constrained OLS for categorical modifiers

from _api_doc_utils import *

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 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.

4 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.

5 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, /)

6 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]

7 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,)