It would be nice to support true multi-level models like those in Bambi. However, this does not fit that nicely in the fit(X,y) framework.
We may need to add additional keyword arguments to fit. This is permitted, as long as they are strictly data dependent
