Up

比率のクロス表のベイズ分析

ロジスティック回帰分析

 

京都大学のホームページに掲載されている平成19年度と平成29年度の一般入試諸統計から合格者が200名以上の学部の合格者数内数（「認」と「他」を除く）を男女別にまとめると図１のようになった。

図１　Excelデータ（data1929.xlsx）

 

女子の比率をロジスティック回帰モデルで分析してみる。ロジスティック回帰モデルについては、＜岡本安晴「心理学データ分析と測定」勁草書房＞第４章で説明している。

表１のように各条件（行i、列j）における注目している事象の確率（図１の場合、女子の確率）をpijで表す。

 

表１　行条件iと列条件jにおける確率pij

 

列１

・・・

列ｊ

・・・

行１

p11

・・・

・・・

行i

・・・

・・・

 

行条件iと列条件jの組み合わせにおける下段の度数（図１では女子の数）が総数（図１では女子の数と男子の数の和）の２項分布における度数として確率pijの推定を行う。確率pijは、比率「女子の数／（女子の数＋男子の数）」に対応するが、ベイズ分析を行うために、総数「女子の数＋男子の数」に対する該当事象の度数が「女子の数」である２項分布を適用してpijの分析を行う。

線形回帰モデルを用いるために、確率pijの値域［０、１］をロジット変換した値gijについて考える（Kruschke,2015）。

である。

ロジット変換値gijを行主効果ai、列主効果bj、交互作用(ab)jkの和に分解する。

ここで、制約条件を

と設定する。

このとき、

である。KRとKCは行条件の数と列条件の数を表す。

ロジスティック回帰モデルでは、式（１）より

となる。

各セルにおける度数は、確率pijの２項分布に従うとして、本ウェブサイト後半に示すリスト１のようにStanスクリプトを用意した。リスト１のStanスクリプトを利用して表１のベイズ分析を行うPythonスクリプトは、本ウェブサイト後半に示すリスト２として用意した。これらのファイルなどは圧縮ファイルfilesLogisticReg.zipとしてまとめたので、ダウンロード解凍すれば自由に利用できる。PythonにおけるStanの使い方については、＜岡本安晴「いまさら聞けない1Pythonでデータ分析」丸善出版＞の第10章～第12章で説明している。

リスト１ではパラメータgijの事後分布が求められるので、リスト２において主効果および交互作用を式（４）～（７）に基づく以下のスクリプトで求めている。

 

n_samples = len(fit['g'])

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 vg in fit['g']:

    gdd = vg.mean()             #   g..

    gdc = vg.mean(axis = 0)     #   g.j

    grd = vg.mean(axis = 1)     #   gi.

    r_e[i] = grd - gdd          #   行主効果

    c_e[i] = gdc - gdd          #   列主効果

    rc_e[i] = vg - grd.reshape((KR, 1)) - gdc + gdd     #   交互作用

    i += 1

 

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

図２

 

入力データファイルは、図１の形式で用意する。第1行目に、列条件のカテゴリ（図１の場合は、学部名）を書き、第１列目に行条件のカテゴリ（図１では年度）、第２列目に度数のカテゴリを書く。度数データは２行目以降、３列目から該当するセルに書く。各行条件・列条件に対応するセルは、上に非該当カテゴリ（図１では「男子」、２項分布の説明では「０」で表されることが多い）の度数、下の段に該当カテゴリ（図１では「女子」、２項分布の説明では「１」で表されることが多い）の度数を書く。列条件あるいは行条件は図１のものより多くてもよい。列および行条件の個数は２以上であればよい。

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

図３

 

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

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

図４

 

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

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

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

図５

 

平成19年度に比べて平成29年度の分布が負の方向に移動していることから、女子の比率に減少の傾向が窺われる。

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

図６

 

列の主効果から女子の比率は、文学部が最も高く、以下、農学部、法学部、理学部、工学部の順に女子の比率が下がっていることが分かる。理学部と工学部の女子の比率が他の学部より低いことが明瞭に表されている。

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

図７

 

交互作用の事後分布の両側に位置しているのは、平成29年×法学部と兵籍19年×法学部の交互作用の事後分布であるが、これは主効果で説明される効果以上に農学部の女子の比率が低い状態から高い状態に増加したことを示している。その次の内側の事後分布の組は平成19年×文学部と平成29年×文学部の交互作用である。これは、主効果で説明される値を基準としたとき、文学部はそれ以上に女子の高い比率から低い比率に変化したことを示している。

図７のウィンドウを閉じると、次にpijの事後分布の中央値のグラフが表示される（図８）。

図８

 

平成19年度と平成29年度を比べたとき、文学部の女子の比率が減少して、法学部の女子の比率が増加していることが分かる。これは、交互作用のグラフ（図７）に対応している。

 

図９

 

図８のウィンドウを閉じると、次に横軸に平成19年と平成29年（行条件）を取ったpijの事後分布の中央値のグラフが表示される。文学部の女子の比率の減少と法学部の女子の比率の増加が表されいる。農学部の女子の比率の減少が見られるが、これは交互作用としては顕著ではない（図７）。

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

図10　実行終了時の端末画面

 

実行終了後、図４で指定したテキストファイルををPythonのIDLEのエディタで開こうとすると、図８のようなダイアログボックスが表示されることがある。

図８

 

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

図９、

 

出力テキストファイルは、以下のような内容である。事後分布の平均値、中央値、四分位数、95％CIが出力されている。

 

入力データファイル = D:/yasuharu/www/LaCoocan/Python/IntroPyStat/BysLogistic/filesLogisticReg/data1929.xlsx

 

行変数：['平成19年', '平成29年']

列変数：['文学部' '法学部' '理学部' '工学部' '農学部']

 

 

P = 女子 / (男子 + 女子)

 

行主効果

 

    平成19年

      平均値 = 0.037

      中央値 = 0.038

      Q1 = 0.007,   Q3 = 0.067

      95%確信区間 = [-0.050, 0.125]

 

    平成29年

      平均値 = -0.037

      中央値 = -0.038

      Q1 = -0.067,   Q3 = -0.007

      95%確信区間 = [-0.125, 0.050]

 

列主効果

 

    文学部

      平均値 = 1.074

      中央値 = 1.076

      Q1 = 1.016,   Q3 = 1.131

      95%確信区間 = [0.909, 1.233]

 

    法学部

      平均値 = 0.169

      中央値 = 0.170

      Q1 = 0.115,   Q3 = 0.223

      95%確信区間 = [0.005, 0.329]

 

    理学部

      平均値 = -0.811

      中央値 = -0.808

      Q1 = -0.881,   Q3 = -0.740

      95%確信区間 = [-1.020, -0.609]

 

    工学部

      平均値 = -1.130

      中央値 = -1.128

      Q1 = -1.182,   Q3 = -1.077

      95%確信区間 = [-1.285, -0.979]

 

    農学部

      平均値 = 0.698

      中央値 = 0.696

      Q1 = 0.645,   Q3 = 0.753

      95%確信区間 = [0.541, 0.861]

 

交互作用

 

    平成19年×文学部

      平均値 = 0.124

      中央値 = 0.125

      Q1 = 0.064,   Q3 = 0.181

      95%確信区間 = [-0.050, 0.300]

 

    平成19年×法学部

      平均値 = -0.198

      中央値 = -0.199

      Q1 = -0.257,   Q3 = -0.141

      95%確信区間 = [-0.365, -0.031]

 

    平成19年×理学部

      平均値 = 0.047

      中央値 = 0.048

      Q1 = -0.024,   Q3 = 0.119

      95%確信区間 = [-0.166, 0.262]

 

    平成19年×工学部

      平均値 = -0.046

      中央値 = -0.045

      Q1 = -0.098,   Q3 = 0.007

      95%確信区間 = [-0.199, 0.105]

 

    平成19年×農学部

      平均値 = 0.074

      中央値 = 0.075

      Q1 = 0.020,   Q3 = 0.126

      95%確信区間 = [-0.081, 0.231]

 

    平成29年×文学部

      平均値 = -0.124

      中央値 = -0.125

      Q1 = -0.181,   Q3 = -0.064

      95%確信区間 = [-0.300, 0.050]

 

    平成29年×法学部

      平均値 = 0.198

      中央値 = 0.199

      Q1 = 0.141,   Q3 = 0.257

      95%確信区間 = [0.031, 0.365]

 

    平成29年×理学部

      平均値 = -0.047

      中央値 = -0.048

      Q1 = -0.119,   Q3 = 0.024

      95%確信区間 = [-0.262, 0.166]

 

    平成29年×工学部

      平均値 = 0.046

      中央値 = 0.045

      Q1 = -0.007,   Q3 = 0.098

      95%確信区間 = [-0.105, 0.199]

 

    平成29年×農学部

      平均値 = -0.074

      中央値 = -0.075

      Q1 = -0.126,   Q3 = -0.020

      95%確信区間 = [-0.231, 0.081]

 

 

引用文献

Kruschke, J. K. (2015). Doing Bayesian Data Analysis, 2nd ed. Academic Press.

 

 

リスト１　ロジスティック回帰分析用Stanスクリプト（LogisModel.stan）

 

data {

    int KR;

    int KC;

    int X[KR, KC, 2];

}

parameters {

    real g[KR, KC];

}

transformed parameters {

    real<lower = 0.0, upper = 1.0> p[KR, KC];

    for (i in 1:KR)

        for (j in 1:KC)

            p[i][j] = inv_logit(g[i][j]);

}

model {

    for (i in 1:KR)

        for (j in 1:KC) {

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

    }

    for (i in 1:KR)

        for (j in 1:KC) {

            X[i][j][2] ~ binomial(X[i][j][1] + X[i][j][2], p[i][j]);       

        }

}

 

 

リスト２　リスト１のStanスクリプトを用いるPythonスクリプト（MainLogistic.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

 

print(data)

 

KC = len(data[0]) - 2

print(KC)

KR = (len(data) - 1) // 2

print(KR)

 

NameC = data[0][2:]

NameR = []

for i in range(KR):

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

print(NameC)

print(NameR)

print('\nP = {0} / ({1} + {0})\n'.format(data[2][1], data[1][1]))

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

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

fout.write('\n\nP = {0} / ({1} + {0})\n'.format(data[2][1], data[1][1]))

 

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

for i in range(KR):

    for j in range(KC):

        X[i][j][0] = data[i*2 + 1][j + 2]

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

 

for i in range(KR):

    for j in range(KC):

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

    print()

    for j in range(KC):

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

    print()

 

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

 

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

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

print(fit)

 

SavedFiles = []

 

#

#       主効果と交互作用

#

n_samples = len(fit['g'])

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 vg in fit['g']:

    gdd = vg.mean()             #   g..

    gdc = vg.mean(axis = 0)     #   g.j

    grd = vg.mean(axis = 1)     #   gi.

    r_e[i] = grd - gdd          #   行主効果

    c_e[i] = gdc - gdd          #   列主効果

    rc_e[i] = vg - grd.reshape((KR, 1)) - gdc + gdd     #   交互作用

    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()

 

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

for i in range(KR):

    for j in range(KC):

        p_med[i][j] = np.percentile(fit['p'][:, i, j], 50)

 

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

for i in range(KR):

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

plt.ylabel('$\hatp_{ij}$')

plt.xticks(range(KC), NameC)

plt.legend()

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

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

plt.show()

 

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

for j in range(KC):

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

plt.ylabel('$\hatp_{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