StanとRでベイズ統計モデリング(アヒル本)をPythonにしてみる - 8.4 ロジスティック回帰の階層モデル

実行環境

インポート

import numpy as np
from scipy import stats
import pandas as pd
import pystan
import matplotlib.pyplot as plt
from matplotlib.figure import figaspect
import seaborn as sns
%matplotlib inline

データ読み込み

attendance41 = pd.read_csv('./data/data-attendance-4-1.txt')
attendance42 = pd.read_csv('./data/data-attendance-4-2.txt')

8.4 ロジスティック回帰の階層モデル

8.4.4 Stanで実装

N = attendance41.index.size
C = attendance42['CourseID'].max()
I = attendance42.index.size
conv = pd.Series([0, 0.2, 1], index=('A', 'B', 'C'))
data = dict(
    N=N,
    C=C,
    I=I,
    A=attendance41['A'],
    Score=attendance41['Score']/200,
    PID=attendance42['PersonID'],
    CID=attendance42['CourseID'],
    W=conv[attendance42['Weather']],
    Y=attendance42['Y']
)
fit = pystan.stan('./stan/model8-8.stan', data=data, pars=('b', 'b_P', 'b_C', 's_P', 's_C', 'q'), seed=1234)

8.4.5 推定結果の解釈

ms = fit.extract()
N_mcmc = ms['lp__'].size

param_names = ['mcmc'] + ['b{}'.format(i+1) for i in range(4)] + ['s_P', 's_C']
d_est = pd.DataFrame(np.hstack([np.arange(N_mcmc).reshape((-1, 1)), ms['b'], ms['s_P'].reshape((-1, 1)), ms['s_C'].reshape((-1, 1))]), columns=param_names)
d_qua = d_est.loc[:, param_names[1:]].quantile((0.025, 0.5, 0.975)).T.reset_index()
d_qua.columns = ('X', 'p2.5', 'p50', 'p97.5')
d_qua['X'] = pd.Categorical(d_qua['X'])
d_melt = pd.melt(d_est, id_vars=('mcmc'), var_name='X')
d_melt['X'] = pd.Categorical(d_melt['X'])

_, (ax1, ax2) = plt.subplots(1, 2, figsize=figaspect(1/2))
sns.violinplot('value', 'X', data=d_melt, order=param_names[1:], scale='width', inner=None, orient='h', color='w', ax=ax1)
ax1.errorbar('p50', 'X', data=d_qua, xerr=[d_qua['p50']-d_qua['p2.5'], d_qua['p97.5']-d_qua['p50']], fmt='o', c='k')
plt.setp(ax1, xlabel='value', ylabel='parameter')

param_names = ['mcmc'] + ['b_C{}'.format(i+1) for i in range(10)]
d_est = pd.DataFrame(np.hstack([np.arange(N_mcmc).reshape((-1, 1)), ms['b_C']]), columns=param_names)

def get_map(col):
    kernel = stats.gaussian_kde(col)
    dens_x = np.linspace(col.min(), col.max(), kernel.n)
    dens_y = kernel.pdf(dens_x)
    mode_i = np.argmax(dens_y)
    mode_x = dens_x[mode_i]
    mode_y = dens_y[mode_i]
    return pd.Series([mode_x, mode_y], index=['X', 'Y'])

d_mode = d_est.loc[:, param_names[1:]].apply(get_map).T

d_est.loc[:, param_names[1:]].apply(lambda s: sns.kdeplot(s, shade=True, legend=False, color='k', alpha=0.15, ax=ax2))
ax2.vlines(d_mode['X'], d_mode['Y'], 0, linestyles='dashed', alpha=0.6)
sns.rugplot(d_mode['X'], color='k', ax=ax2)
plt.setp(ax2, xlabel='value', ylabel='density', xticks=np.arange(-4, 5, 2))
plt.tight_layout()

plt.show()