見出し画像

第18章「"覚えやすさ"の判断メカニズム」のベイズモデリングをPyMC Ver.5 で

ネイピア DS

この記事は、テキスト「たのしいベイズモデリング２」の第18章「"覚えやすさ"の判断メカニズム」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。

今回も心理学術的な色合いの濃いベイズ分析です。
この章は言葉・用語の覚えやすさを４段階（簡単～難しい）に当てはめる際に、どのようなロジック（経路）で判断するかを項目反応ツリーモデルで分析します。

日本語を勉強する外国人のイラスト：「いらすとや」さんより

テキストは brms を用いています。
R の回帰式記述法である formula を用いて Stan モデルを描けるライブラリだそうです。
そこで今回の PyMC 化は formula で記述できる「Bambi」を使っていきます！

Bambiの推論結果はテキストとまあまあ近い結果になりました。
ではでは PyMCのベイズモデリングの世界を楽しんでまいりましょう！

テキストの紹介、引用表記、シリーズまえがき、PyMC等のバージョン情報は、このリンクの記事をご参照ください。

テキストで使用するデータは、R・Stan等のサンプルスクリプトとともに、出版社サイトからダウンロードして取得できます。


サマリー

テキストの概要

執筆者　　　：　山根嵩史 先生
モデル難易度：　★★★・・　（ふつう）

自己評価

評点

$$
\begin{array}{c:c:c}
実装精度 & ★★★★・ & ほぼ納得 \\
結果再現度 & ★★★★・ & あと一歩 \\
楽しさ & ★★★★★ & 楽しい！ \\
\end{array}
$$

花型の評価印のイラスト：「いらすとや」さんより

評価ポイント

  • 項目反応ツリーモデル自体を理解するのに時間がかかったので、モデル難易度を少し高めにしました。
    brms の formula に記述された固定効果・変量効果を織り込む回帰式は平易なのですが、「ツリー構造＋マッピング行列」と「データの中身・変数の値の持ち方」が頭の中で繋がらなかったのです。
    R 環境で R スクリプトを動かして、変数の値の変化を確認しながら、ツリー構造とマッピング行列との関係性を把握するように努めました。

工夫・喜び・反省

  • 今回は pymc をインポートしていません！
    Bambi と arviz だけで PyMC モデリング＆推論を実践しました。 

子鹿・バンビのイラスト：「いらすとや」さんより

モデルの概要

テキストの調査・実験の概要

■ 実験の概要
実験参加者が提示された「ひらがな３文字」の文字列の覚えやすさを「１：簡単～４：難しい」の４段階で回答する実験です。

  • 実験に用いられる文字列には次の２種類があります。

    • 有意味条件（M）
      「すいか」「さくら」のように意味を持つ「有意味語」を使用

    • 無意味条件（N）
      「とぱる」「しれぎ」のように意味無しの「無意味綴り」を使用

  • 実験参加者数は有意味条件10名、無意味条件10名です。

  • 実験参加者はアサインされた条件で20項目に回答します。

データの数は有意味条件・無意味条件それぞれで、実験参加者数１０×回答項目数２０＝２００ずつあります。

実験の詳細を知りたい方はこちらの先生の論文（PDFダウンロード）をお読みください！

すいかのイラスト「丸々すいかとカットすいか」：「いらすとや」さんより

■ 分析の概要
実験参加者が「提示された文字列が覚えやすさ４段階のどれに該当するか判断」するときに「４段階の選択にどのように反応するか」に興味があります。
テキストは反応の仕方に関して、２つの項目反応ツリーモデルを検討します。

① 線形反応ツリーモデル（LRT）
１段階目の「簡単」から順を追って当てはまる段階を吟味・選択するモデルです。
楕円で表した要素を「ノード」と呼びます。

② 入れ子反応ツリーモデル（NRT）
おおむね簡単か難しいかを先に吟味・判断して、その後、簡単の度合い（簡単かやや簡単）・難しいの度合い（やや難しいか難しい）を吟味・判断する２ステップのモデルです。

テキストはこの２つのツリーモデルの判断ステップを仮定して、覚えやすさの判断メカニズムを探ります。
またモデリングにあたって、ツリーモデルはデータの持ち方で表現されます。

一旦まとめます。

興味のあること
覚えやすさ＝学習のしやすさ＝学習容易性の判断メカニズムをモデル化してパラメータを推定します。

「覚えやすさ」に用いる材料
意味を持つ単語の覚えやすさを判断する「有意味条件」と意味を持たない単語の覚えやすさを判断する「無意味条件」の２つの条件を取り扱います。

想定するメカニズム
判断メカニズムには項目反応ツリーモデルである「線形反応ツリーモデル」と「入れ子反応ツリーモデル」の２つのモデルを「データの持ち方」にて扱います。

■ マッピング行列
このあたりからややこしくなります。
項目反応ツリーモデルの樹形図を行列形式に変換します。
選択した回答を行に、楕円形のノードを列に、樹形図の線に付した選択肢を値に持つ行列です。

先に掲載した２つのツリーモデルを一般化して再掲し、マッピング行列に変換します。
テキストの図18.1、図18.2に相当します。

① 線形反応ツリーモデルとマッピング行列
楕円形のノードの中には「簡単？」などの選択肢が入っていると想像しましょう。
矢印の最終的な行先である「Y=…」の部分が「１：簡単～４：難しい」から選んだ回答値です。
矢印に付した値０・１は各ノードで判断された値に相当します。
各ノードにおいて、０は易しい方への判断、１は難しい方への判断と言えるようです。

見方の例です。
ツリーの回答値「Y=2」（やや簡単）の経路に注目します。
最初のノード「Y1*」の判断は１（難しい方へ進む）、続くノード「Y2*」の判断は０（易しい方へ進む）です。
マッピング行列に視線を移します。
「Y=2」について、「Y1*」の値は１、「Y2*」の値は０です。
マッピング行列はツリーの経路情報を持っているのです。
列の左側から順に値を見ていく感じです。
なお、Y3* へは行かないので値は「－」です。

このマッピング行列の値 0, 1 がベイズモデリングの目的変数になります。

② 入れ子反応ツリーモデルとマッピング行列

分析の概要を図示します。

テキストのモデリング

テキスト記載のモデル数式と brms の formula 記述の数式をピタリと合わせて説明することが困難なので、ここでは一旦、テキスト掲載のモデル数式を書きます。

■ 目的変数と関心のあるパラメータ
目的変数は、実験参加者$${p}$$が文字列の項目番号$${i}$$に対する回答をする際、ノード$${n}$$で選んだ経路の値$${Y^*_{pin}}$$です。
$${Y^*_{pin}}$$にはマッピング行列の値がセットされます。
例えば、線形反応ツリーモデルを仮定する場合、ID＝１の実験参加者が項目番号３の文字列を見て、番号２のノード「Y3*」で経路０を選んだ場合、$${Y^*_{1,3,2}=0}$$です。

関心のあるパラメータは項目$${Item}$$の固定効果を示す回帰係数$${\beta_i}$$、ノード$${Node}$$の固定効果を示す回帰係数$${\beta_n}$$などです。

■ モデル数式
１行目の尤度関数はベルヌーイ分布を用いています。
２行目のベルヌーイ分布の確率パラメータは、ノード、項目、実験参加者を説明変数とするロジスティック回帰モデルを用いています。

$$
\begin{align*}
Y^*_{pin} &\sim \text{Bernoulli}\ (\theta_{pin}) \\
\log \left(\cfrac{\theta_{pin}}{1 - \theta_{pin}} \right) &= \beta_0 + \beta_n \times Node + \beta_i \times Item + \beta_p \times Person \\
\end{align*}
$$

テキストより引用

テキストは「典型的に項目$${Item}$$は固定効果、実験参加者$${Person}$$は変量効果としてモデルに投入される」と補足しています。

ファッションモデルのイラスト（女性）：「いらすとや」さんより

■ ３つのベイズモデル
テキストは各ノードにおける判断が似ているのか異なるのかを確認する目的で、３つのベイズモデルを構築します。

$$
\begin{array}{l:l}
モデル1 & 切片のみ変量効果を指定 \\\\&全てのノードにおいて\\&共通の判断が行われると仮定 \\\\
\hdashline
モデル2 & ノードごとのランダム傾きを投入 \\\\ & 各ノードにおいて\\&異なる判断が行われると仮定 \\\\
\hdashline
モデル3 & ノードY^*_1とノードY^*_2およびノードY^*_3で \\ &異なるランダム傾きを投入  \\\\ &ノードY^*_1とノードY^*_2およびノードY^*_3において\\
&異なる判断が行われると仮定 \\\\
\end{array}
$$

モデル１は、２ツリーモデル × ２条件＝４ベイズモデルを扱います。
モデル２と３は、入れ子反応ツリーモデル × ２条件 × ２ベイズモデル＝４ベイズモデルを扱います。
全部で８つのベイズモデリングを構築します！大量です！

こぼれいくら丼のイラスト：「いらすとや」さんより

■分析・分析結果
分析の仕方や分析数値はテキストの記述が正確だと思いますので、テキストの読み込みをおすすめします。
PyMCの自己流モデルの推論値を利用した分析は「PyMC実装」章をご覧ください。

PyMC実装

Let's enjoy PyMC & Python !

準備・データ確認

１．インポート

### インポート

# 数値・確率計算
import pandas as pd
import numpy as np

# PyMC
import bambi as bmb
import arviz as az

# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'

# グラフ描画
from graphviz import Digraph

# ユーティリティ
import pickle

# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')

２．説明用図表の作成
「モデルの概要」章で用いる図表を作成します。
graphviz ライブラリの有向グラフ Digraph を使用します。

① 線形反応ツリーモデルのイメージ図

### 線形反応ツリーモデルのイメージ by graphviz

## 設定 neatoは頂点の位置を指定できる
g = Digraph(engine='neato')

## node：頂点の作成, color='white'で輪郭線を消す
g.node('簡単？', pos='2, 3!')
g.node('やや簡単？', pos='3, 2!')
g.node('やや難しい？', pos='4, 1!')
g.node('簡単', pos='0, 0!', shape='box')
g.node('やや簡単', pos='1.5, 0!', shape='box')
g.node('やや難しい', pos='3, 0!', shape='box')
g.node('難しい', pos='5, 0!', shape='box')

# ## edge：辺の作成, labelで辺にラベル値を表示
g.edge('簡単？', '簡単', label='易')
g.edge('簡単？', 'やや簡単？', label='難')
g.edge('やや簡単？', 'やや簡単', label='易')
g.edge('やや簡単？', 'やや難しい？', label='難 ')
g.edge('やや難しい？', 'やや難しい', label='易')
g.edge('やや難しい？', '難しい', label='難 ')

## グラフの描画
g

【実行結果】

② 入れ子反応ツリーモデルのイメージ

### 入れ子反応ツリーモデルのイメージ by graphviz

## 設定 neatoは頂点の位置を指定できる
g = Digraph(engine='neato')

## node：頂点の作成, color='white'で輪郭線を消す
g.node('簡単？難しい？', pos='3, 3!')
g.node('簡単さ？', pos='1, 1.5!')
g.node('難しさ？', pos='5, 1.5!')
g.node('簡単', pos='0, 0!', shape='box')
g.node('やや簡単', pos='2, 0!', shape='box')
g.node('やや難しい', pos='4, 0!', shape='box')
g.node('難しい', pos='6, 0!', shape='box')

# ## edge：辺の作成, labelで辺にラベル値を表示
g.edge('簡単？難しい？', '簡単さ？', label='易')
g.edge('簡単？難しい？', '難しさ？', label='難      ')
g.edge('簡単さ？', '簡単', label='易')
g.edge('簡単さ？', 'やや簡単', label='難   ')
g.edge('難しさ？', 'やや難しい', label='易')
g.edge('難しさ？', '難しい', label='難   ')

## グラフの描画
g

【実行結果】

③ 図18.1の線形反応ツリーモデルとマッピング行列

### 線形反応ツリーモデルの可視化 by graphviz ※図18.1

## 設定 neatoは頂点の位置を指定できる
g = Digraph(engine='neato')

## node：頂点の作成, color='white'で輪郭線を消す
g.node('Y1*', pos='2, 3!')
g.node('Y2*', pos='3, 2!')
g.node('Y3*', pos='4, 1!')
g.node('Y=1', pos='0, 0!', color='white')
g.node('Y=2', pos='1.5, 0!', color='white')
g.node('Y=3', pos='3, 0!', color='white')
g.node('Y=4', pos='5, 0!', color='white')

# ## edge：辺の作成, labelで辺にラベル値を表示
g.edge('Y1*', 'Y=1', label='0')
g.edge('Y1*', 'Y2*', label='1')
g.edge('Y2*', 'Y=2', label='0')
g.edge('Y2*', 'Y3*', label='1')
g.edge('Y3*', 'Y=3', label='0')
g.edge('Y3*', 'Y=4', label='1')

## グラフの描画
g

【実行結果】

### 線形反応ツリーモデルのマッピング行列のイメージ ※図18.1
vals = [[0, '-', '-'], [1, 0, '-'], [1, 1, 0], [1, 1, 1]]
index = ['Y=1', 'Y=2', 'Y=3', 'Y=4']
columns = ['Y1*', 'Y2*', 'Y3*']
pd.DataFrame(vals, index=index, columns=columns)

【実行結果】

③ 図18.2の入れ子反応ツリーモデルとマッピング行列

### 入れ子反応ツリーモデルの可視化 by graphviz ※図18.2

## 設定 neatoは頂点の位置を指定できる
g = Digraph(engine='neato')

## node：頂点の作成, color='white'で輪郭線を消す
g.node('Y1*', pos='2, 3!')
g.node('Y2*', pos='1, 1.5!')
g.node('Y3*', pos='3, 1.5!')
g.node('Y=1', pos='0, 0!', color='white')
g.node('Y=2', pos='1.5, 0!', color='white')
g.node('Y=3', pos='2.5, 0!', color='white')
g.node('Y=4', pos='4, 0!', color='white')

# ## edge：辺の作成, labelで辺にラベル値を表示
g.edge('Y1*', 'Y2*', label='0')
g.edge('Y1*', 'Y3*', label='1   ')
g.edge('Y2*', 'Y=1', label='0')
g.edge('Y2*', 'Y=2', label='1  ')
g.edge('Y3*', 'Y=3', label='0')
g.edge('Y3*', 'Y=4', label='1   ')

## グラフの描画
g

【実行結果】

### 線形反応ツリーモデルのマッピング行列のイメージ ※図18.2
vals = [[0, 0, '-'], [0, 1, '-'], [1, '-', 0], [1, '-', 1]]
index = ['Y=1', 'Y=2', 'Y=3', 'Y=4']
columns = ['Y1*', 'Y2*', 'Y3*']
pd.DataFrame(vals, index=index, columns=columns)

【実行結果】

３．データの読み込みと前処理
２つのcsvファイルををpandasのデータフレームに読み込みます。
・有意味条件の回答データ「EOLdatM.csv」
・無意味条件の回答データ「EOLdatN.csv」

### データの読み込み
data_m_orgn = pd.read_csv('EOLdatM.csv')  # 有意味条件
data_n_orgn = pd.read_csv('EOLdatN.csv')  # 無意味条件

print('【有意味条件】')
display(data_m_orgn)
print('\n【無意味条件】')
display(data_n_orgn)

【実行結果】
行が実験参加者ID単位、列が回答項目単位、値は「１：簡単～４：難しい」の中から選択した回答値です。

マッピング行列を用いたデータ前処理に進みます。
最初にマッピング行列を定義します。

### マッピング行列の作成
# 線形反応ツリーモデルのマッピング行列
lrt_map = np.array([[0, np.nan, np.nan],
                    [1, 0, np.nan],
                    [1, 1, 0],
                    [1, 1, 1]])
# 入れ子反応ツリーモデルのマッピング行列
nrt_map = np.array([[0, 0, np.nan],
                    [0, 1, np.nan],
                    [1, np.nan, 0],
                    [1, np.nan, 1]])
# 表示
print('線形反応ツリーモデルのマッピング行列')
print(lrt_map)
print('\n入れ子反応ツリーモデルのマッピング行列')
print(nrt_map)

【実行結果】

続いて、データ変換関数の定義です。
R の Dendrify 関数に相当する Python のライブラリ・関数を見つけられなかったので、自作しました。
引数は、変換対象データ（横持ちデータ）：df、マッピング行列：map_mtx です。
戻り値は、マッピング行列の値を設定した縦持ちデータです。

### データの前処理：マッピング行列を用いたデータ変換処理を関数化

# 引数 df:csvを読み込んだデータフレーム, map_mtx:マッピング行列
def make_data(df, map_mtx):
    # 変換後データを格納する一時リストの初期化
    data_list = []
    # 項目i（列）ごとに繰り返し処理
    for i in range(df.iloc[:, 1:].shape[1]):
        # 回答者p（行）ごとに繰り返し処理
        for p in range(df.iloc[:, 1:].shape[0]):
            # 回答値に合致するマッピング行列の列ごとに繰り返し処理
            for j, n in enumerate(map_mtx[df.iloc[p, i + 1] - 1]):
                # マッピング行列の値が0か1の時に変換後データを作成
                if (n==0) or (n==1):
                    # 変換後データを一時リストに追加
                    # [value, item, person, node, sub]
                    data_list.append([int(n), f'item{i+1:02}', f'person{p+1:02}',
                                      f'node{j+1}', f'item{i+1:02}:node{j+1}'])
    # 一時リストを戻り値用のデータフレームに設定
    res_df = pd.DataFrame(np.array(data_list))
    # データフレームの列名の変更
    res_df.columns = ['value', 'item', 'person', 'node', 'sub']
    # valueを整数型に変換 <---- ★★１日悩んだ結果の解決策コード：整数化★★
    res_df['value'] = res_df['value'].astype('int')
    
    return res_df

データ変換に進みます。
有意味条件・無意味条件 × 線形反応ツリーモデル・入れ子反応ツリーモデルの４つのデータを変換・作成します。

【ツリーモデルを仮定したデータ】
２つのツリーモデルの構造を踏襲したマッピング行列を用いてデータ作成します。
したがいまして、特定のツリーモデルを仮定してデータ作成することになります。
データがいずれかのツリーモデルを織り込み済みであるということです。

① 有意味条件 × 線形反応ツリーモデルLRT

### データの前処理：マッピング行列を用いたデータ変換　有意味条件・LRT
data_m_lrt = make_data(data_m_orgn, lrt_map)
display(data_m_lrt)

【実行結果】
データは、実験参加者 person が文字列 item の回答の際にノード番号 node で選択した値 value（0 か 1）を示しています。
表の見方です。
・value：マッピング行列の値であり、目的変数
・item：回答項目（20個の３文字の文字列の識別子）
・person：実験参加者ID
・node：反応ツリーモデルのノード番号
item、person、node はカテゴリ変数です。

② 無意味条件 × 線形反応ツリーモデルLRT

### データの前処理：マッピング行列を用いたデータ変換　無意味条件・LRT
data_n_lrt = make_data(data_n_orgn, lrt_map)
display(data_n_lrt)

【実行結果】

③ 有意味条件 × 入れ子反応ツリーモデルNRT

### データの前処理：マッピング行列を用いたデータ変換　有意味条件・NRT
data_m_nrt = make_data(data_m_orgn, nrt_map)
display(data_m_nrt)

【実行結果】

④ 無意味条件 × 入れ子反応ツリーモデルNRT

### データの前処理：マッピング行列を用いたデータ変換　無意味条件・NRT
data_n_nrt = make_data(data_n_orgn, nrt_map)
display(data_n_nrt)

【実行結果】

４パターンの繰り返し処理が続きました。お疲れ様です。
これからもこの複数パターンによる繰り返し処理が続きますので、よろしくお願いいたします！

モデル１の構築

テキストの 18.3 節「データとモデルの当てはめ」、18.4 節「結果と考察」で用いる「モデル１：切片のみ変量効果を指定したモデル」の４つのモデリングを行います。
分析対象データが異なる点を除くと、モデルの内容は全て同じです。

【分析の視点】
パラメータの事後分布サンプルデータを用いて分析することは、「項目の固定効果」と「ノードの固定効果」です。
項目の固定効果は、個々の文字列の項目ごとの係数を用いて、有意味条件（意味のある語）と無意味条件（意味のない綴り）の違いを分析します。
ノードの固定効果は、ノード単位の係数の要約統計量を用いて、線形反応ツリーモデルと入れ子反応ツリーモデルの違い、有意味条件と無意味条件の違いを分析します。

モデルの数式表現
目指したいBambiのモデルの「なんちゃって数式」表記です。
formula 形式で書きます。

$$
value \sim 0 + node + item + (1\ |\ person)
$$

実験参加者$${person}$$のランダム切片モデルです。
目的変数$${value}$$はベルヌーイ分布に従います。

では４つのベイズモデリングを開始します！

１．有意味条件・線形反応ツリーモデル

(1) モデルの定義

### モデルの定義
model_m_lrt = bmb.Model(
        formula='value ~ 0 + node + item +  (1 | person)',    # フォーミュラ式
        data=data_m_lrt,                                 # データ
        family='bernoulli',                              # 目的変数の誤差分布
        # link='logit',                                  # （省略可）リンク関数
        # categorical=['item', 'person', 'node'],        # （省略可）カテゴリ変数
)

【モデル注釈】
family に目的変数の誤差分布であるベルヌーイ分布を指定します。
リンク関数 link はベルヌーイ分布のデフォルト値 logit を使用します。
以後の３モデルでは 使用するデータ data の指定が変わるだけです。

(2) モデルの外観の確認

### モデルの表示
print(model_m_lrt)

【実行結果】
ノード・項目の固定効果（回帰係数）が従う事前分布、実験参加者の変量効果が従う事前分布などは、Bambiがデータを解析して自動設定してくれます。

### モデルの可視化
model_m_lrt.build()
model_m_lrt.graph()

【実行結果】
Bambiのモデルをグラフ形式で表示するには、事前にモデルの構築「build」かモデルの学習「fit」をしておく必要があります。

(3) 事後分布からのサンプリング
Bambi では、モデルに対して fit を行うことでMCMCを実行できます。
分析時にモデルの評価指標 WAIC と loo を計算するので、MCMC実行時に対数尤度を算出するように「idata_kwargs={'log_likelihood': True}」を指定します。
乱数生成数はテキストよりも少ないです。
nuts_sampler='numpyro'とすることで、numpyroをNUTSサンプラーに利用できます。
処理時間はおよそ1分でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　1分
# テキスト: iter=100000, warmup=50000, chains=4, thin=10
idata_m_lrt = model_m_lrt.fit(
    draws=1000, tune=1000, chains=4, target_accept=0.85, nuts_sampler='numpyro', 
    idata_kwargs={'log_likelihood': True}, random_seed=1969)

【実行結果】省略

(4) サンプリングデータの確認
$${\hat{R}}$$、トレースプロットを確認します。
事後分布の収束確認は$${\hat{R} \leq 1.1}$$とします。
ここでは$${\hat{R} > 1.01}$$のパラメータが無いことを確認します。

### r_hat>1.1の確認
# 設定
idata_in = idata_m_lrt   # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
$${\hat{R} > 1.01}$$のパラメータは０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

### 推論データの要約統計情報の表示
az.summary(idata_m_lrt, hdi_prob=0.95, round_to=3)

【実行結果】

トレースプロットで事後分布サンプリングデータの様子を確認します。

### トレースプロットの表示
az.plot_trace(idata_m_lrt, compact=True)
plt.tight_layout();

【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。

(5) 推論データの保存
推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
pickleで保存します。

### idataの保存 pickle
file = r'idata_m_lrt_ch18.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata_m_lrt, f)

読み込みコードは次のとおりです。

### idataの読み込み pickle
file = r'idata_m_lrt_ch18.pkl'
with open(file, 'rb') as f:
    idata_m_lrt_load = pickle.load(f)

(6) 予習
分析時に用いる項目のパラメータ推定値のフォレストプロットの描画、 WAIC、loo の算出を練習しましょう。

① フォレストプロットの描画
arviz の plot_forest 関数を利用します。
予めデータの平均値で昇順ソートしておきます。

### フォレストプロットの描画 ※図18.4に相当

# 平均値で降順ソートするためのデータ準備
item_means = idata_m_lrt.posterior['item'].mean(('chain', 'draw'))
sorted_items = (idata_m_lrt.posterior['item_dim']
                .sortby(item_means, ascending=False))
# フォレストプロットの描画
az.plot_forest(idata_m_lrt, var_names=['item'], coords={'item_dim': sorted_items},
               combined=True, hdi_prob=0.95,figsize=(4, 4))
plt.axvline(0, color='red', ls='--')
plt.xlim([-8, 8])
plt.grid(lw=0.5);

【実行結果】

② WAIC の算出
arviz の waic 関数を利用します。

### WAICの算出
az.waic(idata_m_lrt, scale='deviance')

【実行結果】

③ loo の算出
arviz の loo 関数を利用します。

### LOOの算出
az.loo(idata_m_lrt, scale='deviance')

【実行結果】

２．有意味条件・入れ子反応ツリーモデル

(1) モデルの定義

### モデルの定義
model_m_nrt = bmb.Model(
        formula='value ~ 0 + node + item +  (1 | person)',    # フォーミュラ式
        data=data_m_nrt,                                 # データ
        family='bernoulli',                              # 目的変数の誤差分布
)

(2) モデルの外観の確認

### モデルの表示
print(model_m_nrt)

【実行結果】

### モデルの可視化
model_m_nrt.build()
model_m_nrt.graph()

【実行結果】

(3) 事後分布からのサンプリング
処理時間はおよそ15秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　15秒
# テキスト: iter=100000, warmup=50000, chains=4, thin=10
idata_m_nrt = model_m_nrt.fit(
    draws=1000, tune=1000, chains=4, target_accept=0.95, nuts_sampler='numpyro', 
    idata_kwargs={'log_likelihood': True}, random_seed=1969)

【実行結果】省略

(4) サンプリングデータの確認
$${\hat{R}}$$、トレースプロットを確認します。
事後分布の収束確認は$${\hat{R} \leq 1.1}$$とします。
ここでは$${\hat{R} > 1.01}$$のパラメータが無いことを確認します。

### r_hat>1.1の確認
# 設定
idata_in = idata_m_nrt   # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
$${\hat{R} > 1.01}$$のパラメータは０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

### 推論データの要約統計情報の表示
az.summary(idata_m_nrt, hdi_prob=0.95, round_to=3)

【実行結果】

トレースプロットで事後分布サンプリングデータの様子を確認します。

### トレースプロットの表示
az.plot_trace(idata_m_nrt, compact=True)
plt.tight_layout();

【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。

(5) 推論データの保存
推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
pickleで保存します。

### idataの保存 pickle
file = r'idata_m_nrt_ch18.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata_m_nrt, f)

読み込みコードは次のとおりです。

### idataの読み込み pickle
file = r'idata_m_nrt_ch18.pkl'
with open(file, 'rb') as f:
    idata_m_nrt_load = pickle.load(f)

３．無意味条件・線形反応ツリーモデル

(1) モデルの定義

### モデルの定義
model_n_lrt = bmb.Model(
        formula='value ~ 0 + node + item +  (1 | person)',    # フォーミュラ式
        data=data_n_lrt,                                 # データ
        family='bernoulli',                              # 目的変数の誤差分布
)

(2) モデルの外観の確認

### モデルの表示
print(model_n_lrt)

【実行結果】

### モデルの可視化
model_n_lrt.build()
model_n_lrt.graph()

【実行結果】

(3) 事後分布からのサンプリング
処理時間はおよそ10秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　10秒
# テキスト: iter=100000, warmup=50000, chains=4, thin=10
idata_n_lrt = model_n_lrt.fit(
    draws=1000, tune=1000, chains=4, target_accept=0.85, nuts_sampler='numpyro', 
    idata_kwargs={'log_likelihood': True}, random_seed=1969)

【実行結果】省略

(4) サンプリングデータの確認
$${\hat{R}}$$、トレースプロットを確認します。
事後分布の収束確認は$${\hat{R} \leq 1.1}$$とします。
ここでは$${\hat{R} > 1.01}$$のパラメータが無いことを確認します。

### r_hat>1.1の確認
# 設定
idata_in = idata_n_lrt   # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
$${\hat{R} > 1.01}$$のパラメータは０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

### 推論データの要約統計情報の表示
az.summary(idata_n_lrt, hdi_prob=0.95, round_to=3)

【実行結果】

トレースプロットで事後分布サンプリングデータの様子を確認します。

### トレースプロットの表示
az.plot_trace(idata_n_lrt, compact=True)
plt.tight_layout();

【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。

(5) 推論データの保存
推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
pickleで保存します。

### idataの保存 pickle
file = r'idata_n_lrt_ch18.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata_n_lrt, f)

読み込みコードは次のとおりです。

### idataの読み込み pickle
file = r'idata_n_lrt_ch18.pkl'
with open(file, 'rb') as f:
    idata_n_lrt_load = pickle.load(f)

４．無意味条件・入れ子反応ツリーモデル

(1) モデルの定義

### モデルの定義
model_n_nrt = bmb.Model(
        formula='value ~ 0 + node + item +  (1 | person)',    # フォーミュラ式
        data=data_n_nrt,                                 # データ
        family='bernoulli',                              # 目的変数の誤差分布
)

(2) モデルの外観の確認

### モデルの表示
print(model_n_nrt)

【実行結果】

### モデルの可視化
model_n_nrt.build()
model_n_nrt.graph()

【実行結果】

(3) 事後分布からのサンプリング
処理時間はおよそ10秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　10秒
# テキスト: iter=100000, warmup=50000, chains=4, thin=10
idata_n_nrt = model_n_nrt.fit(
    draws=1000, tune=1000, chains=4, target_accept=0.8, nuts_sampler='numpyro', 
    idata_kwargs={'log_likelihood': True}, random_seed=1969)

【実行結果】省略

(4) サンプリングデータの確認
$${\hat{R}}$$、トレースプロットを確認します。
事後分布の収束確認は$${\hat{R} \leq 1.1}$$とします。
ここでは$${\hat{R} > 1.01}$$のパラメータが無いことを確認します。

### r_hat>1.1の確認
# 設定
idata_in = idata_n_nrt   # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
$${\hat{R} > 1.01}$$のパラメータは０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

### 推論データの要約統計情報の表示
az.summary(idata_n_nrt, hdi_prob=0.95, round_to=3)

【実行結果】

トレースプロットで事後分布サンプリングデータの様子を確認します。

### トレースプロットの表示
az.plot_trace(idata_n_nrt, compact=True)
plt.tight_layout();

【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。

(5) 推論データの保存
推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
pickleで保存します。

### idataの保存 pickle
file = r'idata_n_nrt_ch18.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata_n_nrt, f)

読み込みコードは次のとおりです。

### idataの読み込み pickle
file = r'idata_n_nrt_ch18.pkl'
with open(file, 'rb') as f:
    idata_n_nrt_load = pickle.load(f)

ちょっと一息つきます。

モデル１の結果・分析

１．項目の固定効果
テキストは項目$${item}$$の固定効果に注目します。
項目の固定効果（係数）の事後分布サンプルデータのフォレストプロットを描画します。
テキストの図18.4～18.7に相当します。

### フォレストプロットの描画 ※図18.4～7に相当

## フォレストプロット描画関数の定義
def forest_plot(idata, ax, model, cond):
    # 平均値で降順ソートするためのデータ準備
    item_means =idata.posterior['item'].mean(('chain', 'draw'))
    sorted_items = (idata.posterior['item_dim']
                    .sortby(item_means, ascending=False))
    # フォレストプロットの描画
    az.plot_forest(idata, var_names=['item'], coords={'item_dim': sorted_items},
                   combined=True, hdi_prob=0.95, ax=ax)
    # 修飾
    ax.axvline(0, color='red', ls='--')
    ax.set(xlim=[-8, 8],
           title=f'{model}反応ツリーモデルの項目の固定効果({cond}条件)')
    ax.grid(lw=0.3)

## 描画
# 描画領域の指定
fig, ax = plt.subplots(2, 2, figsize=(10, 10))
# ４つの推論データからフォレストプロットを描画
forest_plot(idata_m_lrt, ax[0, 0], '線形', '有意味')
forest_plot(idata_m_nrt, ax[1, 0], '入れ子', '有意味')
forest_plot(idata_n_lrt, ax[0, 1], '線形', '無意味')
forest_plot(idata_n_nrt, ax[1, 1], '入れ子', '無意味')
# 表示関連
plt.tight_layout()
plt.show()

【実行結果】

【分析～テキストにならって】
線形反応ツリーモデルも入れ子反応ツリーモデルも有意味条件の固定効果のほうがバラツキが大きく、無意味条件の効果のバラツキは小さいです。
テキストいわく「参加者は有意味語と無意味綴りでは異なる反応をしていることがわかる」です！

以下、個人的な感想です。

有意味語の場合、一般的な認知度合いやイメージの残りやすさなどが覚えやすさを決定づけるのかもです。
具体的には、よく使う単語は「よく経験しているから」覚えやすい、滅多に目にしないものごとに関する単語は「経験が少なくて」覚えにくい、と判断される可能性があるのかもと思った次第です。

一方、無意味のひらがな並びの「記号」的な言葉の場合、言葉や言葉が意味するものごとの使用経験は役に立たず、覚えやすさが一律に定まる可能性があるのかもです。

２．ノードの固定効果
続いてノード$${node}$$の固定効果の確認に進みます。
テキストでは固定効果の要約統計量と、推定結果をロジットリンク関数で０～１の範囲に収めた「主効果」の推定値を吟味します。

では要約統計量を算出しましょう。
テキストの表18.1、18.2に相当します。

### 線形反応ツリーモデルにおけるノードの固定効果の事後統計量 ※表18.1に相当

# 事後統計量算出関数の定義
def calc_stats(x, cond, node_no):
    ci = np.quantile(x, q=[0.025, 0.975])
    return cond, node_no, np.mean(x), np.std(x), ci[0], ci[1]

# 推論データからノードのサンプリングデータを抽出
m_lrt = idata_m_lrt.posterior.node.stack(sample=('chain', 'draw')).data
n_lrt = idata_n_lrt.posterior.node.stack(sample=('chain', 'draw')).data

# 事後統計量を算出してデータフレーム化
stats_df1 = pd.DataFrame(
            {calc_stats(m_lrt[0], '有意味', 'node1'),
             calc_stats(m_lrt[1], '有意味', 'node2'),
             calc_stats(m_lrt[2], '有意味', 'node3'),
             calc_stats(n_lrt[0], '無意味', 'node1'),
             calc_stats(n_lrt[1], '無意味', 'node2'),
             calc_stats(n_lrt[2], '無意味', 'node3'),},
            columns=['条件', 'ノード', 'EAP', 'posd.sd', '2.5%CI', '97.5%CI'])
stats_df1 = stats_df1.sort_values(['条件', 'ノード']).reset_index(drop=True)
# データフレームの表示
display(stats_df1.round(2))

【実行結果】
若干テキストの値と相違しますが、傾向的には似ているでしょう。

### 線形反応ツリーモデルにおけるノードの固定効果の事後統計量 ※表18.2に相当

# 推論データからノードのサンプリングデータを抽出
m_nrt = idata_m_nrt.posterior.node.stack(sample=('chain', 'draw')).data
n_nrt = idata_n_nrt.posterior.node.stack(sample=('chain', 'draw')).data

# 事後統計量を算出してデータフレーム化
stats_df2 = pd.DataFrame(
            {calc_stats(m_nrt[0], '有意味', 'node1'),
             calc_stats(m_nrt[1], '有意味', 'node2'),
             calc_stats(m_nrt[2], '有意味', 'node3'),
             calc_stats(n_nrt[0], '無意味', 'node1'),
             calc_stats(n_nrt[1], '無意味', 'node2'),
             calc_stats(n_nrt[2], '無意味', 'node3'),},
            columns=['条件', 'ノード', 'EAP', 'posd.sd', '2.5%CI', '97.5%CI'])
stats_df2 = stats_df2.sort_values(['条件', 'ノード']).reset_index(drop=True)
# データフレームの表示
display(stats_df2.round(2))

【実行結果】
こちらも若干テキストの値と異なります。

続いてノードの主効果です。
まずロジットリンク関数で0～1の範囲に収めたデータを作成します。
線形反応ツリーモデルと入れ子反応ツリーモデルの別に作成します。

### 事後分布サンプリングデータをロジットリンク関数で0～1の範囲に収める
# 線形反応ツリーモデル

## ロジットリンク関数（＝標準シグモイド関数）の定義
def sigmoid(x):
    return 1 / (1 + (np.exp(-x)))

## 初期値設定
# サンプリングデータの個数
num_sample = m_lrt.shape[1]

## データフレームに格納するデータのリストの準備
# 条件
cond_list = ['有意味']*3 + ['無意味']*3
# ノード
node_list = ['node1', 'node2', 'node3']*2
# サンプリングデータ
est_list = [m_lrt[0], m_lrt[1], m_lrt[2], n_lrt[0], n_lrt[1], n_lrt[2]]
# 格納するデータフレームの初期化
stats_df3 = pd.DataFrame()

## ロジットリンク関数でサンプリングデータを変換してデータフレームに格納
for (cond, node, est) in zip(cond_list, node_list, est_list):
    tmp = pd.DataFrame({'condition': [cond] * num_sample,
                        'node': [node] * num_sample,
                        'estimate': sigmoid(est)})
    stats_df3 = pd.concat([stats_df3, tmp], axis=0)

## データフレームの表示
display(stats_df3)

【実行結果】
線形反応ツリーモデルのノードの主効果です。

### 事後分布サンプリングデータをロジットリンク関数で0～1の範囲に収める
# 入れ子反応ツリーモデル

## データフレームに格納するデータのリストの準備
# サンプリングデータ
est_list = [m_nrt[0], m_nrt[1], m_nrt[2], n_nrt[0], n_nrt[1], n_nrt[2]]
# 格納するデータフレームの初期化
stats_df4 = pd.DataFrame()

## ロジットリンク関数でサンプリングデータを変換してデータフレームに格納
for (cond, node, est) in zip(cond_list, node_list, est_list):
    tmp = pd.DataFrame({'condition': [cond] * num_sample,
                        'node': [node] * num_sample,
                        'estimate': sigmoid(est)})
    stats_df4 = pd.concat([stats_df4, tmp], axis=0)

## データフレームの表示
display(stats_df4)

【実行結果】
入れ子反応ツリーモデルのノードの主効果です。

では描画に進みます。
描画しやすいデータに整えて描画します。
テキストの図18.8、18.9に相当します。

線形反応ツリーモデルから。

### ノードの主効果の描画のためのデータの作成  線形反応ツリーモデル

# 2.5%分位数関数・97.5%分位数関数の定義
def q025(x):
    return np.quantile(x, q=0.025)
def q975(x):
    return np.quantile(x, q=0.975)

# 条件・ノードごとの平均・95%CI区間を算出してデータフレーム化
stats_df3_plot = (stats_df3.groupby(['condition', 'node'])['estimate']
                  .agg(['mean', q025, q975]).reset_index())

# データフレームの表示
display(stats_df3_plot)

【実行結果】

### 線形反応ツリーモデルにおけるノード主効果の描画 ※図18.8に相当

# 有意味条件と無意味条件のデータを取り出し
plot_m = stats_df3_plot[:3]  # 有意味条件
plot_n = stats_df3_plot[3:]  # 無意味条件

plt.figure(figsize=(5, 4))
ax = plt.subplot()
# 有意味条件のエラーバープロットの描画
ax.errorbar(x=np.arange(1, 4) - 0.2,
            y=plot_m['mean'],
            yerr=np.vstack([abs(plot_m['q025'] - plot_m['mean']).values,
                            abs(plot_m['q975'] - plot_m['mean']).values]),
            capsize=8, fmt='d', label='有意味条件')
# 無意味条件のエラーバープロットの描画
ax.errorbar(x=np.arange(1, 4) + 0.2,
            y=plot_n['mean'],
            yerr=np.vstack([abs(plot_n['q025'] - plot_n['mean']).values,
                            abs(plot_n['q975'] - plot_n['mean']).values]),
            capsize=8, fmt='o', label='無意味条件')
# 修飾
ax.grid(lw=0.5)
ax.set(xticks=[0, 1, 2, 3, 4], xticklabels=['', 'node1', 'node2', 'node3', ''],
       xlabel='node', ylabel='estimate', ylim=(-0.1, 1.1))
plt.legend(title='条件');

【実行結果】
線形反応ツリーモデルのプロットです。
テキストのプロットとよく似ていると思います！

続いて入れ子反応ツリーモデルです。

### ノードの主効果の描画のためのデータの作成  入れ子反応ツリーモデル

# 条件・ノードごとの平均・95%CI区間を算出してデータフレーム化
stats_df4_plot = (stats_df4.groupby(['condition', 'node'])['estimate']
                  .agg(['mean', q025, q975]).reset_index())

# データフレームの表示
display(stats_df4_plot)

【実行結果】

### 入れ子反応ツリーモデルにおけるノード主効果の描画 ※図18.9に相当

# 有意味条件と無意味条件のデータを取り出し
plot_m = stats_df4_plot[:3]  # 有意味条件
plot_n = stats_df4_plot[3:]  # 無意味条件

plt.figure(figsize=(5, 4))
ax = plt.subplot()
# 有意味条件のエラーバープロットの描画
ax.errorbar(x=np.arange(1, 4) - 0.2,
            y=plot_m['mean'],
            yerr=np.vstack([abs(plot_m['q025'] - plot_m['mean']).values,
                            abs(plot_m['q975'] - plot_m['mean']).values]),
            capsize=8, fmt='d', label='有意味条件')
# 無意味条件のエラーバープロットの描画
ax.errorbar(x=np.arange(1, 4) + 0.2,
            y=plot_n['mean'],
            yerr=np.vstack([abs(plot_n['q025'] - plot_n['mean']).values,
                            abs(plot_n['q975'] - plot_n['mean']).values]),
            capsize=8, fmt='o', label='無意味条件')
# 修飾
ax.grid(lw=0.5)
ax.set(xticks=[0, 1, 2, 3, 4], xticklabels=['', 'node1', 'node2', 'node3', ''],
       xlabel='node', ylabel='estimate', ylim=(-0.1, 1.1))
plt.legend(title='条件');

【実行結果】
入れ子反応ツリーモデルのプロットです。
こちらもテキストのプロットに近いです！

【分析～テキストにならって】
テキストは次の２点を指摘しています。

１．線形反応ツリーモデルでは条件にかかわらずノード$${Y^*_1}$$からノード$${Y^*_3}$$にかけて０の分岐が選択されがち
→ 参加者が段階的に項目の覚えやすさを判断していると考えられる

２．どちらのツリーモデルも無意味綴りの方が１（難しい方向）への分岐が選ばれる可能性が高い

２番めの指摘は納得しました！
意味あり言葉の方が覚え易くて、意味無し言葉の方が覚えにくいのは、感覚的に理解しやすいです。
１番めの指摘は、データから指摘の内容を読み取ることができませんでした（汗）

そしてテキストは、いったん「入れ子反応ツリーモデル」を選択して、次の分析に進みます。

モデル２・３の構築

テキストは３つのノードでなされる判断が「似ている」のか「質的に異なる」のかを検証するために、２つのベイズモデルを追加します。
用いるデータは「入れ子反応ツリーモデル」です。

【モデル２】
モデル２は各ノードで異なる判断が行われると仮定するモデルであり、ノードごとのランダム傾きをベイズモデルに追加しています。

【モデル３】
モデル３は「最初のノード$${Y^*_1}$$」と「残りのノード$${Y^*_2,\ Y^*_3}$$」とで異なる判断が行われると仮定するモデルであり、ノード$${Y^*_1}$$とノード$${Y^*_2,\ Y^*_3}$$で異なるランダム傾きをベイズモデルに追加しています。

【分析の視点】
モデルの評価指標 WAIC・looを用いて、モデル１・２・３の中から「良いモデル」を選ぶことで「ノードでなされる判断」（似てる／異なる）の検証を行います。

ではモデル２から構築を開始しましょう。

モデル２の数式表現
目指したいBambiのモデルの「なんちゃって数式」表記です。
formula 形式で書きます。

$$
value \sim 0 + node + item + (0 + node\ |\ person)
$$

ノード$${node}$$のランダム傾き「$${(0 + node\ |\ person)}$$」を追加したモデルです。
目的変数$${value}$$はベルヌーイ分布に従います。

モデル２・有意味条件

(1) モデルの定義

### モデルの定義
model_m_nrt2 = bmb.Model(
        formula='value ~ 0 + node + item + (0 + node | person)', # フォーミュラ式
        data=data_m_nrt,                                 # データ
        family='bernoulli',                              # 目的変数の誤差分布
)

(2) モデルの外観の確認

### モデルの表示
model_m_nrt2

【実行結果】
各パラメータが従う事前分布は、Bambiがデータを解析して自動設定してくれます。

### モデルの可視化
model_m_nrt2.build()
model_m_nrt2.graph()

【実行結果】

(3) 事後分布からのサンプリング
処理時間はおよそ10秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　10秒
# テキスト: iter=100000, warmup=50000, chains=4, thin=10
idata_m_nrt2 = model_m_nrt2.fit(
    draws=1000, tune=1000, chains=4, target_accept=0.9, nuts_sampler='numpyro', 
    idata_kwargs={'log_likelihood': True}, random_seed=1969)

【実行結果】省略

(4) サンプリングデータの確認
$${\hat{R}}$$、トレースプロットを確認します。
事後分布の収束確認は$${\hat{R} \leq 1.1}$$とします。
ここでは$${\hat{R} > 1.01}$$のパラメータが無いことを確認します。

### r_hat>1.1の確認
# 設定
idata_in = idata_m_nrt2  # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
$${\hat{R} > 1.01}$$のパラメータは０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

### 推論データの要約統計情報の表示
az.summary(idata_m_nrt2, hdi_prob=0.95, round_to=3)

【実行結果】

トレースプロットで事後分布サンプリングデータの様子を確認します。

### トレースプロットの表示
az.plot_trace(idata_m_nrt2, compact=True)
plt.tight_layout();

【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。

(5) 推論データの保存
推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
pickleで保存します。

### idataの保存 pickle
file = r'idata_m_nrt2_ch18.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata_m_nrt2, f)

読み込みコードは次のとおりです。

### idataの読み込み pickle
file = r'idata_m_nrt2_ch18.pkl'
with open(file, 'rb') as f:
    idata_m_nrt2_load = pickle.load(f)

モデル２・無意味条件

(1) モデルの定義

### モデルの定義
model_n_nrt2 = bmb.Model(
        formula='value ~ 0 + node + item + (0 + node | person)', # フォーミュラ式
        data=data_n_nrt,                                 # データ
        family='bernoulli',                              # 目的変数の誤差分布
)

(2) モデルの外観の確認

### モデルの表示
model_n_nrt2

【実行結果】
各パラメータが従う事前分布は、Bambiがデータを解析して自動設定してくれます。

### モデルの可視化
model_n_nrt2.build()
model_n_nrt2.graph()

【実行結果】

(3) 事後分布からのサンプリング
処理時間はおよそ10秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　10秒
# テキスト: iter=100000, warmup=50000, chains=4, thin=10
idata_n_nrt2 = model_n_nrt2.fit(
    draws=1000, tune=1000, chains=4, target_accept=0.8, nuts_sampler='numpyro', 
    idata_kwargs={'log_likelihood': True}, random_seed=1969)

【実行結果】省略

(4) サンプリングデータの確認
$${\hat{R}}$$、トレースプロットを確認します。
事後分布の収束確認は$${\hat{R} \leq 1.1}$$とします。
ここでは$${\hat{R} > 1.01}$$のパラメータが無いことを確認します。

### r_hat>1.1の確認
# 設定
idata_in = idata_n_nrt2  # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
$${\hat{R} > 1.01}$$のパラメータは０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

### 推論データの要約統計情報の表示
az.summary(idata_n_nrt2, hdi_prob=0.95, round_to=3)

【実行結果】

トレースプロットで事後分布サンプリングデータの様子を確認します。

### トレースプロットの表示
az.plot_trace(idata_n_nrt2, compact=True)
plt.tight_layout();

【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。

(5) 推論データの保存
推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
pickleで保存します。

### idataの保存 pickle
file = r'idata_n_nrt2_ch18.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata_n_nrt2, f)

読み込みコードは次のとおりです。

### idataの読み込み pickle
file = r'idata_n_nrt2_ch18.pkl'
with open(file, 'rb') as f:
    idata_n_nrt2_load = pickle.load(f)

モデル３の構築に移ります。

モデル３の数式表現
目指したいBambiのモデルの「なんちゃって数式」表記です。
formula 形式で書きます。

$$
value \sim 0 + node + item + (0 + node2\ |\ person)
$$

$${node2}$$はノード$${Y^*_2, Y^*_3}$$のときに True、ノード$${Y^*_1}$$のときに False が設定されたカテゴリ変数です。
$${node2}$$のランダム傾き「$${(0 + node2\ |\ person)}$$」を追加したモデルです。
目的変数$${value}$$はベルヌーイ分布に従います。

$${node2}$$をデータに設定します。

### データの前処理
# node2,3フラグを追加
data_m_nrt['node2'] = (data_m_nrt['node'] != 'node1').astype('category')
data_n_nrt['node2'] = (data_n_nrt['node'] != 'node1').astype('category')

モデル３・有意味条件

(1) モデルの定義

### モデルの定義
model_m_nrt3 = bmb.Model(
        formula='value ~ 0 + node + item + (0 + node2 | person)', # フォーミュラ式
        data=data_m_nrt,                                 # データ
        family='bernoulli',                              # 目的変数の誤差分布
)

(2) モデルの外観の確認

### モデルの表示
model_m_nrt3

【実行結果】
各パラメータが従う事前分布は、Bambiがデータを解析して自動設定してくれます。

### モデルの可視化
model_m_nrt3.build()
model_m_nrt3.graph()

【実行結果】

(3) 事後分布からのサンプリング
処理時間はおよそ10秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　10秒
# テキスト: iter=100000, warmup=50000, chains=4, thin=10
idata_m_nrt3 = model_m_nrt3.fit(
    draws=1000, tune=1000, chains=4, target_accept=0.85, nuts_sampler='numpyro',
    idata_kwargs={'log_likelihood': True}, random_seed=1969)

【実行結果】省略

(4) サンプリングデータの確認
$${\hat{R}}$$、トレースプロットを確認します。
事後分布の収束確認は$${\hat{R} \leq 1.1}$$とします。
ここでは$${\hat{R} > 1.01}$$のパラメータが無いことを確認します。

### r_hat>1.1の確認
# 設定
idata_in = idata_m_nrt3  # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
$${\hat{R} > 1.01}$$のパラメータは０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

### 推論データの要約統計情報の表示
az.summary(idata_m_nrt3, hdi_prob=0.95, round_to=3)

【実行結果】

トレースプロットで事後分布サンプリングデータの様子を確認します。

### トレースプロットの表示
az.plot_trace(idata_m_nrt3, compact=True)
plt.tight_layout();

【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。

(5) 推論データの保存
推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
pickleで保存します。

### idataの保存 pickle
file = r'idata_m_nrt3_ch18.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata_m_nrt3, f)

読み込みコードは次のとおりです。

### idataの読み込み pickle
file = r'idata_m_nrt3_ch18.pkl'
with open(file, 'rb') as f:
    idata_m_nrt3_load = pickle.load(f)

モデル３・無意味条件

(1) モデルの定義

### モデルの定義
model_n_nrt3 = bmb.Model(
        formula='value ~ 0 + node + item + (0 + node2 | person)', # フォーミュラ式
        data=data_n_nrt,                                 # データ
        family='bernoulli',                              # 目的変数の誤差分布
)

(2) モデルの外観の確認

### モデルの表示
model_n_nrt3

【実行結果】
各パラメータが従う事前分布は、Bambiがデータを解析して自動設定してくれます。

### モデルの可視化
model_n_nrt3.build()
model_n_nrt3.graph()

【実行結果】

(3) 事後分布からのサンプリング
処理時間はおよそ10秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　10秒
# テキスト: iter=100000, warmup=50000, chains=4, thin=10
idata_n_nrt3 = model_n_nrt3.fit(
    draws=1000, tune=1000, chains=4, target_accept=0.85, nuts_sampler='numpyro',
    idata_kwargs={'log_likelihood': True}, random_seed=1969)

【実行結果】省略

(4) サンプリングデータの確認
$${\hat{R}}$$、トレースプロットを確認します。
事後分布の収束確認は$${\hat{R} \leq 1.1}$$とします。
ここでは$${\hat{R} > 1.01}$$のパラメータが無いことを確認します。

### r_hat>1.1の確認
# 設定
idata_in = idata_n_nrt3  # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
$${\hat{R} > 1.01}$$のパラメータは０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

### 推論データの要約統計情報の表示
az.summary(idata_n_nrt3, hdi_prob=0.95, round_to=3)

【実行結果】

トレースプロットで事後分布サンプリングデータの様子を確認します。

### トレースプロットの表示
az.plot_trace(idata_n_nrt3, compact=True)
plt.tight_layout();

【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。

(5) 推論データの保存
推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
pickleで保存します。

### idataの保存 pickle
file = r'idata_n_nrt3_ch18.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata_n_nrt3, f)

読み込みコードは次のとおりです。

### idataの読み込み pickle
file = r'idata_n_nrt3_ch18.pkl'
with open(file, 'rb') as f:
    idata_n_nrt3_load = pickle.load(f)

一息つきます。ふぅー。

モデル１・２・３の比較

WAICとlooを用いて３つのモデルの中から良いモデルを選択します。
テキストの表18.3に相当する表を作成します。

### 入れ子反応ツリーモデルにおけるモデル１～３の適合度指標 ※表18.3に相当

## 初期値設定
# WAIC,loo計算値を格納する一時配列
stats5 = np.zeros((6, 4))
# forループで回す推論データのリスト
idata_list = [idata_m_nrt, idata_m_nrt2, idata_m_nrt3,
              idata_n_nrt, idata_n_nrt2, idata_n_nrt3]

## 推論データごとにWAIC,looの算出を繰り返し処理
for i, idata in enumerate(idata_list):
    stats5[i, 0] = az.waic(idata, scale='deviance')[0]  # WAIC 推定値
    stats5[i, 1] = az.waic(idata, scale='deviance')[1]  # WAIC SE
    stats5[i, 2] = az.loo(idata, scale='deviance')[0]   # loo 推定値
    stats5[i, 3] = az.loo(idata, scale='deviance')[1]   # loo 推定値

## データフレーム化
# 行名のマルチインデックスの設定
multi_idx_names1 = ['有意味']*3 + ['無意味']*3
multi_idx_names2 = ['model1', 'model2', 'model3']*2
multi_idx = pd.MultiIndex.from_arrays([multi_idx_names1, multi_idx_names2])
# 列名のマルチインデックスの設定
multi_col_names1 = ['WAIC']*2 + ['loo']*2
multi_col_names2 = ['推定値', 'SE']*2
multi_col = pd.MultiIndex.from_arrays([multi_col_names1, multi_col_names2])
# データフレーム化
stats_df5 = pd.DataFrame(stats5, index=multi_idx, columns=multi_col)
# データフレームの表示
display(stats_df5
        .style.set_properties(**{'background-color': 'yellow'}, 
                            subset=pd.IndexSlice[('有意味', 'model1'), :])
        .set_properties(**{'background-color': 'yellow'}, 
                            subset=pd.IndexSlice[('無意味', 'model2'), :])
        .format('{:.1f}'))

【実行結果】
WAIC・looともに、推定値が小さいほうが良いモデルと判断されます。

【分析～テキストにならって】
■ 有意味条件の場合
ノードごとの変量効果を仮定しない model1 がもっとも良いモデルです。
テキストは「有意味条件の覚えやすさの判断はノードごとに異なるものではなく、一元的である可能性がある」としています。
■ 無意味条件の場合
ノードごとのランダム傾きを加味した model2 がもっとも良いモデルです。
テキストは「無意味条件の覚えやすさの判断は一元的ではなく、各ノードで質的に異なる判断が行われていることが推察される」としています。

ちなみに、線形反応ツリーモデル・有意味条件は、以下のWAICとlooの値が示すとおり、model1 ~ 3よりも「良いモデル」かもです・・・

$$
\begin{array}{}
\text{WAIC} & \text{loo} \\
\hline385.88 & 387.33
\end{array}
$$

以上で第18章は終了です。

おわりに

大量のモデル群

今回の章は８つのベイズモデルを構築・比較しました。
たくさんのモデルが出現したことで何かと混乱気味です。
しかも「条件」 ×「項目反応ツリー」 ×「ランダム効果」といったモデルの要素が多くて、混乱度がマシマシです！
複雑な注文呪文を思い出しました。

もやしがたくさん乗ったラーメンのイラスト：「いらすとや」さんより

モデルの覚えやすさにも心理学的な判断メカニズムが存在するのでしょうか？
先生、ぜひ論文を。

おわり

シリーズの記事

次の記事

前の記事

目次


ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。

いいなと思ったら応援しよう！

この記事が参加している募集

#新生活をたのしく

47,054件