pythonのstatsmodelsによる階層線形モデル（混合線形モデル)¶

Rの階層線形モデルlme4のlmerと同等のモデルがstatsmodelsにあるのでそれを利用してみる。

pythonのデータ処理になれている場合はこちらの方が良いか。

import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline

import pandas as pd
test_data = pd.read_csv("dat/nine-school-testdata.csv", sep=",", index_col=0)

test_data.head()

単回帰¶

import statsmodels.api as sm
Y = test_data['normexam'].values
df_X = test_data[["standLRT"]]
X = df_X.values
X1 = sm.add_constant(X) # 1列めに1.0を入れるだけ
model = sm.OLS(Y,X1)
result = model.fit()
result.summary()

# pandasのDataFrameを利用するにはstatsmodels.formula.apiを利用する
import statsmodels.formula.api as smf
result2 = smf.ols(formula="normexam ~ standLRT", data=test_data).fit()
result2.summary()
result2.fittedvalues[:5]
test_data['fitted_olm'] = result2.fittedvalues

混合線形回帰¶

statsmodelsのmixedlmを用いる。 https://www.statsmodels.org/stable/mixed_linear.html

結果はRのlmerを利用したときと同じ

利用例のページ http://www.statsmodels.org/0.6.1/examples/notebooks/generated/mixed_lm_example.html

import statsmodels.formula.api as smf
md = smf.mixedlm("normexam ~ standLRT", test_data, groups=test_data["school"], re_formula="~standLRT")
mdf = md.fit()
print(mdf.summary())

            Mixed Linear Model Regression Results
=============================================================
Model:              MixedLM   Dependent Variable:   normexam 
No. Observations:   514       Method:               REML     
No. Groups:         9         Scale:                0.5336   
Min. group size:    29        Likelihood:           -586.9361
Max. group size:    80        Converged:            Yes      
Mean group size:    57.1                                     
-------------------------------------------------------------
                     Coef. Std.Err.   z   P>|z| [0.025 0.975]
-------------------------------------------------------------
Intercept            0.117    0.163 0.716 0.474 -0.203  0.437
standLRT             0.618    0.076 8.127 0.000  0.469  0.767
Group Var            0.229    0.166                          
Group x standLRT Cov 0.086    0.069                          
standLRT Var         0.040    0.035                          
=============================================================

# ランダム効果の係数部分

mdf.random_effects

{'学校1': Group      -0.371025
 standLRT   -0.103408
 dtype: float64, '学校2': Group       0.127820
 standLRT    0.034374
 dtype: float64, '学校3': Group       0.554811
 standLRT    0.316915
 dtype: float64, '学校4': Group       0.367999
 standLRT    0.089224
 dtype: float64, '学校5': Group      -0.849615
 standLRT   -0.319189
 dtype: float64, '学校6': Group       0.408482
 standLRT    0.076501
 dtype: float64, '学校7': Group      -0.435621
 standLRT   -0.188390
 dtype: float64, '学校8': Group       0.267753
 standLRT    0.111118
 dtype: float64, '学校9': Group      -0.070604
 standLRT   -0.017145
 dtype: float64}

# 各データに対するフィッティング値
mdf.fittedvalues[:5]

1767   -0.105544
1768   -1.083606
1769   -0.445739
1770    0.064554
1771   -0.785935
dtype: float64

test_data['fitted_lmer'] = mdf.fittedvalues

test_data.head()

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
test_data['color'] = test_data['school'].str.replace("学校","").astype(int) # 色は整数値で指定
test_data.plot(kind="scatter", ax=axes[0], x="fitted_olm", y="normexam", c="color", colormap='Accent', colorbar=False)
test_data.plot(kind="scatter", ax=axes[1], x="fitted_lmer", y="normexam", c="color", colormap='Accent', colorbar=False)
print(test_data[['fitted_olm','normexam']].corr())
print(test_data[['fitted_lmer','normexam']].corr())

            fitted_olm  normexam
fitted_olm    1.000000  0.660071
normexam      0.660071  1.000000
             fitted_lmer  normexam
fitted_lmer     1.000000  0.774337
normexam        0.774337  1.000000

test_data.plot(kind="scatter", y="fitted", x="standLRT")

<matplotlib.axes._subplots.AxesSubplot at 0x7f44e9119048>

	school	sex	normexam	standLRT	type
1767	学校1	F	-0.852670	0.288453	Sngl
1768	学校1	F	-0.338842	-1.612530	Sngl
1769	学校1	F	0.194149	-0.372758	Sngl
1770	学校1	F	-0.129085	0.619059	Sngl
1771	学校1	F	0.004322	-1.033970	Sngl

Dep. Variable:	y	R-squared:	0.436
Model:	OLS	Adj. R-squared:	0.435
Method:	Least Squares	F-statistic:	395.3
Date:	Mon, 22 Apr 2019	Prob (F-statistic):	1.30e-65
Time:	11:46:39	Log-Likelihood:	-649.55
No. Observations:	514	AIC:	1303.
Df Residuals:	512	BIC:	1312.
Df Model:	1
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
const	0.1780	0.038	4.696	0.000	0.104	0.252
x1	0.7444	0.037	19.882	0.000	0.671	0.818

Omnibus:	0.186	Durbin-Watson:	1.460
Prob(Omnibus):	0.911	Jarque-Bera (JB):	0.153
Skew:	-0.042	Prob(JB):	0.926
Kurtosis:	3.003	Cond. No.	1.06