Up

新型コロナウィルスワクチンの効果

２項分布によるベイズ分析の試み

Microsoft Edgeで表示したときに不具合があれば、Edgeのメニュ「履歴｜閲覧データをクリア」を選んで表示される「閲覧データをクリア」フォームで「今すぐクリア」ボタンを押してクリアしてからEdgeで再表示して下さい。

PCR検査結果の統計学的信頼性と感染者比率の関係を扱ったウェブサイトも用意しました。

新型コロナウィルスのワクチンの効果について、池上彰（週刊文春、2020.12.17、pp.46−47）は

ファイザーによると、対象者41,135人をワクチン群と偽薬群に分けて効果を調べたとき、偽薬群の発症者は162人で、ワクチン群は8人であった。

と説明している。これらの数値の統計学的分析をベイズ法によって試みる。ベイズ統計分析によって分かり易い分析結果が得られると考えた。

いま、発症したか否かのデータであるので２項分布を用いて効果のベイズ分析を行ってみる。２項分布については、岡本安晴（2019）、pp．169−170に簡単な説明がある。

発症確率をとし、N人中、発症者がｋ人である確率をP(k,N,p)で表せば、この確率が２項分布に従うとすると次式が成り立つ。

２群の各人数は示されていないが、仮に半分ずつに近い偽薬群が20,568人、ワクチン群が20,567人であったとする。総数が奇数で等分に分けられないので、偽薬群の発症者数の方が多いことに合わせて偽薬群の総数を１人多い分け方にした。

偽薬群の発症確率をp0、ワクチン群の発症確率をp1としたとき、ワクチンの効果はp1の方がp0よりどれだけ小さくなったかによって表される。この効果を次式で与えられるfで表す。すなわち、

である。P1/p0はワクチン群の発症率の偽薬群の発症率に対する比であり、１から引くことにより発症率の減少率を表している。(p0-p1)/p0は、発症率の減少量のp0に対する比であるので減少率を表す。

上式より、

を得る。発症率p0は、１から効果量f/100を減じた比率分(1-f/100)に減少してp1となる。

偽薬群の人数をN0、ワクチン群の人数をN1、偽薬群での発症者数をk0、ワクチン群での発症者数をk1と置く。偽薬群での発症確率をp0、ワクチン群での発症確率をp1と置き、２項分布によるベイズ分析を行うスクリプトを以下のように書いてみた。PyMC3によるベイズ分析は、このウェブサイトの後半に示した。

リスト１　Stanによるスクリプト例

import pystan

import matplotlib.pyplot as plt

import seaborn as sb

import numpy as np

ModelScript = """

data {

int N0;

int N1;

int k0;

int k1;

}

parameters {

real<lower = 0.0, upper = 1.0> p0;

real<lower = 0.0, upper = 100.0> f;

}

transformed parameters {

real<lower = 0.0, upper = 1.0> p1;

p1 = p0 * (1.0 - f / 100.0);

}

model {

k0 ~ binomial(N0, p0);

k1 ~ binomial(N1, p1);

}

"""

Data = {'N0': 20568, 'N1': 20567,

'k0': 162, 'k1': 8}

sm = pystan.StanModel(model_code = ModelScript)

fit = sm.sampling(data = Data, n_jobs = 1)

print(fit)

plt.rcParams['font.family'] = ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']

p0 = np.percentile(fit['p0'], 50.0)

p1 = np.percentile(fit['p1'], 50.0)

sb.kdeplot(fit['p0'], label = '偽 薬')

sb.kdeplot(fit['p1'], label = 'ワクチン')

plt.title('ワクチン: {0:.2e} 偽 薬: {1:.2e}'.format(p1, p0), fontsize = 16)

plt.yticks([])

plt.xlabel('確 率', fontsize = 16)

plt.legend(fontsize = 14)

plt.tight_layout()

f_2p5, f_50, f_97p5 = np.percentile(fit['f'], [2.5, 50.0, 97.5])

plt.figure()

sb.kdeplot(fit['f'])

plt.title('中央値 = {0:.1f}\n95%CI = [{1:.1f}, {2:.1f}]'.format(f_50, f_2p5, f_97p5),

fontsize = 16)

plt.yticks([])

plt.xlabel('効 果', fontsize = 14)

plt.tight_layout()

plt.show()

上のスクリプトを実行すると下図のグラフが表示される。

図１

2枚のグラフが重なって表示された場合は、ドラッグして移動すればよい。ベイズ分析のStanスクリプトの実行は、Python用のインターフェースであるPyStanにより行っているが、これらのことは＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞を参照されたい。

図１の左側のグラフが、それぞれの群での発症確率の事後分布を示すものである。偽薬群より左側にワクチン群の分布が位置していて、２つの分布は明瞭に分かれている。すなわち、ワクチン群の発症率は、偽薬群の発症率より低いことが明確に分かる。

右側のグラフは、効果の事後分布である。中央値が約95％であり、95％確信区間は「90.0、97.4」となっている。90％以上の効果は期待できそうである。効果fが90％であれば、発症率は、すなわちワクチン群の発症率は偽薬群の発症率の1/10に減少していることになる。効果の95％確信区間が「90.0、97.4」であることから、効果は90％以上であることがほぼ確実であると考えられるので、ワクチン群の発症率は偽薬群の発症率の1／10以下であると言える。

なお、Jupyter Notebook用のファイルSampleScript.ipynbも用意した。ファイル名SampleScript.ipynbをマウスの右ボタンでクリックしてダウンロードしたファイルを、Jupyter Notebookから開けば実行できる。

PyMC3による分析

下記ウェブサイト

https://pymc-devs.medium.com/the-future-of-pymc3-or-theano-is-dead-long-live-theano-d8005f8a0e9b

https://pypi.org/project/pymc3/

によれば、PyMC4の開発は中止となり、PyMC3の開月が続けられることになったようである（2021.02）。

PyMC3に将来性があると考え、PyMC3によるスクリプト例も用意してみた。

スクリプト例をリスト２に示す。

リスト２　PyMC3によるスクリプト例

import pymc3 as pm

import numpy as np

import theano.tensor as tt

import matplotlib.pyplot as plt

import seaborn as sb

idx = np.array([0, 1], dtype = int)

Data = np.array([162, 8]) # 発症者数

tn = np.array([20568, 20567]) # 各群の総数

with pm.Model() as c_b:

p0 = pm.Uniform('p0', lower = 0.0, upper = 1.0) # 偽薬群の発症確率

f = pm.Uniform('f', lower = 0.0, upper = 1.0) # 効果

p1 = pm.model.Deterministic('p1', p0 * (1 - f)) # ワクチン群の発症確率

ps = tt.stack(p0, p1)

y = pm.Binomial('y', n = tn[idx], p = ps[idx], observed = Data)

trace = pm.sample(3000)

plt.rcParams['font.family'] = ['Yu Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']

plt.figure(figsize = (10, 6))

sb.kdeplot(trace['p0'], label = '偽薬')

sb.kdeplot(trace['p1'], label = 'ワクチン')

plt.title('ワクチン: {0:.2e} 偽薬: {1:.2e}'.format(

trace['p1'].mean(), trace['p0'].mean()),

fontsize = 24)

plt.xlim([0.0, 0.01])

plt.xticks(fontsize = 20)

plt.xlabel('発症確率', fontsize = 20)

plt.yticks([])

plt.ylabel('')

plt.legend(fontsize = 20)

plt.savefig('Fig1.png')

plt.show()

t_f = [v*100 for v in trace['f']] # ％に変換

f025p, fmed, f975p = np.percentile(t_f, [2.5, 50, 97.5])

plt.figure(figsize = (10, 6))

plt.title('中央値 = {0:.1f}\n95% CI = [{1:.1f}, {2:.1f}]'.format(

fmed, f025p, f975p), fontsize = 20)

sb.kdeplot(t_f)

plt.xlim([85, 100])

plt.xticks(fontsize = 20)

plt.xlabel('効果 f (%)', fontsize = 20)

plt.ylabel('')

plt.yticks([])

plt.savefig('Fig2.png')

plt.show()

リスト２をWindowsで実行するときは、Jupyter Notebook上で実行するとよい。

リスト２を実行すると、発症確率の事後分布として図2.1のグラフを得る。

図2.1

効果の事後分布は、図2.2のようにグラフになる。

図2.2

サンプリングによるランダムな変動を考慮すれば、PyMC3による事後分布の推定とStanによる事後分布の推定が一致していることが確認できたと言える。

PyMC3についても、拙著＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で説明している。

Up