Вы можете попробовать модули в sklearn, переменная отклика может быть 2 или более размерной, и я думаю, что она работает для OLS (линейная регрессия), лассо, гребень. Модели в statsmodels могут сделать только 1 ответ (только что проверил) ).
Пример набора данных:
from sklearn.datasets import load_iris
iris = load_iris()
df = pd.DataFrame(data= iris['data'],
columns= iris['feature_names'] )
df.shape
(150, 4)
Теперь мы подходим:
from sklearn import linear_model
clf = linear_model.LinearRegression()
clf.fit(df[['sepal length (cm)']],df[['petal length (cm)','petal width (cm)']])
clf.coef_
array([[1.85843298],
[0.75291757]])
Вы можете видеть, что коэффициенты такие же, как при подборе одного ответа в этот случай:
clf.fit(df[['sepal length (cm)']],df[['petal width (cm)']])
clf.coef_
array([[0.75291757]])