Add support to `ModelSpec` for extracting metadata about factor levels used in specific columns.

[s3alfisc]

Hi @matthewwardrop ,

I noticed a change in the naming of treatment encoded categoricals shipped with formulaic 1.1.0. Is this a desired change?
In some circumstances, treatment encoded variables are now called C(state)[2] instead of previously C(state)[T.2].

Is this an intended change / intended behavior?

import pandas as pd 
import formulaic
from formulaic import Formula

# Formulaic >= 1.1.0 
print(formulaic.__version__)
# 1.1.0
data = pd.read_csv("https://raw.githubusercontent.com/py-econometrics/pyfixest/refs/heads/master/pyfixest/did/data/df_het.csv")
data[["state", "unit"]].dtypes
# np.int64, np.int64

mm = Formula("~ -1 +  C(state) + C(unit)").get_model_matrix(data)
mm.model_spec.column_names[0:5]
# ('C(state)[1]', 'C(state)[2]', 'C(state)[3]', 'C(state)[4]', 'C(state)[5]')

# Formulaic < 1.1.0 
print(formulaic.__version__)
# 1.0.2
data = pd.read_csv("https://raw.githubusercontent.com/py-econometrics/pyfixest/refs/heads/master/pyfixest/did/data/df_het.csv")

mm = Formula("~ -1 +  C(state) + C(unit)").get_model_matrix(data)
mm.model_spec.column_names[0:5]
# ('C(state)[T.1]','C(state)[T.2]','C(state)[T.3]',

Best, Alex

[matthewwardrop]

Hi @s3alfisc ,

Apologies for breaking your code!

Yes, this was intentional; see: #220 . In short, formulaic's previous behaviour was different than patsy's; whereby the T. was only prepended if contrasts were being created, not when the full rank dummies were being used. I did not think anyone was hard-coding against category column names, otherwise I would have made a bigger deal of this in the release notes (where it is mentioned, but not as a breaking change).

Do you need additional metadata to be propagated into the model spec so you do not need to do a lookup with column strings?

Best,
Matthew

[s3alfisc]

Great, then everything makes sense!

Apologies for breaking your code!
No worries about that at all - already when writing my code I thought that I was inviting you to break it eventually!

Metadata that would help me would be the mapping from a column / coefficient name to a categorical variable and level:
m(column name) -> (variable name, level).

I.e. if we have a column named C(f1)[T.1] , we'd have m(C(f1)[T.1]) -> (C(f1), 1).

Maybe there is already an option available to do this programmatically? (I suppose so as formulaic must create the column names somewhere 😄 ?)

If you want to, you can stop reading here, below I just add some detail / context on what I would use this mapping for.

Details

The code section where I use the column name parser estimates the fixed effects coefficients "post hoc" to allow for predictions with fixed effects. More formally, we estimate a regression model

$$ Y_i = X_i \beta + \alpha D_i + u_i $$

where D_i is a matrix of dummies, via FWL and retrieve \beta. To make a prediction, we need an estimate of \alpha, which we obtain by solving a sparse system of equations.

If we make a prediction on the initial data set, everything is easy: we can simply add $X_i \hat{\beta}$ and $\hat{\alpha} D_i$.

If we want to make a prediction on new data, things are more difficult because the new data might contain new levels for the fixed effects, for which we do not have any estimates. Predictions with such fixed effect levels will get a value of np.nan.

To handle this, we store the estimated alpha in a dictionary with keys fixed effect x level and do a lookup at prediction time. To create the dictionary levels and keys, we currently deparse the model column name.

For example, for a coefficient named C(f1)[T.1], we set a column name of C(f1) and a level of 1, and we get both by using a regex on the column name. The mapping of coefficient name -> (variable name, level) could then replace the string parsing.

[matthewwardrop]

That makes sense (and thanks for the context). How do you handle interactions in the fixed effects? Or is that unsupported? A generic column_name to (variable_name, level) on Formulaic's doesn't work because of such interaction terms.

What we could do is add something like:

...
model_spec.factors_by_column -> {
    "Intercept": {1: None},
    "C(f1)[T.1]:a": {C(f1): (1,), a: None},
}

where the output here is a mapping from the column name to another mapping from the factor object to the level specification, which could be nested (hence the tuple specification, which is also nice because it clearly distinguishes it from None, which indicates that no level was selected).

Would that be useful?

[s3alfisc]

How do you handle interactions in the fixed effects?

Through a hack: the only way to interact fixed effects in the demeaned part of the formula is via a special syntax, the ^ operator, which simply glues var1 and var2 as strings. The resulting variable would then be called var1_var2 and formulaic would name it something as C(var1_var2)[T.level{var1}_level{var2}], and as a result, the pattern matching would work =) I should probably define a custom transform for this as described here? Generally I am still far away from leveraging all of formulaic's power. For example, there's another custom operator, i(var1, ref1), which is simply a fixest untility for writing C(var1, contr.treatment(ref1)). By now, it might not surprise you that I am also using string matching to translate one syntax into the other ...

Yes, I think your suggestion above would indeed help!