Bayesian Modeling of Chapter 18 "Mechanism of Judging 'Ease of Memorization'" with PyMC Ver. 5
This article is a statistical documentary depicting the process of "experimentally" implementing the Bayesian model from Chapter 18 "Mechanism of Judging 'Ease of Memorization'" of the text "Fun Bayesian Modeling 2" using PyMC Ver. 5.
This time, it is another Bayesian analysis with a strong academic psychological flavor.
This chapter analyzes the logic (pathway) used to judge the ease of memorizing words and terms when classifying them into four levels (easy to difficult) using an item response tree model.
The text uses brms.
It is said to be a library that allows you to write Stan models using the formula syntax, which is the regression formula notation in R.
So, for this PyMC conversion, I will be using "Bambi," which can be written with formulas!
The inference results from Bambi were fairly close to the results in the text.
Well then, let's enjoy the world of PyMC Bayesian modeling!
Please refer to the article at this link for the introduction of the text, citation notation, series preface, and version information for PyMC, etc.
The data used in the text, along with sample scripts for R, Stan, etc., can be downloaded from the publisher's website.
Summary
Overview of the text
Author: Dr. Takashi Yamane
Model Difficulty: ★★★.. (Average)
Self-evaluation
Rating
$$
\begin{array}{c:c:c}
Implementation Accuracy & ★★★★. & Almost satisfied \\
Result Reproducibility & ★★★★. & One step away \\
Fun & ★★★★★ & Fun! \\
\end{array}
$$
Evaluation Points
It took time to understand the item response tree model itself, so I set the model difficulty a bit higher.
The regression formula incorporating fixed effects and random effects written in the brms formula is simple, but the "tree structure + mapping matrix" and "the contents of the data/how variable values are held" did not connect in my head.
I tried to grasp the relationship between the tree structure and the mapping matrix by running R scripts in an R environment and checking the changes in variable values.
Ingenuity, Joy, and Reflection
I did not import pymc this time!
I practiced PyMC modeling & inference using only Bambi and arviz.
Model Overview
Overview of the text's research and experiments
■ Overview of the experiment
This is an experiment where participants rate the ease of memorizing presented "three-character hiragana" strings on a 4-point scale from "1: Easy" to "4: Difficult".
-
There are two types of strings used in the experiment:
Meaningful Condition (M)
Uses "meaningful words" that have meaning, such as "suika" (watermelon) or "sakura" (cherry blossom).Meaningless Condition (N)
Uses "meaningless strings" that have no meaning, such as "toparu" or "shiregi".
The number of participants is 10 for the meaningful condition and 10 for the meaningless condition.
Participants answer 20 items in their assigned condition.
The number of data points is 200 for each condition, calculated as 10 participants × 20 response items.
If you would like to know the details of the experiment, please readthe professor's paper here (PDF download)!
■ Overview of the analysis
We are interested in "how participants react when choosing from the 4-point scale" when they "judge which of the 4 levels of ease of memorization the presented string corresponds to."
The text examines two item response tree models regarding how these reactions occur.
① Linear Response Tree Model (LRT)
This is a model where one examines and selects the applicable level in order, starting from the first level of "Easy."
Elements represented by ellipses are called "nodes."
② Nested Response Tree Model (NRT)
This is a two-step model where one first examines and judges whether it is generally easy or difficult, and then examines and judges the degree of ease (easy or somewhat easy) or the degree of difficulty (somewhat difficult or difficult).
The text assumes the judgment steps of thesetwo tree models to explore the judgment mechanism of ease of memorization.
Also, in modeling, tree models are represented by how the data is structured.
Let's summarize for now.
Points of interest
We model the judgment mechanism of ease of memorization—that is, ease of learning—and estimate the parameters.
Materials used for "ease of memorization"
We handle two conditions: the "meaningful condition," which judges the ease of memorizing words with meaning, and the "meaningless condition," which judges the ease of memorizing words without meaning.
Assumed mechanism
For the judgment mechanism, we handle two item response tree models, the "Linear Response Tree Model" and the "Nested Response Tree Model," through the "way data is held."
■ Mapping Matrix
Things get a bit complicated from here.
We convert the tree diagrams of the item response tree models into a matrix format.
It is a matrix where the selected answers are rows, the elliptical nodes are columns, and the choices attached to the tree diagram lines are the values.
We will generalize and re-present the two tree models mentioned earlier and convert them into mapping matrices.
These correspond to Figures 18.1 and 18.2 in the text.
① Linear Response Tree Model and Mapping Matrix
Let's imagine that choices like "Easy?" are inside the elliptical nodes.
The "Y=..." part, which is the final destination of the arrows, is the response value chosen from "1: Easy" to "4: Difficult."
The values 0 and 1 attached to the arrows correspond to the values judged at each node.
At each node, 0 seems to represent a judgment toward the easier side, and 1 toward the more difficult side.
Here is an example of how to read it.
Focus on the path for the tree response value "Y=2" (somewhat easy).
The judgment at the first node "Y1*" is 1 (proceed toward the difficult side), and the subsequent judgment at node "Y2*" is 0 (proceed toward the easy side).
Now look at the mapping matrix.
For "Y=2," the value of "Y1*" is 1, and the value of "Y2*" is 0.
The mapping matrix holds the path information of the tree.
It's like looking at the values in order from the left side of the columns.
Note that since it does not go to Y3*, the value is "-".
The values 0 and 1 in thismapping matrix become the objective variables for Bayesian modeling.
② Nested Response Tree Model and Mapping Matrix
I will illustrate the overview of the analysis.
Modeling the text
Since it is difficult to explain by perfectly matching the model formulas in the text with the formula descriptions in brms, I will first write down the model formulas as they appear in the text.
■ Objective variables and parameters of interest
The objective variable is the value $${Y^*_{pin}}$$ of the path chosen at node $${n}$$ when experimental participant $${p}$$ responds to string item number $${i}$$.
The value of the mapping matrix is set to $${Y^*_{pin}}$$.
For example, assuming a linear response tree model, if an experimental participant with ID=1 looks at the string for item number 3 and chooses path 0 at node "Y3*" (number 2), then $${Y^*_{1,3,2}=0}$$.
Parameters of interest include the regression coefficient $${\beta_i}$$ indicating the fixed effect of the item $${Item}$$, and the regression coefficient $${\beta_n}$$ indicating the fixed effect of the node $${Node}$$, among others.
■ Model formulas
The likelihood function in the first line uses a Bernoulli distribution.
The probability parameter of the Bernoulli distribution in the second line uses a logistic regression model with nodes, items, and experimental participants as explanatory variables.
$$
\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*}
$$
The text adds that "typically, $${Item}$$ is entered into the model as a fixed effect, and $${Person}$$ as a random effect."
■ Three Bayesian models
The text constructs three Bayesian models for the purpose of confirming whether judgments at each node are similar or different.
$$
\begin{array}{l:l}
Model 1 & Random effect specified for intercept only \\\\&Assuming common judgment \\\\&is made at all nodes \\\\
\hdashline
Model 2 & Random slopes for each node included \\\\ & Assuming different judgments \\\\&are made at each node \\\\
\hdashline
Model 3 & Different random slopes included for node Y^*_1, node Y^*_2, and node Y^*_3 \\ & Assuming different judgments \\\\ & are made at node Y^*_1, node Y^*_2, and node Y^*_3 \\\\
\end{array}
$$
Model 1 handles 2 tree models × 2 conditions = 4 Bayesian models.
Models 2 and 3 handle nested response tree models × 2 conditions × 2 Bayesian models = 4 Bayesian models.
A total of 8 Bayesian models are constructed! That's a lot!
■ Analysis and analysis results
I believe the text's descriptions of the analysis methods and numerical results are accurate, so I recommend reading the text.
Please see the "PyMC Implementation" chapter for analysis using inferred values from my own PyMC models.
PyMC Implementation
Let's enjoy PyMC & Python !
Preparation and data verification
1. Import
### インポート
# 数値・確率計算
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')2. Creating explanatory charts
I will create the charts used in the "Model Overview" chapter.
I will use the directed graph Digraph from the graphviz library.
① Image diagram of the linear response tree model
### 線形反応ツリーモデルのイメージ 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[Execution Results]
(2) Image of the nested reaction tree model
### 入れ子反応ツリーモデルのイメージ 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[Execution Results]
(3) Linear reaction tree model and mapping matrix in Figure 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[Execution Results]
### 線形反応ツリーモデルのマッピング行列のイメージ ※図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)[Execution Results]
(3) Nested reaction tree model and mapping matrix in Figure 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[Execution Results]
### 線形反応ツリーモデルのマッピング行列のイメージ ※図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)[Execution Results]
3. Data loading and preprocessing
Load the two csv files into pandas dataframes.
- Response data for the meaningful condition "EOLdatM.csv"
- Response data for the meaningless condition "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)[Execution Results]
Rows represent experimental participant IDs, columns represent response items, and values are the response values selected from "1: Easy to 4: Difficult".
Proceed to data preprocessing using the mapping matrix.
First, define the mapping matrix.
### マッピング行列の作成
# 線形反応ツリーモデルのマッピング行列
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)[Execution Results]
Next is the definition of the data conversion function.
Since I could not find a Python library or function equivalent to R's Dendrify function, I created my own.
The arguments are the data to be converted (wide-format data): df, and the mapping matrix: map_mtx.
The return value is long-format data with the mapping matrix values set.
### データの前処理：マッピング行列を用いたデータ変換処理を関数化
# 引数 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_dfProceed to data conversion.
Convert and create four datasets: meaningful condition/meaningless condition × linear reaction tree model/nested reaction tree model.
[Data assuming a tree model]
Create data using a mapping matrix that follows the structure of the two tree models.
Therefore, data is created assuming a specific tree model.
This means the data already incorporates one of the tree models.
(1) Meaningful condition × Linear reaction tree model LRT
### データの前処理：マッピング行列を用いたデータ変換 有意味条件・LRT
data_m_lrt = make_data(data_m_orgn, lrt_map)
display(data_m_lrt)[Execution Results]
The data shows the value (0 or 1) selected by the experimental participant 'person' at node number 'node' when responding to the string 'item'.
How to read the table:
- value: The value of the mapping matrix, which is the objective variable
- item: Response item (identifier for 20 three-character strings)
- person: Experimental participant ID
- node: Node number of the reaction tree model
item, person, and node are categorical variables.
(2) Meaningless condition × Linear reaction tree model LRT
### データの前処理：マッピング行列を用いたデータ変換 無意味条件・LRT
data_n_lrt = make_data(data_n_orgn, lrt_map)
display(data_n_lrt)[Execution Results]
(3) Meaningful Condition × Nested Reaction Tree Model NRT
### データの前処理：マッピング行列を用いたデータ変換 有意味条件・NRT
data_m_nrt = make_data(data_m_orgn, nrt_map)
display(data_m_nrt)[Execution Result]
(4) Meaningless Condition × Nested Reaction Tree Model NRT
### データの前処理：マッピング行列を用いたデータ変換 無意味条件・NRT
data_n_nrt = make_data(data_n_orgn, nrt_map)
display(data_n_nrt)[Execution Result]
The repetition process for the 4 patterns has continued. Thank you for your hard work.
This repetition process involving multiple patterns will continue, so thank you for your cooperation!
Construction of Model 1
We will perform four modeling tasks for "Model 1: A model with random effects specified only for the intercept," which is used in Section 18.3 "Data and Model Fitting" and Section 18.4 "Results and Discussion" of the text.
Except for the fact that the data to be analyzed is different, the content of the models is identical.
[Perspective of Analysis]
Analyzing using the posterior distribution sample data of the parameters involves "fixed effects of items" and "fixed effects of nodes."
The fixed effects of items analyze the difference between the meaningful condition (meaningful words) and the meaningless condition (meaningless spellings) using the coefficients for each individual string item.
The fixed effects of nodes analyze the difference between the linear reaction tree model and the nested reaction tree model, as well as the difference between the meaningful and meaningless conditions, using the summary statistics of the coefficients for each node.
Mathematical Expression of the Model
This is the "pseudo-mathematical" notation for the Bambi model we aim to create.
We will write it in formula format.
$$
value \sim 0 + node + item + (1\ |\ person)
$$
This is a random intercept model for experimental participant $${person}$$.
The objective variable $${value}$$ follows a Bernoulli distribution.
Now, let's begin the four Bayesian modeling tasks!
1. Meaningful Condition / Linear Reaction Tree Model
(1) Definition of the Model
### モデルの定義
model_m_lrt = bmb.Model(
formula='value ~ 0 + node + item + (1 | person)', # フォーミュラ式
data=data_m_lrt, # データ
family='bernoulli', # 目的変数の誤差分布
# link='logit', # （省略可）リンク関数
# categorical=['item', 'person', 'node'], # （省略可）カテゴリ変数
)[Model Annotation]
Specify the Bernoulli distribution, which is the error distribution of the objective variable, for the family.
The link function 'link' uses the default value 'logit' for the Bernoulli distribution.
In the subsequent three models, only the specification of the data used will change.
(2) Confirmation of Model Appearance
### モデルの表示
print(model_m_lrt)[Execution Result]
Bambi automatically analyzes the data and sets the prior distributions followed by the fixed effects (regression coefficients) of nodes and items, as well as the prior distributions followed by the random effects of experimental participants.
### モデルの可視化
model_m_lrt.build()
model_m_lrt.graph()[Execution Result]
To display the Bambi model in graph format, you must first perform model construction ("build") or model training ("fit").
(3) Sampling from the Posterior Distribution
In Bambi, you can execute MCMC by performing a 'fit' on the model.
Since we calculate the model evaluation indices WAIC and loo during analysis, we specify 'idata_kwargs={'log_likelihood': True}' so that the log-likelihood is calculated during MCMC execution.
The number of random number generations is smaller than in the text.
By setting 'nuts_sampler='numpyro'', you can use numpyro as the NUTS sampler.
The processing time was approximately 1 minute.
### 事後分布からのサンプリング ※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)[Execution Result] Omitted
(4) Checking Sampling Data
Check $$\hat{R}$$ and the trace plots.
We will use $$\hat{R} \leq 1.1$$ to confirm convergence of the posterior distribution.
Here, we confirm that there are no parameters with $$\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())[Execution Result]
There were 0 parameters with $$\hat{R} > 1.01$$.
We were able to confirm that all parameters satisfy $$\hat{R} \leq 1.1$$.
We will roughly check the summary statistics and trace plots of the posterior distribution sampling data.
### 推論データの要約統計情報の表示
az.summary(idata_m_lrt, hdi_prob=0.95, round_to=3)[Execution Result]
We will check the state of the posterior distribution sampling data using trace plots.
### トレースプロットの表示
az.plot_trace(idata_m_lrt, compact=True)
plt.tight_layout();[Execution Result]
From the graph on the left, we can see that the four Markov chains have almost the same distribution.
In the graph on the right, the lines are drawn evenly.
It appears to have converged.
(5) Saving Inference Data
Let's save the inference data to a file in case we (might) reuse it.
We will save it using pickle.
### idataの保存 pickle
file = r'idata_m_lrt_ch18.pkl'
with open(file, 'wb') as f:
pickle.dump(idata_m_lrt, f)The code for loading is as follows.
### idataの読み込み pickle
file = r'idata_m_lrt_ch18.pkl'
with open(file, 'rb') as f:
idata_m_lrt_load = pickle.load(f)(6) Preparation
Let's practice drawing forest plots of parameter estimates for items used in the analysis, and calculating WAIC and loo.
① Drawing a Forest Plot
We will use the plot_forest function from arviz.
We will sort the data in ascending order by the mean value in advance.
### フォレストプロットの描画 ※図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);[Execution Result]
② Calculating WAIC
We will use the waic function from arviz.
### WAICの算出
az.waic(idata_m_lrt, scale='deviance')[Execution Result]
③ Calculating loo
We will use the loo function from arviz.
### LOOの算出
az.loo(idata_m_lrt, scale='deviance')[Execution Result]
2. Meaningful Condition / Nested Reaction Tree Model
(1) Defining the Model
### モデルの定義
model_m_nrt = bmb.Model(
formula='value ~ 0 + node + item + (1 | person)', # フォーミュラ式
data=data_m_nrt, # データ
family='bernoulli', # 目的変数の誤差分布
)(2) Checking the Model Structure
### モデルの表示
print(model_m_nrt)[Execution Result]
### モデルの可視化
model_m_nrt.build()
model_m_nrt.graph()[Execution Result]
(3) Sampling from the posterior distribution
The processing time was approximately 15 seconds.
### 事後分布からのサンプリング ※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)[Execution result] Omitted
(4) Confirmation of sampling data
Check $$\hat{R}$$ and the trace plot.
The convergence of the posterior distribution is confirmed with $$\hat{R} \leq 1.1$$.
Here, we confirm that there are no parameters with $$\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())[Execution result]
There were 0 parameters with $$\hat{R} > 1.01$$.
We were able to confirm that all parameters satisfy $$\hat{R} \leq 1.1$$.
We will roughly check the summary statistics and trace plot of the posterior distribution sampling data.
### 推論データの要約統計情報の表示
az.summary(idata_m_nrt, hdi_prob=0.95, round_to=3)[Execution result]
We will check the state of the posterior distribution sampling data using a trace plot.
### トレースプロットの表示
az.plot_trace(idata_m_nrt, compact=True)
plt.tight_layout();[Execution result]
From the graph on the left, it can be seen that the four Markov chains have almost the same distribution.
In the graph on the right, the lines are drawn evenly.
It appears to have converged.
(5) Saving inference data
Let's save the inference data to a file in case we (might) reuse it later.
We will save it using pickle.
### idataの保存 pickle
file = r'idata_m_nrt_ch18.pkl'
with open(file, 'wb') as f:
pickle.dump(idata_m_nrt, f)The code for loading is as follows.
### idataの読み込み pickle
file = r'idata_m_nrt_ch18.pkl'
with open(file, 'rb') as f:
idata_m_nrt_load = pickle.load(f)
3. Meaningless condition/linear reaction tree model
(1) Model definition
### モデルの定義
model_n_lrt = bmb.Model(
formula='value ~ 0 + node + item + (1 | person)', # フォーミュラ式
data=data_n_lrt, # データ
family='bernoulli', # 目的変数の誤差分布
)(2) Confirmation of model appearance
### モデルの表示
print(model_n_lrt)[Execution result]
### モデルの可視化
model_n_lrt.build()
model_n_lrt.graph()[Execution result]
(3) Sampling from the posterior distribution
The processing time was approximately 10 seconds.
### 事後分布からのサンプリング ※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)[Execution result] Omitted
(4) Confirmation of sampling data
Check $$\hat{R}$$ and the trace plot.
The convergence of the posterior distribution is confirmed with $$\hat{R} \leq 1.1$$.
Here, we confirm that there are no parameters with $$\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())[Execution result]
There were 0 parameters with $$\hat{R} > 1.01$$.
We were able to confirm that all parameters satisfy $$\hat{R} \leq 1.1$$.
We will roughly check the summary statistics and trace plot of the posterior distribution sampling data.
### 推論データの要約統計情報の表示
az.summary(idata_n_lrt, hdi_prob=0.95, round_to=3)[Execution Results]
Check the state of the posterior distribution sampling data using trace plots.
### トレースプロットの表示
az.plot_trace(idata_n_lrt, compact=True)
plt.tight_layout();[Execution Results]
From the graph on the left, it can be seen that the four Markov chains have almost the same distribution.
In the graph on the right, the lines are drawn evenly.
It is considered to have converged.
(5) Saving inference data
Let's save the inference data to a file in case we (might) want to reuse it later.
We will save it using pickle.
### idataの保存 pickle
file = r'idata_n_lrt_ch18.pkl'
with open(file, 'wb') as f:
pickle.dump(idata_n_lrt, f)The code for loading is as follows.
### idataの読み込み pickle
file = r'idata_n_lrt_ch18.pkl'
with open(file, 'rb') as f:
idata_n_lrt_load = pickle.load(f)
4. Meaningless condition/nested reaction tree model
(1) Model definition
### モデルの定義
model_n_nrt = bmb.Model(
formula='value ~ 0 + node + item + (1 | person)', # フォーミュラ式
data=data_n_nrt, # データ
family='bernoulli', # 目的変数の誤差分布
)(2) Checking the model structure
### モデルの表示
print(model_n_nrt)[Execution Results]
### モデルの可視化
model_n_nrt.build()
model_n_nrt.graph()[Execution Results]
(3) Sampling from the posterior distribution
The processing time was approximately 10 seconds.
### 事後分布からのサンプリング ※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)[Execution Results] Omitted
(4) Checking sampling data
Check the $$\hat{R}$$ and trace plots.
The convergence check for the posterior distribution is set to $$\hat{R} \leq 1.1$$.
Here, we confirm that there are no parameters with $$\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())[Execution Results]
There were 0 parameters with $$\hat{R} > 1.01$$.
We were able to confirm that all parameters satisfy $$\hat{R} \leq 1.1$$.
Let's roughly check the summary statistics and trace plots of the posterior distribution sampling data.
### 推論データの要約統計情報の表示
az.summary(idata_n_nrt, hdi_prob=0.95, round_to=3)[Execution Results]
Check the state of the posterior distribution sampling data using trace plots.
### トレースプロットの表示
az.plot_trace(idata_n_nrt, compact=True)
plt.tight_layout();[Execution Results]
From the graph on the left, it can be seen that the four Markov chains have almost the same distribution.
In the graph on the right, the lines are drawn evenly.
It is considered to have converged.
(5) Saving inference data
Let's save the inference data to a file in case we (might) want to reuse it later.
We will save it using pickle.
### idataの保存 pickle
file = r'idata_n_nrt_ch18.pkl'
with open(file, 'wb') as f:
pickle.dump(idata_n_nrt, f)The code for loading is as follows.
### idataの読み込み pickle
file = r'idata_n_nrt_ch18.pkl'
with open(file, 'rb') as f:
idata_n_nrt_load = pickle.load(f)Let's take a short break.
Model 1 Results and Analysis
1. Fixed effects of items
The text focuses on the fixed effects of items $${item}$$.
We will draw a forest plot of the posterior distribution sample data for the fixed effects (coefficients) of the items.
This corresponds to Figures 18.4 to 18.7 in the text.
### フォレストプロットの描画 ※図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()[Execution Results]
[Analysis - Following the text]
In both the linear response tree model and the nested response tree model, the fixed effects of the meaningful condition show greater variance, while the variance of the effects in the meaningless condition is smaller.
As the text says, "It can be seen that participants react differently to meaningful words and meaningless strings"!
The following are my personal thoughts.
In the case of meaningful words, general cognitive familiarity or how easily an image remains might determine memorability.
Specifically, I thought it possible that words used frequently are judged as easier to remember because "they are experienced often," while words related to things rarely seen are judged as harder to remember because "experience is limited."
On the other hand, in the case of "symbolic" words made of meaningless hiragana sequences, experience with the words or the things they represent is not useful, so it is possible that memorability is determined uniformly.
2. Fixed effects of nodes
Next, we proceed to check the fixed effects of nodes $${node}$$.
In the text, we examine the summary statistics of the fixed effects and the estimated values of the "main effects" obtained by fitting the estimation results into the range of 0 to 1 using a logit link function.
Now, let's calculate the summary statistics.
This corresponds to Tables 18.1 and 18.2 in the text.
### 線形反応ツリーモデルにおけるノードの固定効果の事後統計量 ※表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))[Execution Results]
Although there are slight differences from the values in the text, the trends are likely similar.
### 線形反応ツリーモデルにおけるノードの固定効果の事後統計量 ※表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))[Execution Results]
These also differ slightly from the values in the text.
Next are the main effects of the nodes.
First, we create data fitted into the range of 0 to 1 using the logit link function.
We create them separately for the linear response tree model and the nested response tree model.
### 事後分布サンプリングデータをロジットリンク関数で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)[Execution Results]
Main effects of nodes in the linear response tree model.
### 事後分布サンプリングデータをロジットリンク関数で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)[Execution Results]
Main effects of nodes in the nested response tree model.
Now, let's proceed to plotting.
We will organize the data to make it easier to plot and then draw it.
This corresponds to Figures 18.8 and 18.9 in the text.
Starting with the linear response tree model.
### ノードの主効果の描画のためのデータの作成 線形反応ツリーモデル
# 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)[Execution Results]
### 線形反応ツリーモデルにおけるノード主効果の描画 ※図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='条件');[Execution Results]
This is the plot for the linear response tree model.
I think it looks very similar to the plot in the text!
Next is the nested response tree model.
### ノードの主効果の描画のためのデータの作成 入れ子反応ツリーモデル
# 条件・ノードごとの平均・95%CI区間を算出してデータフレーム化
stats_df4_plot = (stats_df4.groupby(['condition', 'node'])['estimate']
.agg(['mean', q025, q975]).reset_index())
# データフレームの表示
display(stats_df4_plot)[Execution Results]
### 入れ子反応ツリーモデルにおけるノード主効果の描画 ※図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='条件');[Execution Results]
This is the plot for the nested response tree model.
This is also close to the plot in the text!
【Analysis - Following the text】
The text points out the following two points:
1. In the linear response tree model, a branch of 0 tends to be selected from node $${Y^*_1}$$ to node $${Y^*_3}$$ regardless of the conditions.
→ It is thought that participants are judging the ease of memorizing items in stages.
2. In both tree models, it is more likely that a branch to 1 (the difficult direction) is chosen for nonsense syllables.
I understood the second point!
It is intuitively easy to understand that meaningful words are easier to memorize and nonsense words are harder to memorize.
Regarding the first point, I could not read the content of the point from the data (sweat).
And the text proceeds to the next analysis by selecting the "nested response tree model" for now.
Construction of Models 2 and 3
The text adds two Bayesian models to verify whether the judgments made at the three nodes are "similar" or "qualitatively different".
The data used is the "nested response tree model".
【Model 2】
Model 2 is a model that assumes different judgments are made at each node, and adds a random slope for each node to the Bayesian model.
【Model 3】
Model 3 is a model that assumes different judgments are made between the "first node $${Y^*_1}$$" and the "remaining nodes $${Y^*_2,\ Y^*_3}$$", and adds different random slopes for node $${Y^*_1}$$ and nodes $${Y^*_2,\ Y^*_3}$$ to the Bayesian model.
【Perspective of Analysis】
Using the model evaluation indices WAIC and loo, we will verify the "judgments made at the nodes" (similar/different) by choosing a "good model" from models 1, 2, and 3.
Now, let's start building from Model 2.
Mathematical expression of Model 2
This is the "pseudo-mathematical" notation for the Bambi model I want to aim for.
I will write it in formula format.
$$
value \sim 0 + node + item + (0 + node\ |\ person)
$$
This is a model with the random slope "$${(0 + node\ |\ person)}$$" for the node $${node}$$ added.
The objective variable $${value}$$ follows a Bernoulli distribution.
Model 2 - Meaningful condition
(1) Definition of the model
### モデルの定義
model_m_nrt2 = bmb.Model(
formula='value ~ 0 + node + item + (0 + node | person)', # フォーミュラ式
data=data_m_nrt, # データ
family='bernoulli', # 目的変数の誤差分布
)(2) Confirmation of the model appearance
### モデルの表示
model_m_nrt2【Execution result】
Bambi automatically sets the prior distribution that each parameter follows by analyzing the data.
### モデルの可視化
model_m_nrt2.build()
model_m_nrt2.graph()【Execution result】
(3) Sampling from the posterior distribution
The processing time was approximately 10 seconds.
### 事後分布からのサンプリング ※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)【Execution result】 Omitted
(4) Confirmation of sampling data
Check $${\hat{R}}$$ and the trace plot.
The convergence check of the posterior distribution is set to $${\hat{R} \leq 1.1}$$.
Here, I will confirm that there are no parameters with $${\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())[Execution Result]
There were 0 parameters with $$\hat{R} > 1.01$$.
We were able to confirm that all parameters satisfy $$\hat{R} \leq 1.1$$.
Let's roughly check the summary statistics and trace plots of the posterior distribution sampling data.
### 推論データの要約統計情報の表示
az.summary(idata_m_nrt2, hdi_prob=0.95, round_to=3)[Execution Result]
Let's check the state of the posterior distribution sampling data using trace plots.
### トレースプロットの表示
az.plot_trace(idata_m_nrt2, compact=True)
plt.tight_layout();[Execution Result]
From the graph on the left, we can see that the four Markov chains have almost the same distribution.
In the graph on the right, the lines are drawn evenly.
It appears to have converged.
(5) Saving Inference Data
Let's save the inference data to a file in case we (might) want to reuse it.
We will save it using pickle.
### idataの保存 pickle
file = r'idata_m_nrt2_ch18.pkl'
with open(file, 'wb') as f:
pickle.dump(idata_m_nrt2, f)The code for loading is as follows.
### idataの読み込み pickle
file = r'idata_m_nrt2_ch18.pkl'
with open(file, 'rb') as f:
idata_m_nrt2_load = pickle.load(f)
Model 2: Meaningless Condition
(1) Model Definition
### モデルの定義
model_n_nrt2 = bmb.Model(
formula='value ~ 0 + node + item + (0 + node | person)', # フォーミュラ式
data=data_n_nrt, # データ
family='bernoulli', # 目的変数の誤差分布
)(2) Checking Model Appearance
### モデルの表示
model_n_nrt2[Execution Result]
Bambi analyzes the data and automatically sets the prior distributions that each parameter follows.
### モデルの可視化
model_n_nrt2.build()
model_n_nrt2.graph()[Execution Result]
(3) Sampling from the Posterior Distribution
The processing time was approximately 10 seconds.
### 事後分布からのサンプリング ※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)[Execution Result] Omitted
(4) Checking Sampling Data
We will check $$\hat{R}$$ and the trace plots.
We will use $$\hat{R} \leq 1.1$$ to confirm convergence of the posterior distribution.
Here, we will confirm that there are no parameters with $$\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())[Execution Result]
There were 0 parameters with $$\hat{R} > 1.01$$.
We were able to confirm that all parameters satisfy $$\hat{R} \leq 1.1$$.
Let's roughly check the summary statistics and trace plots of the posterior distribution sampling data.
### 推論データの要約統計情報の表示
az.summary(idata_n_nrt2, hdi_prob=0.95, round_to=3)[Execution Result]
Let's check the state of the posterior distribution sampling data using trace plots.
### トレースプロットの表示
az.plot_trace(idata_n_nrt2, compact=True)
plt.tight_layout();[Execution Result]
From the graph on the left, we can see that the four Markov chains have almost the same distribution.
In the graph on the right, the lines are drawn evenly.
It appears to have converged.
(5) Saving inference data
Let's save the inference data to a file in case we (might) want to reuse it later.
We will save it using pickle.
### idataの保存 pickle
file = r'idata_n_nrt2_ch18.pkl'
with open(file, 'wb') as f:
pickle.dump(idata_n_nrt2, f)The code for loading it is as follows.
### idataの読み込み pickle
file = r'idata_n_nrt2_ch18.pkl'
with open(file, 'rb') as f:
idata_n_nrt2_load = pickle.load(f)Moving on to building Model 3.
Mathematical expression for Model 3
This is a "pseudo-mathematical" representation of the Bambi model we want to aim for.
We will write it in formula format.
$$
value \sim 0 + node + item + (0 + node2\ |\ person)
$$
$${node2}$$ is a categorical variable set to True for nodes $${Y^*_2, Y^*_3}$$ and False for node $${Y^*_1}$$.
This is a model with the random slope "$${(0 + node2\ |\ person)}$$" for $${node2}$$ added.
The objective variable $${value}$$ follows a Bernoulli distribution.
We will set $${node2}$$ in the data.
### データの前処理
# node2,3フラグを追加
data_m_nrt['node2'] = (data_m_nrt['node'] != 'node1').astype('category')
data_n_nrt['node2'] = (data_n_nrt['node'] != 'node1').astype('category')
Model 3: Meaningful condition
(1) Model definition
### モデルの定義
model_m_nrt3 = bmb.Model(
formula='value ~ 0 + node + item + (0 + node2 | person)', # フォーミュラ式
data=data_m_nrt, # データ
family='bernoulli', # 目的変数の誤差分布
)(2) Checking the model structure
### モデルの表示
model_m_nrt3[Execution Result]
Bambi will analyze the data and automatically set the prior distributions that each parameter follows.
### モデルの可視化
model_m_nrt3.build()
model_m_nrt3.graph()[Execution Result]
(3) Sampling from the posterior distribution
The processing time was approximately 10 seconds.
### 事後分布からのサンプリング ※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)[Execution Result] Omitted
(4) Checking the sampling data
We will check $${\hat{R}}$$ and the trace plots.
We will use $${\hat{R} \leq 1.1}$$ to confirm the convergence of the posterior distribution.
Here, we will confirm that there are no parameters with $${\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())[Execution Result]
There were 0 parameters with $${\hat{R} > 1.01}$$.
We were able to confirm that all parameters satisfy $${\hat{R} \leq 1.1}$$.
We will roughly check the summary statistics and trace plots of the posterior distribution sampling data.
### 推論データの要約統計情報の表示
az.summary(idata_m_nrt3, hdi_prob=0.95, round_to=3)[Execution Result]
We will check the state of the posterior distribution sampling data using trace plots.
### トレースプロットの表示
az.plot_trace(idata_m_nrt3, compact=True)
plt.tight_layout();[Execution Result]
From the graph on the left, we can see that the four Markov chains have almost the same distribution.
In the graph on the right, the lines are drawn evenly.
It appears to have converged.
(5) Saving inference data
Let's save the inference data to a file in case we (might) want to reuse it later.
We will save it using pickle.
### idataの保存 pickle
file = r'idata_m_nrt3_ch18.pkl'
with open(file, 'wb') as f:
pickle.dump(idata_m_nrt3, f)The code for loading it is as follows.
### idataの読み込み pickle
file = r'idata_m_nrt3_ch18.pkl'
with open(file, 'rb') as f:
idata_m_nrt3_load = pickle.load(f)
Model 3: Nonsense condition
(1) Model definition
### モデルの定義
model_n_nrt3 = bmb.Model(
formula='value ~ 0 + node + item + (0 + node2 | person)', # フォーミュラ式
data=data_n_nrt, # データ
family='bernoulli', # 目的変数の誤差分布
)(2) Checking the model structure
### モデルの表示
model_n_nrt3[Execution result]
Bambi analyzes the data and automatically sets the prior distributions that each parameter follows.
### モデルの可視化
model_n_nrt3.build()
model_n_nrt3.graph()[Execution result]
(3) Sampling from the posterior distribution
The processing time was approximately 10 seconds.
### 事後分布からのサンプリング ※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)[Execution result] Omitted
(4) Checking sampling data
We will check the R-hat and trace plots.
We will use R-hat <= 1.1 to confirm convergence of the posterior distribution.
Here, we confirm that there are no parameters with R-hat > 1.01.
### r_hat>1.1の確認
# 設定
idata_in = idata_n_nrt3 # idata名
threshold = 1.01 # しきい値
# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())[Execution result]
There were 0 parameters with R-hat > 1.01.
We were able to confirm that all parameters have R-hat <= 1.1.
We will roughly check the summary statistics and trace plots of the posterior distribution sampling data.
### 推論データの要約統計情報の表示
az.summary(idata_n_nrt3, hdi_prob=0.95, round_to=3)[Execution result]
We will check the state of the posterior distribution sampling data using trace plots.
### トレースプロットの表示
az.plot_trace(idata_n_nrt3, compact=True)
plt.tight_layout();[Execution result]
From the graph on the left, we can see that the four Markov chains are following almost the same distribution.
In the graph on the right, the lines are drawn evenly.
It appears to have converged.
(5) Saving inference data
Let's save the inference data to a file in case we (might) want to reuse it later.
We will save it using pickle.
### idataの保存 pickle
file = r'idata_n_nrt3_ch18.pkl'
with open(file, 'wb') as f:
pickle.dump(idata_n_nrt3, f)The code for loading it is as follows.
### idataの読み込み pickle
file = r'idata_n_nrt3_ch18.pkl'
with open(file, 'rb') as f:
idata_n_nrt3_load = pickle.load(f)Taking a break. Phew.
Comparison of Models 1, 2, and 3
We will select a good model from the three using WAIC and loo.
We will create a table corresponding to Table 18.3 in the text.
### 入れ子反応ツリーモデルにおけるモデル１～３の適合度指標 ※表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}'))[Execution Results]
For both WAIC and loo, a smaller estimated value indicates a better model.
[Analysis - Following the Text]
■ In the case of the meaningful condition
Model 1, which does not assume random effects for each node, is the best model.
The text states that "judgments of memorability in the meaningful condition may be unified rather than differing by node."
■ In the case of the meaningless condition
Model 2, which incorporates random slopes for each node, is the best model.
The text states that "it is inferred that judgments of memorability in the meaningless condition are not unified, and that qualitatively different judgments are made at each node."
By the way, as the WAIC and loo values below indicate, the linear response tree model for the meaningful condition might be a "better model" than models 1-3...
$$
\begin{array}{}
\text{WAIC} & \text{loo} \\
\hline385.88 & 387.33
\end{array}
$$
This concludes Chapter 18.
Conclusion
A large number of models
In this chapter, we built and compared eight Bayesian models.
I am feeling somewhat confused by the appearance of so many models.
Moreover, with many model elements such as "conditions" × "item response trees" × "random effects," the level of confusion is increasing!
It reminded me of a complex order spell.
Is there also a psychological judgment mechanism for the memorability of models?
Professor, please write a paper on it.
The End
Series Articles
Next Article
Previous Article
Table of Contents
Blog Introduction
I am writing seven series articles on note.
Please come and take a look!
1. Relaxed Statistics
This is a blog that roughly explores probability and statistics using the Statistical Test Grade 2 problem collection as a guide.
Feel free to read it like casual conversation. Please come and take a look.
It corresponds to the Statistical Test Grade 2 Official Problem Collection CBT Edition.
Sample code for Python and EXCEL is also available.
2. Experiment! Fun Bayesian Modeling 1 & 2 with PyMC Ver. 5
I will draw and analyze the Bayesian models used in psychological research from the books "Fun Bayesian Modeling" and "Fun Bayesian Modeling 2" using PyMC Ver. 5.
Like this book, many Bayesian models are written in R + Stan.
I strive to explore the possibilities of PyMC and make Bayesian modeling easy to practice.
Since these are familiar and easy-to-visualize themes, please try running them with PyMC and let's enjoy it together!
3. Experiment! Bayesian Modeling from Iwanami Data Science 1 with PyMC Ver. 5
I will draw and analyze the Bayesian models by four Bayesians from the book "Experiment! Iwanami Data Science Vol. 1" using PyMC Ver. 5.
This book is a great resource for learning the basics of Bayesian programming.
I feel like I've become friends with Bayesian methods by happily running PyMC models.
Everyone, please try running them with PyMC and let's play and learn together!
4. Fun Copy-Coding: Bayesian, Python, etc.
I will blog about the results of my "book copy-coding activities" for Bayesian, Python, and others.
I am mainly working on translations into Python.
I hope this serves as sample code for fellow copy-coders 🍀
5. Introduction to Time Series Analysis for Psychology with R and Stan, using Python and PyMC Ver. 5
I will practice the time series analysis from the book "Introduction to Time Series Analysis for Psychology with R and Stan" using Python and PyMC Ver. 5.
This book is packed with time series analysis themes!
I realized the depth of time series analysis.
I will enjoy learning time series analysis with my favorite language, Python.
6. Writing about Data Science-like things
I write columns on statistics, data analysis, AI, machine learning, and Python on an irregular basis.
There are many articles related to statistics and data science books.
Series on "Statistics," "Python," "Mathematics and Python," and "R" have been created.
7. Python Machine Learning Programming Practice Log
I wrote articles about my various thoughts when studying the book "Python Machine Learning Programming: PyTorch & scikit-learn Edition."
This book is a textbook for scikit-learn and PyTorch.
Please feel free to give it a try.
Thank you very much for reading until the end.
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！