こんな記事を観測した。
新型コロナウイルス感染症に関する検査について｜厚生労働省
Roche社とAbbott社が売っている抗体検査キットを使って、東京、大阪、宮城の住民を無作為に抽出した結果、各社での陽性陰性結果は以下のようになった、という。Roche社での陽性を, Abbott社での陽性を
とする。
| 東京 | 合計 | ||
| 2 | 4 | 6 | |
| 2 | 1963 | 1965 | |
| 合計 | 4 | 1967 | 1971 |
| 大阪 | 合計 | ||
| 5 | 5 | 10 | |
| 11 | 2949 | 2960 | |
| 合計 | 16 | 2954 | 2970 |
| 宮城 | 合計 | ||
| 1 | 6 | 7 | |
| 2 | 3000 | 3002 | |
| 合計 | 3 | 3006 | 3009 |
この結果を持って、例えば東京で安直に2/1971=0.10% が抗体を持っている、と言っている。
これは感度100%、特異度100%のときのみにそう言えるのであって、偽陽性を考慮していないのはこれまでに何回もやったとおりである。
mikuhatsune.hatenadiary.com
mikuhatsune.hatenadiary.com
mikuhatsune.hatenadiary.com
mikuhatsune.hatenadiary.com
ではRoche社とAbbott社の抗体検査キットはいったいどれくらいの感度、特異度なのだろうと思って調べてみると、例えばRoche社ではElecsys Anti-SARS-CoV-2 が感度29/29 （PCRをして確定診断2週間後の患者で）、特異度が5262/5272 （10例がreactive つまり偽陽性）と言っている。
Elecsys® Anti-SARS-CoV-2
Abbott社は論文が出ていて、感度は13/16 （Roche社と同様にPCR確定して2周間以上経過した患者）、特異度 152/153 と言っている。
Clinical Performance of Two SARS-CoV-2 Serologic Assays | Clinical Chemistry | Oxford Academic
抗体を持っている人（有病率）の推定にはこれらの感度、特異度のゆらぎも考慮しよう。
さて、2種類の検査を行なって、それぞれの陽性・陰性の分割表になっているが、これらの分割表が出来る背景には、Roche社の感度、特異度
、Abbott社の感度
、特異度
として
真陽性 の人が
かつ
：
真陽性 の人が
かつ
：
真陽性 の人が
かつ
：
真陽性 の人が
かつ
：
真陰性 の人が
かつ
：
真陰性 の人が
かつ
：
真陰性 の人が
かつ
：
真陰性 の人が
かつ
：
となる。ただし、Roche社とAbbott社の検査は独立としている。これを分割表にすると
となる。こうした上で、人の集団を検査したと考える。
結果としては東京 0.14%（0.033-0.35）、大阪 0.29%（0.13-0.51）、宮城 0.055%（0.0078-0.18）となった。
解析ではRoche社、Abbott社の両方とも陽性となったときに「陽性」としているようなので、上記8パターンのうち両方とも陽性（真陽性 の人が
かつ
：
と、真陰性
の人が
かつ
：
のふたつの場合）を取り出すと
prevalence 0 1 2.5% 0.0003344003 0.9974306 0.0002557886 50% 0.0014207519 0.9989492 0.0010508271 97.5% 0.0034640029 0.9997442 0.0025693801
となり、事前確率（0.0014）より低くなる（1が陽性で 0.00105）。感度が29/29 や13/16 であっても、区間推定とともに考えれば
binom.test(29, 29)$conf.int
[1] 0.8805551 1.0000000 attr(,"conf.level") [1] 0.95
と90%を下回ることがあるので、90%弱のものを2回行なって陽性を判定すると、そもそもの事前確率が0.001と非常に小さいので、陽性を拾ってくるのが困難になる。
ではいずれかひとつの検査が陽性になった場合に「陽性」と判定する場合は
prevalence 0 1 2.5% 0.0003344003 0.9929788 0.003520887 50% 0.0014207519 0.9950443 0.004955743 97.5% 0.0034640029 0.9964791 0.007021175
となり、0.0014 から0.00496 に上がる。がしかし、もともと低い事前確率のものを検査して事後確率が0.005になったからといって、では本当に感染（今回は抗体検査なので既感染だが）しているかどうかを判断するとなると、たぶん違うと思う。
library(vioplot) library(rstan) rstan_options(auto_write=TRUE) options(mc.cores=parallel::detectCores()) code <- " data{ int N[3]; int P[4, 3]; // R+A+ R+A- R-A+ R-A- int<lower=0> Sn_N[2]; // 感度検査の人数 int<lower=0> Sp_N[2]; // 特異度検査の人数 int<lower=0> Sn_P[2]; // 感度検査の陽性人数 int<lower=0> Sp_P[2]; // 特異度検査の陰性人数 } parameters{ real<lower=0, upper=1> p[3]; // prior probability real<lower=0, upper=1> Sn[2]; real<lower=0, upper=1> Sp[2]; } transformed parameters{ real<lower=0, upper=1> theta[4, 3]; for(i in 1:3){ theta[1, i] = p[i]*Sn[1]*Sn[2] + (1-p[i])*(1-Sp[1])*(1-Sp[2]); // Roche + Abbott + theta[2, i] = p[i]*Sn[1]*(1-Sn[2]) + (1-p[i])*(1-Sp[1])*Sp[2]; // Roche + Abbott - theta[3, i] = p[i]*(1-Sn[1])*Sn[2] + (1-p[i])*Sp[1]*(1-Sp[2]); // Roche - Abbott + theta[4, i] = p[i]*(1-Sn[1])*(1-Sn[2]) + (1-p[i])*Sp[1]*Sp[2]; // Roche - Abbott - } } model{ for(i in 1:2){ Sn_P[i] ~ binomial(Sn_N[i], Sn[i]); Sp_P[i] ~ binomial(Sp_N[i], Sp[i]); } for(j in 1:3){ for(i in 1:4){ P[i, j] ~ binomial(N[j], theta[i, j]); } } } //generated quantities{ // int T[N[1]]; // Tokyo // int O[N[2]]; // Osaka // int M[N[3]]; // Miyagi // for(i in 1:N[1]){ // T[i] = categorical_rng(to_vector(theta[,1])); // } // for(i in 1:N[2]){ // O[i] = categorical_rng(to_vector(theta[,2])); // } // for(i in 1:N[3]){ // M[i] = categorical_rng(to_vector(theta[,3])); // } //} " m0 <- stan_model(model_code=code) # Roche社、Abbott社の検査性能（この順で並んでいる） Sn_N <- c(29, 16) Sp_N <- c(5272, 153) Sn_P <- c(29, 13) Sp_P <- c(5262, 152) # 各都市の結果 R+A+ R+A- R-A+ R-A- で並んでいる dat <- list(Tokyo=c(2, 4, 2, 1963), Osaka=c(5, 5, 11, 2949), Miyagi=c(1, 6, 2, 3000)) standata <- list(N=sapply(dat, sum), P=do.call(cbind, dat), Sn_N=Sn_N, Sp_N=Sp_N, Sn_P=Sn_P, Sp_P=Sp_P) fit <- sampling(m0, standata, iter=3000, warmup=1000, chain=4) ex <- extract(fit, pars=head(fit@model_pars, -1)) p <- ex$p alpha <- 0.05 cia <- c(alpha/2, 0.5, 1-alpha/2) pcia <- apply(p, 2, quantile, cia)*100 cols <- c("pink", "orange", "green") par(mar=c(4.5, 4, 2, 2), cex.lab=1.5) vioplot(as.data.frame(p), horizontal=TRUE, side="right", rectCol=NA, drawRect=FALSE, col=cols, xlab="Prevalence [%]", ylab="", yaxt="n") axis(1, cex.axis=1.6) axis(2, at=seq(dat), labels=names(dat), cex.axis=2) abline(h=seq(dat), lty=3) for(i in seq(dat)){ txt <- sprintf("%.2f%s [%.2f, %.2f]", pcia[2,i], "%", pcia[1,i], pcia[3,i]) text(par()$usr[1], i, txt, adj=c(-0.1, 1.5), cex=1.5) } # Tokyo について考える i <- 1 tmp <- cbind("R+A+"=ex$p[,i]*ex$Sn[,1]*ex$Sn[,2], "R+A-"=ex$p[,i]*ex$Sn[,1]*(1-ex$Sn[,2]), "R-A+"=ex$p[,i]*(1-ex$Sn[,1])*ex$Sn[,2], "R+A-"=ex$p[,i]*(1-ex$Sn[,1])*(1-ex$Sn[,2]), "R+A+"=(1-ex$p[,i])*(1-ex$Sp[,1])*(1-ex$Sp[,2]), "R+A-"=(1-ex$p[,i])*(1-ex$Sp[,1])*ex$Sp[,2], "R-A+"=(1-ex$p[,i])*ex$Sp[,1]*(1-ex$Sp[,2]), "R-A-"=(1-ex$p[,i])*ex$Sp[,1]*ex$Sp[,2]) idx <- mapply(function(z) replace(rep(0, 8), z, 1), list(both=grep(".\\+.\\+", colnames(tmp)), any=grep("\\+", colnames(tmp))), SIMPLIFY=FALSE) # ふたつの検査の一致具合で陽性を判断する apply(rbind(prevalence=ex$p[,i], apply(tmp, 1, tapply, idx$both, sum)), 1, quantile, cia) apply(rbind(prevalence=ex$p[,i], apply(tmp, 1, tapply, idx$any, sum)), 1, quantile, cia) # 検査を複数回したときの事後確率を考える Ppostnegative <- function(pre, Sn, Sp){ pre*(1-Sn)/(pre*(1-Sn) + (1-pre)*Sp) } Ppostpositive <- function(pre, Sn, Sp){ pre*Sn/(pre*Sn + (1-pre)*(1-Sp)) } Sn <- c(0.8, 0.9) Sp <- c(0.95, 0.98) pre <- 0.001 Ppostnegative(pre, Sn[1], Sp[1]) Ppostnegative(Ppostnegative(pre, Sn[1], Sp[1]), Sn[2], Sp[2]) Ppostnegative(pre, Sn[2], Sp[2]) Ppostnegative(Ppostnegative(pre, Sn[2], Sp[2]), Sn[1], Sp[1]) Ppostpositive(Ppostnegative(pre, Sn[1], Sp[1]), Sn[2], Sp[2]) Ppostpositive(pre, Sn[1], Sp[1]) Ppostpositive(Ppostpositive(pre, Sn[1], Sp[1]), Sn[2], Sp[2]) Ppostpositive(pre, Sn[2], Sp[2]) Ppostpositive(Ppostpositive(pre, Sn[2], Sp[2]), Sn[1], Sp[1])
