新型肺炎COVID-19の無症状感染者の割合を慶応大学病院のPCR検査の結果からrstanで推定する

慶応大学病院で、入院中で（新型肺炎COVID-19について）無症状の患者67人をPCR検査すると、6%（=4人）の患者で陽性だった、という話。
某感染症の先生は、めちゃめちゃ単純な推定を使って、一番楽観的な推定で47万人、一番悲観的な推定で430万人、東京に無症状の感染者がいるのではないか、と言っている。

慶応のPCR6％の意味 https://t.co/1TIY5jCUrz
— 岩田健太郎 Kentaro Iwata (@georgebest1969) April 22, 2020

慶応のPCR6％の意味 - 楽園はこちら側

この推定が雑なので、もう少し真面目に考えてみる。
PCR検査を受けたのが67人、PCR陽性だったのが4人、であることは確定したデータであるので、これを使って、無症状（67人は全員無症状である）なのに感染している人の割合を $p$ 、PCR検査の感度を $Sn$ 、特異度を $Sp$ とする。
$N=$ 67人の運命は

	真に感染( $p$ )	真に非感染( $1-p$ )	合計
検査で陽性	$NpS_n$	$N(1-p)(1-S_p)$	$4=N\{pS_n+(1-p)(1-S_p)\}$
検査で陰性	$Np(1-S_n)$	$N(1-p)S_p$	$63=N\{p(1-S_n)+(1-p)S_p\}$
	$Np$	$N(1-p)$	$N=67$

となる。
全体の $N=67$ 人から確率 $\theta=pS_n+(1-p)(1-S_p)$ で4人の陽性者が出る（真陽性と偽陽性も含む）とすればよい。ここで、PCRの特性として、感度は30-70%、特異度は99%以上、というのが有識者の認識のようなので、感度は $S_n\sim U(0.3, 0.7)$ 、特異度は $S_p\sim U(0.99, 1)$ からサンプリングされるとする。これで区間付き推定がそれなりにできるはず。

結果は $p=0.14 [0.04, 0.37]$ となった。およそ14%の人が無症状で感染しているようである。
f:id:MikuHatsune:20200424004911p:plain

 2.5%   50% 97.5%
 73 0.040 0.141 0.366

この結果を単純に東京都の人口1395万人に当てはめると、196万人が無症状で感染しているようである。ただし、慶応大学病院に入院している人が、東京都の人口のランダムサンプリングとはそう簡単には言えないので、もっと多いか、もっと少ないかすらもよくわからない。

   2.5%     50%   97.5% 
 560364 1965582 5102252

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
code <- "
data{
  int N_keio;
  int P_keio;
  int Pop; // 対象としている人口
  real<lower=0, upper=1> Sn[2];  // 感度の範囲
  real<lower=0, upper=1> Sp[2];  // 特異度の範囲
}
parameters{
  real<lower=0, upper=1> p; // prior probability
  real<lower=Sn[1], upper=Sn[2]> Sensitivity;
  real<lower=Sp[1], upper=Sp[2]> Specificity;
}
transformed parameters{
  real<lower=0, upper=1> theta;
  theta = p*Sensitivity + (1-p)*(1-Specificity);
}
model{
  P_keio ~ binomial(N_keio, theta);
}
generated quantities{
  int e;
  e = binomial_rng(Pop, p);
}
"
m0 <- stan_model(model_code=code)

N_keio <- 67
P_keio <- 4

standata <- list(N_keio=N_keio, P_keio=P_keio, Sn=c(0.3, 0.7), Sp=c(0.99, 1), Pop=13950000)
fit <- sampling(m0, standata, iter=1000, warmup=500, chain=4)
ex <- extract(fit, pars=head(fit@model_pars, -1))

alpha <- 0.05
cia <- c(alpha/2, 0.5, 1-alpha/2)
hist(ex$p, xlab="無症状性有病者の割合", ylab="頻度", freq=FALSE, las=1, main="", col="yellow")

round(quantile(ex$p, cia), 3) # 無症状感染者の割合
round(quantile(ex$e, cia))     # 東京都の人口で考えたときの無症状感染者の人数