Bayesian Modeling of Chapter 13 "Do Negative Comments Lower the Productivity of Other Members?" with PyMC Ver. 5

This article is a statistical documentary depicting the process of "experimentally" implementing Chapter 13 "Do Negative Comments from Group Members Lower the Performance of Others?" from the text "Fun Bayesian Modeling" using a Bayesian model in PyMC Ver. 5. In this chapter, we investigate whether the performance of a group task decreases when one of the members makes

negative comments such as "this is a boring task" during group work.

This time, I was able to obtain inference results close to those in the text!
I felt frustrated last time because I had to "abstain," so I am very happy with this result (tears).

Do negative comments really lower the motivation of members? Let's
model it with PyMC and enjoy Bayesian statistics!

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 obtained by downloading them from the publisher's website.

Summary

Overview of the text

Author: Dr. Takashi Goto
Model Difficulty: ★★・・・ (Easy)

Self-Evaluation

Rating

$$
\begin{array}{c:c:c}
Implementation Accuracy & ★★★★★& GoooD!! \\
Result Reproduction & ★★★★★& Best✨ \\
Fun & ★★★★★& Fun! \\
\end{array}
$$

Evaluation Points

• Actually, when I first started writing this article, the inference results of my own model were quite different from those in the text. By persistently reviewing the model, I was finally able to reach results close to those in the text.
  

Efforts, Joy, and Reflections

• I could not easily read whether the units of the intercepts and slopes in the model equations of the text corresponded to "overall units," "group units," or "individual units," so I digested the Stan model descriptions to understand the units. My ability to read R and Stan code is low, but I would like to try to understand the code carefully, even if it is at a low level, in the future.
  

Overview of the model

Overview of the survey and experiment in the text

■ Social Loafing
In group work where individual contributions are difficult to observe, we conduct an experiment to see if negative comments from one member cause other members to lose motivation, leading to collective social loafing, and collect data.
Using the collected data, we perform an analysis of two questions using Bayesian modeling that assumes a "group level" to "individual level" hierarchy.

■ Pre-experiment Planning and Procedures
The text describes the raw reality of the experiment, such as setting hypotheses, creating an analysis plan, and revising the plan to handle unforeseen events on the day of the experiment (e.g., participants not showing up).
The details of the experiment are carefully explained using diagrams.
It is very helpful!

Modeling in the Text

■ Objective Variable and Parameters of Interest
The objective variable $${post_i}$$ is the amount of work done in the second half.
The parameter of most interest is $${\gamma_{03}}$$, which is the coefficient for the experimental manipulation (presence or absence of negative comments).

■ Model
The first line of the equation shows the "individual level."
$${b_0}$$ in the second line: individual-level intercept, and $${b_1}$$ in the third line: individual-level slope, indicate the "group level."
The subscript $${i}$$ is the participant identification number, and $${g[i]}$$ is the identification number of the group $${g}$$ to which participant $${i}$$ belongs.

The following is the model that integrates the above equations for implementation.

【Variable Definitions】
$${post_i}$$: Amount of work done by participant $${i}$$ in the second half
$${\mu_i}$$: Mean value of the amount of work done by participant $${i}$$ in the second half
$${\sigma}$$: Standard deviation of the work done in the second half
$${pre(gc)_i}$$: Amount of work done by participant $${i}$$ in the first half - Mean amount of work done in the first half by the group to which the participant belongs (individual-level first-half work)
$${pre(gm)_{g[i]}}$$: Mean amount of work done in the first half for each group (group-level first-half work)
$${peerN_{g[i]}}$$: Number of people in the group
$${cond_{g[i]}}$$: Presence or absence of experimental manipulation
$${\gamma_{00}}$$: Intercept (Note)
$${\gamma_{10}}$$: Coefficient for individual-level first-half work (Note)
$${\gamma_{01}}$$: Coefficient for group-level first-half work
$${\gamma_{02}}$$: Coefficient for group size
$${\gamma_{03}}$$: Coefficient for experimental manipulation
$${\tau_{00}}$$: Standard deviation of $${\gamma_{00}}$$
$${\tau_{10}}$$: Standard deviation of $${\gamma_{10}}$$

(Note) Due to Stan implementation, $${\gamma_{00}}$$ and $${\gamma_{10}}$$ are group-level variables.

■ Analysis and Results
I believe the text's description of the analysis methods and numerical results is accurate, so I recommend reading the text.
For the analysis results using PyMC modeling, please see the "[5. Analysis]" section in the "PyMC Implementation" chapter.

PyMC Implementation

Let's enjoy PyMC & Python !

Preparation and Data Verification

1. Import

### インポート

# ユーティリティ
import pickle

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

# PyMC
import pymc as pm
import arviz as az

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

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

2. Data Loading and Preprocessing
Load the csv file into a pandas DataFrame.

### データの読み込み

# データの読み込み
data_orgn = pd.read_csv('loafing_data.csv')

# データの表示
print('data_orgn.shape:', data_orgn.shape)
display(data_orgn.head())

[Execution Result]
There are 39 rows and 6 columns in total.
The data items are as follows:
id: Participant identification number (39 participants in total)
group: Group identification number (17 groups in total)
cond: Presence or absence of negative speech (0: No speech/control, 1: Speech present/experimental)
pre.work: Work volume in the first half
post.work: Work volume in the second half (objective variable)
peerN: Number of people in the group

3. Data Overview and Statistics

First, let's check the summary statistics.

### 要約統計量の表示
data_orgn.describe().round(2)

[Execution Result]
Overall, it seems that 'the work volume in the second half is larger'.
Perhaps proficiency in the work increased, or the motivation to produce more than in the first half increased?

Following the text, I will draw a box plot.
I will use matplotlib's boxplot.
This corresponds to Figure 13.2 in the text.

### 条件ごとの作業量の箱ひげ図の描画 ★図13.2に相当

fig, ax = plt.subplots(figsize=(6, 3))
# 箱ひげ図の描画
ax.boxplot([data_orgn[data_orgn['cond']==0]['pre.work'],
            data_orgn[data_orgn['cond']==1]['pre.work'],
            data_orgn[data_orgn['cond']==0]['post.work'],
            data_orgn[data_orgn['cond']==1]['post.work']])
# 修飾
ax.set(xticklabels=['前半－統制群', '前半－実験群', '後半－統制群', '後半－実験群'],
       ylabel='作業量')
plt.show()

[Execution Result]
It feels like the experimental group with negative speech has greater variation in both the first and second halves.
Looking at the second-half work, it seems that the average work volume is smaller in the experimental group with negative speech.

I will perform data preprocessing.
I will calculate the first-half work volume for each group (the average value of the work volume of the individuals belonging to it).
Also, I will center the first-half work volume and the number of people in the group (subtract the average value).

### データの前処理

# 集団ごとの前半作業量の平均preGmeanの算出
tmp = (data_orgn.groupby(['group'])[['pre.work']].mean().reset_index()
       .rename(columns={'pre.work': 'preGmean'}))
data = pd.merge(data_orgn, tmp, on='group', how='left')

# pre.workの集団中心化
data['pre.work'] = data['pre.work'] - data['preGmean']

# peerNの全体中心化
data['peerN'] = data['peerN'] - data['peerN'].mean()

# preGmeanの全体中心化
data['preGmean'] = data['preGmean'] - data['preGmean'].mean()

# 加工後のデータの表示
print('data.shape:', data.shape)
display(data.head())

[Execution Result]

Model Construction

Mathematical expression of the model
This is a 'pseudo-mathematical' notation that mixes the atmosphere of the PyMC model I want to aim for.
The units of $$\gamma_{00}$$ and $$\gamma_{10}$$ are 'group units'.

$$
\begin{align*}
\gamma_{01}, \gamma_{02}, \gamma_{03}, \mu_{\gamma_{00}}, \mu_{\gamma_{10}} &\sim \text{Normal}\ (\text{mu}=0,\ \text{sigma}=100)\\
\tau_{00}, \tau_{00} &\sim \text{HalfCauchy}\ (\text{beta}=100)\\
\gamma_{00}[group] &\sim \text{Normal}\ (\text{mu}=\mu_{\gamma_{00}},\ \text{sigma}=\tau_{00}) \\
\gamma_{10}[group] &\sim \text{Normal}\ (\text{mu}=\mu_{\gamma_{10}},\ \text{sigma}=\tau_{10}) \\
\mu &= \gamma_{00}[group] + \gamma_{10}[group] \times pre(gc) + \gamma_{01} \times pre(gm) + \gamma_{02} \times peerN + \gamma_{03} \times cond\\
\sigma &\sim \text{HalfCaushy}\ (\text{beta}=100) \\
likelihood &\sim \text{Normal}\ (\text{mu}=\mu,\ \text{sigma}=\sigma) \\
\end{align*}
$$

1. Model Definition
In the initial value setting, I am indexing the group values by subtracting 1.
I use this for the indices of $$\gamma_{00}, \gamma_{10}$$.
Also, since this model has many explanatory variables, the number of 'data' definitions is also large.

### モデルの定義

# 初期値設定：集団groupの番号とインデックスを作成
group_cat = sorted(data['group'].unique())  # 集団の番号を1から昇順で取得
group_idx = list(data['group'].values - 1)  # 列group値を0始まりに変更

with pm.Model() as model:
    
   ### データ関連定義
   ## coordの定義
   model.add_coord('id', values=data['id'].values, mutable=True)
   model.add_coord('group', values=group_cat, mutable=True)
   
   ## dataの定義
   # 目的変数：個人レベルの後半の作業量post.work
   y = pm.ConstantData('y', value=data['post.work'], dims='id')
   # 個人が所属する集団groupのインデックス
   groupIdx = pm.ConstantData('groupIdx', value=group_idx, dims='id')
   # 個人レベルの前半の作業量pre.work_c
   preGc = pm.ConstantData('preGc', value=data['pre.work'], dims='id')
   # 集団レベルの前半の作業量pre.work_m_gm
   preGm = pm.ConstantData('preGm', value=data['preGmean'], dims='id')
   # 集団の人数peerN_g
   peerN = pm.ConstantData('peerN', value=data['peerN'], dims='id')
   # 実験操作の有無(0:統制、１：実験)cond
   cond = pm.ConstantData('cond', value=data['cond'], dims='id')
   
   ### 事前分布
   # 集団レベルの前半の作業量の係数
   gamma01 = pm.Normal('gamma01', mu=0, sigma=100)
   # 集団の人数の係数
   gamma02 = pm.Normal('gamma02', mu=0, sigma=100)
   # 実験操作（発言の有無）の係数
   gamma03 = pm.Normal('gamma03', mu=0, sigma=100)
   # 切片γ00のmu
   muGamma00 = pm.Normal('muGamma00', mu=0, sigma=100)
   # 個人レベルの前半の作業量γ10のmu
   muGamma10 = pm.Normal('muGamma10', mu=0, sigma=100)
   # 切片γ00の標準偏差
   tau00 = pm.HalfCauchy('tau00', beta=100)
   # 個人レベルの前半の作業量の係数γ10の標準偏差
   tau10 = pm.HalfCauchy('tau10', beta=100)
   # 切片
   gamma00 = pm.Normal('gamma00', mu=muGamma00, sigma=tau00, dims='group')
   # 個人レベルの前半の作業量の係数
   gamma10 = pm.Normal('gamma10', mu=muGamma10, sigma=tau10, dims='group')
   # 後半の作業量の標準偏差
   sigma = pm.HalfCauchy('sigma', beta=100)
   
   ### 後半の作業量の期待値muの計算
   mu = pm.Deterministic('mu',
      (gamma00[groupIdx] + gamma10[groupIdx] * preGc + gamma01 * preGm
       + gamma02 * peerN + gamma03 * cond),
      dims='id')
   
   ### 尤度
   likelihood = pm.Normal('likelihood', mu=mu, sigma=sigma, observed=y, 
                          dims='id')
   
   ### 計算値
   # 効果量
   es = pm.Deterministic('es', gamma03 / sigma)
   # 切片の平均値に対する平均値差の割合
   rate = pm.Deterministic('rate', gamma03 / muGamma00)

[Model Annotation]

• Definition of coord
  You can name coordinates and set the values that those coordinates can take.
  This time, I set the following two:
  - Row coordinates: Name 'data', value 'row index'
  - Group coordinates: Name 'group', value 'integer from 1 to 17'

• Definition of data
  I set the objective variable $${y}$$, the group index $${groupIdx}$$, and four explanatory variables.

• Prior distribution of parameters
  I assume a normal distribution for the coefficients of the explanatory variables, etc., and a half-Cauchy distribution for the standard deviation.

• Calculation of mu
  It is as per the calculation of the expected value mu of the second-half work volume defined in the model.

• Likelihood
  It follows a normal distribution with mu and sigma as parameters.

• Calculated values
  I will calculate the 'effect size' and 'ratio' used in '13.6 Results' of the text.

2. Checking the model structure

### モデルの表示
model

[Execution result]
In the Stan code, multiple vectors are placed within variables, so the number of variables is small.
On the other hand, in this model, I have increased the number of variables to stay as close as possible to the variables in the text's mathematical formulas and the variables representing the inference results.

### モデルの可視化
pm.model_to_graphviz(model)

[Execution result]
Various random variables and explanatory variables are tightly gathered around "mu" as a hub.

3. Sampling from the posterior distribution
The number of random number generations (draws, tune) is the same as in the text.
By setting nuts_sampler='numpyro', you can use numpyro as the NUTS sampler.
The processing time was approximately 3 minutes and 30 seconds.

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを未使用　3分30秒
with model:
    idata = pm.sample(draws=200000, tune=200000, chains=4, target_accept=0.95,
                      nuts_sampler='numpyro', random_seed=1234)

[Execution result] Omitted

4. Checking sampling data
Check the R-hat and trace plots.
The convergence check for the posterior distribution is set to R-hat <= 1.1.

### r_hat>1.1の確認
rhat_idata = az.rhat(idata)
(rhat_idata > 1.1).sum()

[Execution result]
The number of parameters with R-hat > 1.1 is "0".
I was able to confirm that all parameters have R-hat <= 1.1.

Summary statistics of the inference data.

### 推論データの要約統計量
pm.summary(idata, hdi_prob=0.95)

[Execution result]
Because there are many variables to be inferred, it was partially omitted (sweat).

Trace plots.

### トレースプロットの描画
pm.plot_trace(idata, var_names=var_names, combined=True, compact=False,
              figsize=(15, 20))
plt.tight_layout();

[Execution result]
There are quite a few barcodes showing divergence...
In the distribution on the left, there are strange, high-density areas that look like protrusions...

5. Analysis
Check the summary statistics of the main variables corresponding to Table 13.1 in the text.

### 推論データの要約統計量の表示 ★表13.1に対応
var_names = ['muGamma00', 'muGamma10', 'gamma01', 'gamma02', 'gamma03',
             'tau00', 'tau10', 'sigma', 'es', 'rate']
pm.summary(idata, hdi_prob=0.95, var_names=var_names)

[Execution result]
The mean (mean, EAP in the text) and standard deviation (sd, post.sd in the text) were close to the values in the text (phew).

[Analysis]

1. Probability that the workload when there is "negative speech" is lower than when there is no speech

The estimated value of gamma_03 is -2.982 (95% HDI = [-9.079, 3.139]).
On average, negative speech likely has a negative impact.
Calculate the probability that gamma_03 is less than 0.

### γ03が0より小さい確率の計算
gamma03_sample = idata.posterior.gamma03.data.flatten()
print('P(γ03<0)=', (gamma03_sample < 0).sum() / gamma03_sample.size)

[Execution result]

The probability of being less than 0 is 85%.
The result is that there is an 85% probability that the average workload when there is negative speech is lower than when there is no negative speech.

2. Decrease in average workload when there is "negative speech"

The estimated value of the effect size (es) is -0.495 (95% HDI = [-1.552, 0.523]).
The result is that the average workload when there is negative speech is 49.5% of the standard deviation smaller than when there is no negative speech.
Also, the estimated value of the rate is -0.048 (95% HDI = [-0.145, 0.051]).
The result is that the average workload when there is negative speech is 4.8% smaller than when there is no negative speech.
However, since both the effect size and the rate include 0 or more in the 95% HDI range, caution is required when interpreting "smaller".

6. Saving Inference Data (idata)
Let's save the inference data to a file in case you need to reuse it.
Save idata using pickle.

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

The code for loading is as follows.

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

This concludes Chapter 13.

Conclusion

Difficulties in PyMC Modeling and Breaking Through Limits

Among the models up to Chapter 13, there were eight chapters where I felt hesitant and thought, "PyMC implementation seems impossible," when I first saw the text.
There was a possibility that I would have given up on about 60% of the chapters.

This Chapter 13 was also one of those "seems impossible" chapters.
There were many mathematical formulas for the model, and at a glance, I assumed it was "complex."
However, practicing each chapter became a form of discipline, and repeating the process of clearing the modeling became training, allowing me to successfully implement Chapter 13 in PyMC (I overcame it).

Series Articles

Next Article

Previous Article

Table of Contents

Blog Introduction

I am writing four series of articles on note.
Please feel free to 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 as if it were casual conversation. Please take a look.
It corresponds to the Statistical Test Grade 2 Official Problem Collection CBT version.

2. Introduction to Time Series Analysis for Psychology with R and Stan, using Python and PyMC Ver. 5
I am working on the time series analysis topics from the book "Introduction to Time Series Analysis for Psychology with R and Stan" using Python and PyMC Ver. 5.
I believe that practicing with a wealth of themes (topics) will lead to building basic strength in Python and PyMC.
Every day, I work hard at web searches, understanding time series models, grasping Python packages, and translating R/Stan code!
I hope this series will be helpful for beginners in Python time series analysis 🍀

3. Practical Record of Python Machine Learning Programming
I wrote articles about my various thoughts when I studied the book "Python Machine Learning Programming: PyTorch & scikit-learn Edition."
This book is a textbook for scikit-learn and PyTorch.
Please feel free to try it out if you like.

4. Writing About Data Science-like Things
I write columns on statistics, data analysis, AI, machine learning, and Python on an irregular basis.
Series on "Statistics," "Python," "Mathematics and Python," and "R" have been created.
Practical records of Bayesian books are also being posted.

Thank you very much for reading until the end.