Up

度数クロス表データのベイズ分析

ポアッソン回帰モデル

 

図１のデータは、＜芝ら（編）「統計用語辞典」新曜社＞の分割表の例として掲載されているもの（p. 230）を入力したものである。

図１　Excelデータ（datap230.xlsx）

 

図１の形式の度数のクロス表データに対してポアッソン回帰モデルによるベイズ分析を行う。度数データをポアッソン分布で表すとき、パラメータの対数に対して線形モデルが適用される（Reich & Ghosh, 2019）。ポアッソン分布のパラメータの対数を、クロス表の主効果と交互作用の和によって表す（Kruschke, 2015）。

クロス表の分析モデルを一般的に説明するために、表１のように各条件（行i、列j）における効果を表す。

 

表１　主効果と交互作用

 

列１

・・・

列ｊ

・・・

行１

・・・

・・・

行i

・・・

・・・

 

表１における効果を行主効果rai、列主効果caj、交互作用rcaijの和として

と表す。ここで、制約条件

をおく。

制約条件（２）より、以下の式が導かれる。

上の効果に対してポアッソン分布のパラメータを対数関数の逆関数である指数関数によって次のようにおく。

各セルにおける度数は、パラメータがであるポアッソン分布に従うとして、本ウェブサイト後半に示すリスト１のようにStanスクリプトを用意した。リスト１のStanスクリプトによりzijの事後分布が得られる。このzijの事後分布から、式（４）～（６）により主効果および効果量の事後分布を本ウェブサイト後半に示すリスト２における次のスクリプトによって求めている。

 

for vz in fit['z']:

    zdd = vz.mean()             #   z..

    zdc = vz.mean(axis = 0)     #   z.j

    zrd = vz.mean(axis = 1)     #   zi.

    r_e[i] = zrd - zdd          #   行主効果: rai

    c_e[i] = zdc - zdd          #   列主効果: caj

    rc_e[i] = vz - zrd.reshape((KR, 1)) - zdc + zdd     #   交互作用: rcaij

 

スクリプトファイルなどのファイルは、圧縮ファイルfilesPoisson.zipとしてまとめたので、ダウンロード解凍すれば自由に利用できる。なお、StanのPythonにおける利用方法は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第10章～第12章で説明している。

リスト２のスクリプトを実行すると、入力データファイル名の設定を求めるダイアログウィンドウが表示される（図２）。

図２

 

入力データファイルは、図１の形式で用意する。第1行目に、列条件のカテゴリ（図１の場合は、「自宅生」、「寮・下宿生」）を書き、第１列目に行条件のカテゴリ（図１では「男」、「女」）を書く。度数データは２行目以降、２列目から該当するセルに書く。列条件あるいは行条件は図１のものより多くてもよい。列および行条件の個数は２以上であればよい。

図２のダイアログウィンドウにファイル名を設定して「開く」ボタンをクリックすると、画像保存用ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。

図３

 

事後分布のグラフを保存するファイル名であるが、ファイル拡張子が「.png」である任意のファイル名を設定する。設定したファイル名の「.png」の部分が「_ra.png」などに変更された名前のファイルに画像が保存される。

図３のダイアログウィンドウにファイル名を設定して「保存」ボタンをクリックすると、計算結果を保存するためのテキストファイル名の設定を求めるダイアログウィンドウが表示される（図４）。

図４

 

ファイル拡張子が「.txt」の任意のテキストファイル名を設定する。

テキストファイル名を設定して「保存」ボタンをクリックすると、データの読み込みが行われ、Stanスクリプトのコンパイル・ビルドの後、MCMCサンプリングが始まる。

サンプリングが終了すると、まず行の主効果を表すパラメータra（r_e）の事後分布が表示される（図５）。

図５

 

行主効果raは、「男」の事後分布が「女」の事後分布よりかなり明瞭に右側に位置しているが、これは図１のデータにおいて「男」の人数が「女」に人数より多いことに対応している。

図５のウィンドウを閉じると、列の主効果を表すパラメータca（c_e）の事後分布のグラフが表示される（図６）。

図６

 

「自宅生」の事後分布が「寮・下宿生」の事後分布より右側に位置しているが、これは図１のデータにおいて「自宅生」の人数が「寮・下宿生」より多いことに対応している。

図６のウィンドウを閉じると、交互作用のパラメータrca（rc_e）の事後分布が表示される（図７）。

図７

 

図７においてグラフは２本しか見えないが、交互作用「男×自宅生」と「女×寮・下宿生」のグラフが重なっており、交互作用「男×寮・下宿生」と「女×自宅生」のグラフが重なっている。これは、制約条件（２）により１つの交互作用（例えば、rca11）が与えられると他の３つの交互作用（rca12、rca21、rca22）が決まるからである。どの交互作用とどの交互作用が重なっているのかは、スクリプトの実行終了後、出力ファイルを見れば確認できる。

図７の交互作用は、主効果で説明される効果より、「男」においては「自宅生」が少なく、「寮・下宿生」が多いが、「女」においては「自宅生」が多く、「寮・下宿生」が少ないことを表している。

図７のウィンドウを閉じると、の事後分布の中央値のグラフが列条件（いまの場合、「自宅生」と「寮・下宿生」）を横軸として描かれたのが表示される（図８）。

図８

 

図８のウィンドウを閉じると、次にの事後分布の中央値のグラフが行条件（いまの場合、「男」と「女」）を横軸として描かれたものが表示される（図８）。

 

図９

 

図９のウィンドウを閉じると、保存されたファイル名が表示されて、スクリプトの実行終了となる（図10）。

図10

 

出力テキストファイルをPythonのIDLEのエディタで開こうとすると、図８のようなダイアログボックスが表示されることがある。

図８

 

このときは、「s-jis」と設定して「OK」ボタンをクリックすればよい（図９）。

図９、

 

実行終了後、図４で指定したテキストファイルを開くと、以下のようになっている。

 

入力データファイル = D:/yasuharu/www/LaCoocan/Python/IntroPyStat/BysPoisson/filesPoisson/datap230.xlsx

 

行変数：['男', '女']

列変数：['自宅生' '寮・下宿生']

 

Inference for Stan model: anon_model_6e38c9fdf24dcd19d74bd3713ed4a92c.

4 chains, each with iter=2000; warmup=1000; thin=1;

post-warmup draws per chain=1000, total post-warmup draws=4000.

 

            mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat

z[1,1]      2.68  3.9e-3   0.25   2.16   2.51   2.69   2.86   3.15   4157    1.0

z[2,1]       2.8  3.7e-3   0.25   2.29   2.63   2.81   2.98   3.26   4578    1.0

z[1,2]      3.01  3.7e-3   0.23   2.54   2.86   3.02   3.17   3.43   3906    1.0

z[2,2]      1.87  6.1e-3   0.39   1.05   1.62   1.89   2.13   2.57   3996    1.0

lmbd[1,1]  15.02    0.06   3.73   8.64  12.34  14.75  17.39  23.33   4365    1.0

lmbd[2,1]   17.0    0.06   4.22   9.87  13.92  16.61   19.7  25.97   4911    1.0

lmbd[1,2]  20.87    0.07   4.73  12.69  17.54  20.49  23.86  30.95   4148    1.0

lmbd[2,2]   6.96    0.04   2.62   2.85   5.05   6.62   8.45  13.05   4548    1.0

lp__       104.3    0.03   1.43 100.74 103.61 104.63 105.35 106.09   1801    1.0

 

Samples were drawn using NUTS at Tue Aug 11 15:00:09 2020.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at

convergence, Rhat=1).

 

行主効果

 

    男

      平均値 = 0.255

      中央値 = 0.255

      Q1 = 0.159,   Q3 = 0.352

      95%確信区間 = [-0.018, 0.529]

 

    女

      平均値 = -0.255

      中央値 = -0.255

      Q1 = -0.352,   Q3 = -0.159

      95%確信区間 = [-0.529, 0.018]

 

列主効果

 

    自宅生

      平均値 = 0.150

      中央値 = 0.148

      Q1 = 0.055,   Q3 = 0.245

      95%確信区間 = [-0.128, 0.437]

 

    寮・下宿生

      平均値 = -0.150

      中央値 = -0.148

      Q1 = -0.245,   Q3 = -0.055

      95%確信区間 = [-0.437, 0.128]

 

交互作用

 

    男×自宅生

      平均値 = -0.317

      中央値 = -0.316

      Q1 = -0.411,   Q3 = -0.220

      95%確信区間 = [-0.601, -0.045]

 

    男×寮・下宿生

      平均値 = 0.317

      中央値 = 0.316

      Q1 = 0.220,   Q3 = 0.411

      95%確信区間 = [0.045, 0.601]

 

    女×自宅生

      平均値 = 0.317

      中央値 = 0.316

      Q1 = 0.220,   Q3 = 0.411

      95%確信区間 = [0.045, 0.601]

 

    女×寮・下宿生

      平均値 = -0.317

      中央値 = -0.316

      Q1 = -0.411,   Q3 = -0.220

      95%確信区間 = [-0.601, -0.045]

 

 

λijの中央値:

 

λ[男, 自宅生] = 14.745

 

λ[男, 寮・下宿生] = 20.485

 

λ[女, 自宅生] = 16.612

 

λ[女, 寮・下宿生] = 6.615

 

交互作用の結果から、交互作用「男×自宅生」と「女×寮・下宿生」、交互作用「男×寮・下宿生」と「女×自宅生」の事後分布が重なってることが判る。

 

 

引用文献

Kruschke, J. K. (2015). Doing Bayesian data analysis, 2nd ed. Academic Press.

Reich, B. J. & Ghosh, S. K. (2019). Beysian statistical methods. CRC Press.

芝 祐順、渡部 洋、石塚智一（編）（1984）．統計用語辞典．新曜社．

 

 

リスト１　ポアッソン回帰分析用Stanスクリプト（PoissonModel.stan）

 

data {

    int KR;

    int KC;

    int X[KR, KC];

}

parameters {

    real z[KR, KC];

}

transformed parameters {

    real<lower = 0.0> lmbd[KR, KC];

    for (i in 1:KR)

        for (j in 1:KC)

            lmbd[i][j] = exp(z[i][j]);

}

model {

    for (i in 1:KR)

        for (j in 1:KC)

            z[i][j] ~ normal(0.0, 1000.0);

    for (i in 1:KR)

        for (j in 1:KC) {

            X[i][j] ~ poisson(lmbd[i][j]);       

        }

}

 

 

リスト２　リスト１のStanスクリプトを用いるPythonスクリプト（MainPoisson.py）

 

#

#           Yasuharu Okamoto, 2020.08

#

import tkinter as tk

import tkinter.filedialog as fd

import pandas as pd

import numpy as np

import pystan

import matplotlib.pyplot as plt

import seaborn as sb

 

def askinfilename(msg = None, types = [('', '*.*')]):

    """

            入力用ファイル名の設定

    """

    rt = tk.Tk()

    rt.withdraw()

    flnm = fd.askopenfilename(title = msg, filetypes = types)

    rt.destroy()

    return flnm

 

def askoutfilename(msg = None, types = [('', '*.*')], deftype = '.txt'):

    """

            出力用ファイル名の設定

    """

    rt = tk.Tk()

    rt.withdraw()

    flnm = fd.asksaveasfilename(title = msg, filetypes = types,

                                           defaultextension = deftype)

    rt.destroy()

    return flnm

 

finnm = askinfilename('入力データExcelファイル', types = [('Excel(*.xlsx)', '*.xlsx')])

foutnm = askoutfilename(msg = '画像保存ファイル', types = [('画像 (*.png)', '*.png')],

                       deftype = '.png')

ftxtoutnm = askoutfilename(msg = 'テキスト出力ファイル', types = [('Text (*.txt)', '*.txt')],

                           deftype = '.txt')

 

fout = open(ftxtoutnm, 'w')

fout.write('入力データファイル = {}\n'.format(finnm))

 

xlsx_fl = pd.ExcelFile(finnm)

data_xlsx = pd.read_excel(xlsx_fl,  header = None)

data = data_xlsx.values

 

KC = len(data[0]) - 1

print('KC = ', KC)

KR = len(data) - 1

print('KR = ', KR)

 

NameC = data[0][1:]

NameR = []

for i in range(KR):

    NameR.append(data[i + 1][0])

print('NameC: ', NameC)

print('NameR: ', NameR)

 

fout.write('\n行変数：{}\n'.format(NameR))

fout.write('列変数：{}\n'.format(NameC))

 

X = np.empty((KR, KC), dtype = int)

for i in range(KR):

    for j in range(KC):

        X[i][j] = data[i + 1][j + 1]

 

for i in range(KR):

    for j in range(KC):

        print('{0:>10d}'.format(X[i][j]), end = '')

    print()

 

StanData = {'KR': KR, 'KC': KC, 'X': X}

 

sm = pystan.StanModel(file = 'PoissonModel.stan')

fit = sm.sampling(data = StanData, n_jobs = 1)

print('\n', fit)

 

fout.write('\n{}\n'.format(fit))

 

SavedFiles = []

 

#

#       主効果と交互作用

#

n_samples = len(fit['z'])

r_e = np.empty((n_samples, KR))

c_e = np.empty((n_samples, KC))

rc_e = np.empty((n_samples, KR, KC))

i = 0

for vz in fit['z']:

    zdd = vz.mean()             #   g..

    zdc = vz.mean(axis = 0)     #   g.j

    zrd = vz.mean(axis = 1)     #   gi.

    r_e[i] = zrd - zdd          #   行主効果

    c_e[i] = zdc - zdd          #   列主効果

    rc_e[i] = vz - zrd.reshape((KR, 1)) - zdc + zdd     #   交互作用

    i += 1

 

#       日本語表記の設定

plt.rcParams['font.family'] = ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']

 

plt.title('行主効果', fontsize = 20)

fout.write('\n行主効果\n')

for i in range(KR):

    sb.kdeplot(r_e[:, i], label = NameR[i])

    ra_mean = r_e[:, i].mean()

    ra_025p, ra_25p, ra_50p, ra_75p, ra_975p = \

        np.percentile(r_e[:, i], [2.5, 25, 50, 75, 97.5])

    fout.write('\n    {}\n'.format(NameR[i]))

    fout.write('      平均値 = {0:>.3f}\n'.format(ra_mean))

    fout.write('      中央値 = {0:>.3f}\n'.format(ra_50p))

    fout.write('      Q1 = {0:>.3f},   Q3 = {1:>.3f}\n'.format(ra_25p, ra_75p))

    fout.write('      95%確信区間 = [{0:>.3f}, {1:>.3f}]\n'.format(ra_025p, ra_975p))

 

plt.yticks([])

plt.legend()

plt.tight_layout()

plt.savefig(foutnm[:-4] + '_ra.png')

SavedFiles.append(foutnm[:-4] + '_ra.png')

plt.show()

 

plt.title('列主効果', fontsize = 20)

fout.write('\n列主効果\n')

for j in range(KC ):

    sb.kdeplot(c_e[:, j], label = NameC[j])

    ca_mean = c_e[:, j].mean()

    ca_025p, ca_25p, ca_50p, ca_75p, ca_975p = \

        np.percentile(c_e[:, j], [2.5, 25, 50, 75, 97.5])

    fout.write('\n    {}\n'.format(NameC[j]))

    fout.write('      平均値 = {0:>.3f}\n'.format(ca_mean))

    fout.write('      中央値 = {0:>.3f}\n'.format(ca_50p))

    fout.write('      Q1 = {0:>.3f},   Q3 = {1:>.3f}\n'.format(ca_25p, ca_75p))

    fout.write('      95%確信区間 = [{0:>.3f}, {1:>.3f}]\n'.format(ca_025p, ca_975p))

plt.yticks([])

plt.legend()

plt.tight_layout()

plt.savefig(foutnm[:-4] + '_ca.png')

SavedFiles.append(foutnm[:-4] + '_ca.png')

plt.show()

 

plt.title('交互作用', fontsize = 20)

fout.write('\n交互作用\n')

for i in range(KR):

    for j in range(KC):

        sb.kdeplot(rc_e[:, i, j], label = r'{0}$\times${1}'.format(

            NameR[i], NameC[j]))

        ga_mean = rc_e[:, i, j].mean()

        ga_025p, ga_25p, ga_50p, ga_75p, ga_975p = \

            np.percentile(rc_e[:, i, j], [2.5, 25, 50, 75, 97.5])

        fout.write('\n    {0}×{1}\n'.format(

                    NameR[i], NameC[j]))

        fout.write('      平均値 = {0:>.3f}\n'.format(ga_mean))

        fout.write('      中央値 = {0:>.3f}\n'.format(ga_50p))

        fout.write('      Q1 = {0:>.3f},   Q3 = {1:>.3f}\n'.format(ga_25p, ga_75p))

        fout.write('      95%確信区間 = [{0:>.3f}, {1:>.3f}]\n'.format(ga_025p, ga_975p))

       

plt.yticks([])

plt.legend()

plt.tight_layout()

plt.savefig(foutnm[:-4] + '_rca.png')

SavedFiles.append(foutnm[:-4] + '_rca.png')

plt.show()

 

fout.write('\n\nλijの中央値:\n')

lmbd_med = np.empty((KR, KC))

for i in range(KR):

    for j in range(KC):

        lmbd_med[i][j] = np.percentile(fit['lmbd'][:, i, j], 50)

        fout.write('\nλ[{0}, {1}] = {2:.3f}\n'.format(NameR[i], NameC[j], lmbd_med[i][j]))

 

plt.title('$\hat\lambda_{ij}$の中央値', fontsize = 20)

for i in range(KR):

    plt.plot(range(KC), lmbd_med[i], label = NameR[i])

plt.ylabel('$\hat\lambda_{ij}$')

plt.xticks(range(KC), NameC)

plt.legend()

plt.savefig(foutnm[:-4] + '.lr.png')

SavedFiles.append(foutnm[:-4] + '_lr.png')

plt.show()

 

plt.title('$\hat\lambda_{ij}$の中央値', fontsize = 20)

for j in range(KC):

    plt.plot(range(KR), lmbd_med[:, j], label = NameC[j])

plt.ylabel('$\hat\lambda_{ij}$')

plt.xticks(range(KR), NameR)

plt.legend()

plt.savefig(foutnm[:-4] + '_pc.png')

SavedFiles.append(foutnm[:-4] + '_pc.png')

plt.show()

 

fout.close()

SavedFiles.append(ftxtoutnm)

 

print('\n保存ファイル：')

for nm in SavedFiles:

    print('    ', nm)

print()

 

 

Up