Осваивать смешанные линейные модели statsmodels - PullRequest
Новая
       25

Осваивать смешанные линейные модели statsmodels

0 голосов
/ 13 декабря 2018

У меня есть некоторые проблемы с пониманием использования смешанных линейных моделей в statsmodels.

Насколько я понимаю, смешанные модели предоставляют инструмент для проведения линейных регрессий в зависимости от некоторых групповых ограничений.В частности, я могу смоделировать что-то вроде:

y = a1 + b1 * X для X в группе "0"

и

y = (a1 + 1 * a2) + (b1 + 1 * b2) * X для X в группе "1"

И, естественно, обе модели можно суммировать как

y = a1 + (a2 | group = 1) +(b1 + (b2 | group = 1)) * X

Я пытаюсь точно смоделировать это в python, используя statsmodels.

Мой код выглядит так: 

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import statsmodels.api as sm
import statsmodels.formula.api as smf
from xgboost import XGBRegressor

a1 = 1
a2 = 3
b1 = 1.2
b2 = 4
g = 1

xx = np.linspace(-4,4,100)

y   = np.concatenate([a1 + b1*xx, a1 + g * a2 + (b1 * g*b2)*xx])
cat = np.concatenate([np.zeros(len(xx)), np.ones(len(xx))])

plt.plot(np.concatenate([xx,xx]),y, 'x')
plt.title('Raw Data')

Raw Data

data = pd.DataFrame(np.array([y,np.concatenate([xx,xx]),cat]).T, columns = ['y', 'Time', 'Group'])
data['Group'] = data['Group'].astype("category")

md = smf.mixedlm("y ~ 1 + Time", data, groups = data["Group"],     re_formula="~Time")

mdf = md.fit()
print(mdf.summary())

yhat = mdf.predict(data)

plt.plot(np.concatenate([xx,xx]),y, 'x')
plt.plot(np.concatenate([xx,xx]),yhat, 'xr')
plt.title('Fit Mixed Model')

Fit of Mixed Model

xgbm = XGBRegressor(max_depth = 2)
XX = np.array([data['Time'], data['Group'].astype(int), data['Time']*data['Group'].astype(int)]).T
xgbm.fit(XX,y)
yhat2 = xgbm.predict(XX)

plt.plot(np.concatenate([xx,xx]),y, 'x')
plt.plot(np.concatenate([xx,xx]),yhat2, 'xr')
plt.title('Fit Costum Model')

Fit of Costum Model

Линейная смешанная модель просто дает линейное соответствие «средних» данных.Я не понимаю, почему, поскольку я определил групповой признак в смешанной модели (groups = data ['Group']).Модель costum (через xgboost; OLS работает аналогично) дает результаты, как и предполагалось.Он «разбивает» данные в зависимости от их групповой характеристики и дает разные оценки для наклона и перехвата для каждой группы.

Что-то не так в моем определении смешанной модели?

...