PCR検査結果の統計学的解釈
岡本安晴
注意： Mcrosoft Edgeによる表示が乱れた場合は、履歴のキャッシュをクリアしてから、再表示して下さい。
新型コロナウィルスワクチンの効果のベイズ分析を試みたウェブサイトも用意しました。
新型コロナウィルスのPCR検査で「陰性です（検出されませんでした）」という結果が出たときに、実は感染しているという確率はどれくらいであるのか、鼻咽頭ぬぐい液と唾液の場合について、感染者の比率の関数として計算してみた。「検査結果が陰性である」ときに「実は感染している」という確率
は、検査対象母集団における感染者の比率
によって大きく変わる。検査対象母集団における感染者の比率（罹患率）については、本ウェブサイトの後半で扱っている。
検査結果が「陰性」であるときの感染確率は条件付確率として次式（岡本安晴「心理学データ分析と測定」勁草書房、第12.3節「診断検査」参照）で計算できる。
確率
は感度（sensitivity）と呼ばれ、確率
は特異度（specificity）と呼ばれている。PCR検査におけるこれらの値を、Yokotaら（2020）に報告されているものとした。
感染率の関数として確率のグラフを鼻咽頭ぬぐい液の場合について描くと図１のようになる。プログラムは、リスト１に示したものである。
図１ 鼻咽頭ぬぐい液
感染率が上がると、確率も上がることが分かる。すなわち、「陰性」の結果が出ても、感染している確率が高くなる。感染率が100％で全員が感染しているときは、検査結果が陰性の人（確率は、１―感度）も全員感染している。感染率が0％のときは、検査結果が陰性の人（確率は、特異度）も全員感染していない。
感染率の具体例として以下の場合を考えてみた。ただし、感染率は大まかな推定である。
ダイアモンドプリンセス号の場合
国立感染症研究所の「現場からの概況：ダイアモンドプリンセス号におけるCOVID-19症例」（2020年2月19日掲載）によれば、乗船者3711人に対して感染確定者（検査陽性者）531人となっている。感染率を
とおいて、確率を求めると
となる。検査で「陰性」と出ても、50人に1人は感染者であるということになる。
東京都の場合
感染者数は、
https://news.yahoo.co.jp/pages/article/covid19tokyo
に掲載されている2021年3月4日のJX通信社/FASTALERT提供のデータでは、2939人である。東京都の人口は、東京都のHPによれば令和3年現在13960236人である。感染率を
とおいて、確率を求めると
となる。検査で「陰性」と出た場合、33333人に1人は感染者であるということになる。
感染率の関数として確率のグラフを唾液（saliva）の場合について描くと図２のようになる。プログラムは、リスト２に示したものである。
図２ 唾液
感染率が上がると、確率も上がることが分かる。すなわち、「陰性」の結果が出ても、感染している確率が高くなる。
感染率の具体例として以下の場合を考えてみた。ただし、感染率は大まかな推定である。
ダイアモンドプリンセス号の場合
国立感染症研究所の「現場からの概況：ダイアモンドプリンセス号におけるCOVID-19症例」（2020年2月19日掲載）によれば、乗船者3711人に対して感染確定者（検査陽性者）531人となっている。感染率を
とおいて、確率を求めると
となる。検査で「陰性」と出ても、77人に1人は感染者であるということになる。
東京都の場合
感染者数は、
https://news.yahoo.co.jp/pages/article/covid19tokyo
に掲載されている2021年3月4日のJX通信社/FASTALERT提供のデータでは、2939人である。東京都の人口は、東京都のHPによれば令和3年現在13960236人である。
を感染率を
とおいて、確率を求めると
となる。検査で「陰性」と出た場合、62500人に1人は感染者であるということになる。
鼻咽頭ぬぐい液（NPS）と唾液（saliva）の比較を容易にするため、図１と図２のグラフを同じグラフ面上に描いてみた。リスト１とリスト２のプログラムを繋いだものをまとめてリスト３として実行すると図３のグラフが得られる。
図３ 鼻咽頭ぬぐい液と唾液の比較
唾液の方が感染率P(infected)全体に亘って確率P(infected|neg.)の値が小さいことが分かる。
感染率の低い範囲０〜0.1についてグラフをリスト４のプログラムで描くと図４を得る。
図４ 感染率０〜1％における比較
唾液の方が確率P(infected|Neg.)の値が小さい。「陰性」の結果が得られたときは、「唾液」による方法の方が安心である。
Yokota, I., Shane, P. Y., Okada, K, Uroki, Y., Yang, Y., Inao, T., Sakamaki, K., Iwasaki, S., Hayasaka, K., Sugita, J., Nishida, M., Fujisawa, S., & Teshima, T. Mass screening of asymptomatic persons for SARS-Cov-2 using saliva. Clinical Infectious Diseases, 2020.9.25, doi: 10.1093/cid/claa1388.
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
リスト１ 鼻咽頭ぬぐい液（nasopharyngeal swabs: NPS）PCR検査のP(infected|negative)の計算
import numpy as np
import matplotlib.pyplot as plt
class Infected_Neg:
"""
For calculation of P(infected|negative)
"""
def __init__(self,
sensitivity, specificity):
self.Pr_pos_infctd = sensitivity
# P(positive|infected)
self.Pr_neg_noninfctd = specificity
#
P(negative|non-infected)
self.Pr_neg_infctd = 1.0 - self.Pr_pos_infctd # P(negative|infected)
self.Pr_pos_noninfctd = 1.0 - self.Pr_neg_noninfctd # P(positive|non-infected)
def value(self, Pr_infctd):
"""
P(infected|negative) for P(infected)
"""
return self.Pr_neg_infctd * Pr_infctd / \
(self.Pr_neg_infctd * Pr_infctd + self.Pr_neg_noninfctd * (1.0 -
Pr_infctd))
Pr_infctd = np.linspace(0, 1, 1000)
#
#
P(infected|negative) as the function of P(infected)
#
for nasopharyngeal swabs: NPS）PCR test
#
Pr_infctd_neg =
Infected_Neg(sensitivity = 0.86, specificity = 0.9993).value
Pr = [Pr_infctd_neg(v) for v in
Pr_infctd]
plt.plot(Pr_infctd, Pr)
plt.xlabel('P(infected)', fontsize =
14)
plt.ylabel('P(infected|neg.)',
fontsize = 14)
plt.title('Nasopharyngeal swabs (NPS)
PCR test', fontsize = 16)
plt.savefig('naso.png')
plt.show()
Pr_infctd = 0.0002
print(f'{Pr_infctd_neg(Pr_infctd):.7f}
for Pr(infected) = {Pr_infctd}')
Pr_infctd = 0.143
print(f'{Pr_infctd_neg(Pr_infctd):.7f}
for Pr(infected) = {Pr_infctd}')
リスト２ 唾液（saliva）PCR検査のP(infected|negative)の計算
import numpy as np
import matplotlib.pyplot as plt
class Infected_Neg:
"""
For calculation of P(infected|negative)
"""
def __init__(self,
sensitivity, specificity):
self.Pr_pos_infctd = sensitivity
# P(positive|infected)
self.Pr_neg_noninfctd = specificity
#
P(negative|non-infected)
self.Pr_neg_infctd = 1.0 - self.Pr_pos_infctd # P(negative|infected)
self.Pr_pos_noninfctd = 1.0 - self.Pr_neg_noninfctd # P(positive|non-infected)
def value(self, Pr_infctd):
"""
P(infected|negative) for P(infected)
"""
return self.Pr_neg_infctd * Pr_infctd / \
(self.Pr_neg_infctd * Pr_infctd + self.Pr_neg_noninfctd * (1.0 -
Pr_infctd))
Pr_infctd = np.linspace(0, 1, 1000)
#
#
P(infected|negative) as the function of P(infected)
#
for saliva PCR test
#
Pr_infctd_neg =
Infected_Neg(sensitivity = 0.92, specificity = 0.9996).value
Pr_infctd = np.linspace(0, 1, 1000)
Pr = [Pr_infctd_neg(v) for v in
Pr_infctd]
plt.plot(Pr_infctd, Pr)
plt.title('Saliva PCR test', fontsize
= 16)
plt.xlabel('P(infected)', fontsize =
14)
plt.ylabel('P(infected|neg.)',
fontsize = 14)
plt.savefig('saliva.png')
plt.show()
Pr_infctd = 0.0002
print(f'{Pr_infctd_neg(Pr_infctd):.7f}
for Pr(infected) = {Pr_infctd}')
Pr_infctd = 0.143
print(f'{Pr_infctd_neg(Pr_infctd):.7f}
for Pr(infected) = {Pr_infctd}')
リスト３ NPSとsalivaのPCR検査のP(infected|negative)のグラフ
import numpy as np
import matplotlib.pyplot as plt
class Infected_Neg:
"""
For calculation of P(infected|negative)
"""
def __init__(self,
sensitivity, specificity):
self.Pr_pos_infctd = sensitivity
# P(positive|infected)
self.Pr_neg_noninfctd = specificity
#
P(negative|non-infected)
self.Pr_neg_infctd = 1.0 - self.Pr_pos_infctd # P(negative|infected)
self.Pr_pos_noninfctd = 1.0 - self.Pr_neg_noninfctd # P(positive|non-infected)
def value(self, Pr_infctd):
"""
P(infected|negative) for P(infected)
"""
return self.Pr_neg_infctd * Pr_infctd / \
(self.Pr_neg_infctd * Pr_infctd + self.Pr_neg_noninfctd * (1.0 -
Pr_infctd))
Pr_infctd = np.linspace(0, 1, 1000)
#
#
P(infected|negative) as the function of P(infected)
#
for nasopharyngeal swabs: NPS）PCR test
#
Pr_infctd_neg_nps =
Infected_Neg(sensitivity = 0.86, specificity = 0.9993).value
Pr = [Pr_infctd_neg_nps(v) for v in
Pr_infctd]
plt.plot(Pr_infctd, Pr, label = 'NPS')
#
#
P(infected|negative) as the function of P(infected)
#
for saliva PCR test
#
Pr_infctd_neg_saliva =
Infected_Neg(sensitivity = 0.92, specificity = 0.9996).value
Pr = [Pr_infctd_neg_saliva(v) for v in
Pr_infctd]
plt.plot(Pr_infctd, Pr, label =
'Saliva')
plt.xlabel('P(infected)', fontsize =
14)
plt.ylabel('P(infected|neg.)',
fontsize = 14)
plt.legend(fontsize = 14)
plt.savefig('Both.png')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
class Infected_Neg:
"""
For calculation of P(infected|negative)
"""
def __init__(self,
sensitivity, specificity):
self.Pr_pos_infctd =
sensitivity
# P(positive|infected)
self.Pr_neg_noninfctd = specificity
#
P(negative|non-infected)
self.Pr_neg_infctd = 1.0 - self.Pr_pos_infctd # P(negative|infected)
self.Pr_pos_noninfctd = 1.0 - self.Pr_neg_noninfctd # P(positive|non-infected)
def value(self, Pr_infctd):
"""
P(infected|negative) for P(infected)
"""
return self.Pr_neg_infctd * Pr_infctd / \
(self.Pr_neg_infctd *
Pr_infctd + self.Pr_neg_noninfctd * (1.0 - Pr_infctd))
Pr_infctd = np.linspace(0, 0.01, 1000)
#
#
P(infected|negative) as the function of P(infected)
#
for nasopharyngeal swabs: NPS）PCR test
#
Pr_infctd_neg_nps =
Infected_Neg(sensitivity = 0.86, specificity = 0.9993).value
Pr = [Pr_infctd_neg_nps(v) for v in
Pr_infctd]
plt.plot(Pr_infctd, Pr, label = 'NPS')
#
#
P(infected|negative) as the function of P(infected)
#
for saliva PCR test
#
Pr_infctd_neg_saliva =
Infected_Neg(sensitivity = 0.92, specificity = 0.9996).value
Pr = [Pr_infctd_neg_saliva(v) for v in
Pr_infctd]
plt.plot(Pr_infctd, Pr, label =
'Saliva')
plt.xlabel('P(infected)', fontsize =
14)
plt.ylabel('P(infected|neg.)',
fontsize = 14)
plt.legend(fontsize = 14)
plt.tight_layout()
plt.savefig('Both_small.png')
plt.show()
まず、代数的計算による方法を説明する。その後、それに基づくベイズ分析例を示す。
検査母集団における罹患率と検査結果の陽性率の関係は次式で表される。
いま、
とおけば、
すなわち、sensitivityとspecificityが判れば、陽性率
から罹患率
を求めることができる。
次の関数は、陽性率から罹患率を求めるものである。
def EstPopInfctn( sensitivity,
specificity, prop_positive):
"""
検査母数団の罹患率の計算
"""
A = np.array([[sensitivity,
1 - specificity], [1 - sensitivity, specificity]])
inv_A = np.linalg.inv(A)
Results =
np.array([prop_positive, 1 - prop_positive]).T
P_inf_non_inf = inv_A @
Results
return P_inf_non_inf
上の関数を用いて罹患率を求めるスクリプト例を、リストA1に示す。
リストA1では、sennsitivityとspecificityを以下のようにおいている。
sensitivity = 0.90
specificity = 0.999
この設定の下で、ダイアモンドプリンセス号の陽性率0.143に対して罹患率を求めると以下のようになる。
ダイヤモンド号
陽性率 = 0.143
罹患率 = 0.15795328142380427
非罹患率 = 0.8420467185761956
東京都の場合、https://stopcovid19.metro.tokyo.lg.jp/monitoringの記載によると、検査人数5441.4人に対する陽性率は3.4％となっていた（2022.01.07）。
この場合の、罹患率は以下のようになる。
Tokyo
陽性率 = 0.034
罹患率 = 0.036707452725250285
非罹患率 = 0.9632925472747496
罹患率は、検査対象母集団に対するものであることに注意。
import numpy as np
def EstPopInfctn( sensitivity,
specificity, prop_positive):
"""
検査母数団の罹患率の計算
"""
A = np.array([[sensitivity,
1 - specificity], [1 - sensitivity, specificity]])
inv_A = np.linalg.inv(A)
Results =
np.array([prop_positive, 1 - prop_positive]).T
P_inf_non_inf = inv_A @
Results
return P_inf_non_inf
sensitivity = 0.90
specificity = 0.999
print('sensi =', sensitivity)
print('specificity =', specificity)
print('\nダイヤモンド号')
p_positive = 0.143
# 陽性率
est_p_inf = EstPopInfctn(sensitivity,
specificity, p_positive)
print('陽性率
=', p_positive)
print('罹患率
=', est_p_inf[0])
print('非罹患率 =', est_p_inf[1])
print('\nTokyo')
p_positive = 0.034 # 東京都2022．0107
est_p_inf = EstPopInfctn(sensitivity,
specificity, p_positive)
print('陽性率
=', p_positive)
print('罹患率
=', est_p_inf[0])
print('非罹患率 =', est_p_inf[1])
罹患率を求めるベイズ分析は、次の式に基づいて行うことができる。
すなわち、sensitivityとspecificityが与えられているとき、陽性率は罹患率
の関数である。陽性率が与えられると、陽性者数は検査者総数における２項分布で表される。このモデルのStanスクリプトは以下のようになる。
data {
int N_Pos;
int N_All;
real sens; // sensitivity
real spec; // specificity
}
parameters {
real<lower = 0.0, upper =
1.0> prob_infect;
}
transformed parameters {
real<lower = 0.0, upper =
1.0> prob_pos; // probability of positive
prob_pos = sens *
prob_infect + (1.0 - spec) * (1.0 - prob_infect);
}
model {
prob_infect ~ uniform(0.0,
1.0);
N_Pos ~ binomial(N_All,
prob_pos);
}
上のStanスクリプトを用いるPythonスクリプトを、リストB1に示す。
ダイアモンドプリンセス号の場合、陽性者531人、検査総数3711人として分析すると、図B1の罹患率の事後分布を得る。
図B1
東京都の場合、https://stopcovid19.metro.tokyo.lg.jp/monitoringの記載によると、検査人数5441.4人に対する陽性率は3.4％となっていた（2022.01.07）ので、陽性者185人、検査者数5441人として分析すると、図B2の事後分布を得る。
図B2
罹患率は、検査対象母集団におけるものであることに注意。
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
import pystan
import pickle
# グラフの日本語表記の設定
plt.rcParams['font.family'] = ['Yu
Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
sm = pystan.StanModel(file =
'bysestpop.stan')
with open('sm.pkl', 'wb') as f:
pickle.dump(sm, f)
Data = {'N_Pos': 531, 'N_All': 3711,
'sens': 0.90, 'spec': 0.999}
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
sb.kdeplot(fit['prob_infect'])
plt.yticks([])
plt.xlabel('Probability', fontsize =
14)
plt.ylabel('Density', fontsize = 14)
plt.title('罹患率：ダイアモンドプリンセス号', fontsize =
18)
plt.savefig('FigDiamond.png')
plt.show()
Data = {'N_Pos': 185, 'N_All': 5441,
'sens': 0.90, 'spec': 0.999}
with open('sm.pkl', 'rb') as f:
sm = pickle.load(f)
fit = sm.sampling(data = Data, n_jobs
= 1)
print(fit)
sb.kdeplot(fit['prob_infect'])
plt.yticks([])
plt.xlabel('Probability', fontsize =
14)
plt.ylabel('Density', fontsize = 14)
plt.title('罹患率：東京都', fontsize = 18)
plt.savefig('FigTokyo.png')
plt.show()
ベイズ分析とPython用StanであるPyStanの解説を＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。