BUG: Formula with categoricals results in singular matrix

[lorentzenchr]

Issue A

When a formula with a categorical "C(..)" is used without penalty, no level is dropped, i.e. one-hot encoding is used instead of treatment coding.

from glum import GeneralizedLinearRegressor
import pandas as pd

df = pd.DataFrame({"y": [0, 1, 1, 2], "cat": pd.Series(["a", "a", "b", "b"], dtype="category")})
model = GeneralizedLinearRegressor(
    alpha=0,
    formula="y ~ C(cat)"
)
model.fit(df)

results in

LinAlgError: Matrix is singular.

model.feature_names_ outputs ['C(cat)[a]', 'C(cat)[b]'] and reveals that no reference level is chosen.

The is no difference in:

• "y ~ C(cat)"
• "y ~ C(cat, spans_intercept=False)"
• "y ~ C(cat, spans_intercept=True)"

However, the error goes away with GeneralizedLinearRegressor(drop_first=True, ..).

Issue B

According to formulaic, it should be possible to call

model = GeneralizedLinearRegressor(
    alpha=0,
    formula="y ~ C(cat, contr.treatment)"
)
model.fit(df)

But this results in

FactorEvaluationError: Unable to evaluate factor `C(cat, contr.treatment)`. [TypeError: _C() takes 1 positional argument but 2 were given]

Formulas are specified as strings and neither formulaic nor glum nor tabmat have an API reference for, e.g., C. As a user, there is no way to infer the API because C is only ever specified as string inside a formula, i.e., "y ~ C(cat)" because the symbol C is not exposed anywhere.

The only way to see what happens is digging into the source code. There I learned that glum delegates formulas to tabmat. But tabmat overrides the C of formulaic, see https://github.com/Quantco/tabmat/blob/0e2608fcaca9f11830bafdf873a700c89de1391f/src/tabmat/formula.py#L665. This "C" (or _C) has a different API with different parameters.

This deficiency of strings for formulas was one reason I proposed a programmatic alternative to strings in #731.

[stanmart]

I would not classify it as a bug. I think it would be more confusing if glum had different drop_first behavior depending on the regularization parameter. I believe the documentation is more-or-less clear

drop_first (bool, optional (default = False)) – If True, drop the first column when encoding categorical variables. Set this to True when alpha=0 and solver='auto' to prevent an error due to a singular feature matrix. In the case of using a formula with interactions, setting this argument to True ensures structural full-rankness (it is equivalent to ensure_full_rank in formulaic and tabmat).

but of course suggestions for improving it are much appreciated.

The is no difference in: [...]

This is also expected behavior. When drop_first=False, glum does not attempt to ensure structural full-rankness, i.e., drop any levels. spans_intercept is used to signal to glum that a certain catagorical would not span the intercept even if all levels were to be expanded (e.g., because of treating missings as all-zero indicators), so it only has an effect with drop_first set to true.

According to formulaic, it should be possible to call [...]

glum only implements a subset of formulaic's C function. The reason is that categoricals are special-cased in tabmat for performance, so we can only support the most basic encoding types. contr.treatment is the only option that we support, so this parameter is not implemented. If performance is not an issue, one can always fall back on formulaic to create a (pandas) design matrix and pass that to glum.

Formulas are specified as strings and neither formulaic nor glum nor tabmat have an API reference for, e.g., C. As a user, there is no way to infer the API because C is only ever specified as string inside a formula, i.e., "y ~ C(cat)" because the symbol C is not exposed anywhere.

Good point, we should extend the description of the formula parameter with this information.

This deficiency of strings for formulas was one reason I proposed a programmatic alternative to strings in #731.

I agree, one could do better than string formulas. We decided to implement this API, because (1) it is very familiar to R and statsmodels users, and (2) we could rely on existing upstream tools. Of course it comes with all the limitations of stringy solutions (no editor support, static analysis, etc.), but I did not view those as a big downside as it is meant for quick prototyping and throwaway scripts. I'd think that for more complex uses one would rather rely on, e.g., scikit-learn to build a preprocessing pipeline. Nevertheless, I do absolutely see value in a better, non-string based interface (I especially share your desire to be able to include per-term penalties). Unfortunately that would be a much larger undertaking, though.

Edit: oops, sorry, closed it by accident.