新型肺炎COVID-19の抗体検査から集団の有病率をrstanで推定する

結論から言うと、抗体検査を受けた202人の集団において、9.4%（95% 信用区間で5.4%〜13.9%）が感染者である。ただし、この抗体検査を受けた202人は、東京都（とか他の都道府県一般）のランダムサンプリングとは言えないので、単純に東京都人口の9.4%が感染している、とは言えないことに注意。

こんな記事を観測した。
https://www.tokyo-np.co.jp/s/article/2020043090070748.htmlwww.tokyo-np.co.jp
解析に必要な数値の事実ベースで抽出すると

一般市民の4.8%、医療従事者の9.1%が陽性（抗体あり）で、過去に感染していたことが分かった。
一般市民の147人の4.8%％にあたる7人が陽性、医療従事者55人のうち9.1%の5人が陽性だった。
ホームページで希望者を募り、20～80歳の男性123人、女性79人を検査した。
★このうち一カ月以内に発熱のあった人は52人
★同居者でコロナウイルス感染者がいる人は2人
★PCR検査を受診したことがある人は9人
★PCR検査で陽性反応だった1人も含む
★市民・医療従事者を合計した202人全体では5.9%の12人（男女とも6人）が陽性だった。
以前のPCR検査で陰性とされたが、抗体検査で陽性だった人もいた。

このうち、★は今回の解析で使用したデータである。

興味があるのは、「202人中、抗体検査で陽性だったのは12人（5.9%）」という結果から、「202人中、感染しているのはX人」を推定することである。検査しているのだから感染者（抗体検査なので既感染を含むが、面倒なので「今、感染している」こととする）検査陽性の人のことだろう、と思うかもしれないが、それは臨床検査学的にはよくある間違い（というかベイズ統計屋的に間違い）で、検査陽性というのは、
（本当に感染していて、検査で正しく陽性）+（本当は感染していないのに、検査で間違って陽性）
というパターンの足し算の結果、「検査陽性」となっているので、（本当は感染していない人）を加算していることに注意しなければならない。
が、検査の特性として、（本当は感染していないのに、検査で間違って陽性）としてしまう確率などがわかれば、ある程度正しく推定できる。これらの性能は感度と特異度と言われ、検査として売り出すなら絶対に評価している。

ここでは、PCR検査は感度 0.1-0.5、特異度0.99程度、とした。一方で、抗体検査は感度0.3-0.7、特異度0.8-0.99程度とした。というのも、インフルエンザの抗体検査がこれくらいだからである。気になる人は適当にパラメータをいじって追試してほしい。

さて、推定したいのは「感染している人の割合」であるが、202人は属性がバラバラである。一ヶ月以内に発熱のあった人、同居者でコロナウイルス感染者がいる人、PCR検査を受信したことがある人（は、つまり陰性だったのだろうと推測）、PCR検査で陽性だった人、を除外すると、202人中138人は、一体どういう人なのかよくわからないが、ホームページで募った希望者、ということから推測するに、無症状の人とする。つまり、
$p_1$ ：無症状の138人のうち本当に感染している人の割合
$p_2$ ：このうち一カ月以内に発熱のあった52人のうち本当に感染している人の割合
$p_3$ ：同居者でコロナウイルス感染者がいる2人のうち本当に感染している人の割合
$p_4$ ：PCR検査を受診したことがある9人（で、PCR検査陰性だった）のうち本当に感染している人の割合
$p_5$ ：PCR検査で陽性反応だった1人のうち本当に感染している人の割合
をそれぞれ推定することになる。ただし、202人が5つのカテゴリに属し、総計で12人検査陽性、という制限がある。
ここで、変数は5つあるが、PCR検査は「怪しそうな人にやる」のが原則である。「一カ月以内に発熱があった」というのは、いま感染しているかどうかにはあまり寄与しない（※抗体検査で既感染を疑う設定なら意味があるが、この解析では現時点の感染を意味しているので注意）。しかし、無症状よりかは感染していそうな感じはするので、 $p_1< p_2$ くらいの緩い制限を持ちつつ、PCR検査は、事前確率 $p_3$ の集団に行われる、とする。
すると、PCR検査陰性（ $p_4$ ）とPCR検査陽性（ $p_5$ ）の運命がふたつともある。ここで、事前確率 $p_3$ のときの陰性事後確率と陽性事後確率は、実はそれぞれ計算できて
$p_4=\frac{p_3 (1-S_n)}{p_3 (1-S_n)+(1-p_3)(1-S_p)}$ （PCR検査の陰性事後確率）
（ $p_4=\frac{p_3 (1-S_n)}{p_3 (1-S_n)+(1-p_3)S_p}$ こっちです）
$p_5=\frac{p_3S_n}{p_3S_n + (1-p_3)(1-S_p)}$ （PCR検査の陽性事後確率）
となる。ただし $S_n$ は感度、 $S_p$ は特異度である。
$p$ は各カテゴリの感染者割合（確率）なので、実際に（抗体検査で）陽性になる確率 $q$ は
$q=pS_n + (1-p)(1-S_p)$
である。ここで $S_n$ 、 $S_p$ は抗体検査の感度と特異度である。
こちらでも似たようなことをやった。
mikuhatsune.hatenadiary.com

5つのカテゴリに $q_{*}=(q_1, q_2,q_3, q_4, q_5)$ と抗体検査で陽性となる確率がついていて、それぞれ138, 52, 2, 9, 1 人がいて、検査陽性となる人が合計12人なので、ひたすらfor loop で回してbinomial_lpmf関数でtarget += する、としてみる。多項分布とかもっといい方法があれば教えてほしい。
最終的にはgenerated quantities ブロックで、202人5カテゴリの人たちが感染者かそうでないかを $p$ でランダムサンプリング（bernoulli_rng でできる）することでシミュレーションし、最終的に「202人のうち感染者の割合」を推定する。

9.4%（95% CI 5.4%〜13.9%）が感染者と推定された。
f:id:MikuHatsune:20200501163624p:plain

$p_1$ ：無症状の138人のうち本当に感染している人の割合は1.67%
$p_2$ ：このうち一カ月以内に発熱のあった52人のうち本当に感染している人の割合は10.7%
$p_3$ ：同居者でコロナウイルス感染者がいる2人のうち本当に感染している人の割合は85.1%
$p_4$ ：PCR検査を受診したことがある9人（で、PCR検査陰性だった）のうち本当に感染している人の割合は80.1%
$p_5$ ：PCR検査で陽性反応だった1人のうち本当に感染している人の割合は99.8%
だった。
f:id:MikuHatsune:20200501163649p:plain
本当は「医療従事者であるか」とか「男性か女性か」の属性も加味したかったが断念した。

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

code <- "
data{
  int Ncat;        //カテゴリの数
  int Ns[Ncat]; //カテゴリごとの人数
  int N_pos;     //陽性患者 12人
  real<lower=0, upper=1> Sn_PCR[2];  // PCR検査の感度の範囲
  real<lower=0, upper=1> Sp_PCR[2];  // PCR検査の特異度の範囲
  real<lower=0, upper=1> Sn_Anti[2];   // 抗体検査の感度の範囲
  real<lower=0, upper=1> Sp_Anti[2];   // 抗体検査の特異度の範囲
  real<lower=0, upper=1> b1[2];           // 無症状の人の感染者割合の範囲
  real<lower=0, upper=1> b2[2];           // 一ヶ月以内に発熱した人の感染者割合の範囲
  real<lower=0, upper=1> b3[2];           // 同居している人がコロナウイルス感染者の人の感染者割合の範囲
  real<lower=0, upper=1> b4[2];           // PCR陰性者の感染者割合の範囲
  real<lower=0, upper=1> b5[2];           // PCR陽性者の感染者割合の範囲
}
parameters{
  real<lower=Sn_PCR[1], upper=Sn_PCR[2]> Sensitivity; // PCR検査の感度
  real<lower=Sp_PCR[1], upper=Sp_PCR[2]> Specificity; // PCR検査の特異度
  real<lower=Sn_Anti[1], upper=Sn_Anti[2]> Sn;                // 抗体検査の感度
  real<lower=Sp_Anti[1], upper=Sp_Anti[2]> Sp;                // 抗体検査の特異度
  real<lower=b1[1], upper=b1[2]> p1;
  real<lower=b2[1], upper=b2[2]> p2;
  real<lower=b3[1], upper=b3[2]> p3;
}
transformed parameters{
  real<lower=b4[1], upper=b4[2]> p4 = p3*(1-Sensitivity)/(p3*(1-Sensitivity) + (1-p3)*Specificity);
  real<lower=b5[1], upper=b5[2]> p5 = p3*Sensitivity/(p3*Sensitivity + (1-p3)*(1-Specificity));
  vector<lower=0, upper=1>[Ncat] theta;
  theta[1] = p1*Sn + (1-p1)*(1-Sp);
  theta[2] = p2*Sn + (1-p2)*(1-Sp);
  theta[3] = p3*Sn + (1-p3)*(1-Sp);
  theta[4] = p4*Sn + (1-p4)*(1-Sp);
  theta[5] = p5*Sn + (1-p5)*(1-Sp);
}
model{ // 合計12人が陽性になるようにひたすら回す
  for(i in 0:Ns[5]){
    target += binomial_lpmf(i | Ns[5], theta[5]);
    for(j in 0:Ns[3]){
      target += binomial_lpmf(j | Ns[3], theta[3]);
      for(k in 0:(Ns[4]-i-j)){
        target += binomial_lpmf(k | Ns[4], theta[4]);
        for(l in 0:(N_pos-i-j-k)){
          target += binomial_lpmf(l | Ns[2], theta[2]);
          target += binomial_lpmf(N_pos-l-k-j-i | Ns[1], theta[1]);
        }
      }
    }
  }
}
generated quantities{
  int e[sum(Ns)];
  for(i in 1:Ns[1]) e[i] = bernoulli_rng(p1);
  for(i in sum(Ns[1:(2-1)]):sum(Ns[1:2])) e[i] = bernoulli_rng(p2);
  for(i in sum(Ns[1:(3-1)]):sum(Ns[1:3])) e[i] = bernoulli_rng(p3);
  for(i in sum(Ns[1:(4-1)]):sum(Ns[1:4])) e[i] = bernoulli_rng(p4);
  for(i in sum(Ns[1:(5-1)]):sum(Ns[1:5])) e[i] = bernoulli_rng(p5);
}
"
m0 <- stan_model(model_code=code)

Ns <- c(138, 52, 2, 9, 1)
standata <- list(Ncat=length(Ns), Ns=Ns, N_pos=12, 
                 Sn_PCR=c(0.1, 0.5), Sp_PCR=c(0.99, 1),
                 Sn_Anti=c(0.3, 0.7), Sp_Anti=c(0.8, 0.99),
                 b1=c(0.0001, 0.3), b2=c(0.05, 0.8), b3=c(0.5, 0.99), b4=c(0, 1), b5=c(0.3, 1)
                 )

fit <- sampling(m0, standata, iter=5000, warmup=1000, chain=4)
ex <- extract(fit, pars=head(fit@model_pars, -1))

labs <- c("Sensitivity of PCR", "Specificity of PCR", "Sensitivity of immunoassay", "Specificity of immunoassay", "Prevalence of asymptomatic", "Prevalence of fever within a month", "Prevalence of housemate positive", "Prevalence of PCR negative", "Prevalence of PCR positive")
alpha <- 0.05
cia <- c(alpha/2, 0.5, 1-alpha/2)

xl <- c(0, 1)
par(mfcol=c(2, 5), mar=c(4.5, 4.5, 2.5, 2))
for(i in seq(labs)){
  CrI <- quantile(ex[[i]], cia)
  # sprintf("%.2f [%.2f, %.2f]", CrI[2], CrI[1], CrI[3])
  hist(ex[[i]], main=sprintf("%s\n%.2f [%.2f, %.2f]", labs[i], CrI[2], CrI[1], CrI[3]), freq=FALSE, xlab="Probability", nclass=50, border=NA, col="violet", xlim=xl)
}

CrI <- quantile(rowMeans(ex$e), cia)
hist(rowMeans(ex$e), main=sprintf("%s\n%.2f [%.2f, %.2f]", "Prevalence", CrI[2], CrI[1], CrI[3]), freq=FALSE, xlab="Probability", nclass=50, border=NA, col="lightgreen", xlim=xl)