実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で⑤直線あてはめ

2024年2月28日 11:46

「階層ベイズ最初の一歩」章

テキスト「階層ベイズ最初の一歩」の執筆者
久保拓弥先生

この記事は、テキストの「階層ベイズ最初の一歩」章のモデル１「直線あてはめ」の実践を取り扱います。
この章は、給食の変更と児童の身長の関係のベイズ分析を紹介しています。
今回の記事範囲では「不備なモデル」を取り扱います。
次回の「階層ベイズ」に向けての第一歩を踏み出しましょう！

はじめに

岩波データサイエンスVol.1の紹介

この記事は書籍「岩波データサイエンス vol.1」（岩波書店、以下「テキスト」と呼びます）の特集記事「ベイズ推論とMCMCのフリーソフト」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。

テキストは、2015年10月に発売され、ベイズモデリングの様々なソフトウェアを用いたモデリング事例を多数掲載し、ベイズモデリングの楽しさを紹介する素晴らしい書籍です。
入門的なモデルから2次階差を取り扱う空間モデルまで、幅広い難易度のモデルを満喫できます！

このシリーズは、テキストのベイズモデルをPyMC Ver.5に書き換えて実践します。

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「岩波データサイエンス vol.1」
第９刷、編者岩波データサイエンス刊行委員会　岩波書店

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

PyMC環境の準備

Anacondaを用いる環境構築とGoogle ColaboratoryでPyMCを動かす方法について、次の記事にまとめています。
「PyMCを動かすまでの準備」章をご覧ください。

PythonとPyMC
テキストで利用するツールは、R、BUGS、JAGS です。
この記事では、PythonとPyMCを用いたコードに変換してベイズモデリングを実践いたします。

モデル１直線あてはめ

1. 分析シナリオ

「給食の種類の変更が児童の身長に影響を与えるか」を分析します。
テキストは10県の男子小学生（12歳）の身長を２回測定した架空の平均身長データを用いています
測定１回目のときは全学校が同じ給食タイプ、１年後の測定２回目はランダムに選ばれた5県が給食タイプを変更し、残りの5県が１回目と同じ給食タイプを継続しているという条件です。

テキストの早い段階で明かされる秘密があります。
それは・・・
「給食タイプの変更は身長の増減に影響がない」
ように架空データが作られたのです。

データの概略を踏まえた上で、テキストは次の２つのベイズモデルを構築します。

モデル１「直線当てはめ」
・線形回帰モデルです。
・間違った結果を出力するモデルです。
モデル２「階層ベイズ」
・「ある工夫」を階層に加えるベイズモデルです。
・モデル１よりも「よりよいモデル」です。

この記事はモデル１を実践します。

2. 準備

「階層ベイズ最初の一歩」章で利用するパッケージをインポートします。

### インポート

# ユーティリティ
import pickle

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

# 確率統計
import statsmodels.formula.api as smf
import scipy.stats as stats

# PyMC
import pymc as pm
import arviz as az

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

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

3. データの登録と前処理

テキストに「岩波データサイエンスのサポートページから例題のデータを格納した，dという名前のdata.frameをダウンロードしてみましょう」と案内があります。
ただ、私には上記データを見つけることができませんでした。
そこで、テキストの表1のデータを引用させていただきます。

### データの登録 ※テキストの表1を引用しました
data_orgn = pd.DataFrame({
    '調査県': ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J'],
    '給食タイプ': ['乙', '乙', '甲', '乙', '乙', '甲', '甲', '甲', '乙', '甲'],
    '標本サイズ1': [55, 53, 55, 53, 58, 55, 58, 59, 56, 56],
    '標本サイズ2': [51, 49, 53, 52, 55, 53, 53, 57, 51, 50],
    '平均身長1': [151.36, 151.56, 152.22, 153.09, 153.22, 153.31, 152.98,
               153.27, 152.67, 155.37],
    '平均身長2': [157.27, 156.83, 157.08, 156.00, 157.24, 157.22, 157.81,
               158.95, 156.82, 161.71],
    '標準偏差1':  [2.94, 3.07, 3.2, 2.65, 3.07, 3.1, 2.49, 3.08, 2.82, 3.1],
    '標準偏差2': [2.98, 3.14, 3.21, 2.64, 3.03, 3.13, 2.45, 3.06, 2.92, 3.21],
})
display(data_orgn)

# 必要に応じてcsvファイルに保存
# data_orgn.to_csv('./data/kyushoku.csv')

【実行結果】
10県のデータです。
給食タイプ甲が１回目の給食と同じ県、乙が２回目に給食の変更があった県です。
標本サイズ・平均身長・標準偏差の末尾1が１回目の測定、末尾2が２回目の測定です。

データの前処理をします。
１回目・２回目の値を縦持ちに変換します。

### データの前処理　1回目・2回目の横持ちから縦持ちへ変換

# コピーの取得
tmp = data_orgn.copy()
# 一部の列名を変更
tmp.rename({'調査県': 'pref', '給食タイプ': 'X'}, axis=1, inplace=True)
# 給食タイプを0:甲、1:乙に変換
tmp['X'] = tmp['X'].apply(lambda x: int(0) if x=='甲' else int(1))
# 標本サイズ、平均身長、標準偏差について、縦持ち変換＆値を01の二値変換
tmp1 = pd.melt(tmp, id_vars=['pref', 'X'], 
               value_vars=['標本サイズ1', '標本サイズ2'],
               var_name='Age', value_name='N')
tmp1['Age'] = tmp1['Age'].apply(lambda x: int(0) if x=='標本サイズ1' else int(1))
tmp2 = pd.melt(tmp, id_vars=['pref', 'X'],
               value_vars=['平均身長1', '平均身長2'],
               var_name='Age', value_name='mean_Y')
tmp3 = pd.melt(tmp, id_vars=['pref', 'X'],
               value_vars=['標準偏差1', '標準偏差2'],
               var_name='Age', value_name='sd_Y')
# ３つのtmpを結合
data = pd.concat([tmp1, tmp2['mean_Y'], tmp3['sd_Y']], axis=1)
# 列順の変更
data = data[['pref', 'N', 'mean_Y', 'sd_Y', 'Age', 'X']]
display(data)

【実行結果】
テキスト２３ページの「表１を並びかえたもの」が出来上がりました。
prefは調査県、Nは標本サイズ、mean_Yは平均身長、sd_Yは標準偏差、Ageは測定時期（0は最初の測定、1は１年後：２回目の測定）、Xは給食の変更有無（0：変更なし甲、1：変更あり乙）です。

4. データの外観の確認

県別に平均身長の時系列（１回目と２回目）の折れ線グラフを可視化します。
テキストの図1に相当します。

### 県別平均身長の推移※図1に相当
# 描画用データの作成
data_plot = data.copy()
data_plot.columns = ['調査県', '標本サイズ', '平均身長[cm]', '標準偏差', '測定回',
                     '給食タイプ']
# 描画領域の指定
plt.figure(figsize=(3, 4))
# 折れ線グラブの描画
ax = sns.lineplot(data=data_plot, x='測定回', y='平均身長[cm]', hue='調査県',
                  palette='turbo', style='給食タイプ', markers=True)
# 修飾
ax.set(xticks=(0, 1))
ax.legend(bbox_to_anchor=(1, 1))
ax.grid(lw=0.5);

【実行結果】
直線は右肩上がりです。全ての県で１年間に身長が伸びたことが分かります。

点線が給食を変更した県です。
軒並み平均身長が低い感じです。このことが、給食の変更で身長の伸びが低くなるという錯覚に陥る原因になります。

5. 線形回帰モデルへの当てはめ

テキスト２１ページからの『2 モデル１－「直線的あてはめ」を使ってかたづけてしまおう』から始まる線形回帰モデルをPythonで実践します。

ここで用いるモデルは次のような式で表されます。

$$
\begin{align*}
Y_{i,j} &\sim \text{Normal}\ (\mu_{i, j},\ \sigma^2) \\
\mu_{i, j} &= \beta_1 + \beta_2 A_i + \beta_3 X_j A_i \\
\end{align*}
$$

テキストより引用

目的変数$${Y_{i,j}}$$は測定$${i}$$回目の調査県$${j}$$の平均身長です。
説明変数は$${A}$$と$${X}$$です。

$${A_i}$$は測定回・測定年0年目か1年目（i=0,1）、$${X_j}$$は給食タイプの変更の有無（0:変更なし、1:変更あり）を要素にしています。

$${\beta_1, \cdots, \beta_3}$$は回帰モデルの係数です。
$${\beta_1}$$は切片であり、１回目の測定の平均値を示します。

$${\beta_2, \beta_3}$$は説明変数にかかっていて、$${A_i=0}$$（1回目の測定）や$${X_j=0}$$（給食タイプ変更なし）の場合は、係数を含め０になります。
したがいまして、$${\beta_2}$$は１回目から２回目の身長の変化量、$${\beta_3}$$は給食タイプ変更の効果量といえます。
測定回・給食タイプと係数の関わりを示すテキスト掲載の数式をお借りします。

１回目の測定　　　　　　　：$${\mu_{0,j} = \beta_1}$$
２回目・給食タイプ変更なし：$${\mu_{1, 0} = \beta_1 + \beta_2}$$
２回目・給食タイプ変更あり：$${\mu_{1, 1} = \beta_1 + \beta_2 + \beta_3}$$

テキストより引用

給食タイプの変更が身長の伸びに影響するかどうかは、係数$${\beta_3}$$を確認すればいいのです。

6. 一般化線形モデル

上記のモデルをテキストと同じように一般化線形モデルに当てはめてみます。
Pythonのパッケージ statsmodels を利用します。
テキストの glm() と同じように formula で線形回帰モデルを指定しました。

### 一般化線形モデル
result = smf.glm(formula='mean_Y ~ 1 + Age + Age:X', data=data).fit()
display(result.summary())

【実行結果】
赤枠は私が付けました。

$${\beta_1}$$はIntercept、$${\beta_2}$$はAge、$${\beta_3}$$はAge:Xの箇所が該当します。
３つの係数の $${p}$$値（P>|z| ）は0.05未満なので、統計的に有意であるように見えます。

給食タイプの変更有無を比較します。
２回目の給食タイプ変更なしの身長への影響は$${\beta_2}$$ですので、$${5.65}$$cmです。
２回目の給食タイプ変更ありの身長への影響は$${\beta_2}$$と$${\beta_3}$$の係数coefの合計ですので、$${5.65-1.72=3.93}$$cmです。
給食タイプを変更することで$${\boldsymbol{1.72}}$$cm（$${\boldsymbol{\beta_3}}$$相当）伸びが悪いように見えてしまいます。

「給食タイプの変更は身長の増減に影響がない」はずのデータですが、直線的あてはめを採用すると、データの成り立ちと異なる結果が出力されてしまうのです。

線形回帰モデルの出力結果に関する原因はテキストの２４～２５ページにかけて久保先生が丁寧に解説されています。
要約すると「同じ児童の身長を測定しているので対応のある分析を行うべきだが、この線形回帰モデルは対応のない分析である」からだそう。

テキストはこの「不備のある」モデル１からベイズモデルを開始します。
たいへん簡単なモデルなので「ベイズモデルとは何か？」という説明が簡単になり、階層モデルへのステップにするようです。

では、モデル１のPyMCモデリングに進みます。

ベイズモデル

1. モデル数式

線形回帰モデルの数式を参照しつつ、テキストの事前分布と同じ結果になるように数式化しました。

$$
\begin{align*}
\beta_1,\ \beta_2,\ \beta_3 &\sim \text{Uniform}\ (\text{lower}=-10000,\ \text{upper}=10000) \\
\sigma &\sim \text{Uniform}| (\text{lower}=0,\ \text{upper}=10000) \\
\mu &= \beta_1 + (\beta_2 + \beta_3 \times X) \times Age \\
y &\sim \text{Normal}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)
\end{align*}
$$

2. モデルの定義

数式モデルどおり実直にPyMCのモデルを定義しましょう。

### モデルの定義
with pm.Model() as model1:
    
    ### データ関連定義
    ## coordの定義
    model1.add_coord('data', values=data.index, mutable=True)
    model1.add_coord('params', values=range(1, 4), mutable=True)
    
    ## dataの定義
    # 目的変数：身長の平均値Y
    Y = pm.ConstantData('Y', value=data['mean_Y'].values, dims='data')
    # 説明変数：給食タイプx 0:甲（変更なし）, 1:乙（変更あり）
    X = pm.ConstantData('X', value=data['X'].values, dims='data')
    # 説明変数：測定回Age 0:１回目, 1:２回目
    Age = pm.ConstantData('Age', value=data['Age'].values, dims='data')

    ### 事前分布
    # 偏回帰係数β
    beta = pm.Uniform('beta', lower=-10000, upper=10000, dims='params')
    # 尤度関数の標準偏差
    sigma = pm.Uniform('sigma', lower=0, upper=10000)

    ### 線形予測子
    mu = pm.Deterministic('mu', beta[0] + (beta[1] + beta[2] * X) * Age,
                          dims='data')
    
    ### 尤度
    y = pm.Normal('y', mu=mu, sigma=sigma, observed=Y, dims='data')

モデルの内容を表示・可視化してみましょう。

### モデルの表示
model1

【実行結果】

続いてモデルを可視化します。

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

【実行結果】
データやパラメータの繋がりが示されます。
丸い形状の変数には「確率分布」が併記されます。
データ値を示す変数はグレー色が付いています。
説明変数$${X,\ Age}$$と３つ係数$${beta}$$が$${mu}$$の計算に繋がり、$${mu}$$を平均パラメータ・$${sigma}$$を標準偏差パラメータにする正規分布に２０個の観測データ$${y}$$が従う様子が分かります。

3. MCMCの実行と収束の確認

■ MCMC
いよいよベイズプログラミングの真骨頂 MCMC の実行です。
マルコフ連鎖=chainsを4本、バーンイン期間=tuneを1000、利用するサンプル=drawを１chainあたり1000に設定して、合計4000個の事後分布からのサンプリングデータを生成します。テキストの指定数とは異なっています。

サンプリングデータは idata1 に格納されます。

サンプル方法に numpyro を指定しています。処理速度が早くなります。
numpyro を使わない場合は「nuts_sampler='numpyro',」を削除します。

### 事後分布からのサンプリング 20秒
with model1:
    idata1 = pm.sample(draws=1000, tune=1000, chains=4, random_seed=123,
                       nuts_sampler='numpyro')

【実行結果】
処理時間は20秒です。

■ $${\hat{R}}$$で収束確認
収束の確認をします。
指標$${\hat{R}}$$（あーるはっと）の値が$${1.1}$$以下のとき、収束したこととします。
次のコードでは$${\hat{R}>1.01}$$のパラメータの個数をカウントします。
個数が０ならば、$${\hat{R} \leq1.1}$$なので収束したとみなします。

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

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

【実行結果】
$${\hat{R}>1.01}$$のパラメータは０個でした。
$${\hat{R} \leq1.1}$$なので収束したとみなします。

■ $${\hat{R}}$$の値把握と事後統計量の表示

### 事後統計量の表示
pm.summary(idata1, hdi_prob=0.95)

【実行結果】
MCMCでサンプリングデータを生成したすべてのパラメータについて、統計量や診断情報が一覧表示されました。

係数$${beta}$$の平均値なども確認できます。
注目の$${beta_3}$$は平均値$${-1.737}$$、95%HDI$${[-3.572, -0.037]}$$です。
給食タイプの変更が身長の伸びに悪化方向の影響を与えているように見えてしまいます。。。

■ トレースプロットで収束確認
トレースプロットを描画します。
テキストの図5に該当します。
こちらの図でも収束の確認を行えます。

### トレースプロットの描画
pm.plot_trace(idata1, compact=False)
plt.tight_layout();

【実行結果】
左の曲線は、４本のマルコフ連鎖の分布です。
４本の曲線がほぼ同じ分布を描いているので、サンプリングデータが同一の分布から生成されたことが推察できます。
右のノイズみたいな図は、均等に乱雑に描画されています。
この状態は収束していると推察できます。
均等・乱雑ではなくて、何らかの傾向を示すような描画の場合は収束を疑います。

$${\hat{R}}$$とトレースプロットから、分布の収束が確認できましたので、MCMCのサンプリングデータを用いて、推論を継続できます！
やったね！

4. 分析

$${\beta_3}$$の事後分布を描画して、直線的あてはめの悪さを再確認しましょう。
テキストの図6に相当します。

### β3の事後分布プロットを描画 ※図6に相当
pm.plot_posterior(idata1, hdi_prob=0.95, var_names=['beta'], round_to=3,
                  coords={'params': [3]}, ref_val=0, figsize=(5, 3))
plt.tight_layout();

【実行結果】
係数の事後平均値は$${-1.74}$$。給食タイプの変更が身長の伸びに悪化方向の影響を与えるという、誤った推論をしてしまっています。

また、オレンジの文字は、$${\beta_3}$$が0超の確率、0未満の確率を示しています。
$${\beta_3>0}$$となる確率は$${2.6\%}$$です。ほぼテキストと同じです。

$${\beta_3>0}$$の計算は、MCMCサンプリングデータを用いて次のように実施することもできます。

### beta3>0の確率
(idata1.posterior.beta[:, :, 2] > 0).sum().data / 4000

【実行結果】

この先、よりよいモデルを目指して階層モデルを構築するのは、次回のお楽しみにいたしましょう。

5. 推論データの保存

pickle で MCMCサンプリングデータ idata1 を保存できます。

### 推論データの保存
file = r'idata1_ch2.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata1, f)

次のコードで保存したファイルを読み込みできます。

### 推論データの読み込み
file = r'idata1_ch2.pkl'
with open(file, 'rb') as f:
    idata1_ch2_load = pickle.load(f)

ベイズ記事は以上です。

次回予告

「階層ベイズ最初の一歩」章
　モデル２「階層ベイズモデル」

結び

久保先生の章を実践していて気づいたことがあります。
ベイズモデリングに物語性を持たせて学習の進行をはかっていらっしゃるのです。
読者の想像力を借りつつ、物語・モデル側に引き込みつつ、分析シナリオを展開する、こんなストーリーテリングに憧れます。
自分のブログもいつか、そうなれたらいいです（無理）。

おわり

シリーズの記事

次の記事

前の記事

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

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

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

#仕事について話そう

115,101件

#今月の振り返り

24,610件

この記事が気に入ったらサポートをしてみませんか？