In 2025, this is not possible with sklearn (I spent a while looking for a solution)

This is a valid task in multi-class logistic regression. You are asking to find a single set of coefficients that simultaneously explain all of the classes. Sklearn (and other packages) find a set of coefficients for each class, which is why it returns a 6x4 matrix - you have 6 features and 4 targets.

The PyLogit package can fit your model. There is an example of performing logistic regression on a dataset with 4 input features and 4 targets - see Specify and Estimate a Multinomial Logit (MNL) Model.