読んだ。
Estimation of seroprevalence of novel coronavirus disease (COVID-19) using preserved serum at an outpatient setting in Kobe, Japan: A cross-sectional study. | medRxiv
神戸市の入院中の患者で、入院中の検査を適当（ランダムサンプリングのつもりで）にIgG抗体検査をすると、1000人中33人が抗体陽性だったということで、seropositive の頻度は3.3%（95%CI 2.3%〜4.6％）、と言っている。これは二項検定で計算できて
round(binom.test(33, 1000)$conf.int, 3)
[1] 0.023 0.046 attr(,"conf.level") [1] 0.95
となるが、これは1000人から陽性患者を検出する確率がたしかに
であるときに、最尤推定値が
で、二項分布から信頼区間を計算していいのである。すなわち、抗体検査の感度と特異度がともに100（のとき、accuracy が1）のときに、「1000人中、抗体検査陽性になる患者が
だけ存在する」=「1000人中、抗体検査陽性になった患者が
だけ存在する」と言っていいのであって、（抗体検査陽性）=（本当に陽性で、抗体検査も陽性）+（本当は陰性で、抗体検査は陽性）である状況が考慮されていない。
最近やった。
mikuhatsune.hatenadiary.com
mikuhatsune.hatenadiary.com
KURABOの抗体検査は、論文によると（SARS-CoV-2 Antibody Detection Kit (IgM/IgG)）、IgG検査は521人に行われ、感度 0.7638、特異度 1、精度 0.9424だった、という。これを一生懸命分割表で逆算すると
| 真に陽性 | 真に陰性 | 合計 | |
| 検査陽性 | 感度 97 | 0 | 97 |
| 検査陰性性 | 30 | 特異度 394 | 424 |
| 127 | 394 |
となり、この検査精度を推定したときの対象集団の事前確率は、精度（
：accuracy）とすれば
で、24.4%のようである。
Sn <- 0.7638 Sp <- 1 Ac <- 0.9424 N_KURABO <- 521 r <- (Ac - Sp)/(Sn - Sp) M0 <- matrix(c(r*Sn, (1-r)*(1-Sp), r*Sn+(1-r)*(1-Sp), r*(1-Sn), (1-r)*Sp, r*(1-Sn)+(1-r)*Sp, r, (1-r), 1), 3, byrow=TRUE) M <- round(M0*N_KURABO) M[1,1]/M[3,1] M[2,2]/M[3,2] (M[1,1]+M[2,2])/M[3,3]
年齢構成を考えずに、とりあえずバルクの集団として、「真の抗体陽性の人の割合」を推定する。抗体検査の感度はほぼ100%なので、陽性になってしまえばその人はほぼ確実に真に陽性の患者（SPNINという）であるが、感度は75%程度なので、真に陽性である人が、検査陰性になってしまい、取りこぼしている可能性がある。というわけで、ただ単純にで二項分布を考えると、過小評価するおそれがある。
というわけでやった。
4.1%（2.6%〜5.9%）が真に抗体陽性の確率、となった。3.3%（2.3%〜4.6％）よりか多くなっており、論文では感度が75%であることを考慮していない分、過小評価となっている。
library(rstan) rstan_options(auto_write=TRUE) options(mc.cores=parallel::detectCores()) code <- " data{ int N_IgG; // 神戸集団でのIgG検査数 1000 int N_IgG_positive; // 神戸集団でのIgG検査陽性数 33 int N_KURABO_P; // KURABO検査での陽性結果 int N_KURABO_TP; // KURABO検査での真の陽性結果 int N_KURABO_N; // KURABO検査での陰性結果 int N_KURABO_TN; // KURABO検査での真の陰性結果 } parameters{ real<lower=0, upper=1> p; // IgG陽性の有病率 real<lower=0, upper=1> Sn; // KURABOの検査の感度 real<lower=0, upper=1> Sp; // KURABOの検査の特異度 } transformed parameters{ real<lower=0, upper=1> theta = p*Sn + (1-p)*(1-Sp); // KURABO検査で陽性になる確率 } model{ N_KURABO_P ~ binomial(N_KURABO_TP, Sn); // KURABO検査の感度 N_KURABO_N ~ binomial(N_KURABO_TN, Sp); // KURABO検査の特異度 N_IgG_positive ~ binomial(N_IgG, theta); // 神戸集団でのIgG陽性 } " m0 <- stan_model(model_code=code) standata <- list(N_IgG=1000, N_IgG_positive=33, N_KURABO_P=M[1,1], N_KURABO_TP=M[3,1], N_KURABO_N=M[2,2], N_KURABO_TN=M[3,2]) fit <- sampling(m0, standata, iter=30000, warmup=10000, chain=4) ex <- extract(fit, pars=head(fit@model_pars, -1)) alpha <- 0.05 cia <- c(alpha/2, 0.5, 1-alpha/2) labs <- c("IgG positive prevalenve in Kobe", "KURABO Sensitivity", "KURABO Specificity", "KURABO Positive result probability") cols <- c("lightgreen", "violet", "violet", "violet") xl <- c(0, 1) par(mfcol=c(4, 1), mar=c(4.5, 4.5, 3.5, 2), cex.axis=1.5) 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%.3f [%.3f, %.3f]", labs[i], CrI[2], CrI[1], CrI[3]), freq=FALSE, xlab="Probability", nclass=50, border=NA, col=cols[i], cex.main=1.5) } cols <- c(rgb(144,238,144,200,maxColorValue=256), rgb(238,130,238,120,maxColorValue=256)) hist(ex$p, xlab="Prevalence", nclass=100, border=NA, col=cols[1], cex.main=1.5, main="") hist(ex$theta, xlab="Probability", nclass=100, border=NA, col=cols[2], cex.main=1.5, add=TRUE) legend("topright", legend=c("真の推定", "素の推定"), col=cols, pch=15, cex=2, bty="n")
次に、年齢構成がデータとしてあるので、グループ別にやった。ただ分けてfor loopを回すだけである。
10歳未満や90歳以上は標本が少なく、信用区間の幅が非常に広くなってしまったが、それ以外の階級では概ね5%〜10%くらいが真に抗体陽性となった。男女別でみても6.5%〜6.7%とほとんど差はなく、男女合わせて年齢構成を考慮した場合、6.7%の抗体陽性となった。
Age_M <- c(7, 17, 17, 41, 64, 84, 87, 85, 81, 6) Age_F <- c(1, 10, 19, 49, 91, 80, 84, 81, 83, 13) Positive_M <- c(0, 0, 0, 1, 2, 4, 2, 6, 1, 0) Positive_F <- c(0, 0, 0, 2, 3, 2, 3, 3, 4, 0) names(Age_M) <- names(Age_F) <- c("10<", apply(matrix(c(1, head(rep(2:9, each=2), -1)), nr=2)*10, 2, paste, collapse="-"), ">90") code <- " data{ int Cat; int Age_M[Cat]; // 年齢区分の数 int Age_F[Cat]; // 年齢区分の数 int Positive_M[Cat]; int Positive_F[Cat]; int N_KURABO_P; // KURABO検査での陽性結果 int N_KURABO_TP; // KURABO検査での真の陽性結果 int N_KURABO_N; // KURABO検査での陰性結果 int N_KURABO_TN; // KURABO検査での真の陰性結果 } parameters{ vector<lower=0, upper=1>[Cat] pM; // IgG陽性の有病率 vector<lower=0, upper=1>[Cat] pF; // IgG陽性の有病率 real<lower=0, upper=1> Sn; // KURABOの検査の感度 real<lower=0, upper=1> Sp; // KURABOの検査の特異度 } transformed parameters{ vector<lower=0, upper=1>[Cat] thetaM = pM*Sn + (1-pM)*(1-Sp); // KURABO検査で陽性になる確率 vector<lower=0, upper=1>[Cat] thetaF = pF*Sn + (1-pF)*(1-Sp); // KURABO検査で陽性になる確率 } model{ N_KURABO_P ~ binomial(N_KURABO_TP, Sn); // KURABO検査の感度 N_KURABO_N ~ binomial(N_KURABO_TN, Sp); // KURABO検査の特異度 Positive_M ~ binomial(Age_M, thetaM); // 神戸集団でのIgG陽性 Positive_F ~ binomial(Age_F, thetaF); // 神戸集団でのIgG陽性 } generated quantities{ int RM[sum(Age_M)]; int RF[sum(Age_F)]; for(i in 1:Cat){ for(j in ((i==1)?1:(sum(Age_M[1:(i-1)])+1)):sum(Age_M[1:i])){ RM[j] = bernoulli_rng(pM[i]); } for(j in ((i==1)?1:(sum(Age_F[1:(i-1)])+1)):sum(Age_F[1:i])){ RF[j] = bernoulli_rng(pF[i]); } } } " m1 <- stan_model(model_code=code) standata <- list(Cat=length(Age_M), Age_M=Age_M, Age_F=Age_F, Positive_M=Positive_M, Positive_F=Positive_F, N_KURABO_P=M[1,1], N_KURABO_TP=M[3,1], N_KURABO_N=M[2,2], N_KURABO_TN=M[3,2]) fit <- sampling(m1, standata, iter=10000, warmup=5000, chain=4) ex <- extract(fit, pars=head(fit@model_pars, -1)) mapply(function(j) (ifelse(j==1, 1, sum(Age_M[1:(j-1)])+1)):sum(Age_M[1:j]), 1:10) library(vioplot) v <- lapply(list(ex$pM, ex$pF), as.data.frame) v <- mapply(function(z1, z2) cbind.data.frame(z1, rowMeans(z2)), list(ex$pM, ex$pF), list(ex$RM, ex$RF), SIMPLIFY=FALSE) cols <- c("skyblue", "pink") Mcols <- c("blue", "red") xl <- c(0, 1.1) xd <- 0.23 P <- lapply(v, apply, 2, quantile, cia) labs <- c(names(Age_M), "All") par(mar=c(4.5, 4.5, 3.5, 9), cex.axis=1.5, cex=1.25) plot(rev(seq(labs)), type="n", xlim=xl, ylim=c(0.5, length(labs)+0.5), yaxt="n", xlab="Prevalence", ylab="") axis(2, at=rev(seq(labs)), labels=labs, las=1) title("Age adjusted seroprevalence in Kobe") pa <- par()$usr for(i in seq(v)){ vioplot(v[[i]], at=rev(seq(labs)), horizontal=TRUE, side=c("right", "left")[i], add=TRUE, rectCol=Mcols[i], col=cols[i], colMed=Mcols[i], lineCol=NA) for(j in seq(labs)){ txt <- sprintf("%.3f [%.3f, %.3f]", P[[i]][2, j], P[[i]][1, j], P[[i]][3, j]) text(pa[2], length(labs)+1-j+xd*c(1, -1)[i], txt, xpd=TRUE, pos=4) text(pa[2], length(labs)+1-j+xd*c(1, -1)[i], c("M", "F")[i], xpd=TRUE, pos=2, font=2) } } PAll <- quantile(rowMeans(cbind(ex$RM, ex$RF)), cia) txt <- sprintf("%.3f [%.3f, %.3f]", PAll[2], PAll[1], PAll[3]) text(pa[2], pa[3]+xd, txt, xpd=TRUE, pos=4) text(pa[2], pa[3]+xd, "All", xpd=TRUE, pos=2, font=2) abline(h=rev(seq(labs)), lty=3)
なんかベイズの定理がまったく誤解されてるリプ多いがフィッシャー先生も否定し、多くの年輩ドクターも理解してない概念だから仕方ないと言えば仕方ない。が、検査学の勉強せずに検査を語るなかれ。
— 岩田健太郎 Kentaro Iwata (@georgebest1969) April 25, 2020
草生える（）
