新型肺炎COVID-19の神戸における真のIgG抗体陽性患者数をrstanで推定する

R Rpackage Rを使いこなす rstan 医学 感染症学 数学いらずの医科統計学 統計

読んだ。
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人から陽性患者を検出する確率pがたしかにpであるときに、最尤推定値が\frac{33}{1000}で、二項分布から信頼区間を計算していいのである。すなわち、抗体検査の感度と特異度がともに100（のとき、accuracy が1）のときに、「1000人中、抗体検査陽性になる患者がpだけ存在する」=「1000人中、抗体検査陽性になった患者がpだけ存在する」と言っていいのであって、（抗体検査陽性）=（本当に陽性で、抗体検査も陽性）+（本当は陰性で、抗体検査は陽性）である状況が考慮されていない。
最近やった。
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 N=521

となり、この検査精度を推定したときの対象集団の事前確率rは、精度（A_c：accuracy）とすれば
A_c=\frac{rNS_n+(1-r)NS_p}{N}=A_c
r=\frac{A_c-S_p}{S_n-S_p}
で、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%程度なので、真に陽性である人が、検査陰性になってしまい、取りこぼしている可能性がある。というわけで、ただ単純に\frac{33}{1000}で二項分布を考えると、過小評価するおそれがある。
というわけでやった。
4.1%（2.6%〜5.9%）が真に抗体陽性の確率、となった。3.3%（2.3%〜4.6％）よりか多くなっており、論文では感度が75%であることを考慮していない分、過小評価となっている。
f:id:MikuHatsune:20200502225559p:plain
f:id:MikuHatsune:20200502230829p:plain

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%の抗体陽性となった。
f:id:MikuHatsune:20200502225419p:plain

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)


草生える（）