交互作用項のある回帰分析で標準化偏回帰係数を求めるのに標準化はいつするべきか

みたいな相談を受けた。
回帰分析は最近ほとんどしないし、そもそもしたとしても交互作用の項は変数が増えるし解釈も面倒になるのでしたことがないのが本音だが、聞かれたので考えた。
結論から言うと、変数は先に標準化して、そして積を取るようである。
irisを使うわけだが、Sepal.LengthSepal.Width, Petal.Length, Petal.Widthで回帰することを考える。ここで、変数を標準化しない場合、標準化しないデータで得た回帰係数に対して、「目的変数の標準偏差/説明変数の標準偏差」、をかけると標準化した場合のデータで得た回帰係数、つまり標準化偏回帰係数が得られる。

dat <- iris[,-5]
dat_scale <- as.data.frame(scale(dat))

lm0 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data=dat)
lm1 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width, data=dat_scale)

standard_coef_coef <- apply(dat, 2, sd)/sd(dat$Sepal.Length)

# no interaction
lm0$coefficients * standard_coef_coef
lm1$coefficients
lm0$coefficients * standard_coef_coef
 (Intercept)  Sepal.Width Petal.Length  Petal.Width 
   1.8559975    0.3425789    1.5117505   -0.5122442 

lm1$coefficients
  (Intercept)   Sepal.Width  Petal.Length   Petal.Width 
-1.176045e-16  3.425789e-01  1.511751e+00 -5.122442e-01 

確かに一致する。

次に、Petal.LengthPetal.Widthが交互作用を持つと仮定した場合を考える。Rのlm関数で自動的にするならば、入力データを標準化しない場合と、標準化した場合はそれぞれlm2, lm3となる。

lm2 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length * Petal.Width, data=dat)
lm3 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length * Petal.Width, data=dat_scale)

ここで、Rのformula内で*を使えば交互作用となるが、これは事前に入力データ内に、積の変数を作っておけばよい。ここで、dat1は、積の変数を作ってから標準化したもので、dat2は、標準化した変数を使って積の変数にしたものである。

dat0 <- cbind.data.frame(dat, Petal.Length_Petal.Width=dat$Petal.Length*dat$Petal.Width)
dat1 <- as.data.frame(scale(dat0))
dat2 <- cbind.data.frame(dat_scale, Petal.Length_Petal.Width=dat_scale$Petal.Length*dat_scale$Petal.Width)

lm内のformula をベタに書くと

lm4 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Petal.Length_Petal.Width, data=dat0)
lm5 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Petal.Length_Petal.Width, data=dat1)
lm6 <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Petal.Length_Petal.Width, data=dat2)

ここで、標準化していない場合のデータで、*を使って交互作用を書いた場合lm2と、積の変数を作って回帰した場合lm4

lm2

Call:
lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length * Petal.Width, 
    data = dat)

Coefficients:
             (Intercept)               Sepal.Width              Petal.Length               Petal.Width  
                 2.24236                   0.57835                   0.65778                  -0.83072  
Petal.Length:Petal.Width  
                 0.06164  

lm4

Call:
lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + 
    Petal.Length_Petal.Width, data = dat0)

Coefficients:
             (Intercept)               Sepal.Width              Petal.Length               Petal.Width  
                 2.24236                   0.57835                   0.65778                  -0.83072  
Petal.Length_Petal.Width  
                 0.06164  

となり一致する。
一方で、標準化した場合のデータだが、

lm3

Call:
lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length * Petal.Width, 
    data = dat_scale)

Coefficients:
             (Intercept)               Sepal.Width              Petal.Length               Petal.Width  
                 -0.0958                    0.3044                    1.5599                   -0.5515  
Petal.Length:Petal.Width  
                  0.1002  

lm5

Call:
lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + 
    Petal.Length_Petal.Width, data = dat1)

Coefficients:
             (Intercept)               Sepal.Width              Petal.Length               Petal.Width  
               1.556e-16                 3.044e-01                 1.402e+00                -7.647e-01  
Petal.Length_Petal.Width  
               3.508e-01  

lm6

Call:
lm(formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + 
    Petal.Length_Petal.Width, data = dat2)

Coefficients:
             (Intercept)               Sepal.Width              Petal.Length               Petal.Width  
                 -0.0958                    0.3044                    1.5599                   -0.5515  
Petal.Length_Petal.Width  
                  0.1002  

となり、*記法に一致するのは、標準化した変数を使って積の変数にした場合のlm6である。
では、標準化しないで交互作用を含む回帰をしたlm2の係数を、標準化して交互作用を含む回帰をしたものに一致させるかは、単純に「目的変数の標準偏差/説明変数の標準偏差」ではうまくいかなかった。おそらく共分散とかそういうのを使えばなんとかなるのだろうと思ったが数式をいじらずにやっていたので結局答えが出ないまま力尽きてしまった。
たぶん線形代数的に真面目にやるとわかりそうな気がする(やるとは言ってない
mikuhatsune.hatenadiary.com

参考
回帰分析による交互作用の検証
心理統計の注意点:重回帰分析についての注意点

Rのプロット領域の角に文字をおきたい

論文のsub figure でAとかBとかを角に置きたいが、角の座標を取得したい。
Rのプロット領域の仕様として、余白 margin と描出されない外領域がある。
box 関数を使うと余白領域または外領域を囲ってくれるが、これは.External.graphics で関数を参照するのを諦めたのでpar 関数で取れる値でゴリ押しする。

まず、plot 関数で呼び出されるプロット領域は、基本的にインチで情報を持っているようである。これは外領域par()$din、余白含めたプロット領域par()$fin、実際にプロットされる領域par()$pin で決まる。
実際にプロットされる領域は、par()$usrで四隅が取得できるが、これはプロット領域の実数、つまりインチではないので、x軸(またはy軸)描出領域を実際にプロットされる領域のインチ表示であるpar()$pinで補正すればなんとかなる。

なんとかした結果がこれである。
f:id:MikuHatsune:20200902151613p:plain

xs <- c("pin"=10, "fin"=80)
pt.cex <- 3
text.cex <- 1.5
par(oma=c(2, 3, 4, 1), font=2, xpd=TRUE)
xl <- c(0, 100)
plot(1, type="n", xlim=xl, xlab="", ylab="mgp[1]")
for(z in list(c(1,3), c(2, 3), c(1, 4), c(2, 4))){
  points(par()$usr[z[1]], par()$usr[z[2]], col="white", pch=15, cex=pt.cex)
  text(par()$usr[z[1]], par()$usr[z[2]], sprintf("usr[%d, %d]", z[1], z[2]), xpd=TRUE, font=2, col="red", cex=text.cex)
}

box("figure", lty=3)
box("outer", lwd=15, col=4)

segments(xs["pin"], par()$usr[3], y1=par()$usr[4], lwd=3)
points(xs["pin"], mean(par()$usr[3:4]), col="white", pch=15, cex=pt.cex)
text(xs["pin"], mean(par()$usr[3:4]), "pin", col="red", cex=text.cex)

y0 <- par()$usr[4]+diff(par()$usr[3:4])/par()$pin[2]*par()$mai[3]
y1 <- par()$usr[4]
segments(xs["fin"], y0, y1=y1, lwd=3)
points(xs["fin"], mean(c(y0, y1)), col="white", pch=15, cex=pt.cex)
text(xs["fin"], mean(c(y0, y1)), "fin", col="red", cex=text.cex)

y0 <- par()$usr[3]-diff(par()$usr[3:4])/par()$pin[2]*par()$mai[1]
y1 <- par()$usr[3]
segments(xs["fin"], y0, y1=y1, lwd=3)
points(xs["fin"], mean(c(y0, y1)), col="white", pch=15, cex=pt.cex)
text(xs["fin"], mean(c(y0, y1)), "fin", col="red", cex=text.cex)

x0 <- par()$usr[1]-diff(par()$usr[1:2])/par()$pin[1]*par()$mai[2]
x1 <- par()$usr[2]+diff(par()$usr[1:2])/par()$pin[1]*par()$mai[4]
segments(x0, y0, x1, lwd=3)
text(mean(c(x0, x1)), y0, "din is ourter and fin is figure inch", pos=3, col="red", cex=text.cex)

y0 <- par()$usr[4]+diff(par()$usr[3:4])/par()$pin[2]*par()$mai[3]
x0 <- par()$usr[1]-diff(par()$usr[1:2])/par()$pin[1]*par()$mai[2]
text(x0, y0, "A", adj=c(0, 1), cex=3)

y0 <- par()$usr[3]-diff(par()$usr[3:4])/par()$pin[2]*par()$mai[1]
x0 <- par()$usr[1]-diff(par()$usr[1:2])/par()$pin[1]*par()$mai[2]
text(x0, y0, "B", adj=c(0, 0), cex=3)

x1 <- par()$usr[2]+diff(par()$usr[1:2])/par()$pin[1]*par()$mai[4]
text(x1, y0, "C", adj=c(1, 0), cex=3)

y0 <- par()$usr[4]+diff(par()$usr[3:4])/par()$pin[2]*par()$mai[3]
text(x1, y0, "D", adj=c(1, 1), cex=3)

for(i in -1:8){
  mtext(sprintf("line=%d", i), side=3, line=i)
  mtext(sprintf("line=%d", i), side=2, line=i, adj=0.2)
}

二次計画法を久しぶりにやる

二次計画法をRでやる必要が出てきたので久しぶりにやってみる。
二次計画法はRはでは


ただし、行列D が正定値行列でないと

matrix D in quadratic function is not positive definite!

と言われてしまうので、逆行列を持つように微妙に変化させた行列を無理やり作る。

R のときの最適化はmin(-d^T b + \frac{1}{2} b^T D b) with the constraints A^T b \geq b_0 を行うので、それに合わせて引数となる行列とベクトルを書き換える。
meq の引数は上からmeq 個の条件式が等号の場合であり、すべてが不等式ならmeq=0 とする。

例題 http://courses.csail.mit.edu/6.867/wiki/images/a/a7/Qp-cvxopt.pdf

\min(\frac{1}{2}x^2+3x+4y
\min\frac{1}{2}\begin{bmatrix}x\\y\end{bmatrix}^T \begin{bmatrix}1&0\\0&0\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix}+\begin{bmatrix}3\\4\end{bmatrix}^T \begin{bmatrix}x\\y\end{bmatrix}
と書けるのでこれに合うように書く。

library(quadprog)
library(Matrix)
min(-d^T b + 1/2 b^T D b) with the constraints A^T b >= b_0

Amat <- rbind(c(1, 0, 1, -2, -3), c(0, 1, 3, -5, -4))
bvec <- c(0, 0, 15, -100, -80)
dvec <- -c(3, 4)
Dmat <- matrix(c(1, 0, 0, 0), 2)
Dmat <- nearPD(Dmat)$mat
solve.QP(Dmat, dvec, Amat, bvec=bvec, meq=0)$solution
$solution
[1] 0 5
from cvxopt import matrix
from cvxopt import solvers
P = matrix([[1.0,0.0],[0.0,0.0]])
q = matrix([3.0,4.0])
G = matrix([[-1.0,0.0,-1.0,2.0,3.0],[0.0,-1.0,-3.0,5.0,4.0]])
h = matrix([0.0,0.0,-15.0,100.0,80.0])
sol = solvers.qp(P,q,G,h)
list(sol["x"])
[7.131816408856328e-07, 5.000001008391809]

賭ケグルイ(双)のスリーヒットダイスをやってみる

賭ケグルイ(双)の1巻に、スリーヒットダイスという賭けがある。
賭ケグルイ双(漫画)- マンガペディア

聚楽幸子が提案したゲーム。ダイスの123を「DOWN」、456を「UP」として扱い、プレイヤーはDOWNとUPの「連続する3つの出目」を予想してカードに書き、ディーラーが1つのダイスを振り続け、予想した3つの出目が先に出たプレイヤーの勝ちというルール。DOWNは「D」、UPは「U」として表記する。例としてプレイヤー1「UDD」対プレイヤー2「DDU」という予想の勝負であれば、ダイスの出目が「1(D)4(U)5(U)2(D)3(D)」という順番で出現するとなると、523の出目でUDDが成立しプレイヤー1の勝利となる。

DかUが出る確率は各々\frac{1}{2}で、それらが3連になるのはいずれも(\frac{1}{2})^3 だろうと思っていたら、「これは典型的な確率の錯誤」とアメリが言っており、かつ、「連続する3つの出目」のほかにも「出目同士の相性」もあるという二段構えである。

シミュレーションではDUU、DDU、UUD、UDDが早めに出やすく、UDU、DUDがやや遅く、UUU、DDDが最も回数を要する。
f:id:MikuHatsune:20200727162232p:plain

次に、出目で勝てるかどうかを確かめてみる。最速で出目が出るからといって常に勝つとは限らない。例えば、DUUは最速(平均 8回、最頻値 7回)で出るが、UDUに勝てるかどうかは五分五分であるし、DDUに対しては\frac{1}{3} しか勝てない。
UUUはDUUに\frac{1}{8} しか勝てない。というのも、UUUがDUUに勝てるのは最初から3回目までUUUの時だからである。
_8C_2通りの出目の賭け方で、平均的に最も勝てるのはDUUかUDD、ということになる。

 →相手

自分

UUU DUU UDU DDU UUD DUD UDD DDD
UUU 0.000 0.124 0.403 0.299 0.499 0.416 0.397 0.503
DUU 0.876 0.000 0.503 0.333 0.753 0.502 0.500 0.601
UDU 0.597 0.497 0.000 0.375 0.334 0.499 0.497 0.583
DDU 0.701 0.667 0.625 0.000 0.501 0.664 0.248 0.497
UUD 0.501 0.247 0.666 0.499 0.000 0.625 0.665 0.701
DUD 0.584 0.498 0.501 0.336 0.375 0.000 0.500 0.603
UDD 0.603 0.500 0.503 0.752 0.335 0.500 0.000 0.874
DDD 0.497 0.399 0.417 0.503 0.299 0.397 0.126 0.000
ud <- c("U", "D")
ud3 <- apply(expand.grid(ud, ud, ud), 1, paste, collapse="")

n <- 30000
s_generator <- function(size=300){
  tmp <- sample(ud, size=size, replace=TRUE)
  s <- apply(mapply(function(z) tmp[z:(length(tmp)-3+z)], 1:3), 1, paste, collapse="")
  return(s)
}
s0 <- replicate(n, s_generator())

s1 <- mapply(function(z) match(ud3, s0[,z]), seq(ncol(s0)))
m1 <- rowMeans(s1, na.rm=TRUE) + 2
m2 <- apply(s1, 1, median, na.rm=T) + 2

rownames(s1) <- ud3
s2 <- mapply(function(z) table(factor(s1[z,], 1:max(s1, na.rm=TRUE))), seq(ud3))

cols <- jet.colors(length(ud3))
txt <- sprintf("%s (mean=%d, median=%d)", ud3, round(m1), m2)
par(mar=c(5, 5, 2, 2))
matplot(seq(nrow(s2))+2, s2/n, type="l", xlim=c(0, 49), col=cols, lty=1, lwd=4, cex.lab=1.3,
        xlab="最初に3連組み合わせが出現するサイコロ施行回数", ylab="density")
# legend("topright", legend=ud3, pch=15, col=cols, ncol=2, cex=2)
legend("topright", legend=txt, pch=15, col=cols, ncol=1, cex=1.1, y.intersp=1.5)

# 勝負
res <- diag(0, length(ud3))
rownames(res) <- colnames(res) <- ud3
for(i in 1:n){
  m <- match(ud3, s0[,i])
  for(j in seq(ud3)){
    res[j, -j] <- res[j, -j] + (m[j] < m[-j])
  }
}

# 勝率
round(res/n, 3)

内科診療ストロングエビデンス

読んだ。

内科診療 ストロング・エビデンス

内科診療 ストロング・エビデンス

COI:中古でワンコインだったので買ったがそもそも3500円だったらしい。

面白そうだったのが貧血の項で、術後Hb が8を下回った人について、宗教上の理由で輸血を拒否した人の予後調査というのがあった。
Mortality and Morbidity in Patients With Very Low Postoperative Hb Levels Who Decline Blood Transfusion - PubMed
Hb 5 でもなんとか生きるけど、やはり死亡率は上がる。

2014年に更に別の研究が出ていた。
An Update on Mortality and Morbidity in Patients With Very Low Postoperative Hemoglobin Levels Who Decline Blood Transfusion (CME) - PubMed
当然のことながらHb 1 下がるごとに死亡オッズは2倍くらいあがっていくらしい。

市中肺炎の項で、de-escalation をしても予後は変わらないらしいのは興味深かった。
Comparison Between Pathogen Directed Antibiotic Treatment and Empirical Broad Spectrum Antibiotic Treatment in Patients With Community Acquired Pneumonia: A Prospective Randomised Study - PubMed

手術に来る人が予定、緊急かぎらず、どんな内科的治療を受けてきているか、というのをおさらいするのにはまあよかった。
敗血症治療など集中治療的なこととか、術前にどう管理してもっていくか(前述の貧血とか)の話もあったが、ザ・内科医って感じで実際の周術期管理のことはあんま詳しくなさそうだった。

エビデンスにはうるさそうで、とにかくエビデンスをならべ、ROCAUC大好きでc-statistics のことばかり書いていた。
ベイズもやったことがあるようで、JAGSの話が1ページだけ書いてあったのは興味深かった。

中古ワンコインなら買いだが、新品ならオススメはしない。

Rのformula をdata.frame の色々な変数の組み合わせで内容を変えながらやりたいのだが

というような質問を受けた。

glm関数とかでみかける、
A ~ B + C
みたいな式の、A,B,Cをfor loopとかで内容変えながら生成する方法はないかとおもって。
A,B,Cはデータの列名やけど、その名前を変数に格納して渡しても無理でさ。
いちいちもとのデータの列名そのものをA,B,Cとかに変えてやるっていうのも思いついたけど、正攻法とは思えず。

いやおそらくその方法のほうが正攻法っぽいと思われる。

やりたいことを素直に表現するなら、評価式はA ~ B + C と書くのを規則に従ってたくさん文字列として生成して、その文字列をeval 関数に渡すのだが文字列なのでparse 関数で渡してゴリ押ししよう。

# 適当なサンプルデータのdataframe
cols <- head(colnames(iris), -1)

# A ~ B + C のBCの組み合わせ
cmb <- combn(length(cols)-1, 2)

# 4番目のPetal.Width を Petal.Width ~ B + C と推定する
# 普通に全部素直に書いたら
glm(Petal.Width ~ Sepal.Length + Sepal.Width, data=iris)
glm(Petal.Width ~ Sepal.Length + Petal.Length, data=iris)
glm(Petal.Width ~ Sepal.Width + Petal.Length, data=iris)

# 上記式を文字列的に書いて
# 強引に評価式としてRに計算させる

for(i in 1:ncol(cmb)){ # すべてのパターンの数
  # 評価したい式を文字列で生成する
  equation <- sprintf("glm(Petal.Width ~ %s + %s, data=iris)", cols[cmb[1,i]], cols[cmb[2,i]])
  # 文字列をRで計算させる
  print(eval(parse(text=equation)))
}
Call:  glm(formula = Petal.Width ~ Sepal.Length + Sepal.Width, data = iris)

Coefficients:
 (Intercept)  Sepal.Length   Sepal.Width  
     -1.5635        0.7233       -0.4787  

Degrees of Freedom: 149 Total (i.e. Null);  147 Residual
Null Deviance:	    86.57 
Residual Deviance: 22.25 	AIC: 147.5

Call:  glm(formula = Petal.Width ~ Sepal.Length + Petal.Length, data = iris)

Coefficients:
 (Intercept)  Sepal.Length  Petal.Length  
   -0.008996     -0.082218      0.449376  

Degrees of Freedom: 149 Total (i.e. Null);  147 Residual
Null Deviance:	    86.57 
Residual Deviance: 6.144 	AIC: -45.58

Call:  glm(formula = Petal.Width ~ Sepal.Width + Petal.Length, data = iris)

Coefficients:
 (Intercept)   Sepal.Width  Petal.Length  
     -0.7065        0.0994        0.4263  

Degrees of Freedom: 149 Total (i.e. Null);  147 Residual
Null Deviance:	    86.57 
Residual Deviance: 6.082 	AIC: -47.12

新型肺炎COVID-19の厚生労働省が行なった抗体検査から集団の有病率をrstanで推定する

こんな記事を観測した。
新型コロナウイルス感染症に関する検査について|厚生労働省
Roche社とAbbott社が売っている抗体検査キットを使って、東京、大阪、宮城の住民を無作為に抽出した結果、各社での陽性陰性結果は以下のようになった、という。Roche社での陽性をR+, Abbott社での陽性をA+ とする。

東京 A+ A- 合計
R+ 2 4 6
R- 2 1963 1965
合計 4 1967 1971

 

大阪 A+ A- 合計
R+ 5 5 10
R- 11 2949 2960
合計 16 2954 2970

 

宮城 A+ A- 合計
R+ 1 6 7
R- 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

抗体を持っている人(有病率p)の推定にはこれらの感度、特異度のゆらぎも考慮しよう。

さて、2種類の検査を行なって、それぞれの陽性・陰性の分割表になっているが、これらの分割表が出来る背景には、Roche社の感度N_R、特異度C_R、Abbott社の感度N_A、特異度C_A として
真陽性p の人がR+ かつA+pN_R N_A
真陽性p の人がR+ かつA-pN_R(1-N_A)
真陽性p の人がR- かつA+p(1-N_R)N_A
真陽性p の人がR- かつA-p(1-N_R)(1-N_A)
真陰性(1-p) の人がR+ かつA+(1-p)(1-C_R)(1-C_A)
真陰性(1-p) の人がR+ かつA-(1-p)(1-C_R)C_A
真陰性(1-p) の人がR- かつA+(1-p)C_R(1-C_A)
真陰性(1-p) の人がR- かつA-(1-p)C_R C_A
となる。ただし、Roche社とAbbott社の検査は独立としている。これを分割表にすると

A+ A-
R+ pN_R N_A+(1-p)(1-C_R)(1-C_A) pN_R(1-N_A)+(1-p)(1-C_R)C_A
R- p(1-N_R)N_A+(1-p)C_R(1-C_A) p(1-N_R)(1-N_A)+(1-p)C_R C_A

となる。こうした上で、N人の集団を検査したと考える。

結果としては東京 0.14%(0.033-0.35)、大阪 0.29%(0.13-0.51)、宮城 0.055%(0.0078-0.18)となった。
f:id:MikuHatsune:20200617212047p:plain

解析ではRoche社、Abbott社の両方とも陽性となったときに「陽性」としているようなので、上記8パターンのうち両方とも陽性(真陽性p の人がR+ かつA+pN_R N_Aと、真陰性(1-p) の人がR+ かつA+(1-p)(1-C_R)(1-C_A)のふたつの場合)を取り出すと

        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])