Model selection for a “set of linear models”?

Question

I would appreciate your input on which predictive ML model(s) could fit our dataset the best.

The primary features of the dataset are x,y (continous) and c (factor with around L=500 unordered levels). The secondary features (~20) are both continous and factors. We know that the lowest-order correlation is x~y. We also know that modeling x~y for each level in c improves the fit significantly.

We could of course create L independent linear models, but would like to explore a single model where these sub-linear relationships are partially correlated. The reason being that some levels of c have relatively low statistics and would benefit from the global x~y relationship. Also, we would like a single model when extending to secondary features.

Any suggestions are welcome. (In case of specific libraries, R is preferred.)

Answer 1

A generalized linear model (glm) appeared to be work reasonable.

glm(y ~ x * c, data=df) and extensions hereof.