スタチンの効果の多項分布によるベイズ分析
Spiegelhalter (2019)のTable 4.1には、スタチン（statin）を処方された人と偽薬（placebo）を処方された人の心臓発作（heart attack）、脳卒中（Stroke）などの発症率が記載されている。発症率から発症数を算出して表１にまとめた。
|
表１ Frequencies calculated from Table 4.1 of Spiegelhalter (2019) |
||
|
Event |
Placebo |
Statin |
|
Heart attack |
1,212 |
893 |
|
Stroke |
585 |
442 |
|
Death from any cause |
1,509 |
1,325 |
|
Others |
6,961 |
7,609 |
発症数に対して多項分布（岡本、2009）を適用してベイズ分布を行ってみた。
Stanスクリプトは、リスト１のように作成した。
Placebo群とstatin群における心臓発作などの発症頻度をデータとしている。
array[K] int placebo; // Placebo群のデータ
array[K] int statin; // Statin群のデータ
上記、配列において、各群
［Heart attack、Stroke、Death from any cause、Others］
の順である。これは、他の配列（多項分布のパラメータなど）においても同様である。
発症頻度が従う多項分布のパラメータをsimplex型の変数として、以下のように設定している。
simplex[K] theta_p; // Placebo群
simplex[K] theta_s; // Statin群
条件間のパラメータ値の差、および比をgenerated quantitiesにおいて以下のように算出している。
generated quantities {
vector[K] diff_p;
vector[K] ratio_p;
diff_p = theta_s - theta_p;
for (k in 1:K)
ratio_p[k] = theta_s[k] / theta_p[k];
}
リスト１のStanスクリプトを用いてベイズ分析を行うPythonスクリプトをリスト2のように作成した。
表１のデータを次のコードで設定して、リスト１のStanスクリプトを実行している。
Data = {'K': 4,
'placebo': [1212, 585, 1509, 6961],
'statin': [ 893, 442, 1325,
7609]}
リスト２を実行すると、MCMCサンプリング後、多項分布のパラメータの事後分布のグラフが表示される（図１）。
図１
タイトルに記されているplacebo:およびstain:の右側の数値は、それぞれの群の人数である。ほぼ、総計２万人のデータである。
図中のレジェンド中の（P）および（S）は、placebo群およびstatin群を表す。３つの発症確率において、statin群の事後分布がplacebo群より左側にあり、群間の重なりはほとんど認められない。すなわちstatin群の発症確率がplacebo群より小さいことが統計学的に明瞭に示されていると言える。
図１のグラフのWindowを閉じると図２のグラフが表示される。
図２
２つの群のパラメータの差
（statin群のパラメータ） ― （placebo群のパラメータ）
の事後分布のグラフである。発症確率の差の事後分布は、いずれも負の領域にあり、差が負であることが統計学的に明確に示されていると言える。すなわち、statin群の発症確率がplacebo群より小さいということである。
図２のWindowを閉じると、図３のグラフが表示される。
図３
パラメータの比
（Statin群のパラメータ値）／（Placebo群のパラメータ値）
の事後分布である。発症確率のstatin群のplacebo群に対する比の事後分布は統計学的に１より小さいことが明瞭に示されている。Heart attackおよびStrokeの発症率において、statin群のplacebo群に対する比の事後分布の中央値は約3／4である。
上の分析データは、サンプルサイズが約２万人である。これを、比率はそのままにして、10分の1のサイズに変えてみる。すなわち、リスト２における次のコード
Data = {'K': 4,
'placebo': [1212, 585, 1509, 6961],
'statin': [ 893, 442, 1325,
7609]}
を
Data = {'K': 4,
'placebo': [121, 59, 151, 696],
'statin': [ 89, 44, 133,
761]}
に変更して実行してみる。このとき、statin群とplacebo群を合わせたサンプルサイズは約２千人である。
まず、各パラメータの事後分布であるが、図４のようになる。
図４
各発症確率は線のスタイルによって区別されているが、それぞれのstatin群（緑）とplacebo群（赤）の事後分布の分離は図１の分布に比べると重なりがかなり認められる。
各発症パラメータの差の事後分布は、図５のようになる。
図５
図２のグラフと比べると、０を跨いでいるのが認められ、特にstrokeの発症確率の事後分布の0からの分離が明瞭ではない。
図３のWindowを閉じると、図６のグラフが表示される。
図６
比の事後分布も、サンプルサイズを小さくした場合は、事後分布と１との関係が図３と比べると明瞭ではない。
リスト１ 多項分布のStanスクリプト（statin_model.stan）
data {
int K;
array[K] int placebo; // Placebo群のデータ
array[K] int statin; // Statin群のデータ
}
transformed data {
//
// Dirichlet分布のパラメータ
//
一様分布
//
vector[K] alpha_p;
vector[K] alpha_s;
for (i in 1:K) {
alpha_p[i] = 1;
alpha_s[i] = 1;
}
}
parameters {
//
// 多項分布のパラメータ
//
simplex[K] theta_p; // Placebo群
simplex[K] theta_s; // Statin群
}
model {
theta_p ~
dirichlet(alpha_p);
theta_s ~
dirichlet(alpha_s);
placebo ~
multinomial(theta_p);
statin ~
multinomial(theta_s);
}
generated quantities {
vector[K] diff_p;
vector[K] ratio_p;
diff_p = theta_s - theta_p;
for (k in 1:K)
ratio_p[k] = theta_s[k] / theta_p[k];
}
リスト２ 多項分布によるベイズ分析用PyStanスクリプト（statin_model.stan リスト１）を用いるPythonスクリプト
import stan
import pandas as pd
import arviz as az
import seaborn as sb
import matplotlib.pyplot as plt
import numpy as np
Data = {'K': 4,
'placebo': [1212, 585, 1509, 6961],
'statin': [ 893, 442, 1325,
7609]}
n_placebo = np.sum(Data['placebo'])
n_statin = np.sum(Data['statin'])
print('n_placebo =', n_placebo)
print('n_statin =', n_statin)
with open('statin_model.stan') as
f_stan:
sm =
stan.build(f_stan.read(), data=Data)
fit = sm.sample()
i_data = az.from_pystan(posterior=fit,
posterior_model=sm)
print(az.summary(i_data))
d_fm = fit.to_frame()
#
#
多項分布のパラメータの事後分布
#
name_lst = ['Heart attack', 'Stroke',
'Any cause', 'Others']
line_lst = ['--', '-.', ':']
plt.title('Posterior distribution' +
f'\nplacebo:{n_placebo},
statin:{n_statin}')
for i in range(3):
med =
np.median(d_fm[f'theta_p.{i+1}'])
sb.kdeplot(d_fm[f'theta_p.{i+1}'], c='r', ls=line_lst[i],
label='(P):'+ name_lst[i] + f'\nmed.={med:.3f}')
med =
np.median(d_fm[f'theta_s.{i+1}'])
sb.kdeplot(d_fm[f'theta_s.{i+1}'],
c='g', ls=line_lst[i],
label='(S):' + name_lst[i] + f'\nmed.={med:.3f}')
plt.xlabel('Probability')
plt.legend(ncols=2)
plt.tight_layout()
plt.show()
#
#
パラメータの差（Statin
- Placebo）の事後分布
#
plt.title('Posterior distribution' +
f'\nplacebo:{n_placebo},
statin:{n_statin}')
for i in range(4):
med =
np.median(d_fm[f'diff_p.{i+1}'])
sb.kdeplot(d_fm[f'diff_p.{i+1}'], label=name_lst[i] +
f'\nmed.={med:.3f}')
plt.xlabel('Difference of
probabilities: Statin - Placebo')
plt.legend()
plt.show()
#
#
パラメータの比（Statin
/ Placebo）の事後分布
#
plt.title('Posterior distribution' +
f'\nplacebo:{n_placebo},
statin:{n_statin}')
for i in range(4):
med = np.median(d_fm[f'ratio_p.{i+1}'])
sb.kdeplot(d_fm[f'ratio_p.{i+1}'], label=name_lst[i] +
f'\nmed.={med:.3f}')
plt.xlabel('Ratio of probabilities:
Statin / Placebo')
plt.legend()
plt.show()
岡本安晴（2009） 統計学を学ぶための数学入門［下］．培風館
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
Spiegelhalter, D. (2019). The art of statistics: Learning from data. Pelican