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

読んだ。

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

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

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

zip になってるshape ファイルをRのみの操作でダウンロードして扱いたい

という相談を受けた。
昔は地図を扱うときにsp というパッケージを使っていたが、最近ではsf が流行りだそうだ。
notchained.hatenablog.com
インタラクティブにしたければleaflet があるが、普通にプレゼンとか論文の図にしたい、ということでsf を使いたいらしい。
そもそもsf を使うのにいろいろ準備しないといけなくて、ubuntu ではgdal が2.0.0 以上でないとだめなようで、apt-get するのにもppa を更新して最新のものがインストールされるようにしないといけない。
How To Install GDAL/OGR Packages on Ubuntu — mothergeo 0.0.1 documentation

質問ではダウンロード先がzipになっている。
統計GISデータダウンロード | 政府統計の総合窓口
Pythonならzipのままデータを扱ってなかのデータ(csv だったり)を扱うことができたはずだが、Rではやったことがない(あったかもしれないが記憶にない)ので、ローカルにダウンロードしたものを解凍してshpファイルだけ読み込んだらRからファイルを削除するコマンドを実行する、というなんとも言えない作業をすることにした。
1. Simple Features for R • sf

library(sf)
library(stringr)

# 命名規則がよくわからなかったので5ファイル分取ってみる
dlfiles <- sprintf("https://www.e-stat.go.jp/gis/statmap-search/data?dlserveyId=A002005212015&code=131%02d&coordSys=1&format=shape&downloadType=5", 1:5)

sh <- vector("list", length(dlfiles))
for(f in seq(dlfiles)){
 # Download .shp data
 u_shp <- dlfiles[f]
 fname <- sprintf("%02d.zip", f)
 download.file(u_shp, fname)
 unzip(fname)
 sh[[f]] <- st_read(list.files(pattern="shp"))
 file.remove(list.files(pattern=gsub(".shp", "", list.files(pattern="shp")))) # ほかを消す
}

f:id:MikuHatsune:20200607003039p:plain

こんなことをしていたら日本地図を扱う必要が出てきたので逆に聞いたら、jpndistrict というパッケージを教えてもらった。
jpndistrict package | R Documentation

新型肺炎COVID-19の日本の実効再生産数を推定したrstanのコードを解説してみる

西浦先生が日本の実効再生産数を推定した。
コードはrstanで下記から取れる。
https://nbviewer.jupyter.org/github/contactmodel/COVID19-Japan-Reff/tree/master/

解説動画を見逃したのでコードと関連論文からのお勉強になるが、肝としては、
・知りたいのは「感染した日」である。
・診断日もしくは報告日は、データを収集して統計を取っているのでわかる。
・診断されるには検査される必要があるから、だいたい症状か接触歴があって、発症日はそこそこデータがある。
・感染した瞬間、はもちろん発病(はほとんど)していないのでわからない。
という前提がある。PDFの「患者」の観測データについて、の項。
(誰からから感染させられる)ー感染日ー発症日ー診断日/報告日という一連の流れについて、まったく情報がないわけではなく、いままでの数理モデルの知見からは
・誰かから感染させられるー感染日:serial interval もしくはgeneration time と書いてある。厳密には違うようだが、使用しているパラメータはワイブル分布から同一のものとしている。PDFでは「世代時間<=明示的にまだ得られていない間、発病間隔と同じと想定」の項。
・感染日ー発症日:incubation period. 対数正規分布から平均5.6日くらいと適当に設定している(適当にというのはランダム、というわけではなく、過去の推定で当てはまりがよさそうな、という意味合いで、対数正規分布ではくワイブル分布やガンマ分布でもやれないことはない)。
・発症日ー報告日:Right-truncated reporting delay. 打ち切りを考慮したワイブル分布を考慮しているようで、尤度関数L
L(\theta|S_k,O_k,T)=\displaystyle\prod_{k=1}^K\frac{h(S_k-O_k|\theta)}{H(T-O_k|\theta)}
S_kkさんの報告日
O_kkさんの発症日
T:カレンダー上の直近の観測日
h:ワイブル分布の確率密度(PDF)
H:ワイブル分布の累積確率密度(CDF)
となる。


さて、たぶんこれが一番大事な話だと思うが、感染日ー発症日、そして発症日ー診断日/報告日というのをなんとか推定するのに、backprojection というものを利用している。報告日はそれなりにデータがあるので、発症日←診断日/報告日のbackprojection をして発症日を推定し、感染日←発症日を再度して感染日を推定している。rstan のコードではこの部分がsurveillance というパッケージからbackprojNPという関数を使って最尤推定法で求めた点推定値を流用して実効再生産数を求めている。backprojection についての説明は皆無なので、むしろこれを理解するのに時間がかかった。


A method of non-parametric back-projection and its application to AIDS data. - PubMed - NCBI
今回のCOVID-19のように、診断日ははっきりしているが感染日がわからない疾患として、HIVを題材にしている。ある日付t, t=\{1,2.\dots T\}HIVの診断数がY_t とすると、d の期間のincubation period はf_d の確率(PDF)に従うとする。incubation period を考慮した真の推定値\lambda_tポアソン分布に従うとして、尤度関数は
\displaystyle\prod_{t=1}^T(\displaystyle\sum_{i=1}^t \lambda_i f_{t-i})^{y_t} \exp(-\displaystyle\sum_{i=1}^t \lambda_i f_{t-i})
を最大化することで得られる、となっているらしい。\displaystyle\sum_{i=1}^t \lambda_i f_{t-i} の部分が一般的なポアソン分布のパラメータ\lambda 部分だがy_t! で割らなくていいのだろうか。
これはGLMで解くのが難しかったらしいので、EMアルゴリズムで解けるようにしたのが論文である。
\lambda_t というのが日付t=\{1,2,\dots, T\} における真の(平均)感染数とすると、EステップとMステップを合体させた式が
\lambda_t^{new}=\frac{\lambda_t^{old}}{F_{T-t}}\displaystyle\sum_{d=0}^{T-t}\frac{Y_{t+d} f_d}{\displaystyle\sum_{i=1}^{t+d} \lambda_t^{old} f_{t+d-i}}
となる。ここで、F_{T-t}=\displaystyle\sum_{d=0}^{T-t}f_d である。\lambda_t が更新されたら差分でもなんでもとって適当な閾値\epsilon を下回れば終了する。
ただし、このまま\lambda_t を使うと結構ギザギザなので、smoothing を行なっている。これがExpectation-Maximization-Smoothing (EMS) と書いてある所以のようで、上で求めたのを\phi_t^{new}\gets \lambda_t^{new} にしておく。やり方は\phi_t^{new} の前後での重み付け平均をしているだけで、重みw_i={k \choose i}/2^k とする。k=2 だと[0.25, 0.5, 0.25] となる。これを使って
\lambda_t^{new}=\displaystyle\sum_{i=0}^k w_i \phi_{t+i-\frac{k}{2}}^{new} とする。k は偶数とし、あるi を中心に前後の重み付け平均をとる。

surveillanceではsurveillance:::backprojNP.fitおよびsurveillance:::em.step.beckerでEMSを行なっている。普通にやると遅いので、RccpC++を呼び出して行なっている。eq3aというmethod 名が論文での式3a(eq3a)である。
こんな感じになる。感染日と報告日は10-14日くらいずれていそうな感じがある。
f:id:MikuHatsune:20200513203350p:plain

入力するデータがそれなりに整ったので、後はrstanで実効再生産数R_t を推定する。平均的な時間t での感染者数は(domestic とimported でデータは分かれているが、簡便のため略)
E(C)=R_t\displaystyle\sum_{\tau =0}^{t-1}C_{total}(t-\tau)g(\tau)\frac{F(T-t)}{F(T-t+\tau)}
ここで
C_{total}(t):感染者数の総計
g:generation interval の確率密度(PDF)
F:感染日から報告日までの遅れの累積確率(CDF)。これはincubation period のPDFと発症日ー報告日の遅れのPDFとの畳み込みになる(感染日→発症日→報告日)。
である。後はポアソン分布を仮定する。

ここまでくるとrstanのコードがだいぶわかるようになってくる。
冒頭に使用するCPUコア数をオプションで指定すれば、利用可能な最大CPUで並列してくれる。コードでは2chain しか使ってないが、収束が非常に速いので困りはしない。

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

functionsブロックでは、好きな関数を定義してコード内で使えるようにする。rstanは離散確率分布は別に定義しないと柔軟に扱えないことが多く、今回のようにある日(離散)での確率密度はPMFになるので、累積確率分布の差分を取って柔軟にrstan内で扱えるようにする。
generation interval にはワイブル分布、incubation period には対数正規分布なので、それぞれ定義する。

functions {
    /* discretesized version of lognormal distribution */
    vector plognormal(real mu, real sigma, int K) {
        vector[K] res; 
        for (k in 1:K)
            res[k] = exp(lognormal_lcdf(k | mu, sigma));

        return append_row(res[1], res[2:K]-res[1:(K-1)]);        
    }

    /* discretesized version of Weibull distribution */
    vector pweibull(real kappa, real theta, int K) {
        vector[K] res; 
        for (k in 1:K)
            res[k] = -expm1(-pow(1.0*k/theta, kappa));

        return append_row(res[1], res[2:K]-res[1:(K-1)]);  
    }

確率分布の畳み込みを行うのも関数を定義する。rstanにはreverse関数がないようなのでYrevというのを後のtransformed dataブロックで作っているが、X[1:k]Y[k:1]でやれば同じである。ただし、rstanは普通のRみたいにY[k:1]とやるのは受け付けてくれなかった。

/* calculating the convolutions */
// X: first function, Yrev: reversed version of the second function
// K: length of X and Yrev
// the result is the vector of length K-1, because the first element equal zero is omitted
vector convolution(vector X, vector Yrev, int K) {
    vector[K-1] res;
    res[1] = X[1]*Yrev[K];
    for (k in 2:K-1) 
        res[k] = dot_product(head(X, k), tail(Yrev, k));

    return res;        
}

// Special form of the convolution with adjusting to the delay
// F: cumulative distribution function of the reporting delay
vector convolution_with_delay(vector X, vector Yrev, vector F, int K) {
    vector[K-1] res;
    vector[K] Z = X ./ F;

    res[1] = F[2]*Z[1]*Yrev[K];
    for (k in 3:K) 
        res[k-1] = F[k]*dot_product(head(Z, k-1), tail(Yrev, k-1));

    return res;        

modelブロック、R_t が2.4 くらいにいるだろうという事前分布はいいとして、scaleパラメータが1のガンマ分布からサンプリングしたのは単にRt .* conv が期待値になるようにしたかったからなんでしょうか。

model {
    Rt ~ normal(2.4, 2.0);
    Rt_adj ~ normal(2.4, 2.0);
    
    target += gamma_lpdf(tail(domestic_backproj, K-1) | Rt .* conv + 1e-13, 1.0)
            + gamma_lpdf(tail(domestic_backproj, K-1) | Rt_adj .* conv_delay_adj + 1e-13, 1.0);
}


backprojection はEMSでなくてもベイズ的に推定できるはずなので、これもrstan に組み込めばよさそうだが、まだこれができていないので保留。

新型肺炎COVID-19の山梨大学医学部附属病院でのPCR検査の結果をrstanで解析する

山梨大学医学部附属病院で、慶応大学病院と同じように入院中の患者にPCR検査を行うと、370人中全員が陰性だったらしい。
山梨大病院、すべての入院患者らにPCR検査 全員陰性 [新型コロナウイルス]:朝日新聞デジタル

山梨大の島田真路学長は「陽性者がいる場合も考え、(入院患者への検査を)先駆的に取り組んだ。全員が陰性となり、ほっとしている。入院患者へのPCR検査をやらないと安心できない。本来はすべての病院で積極的にやるべきだ」と語った。

慶応大学病院のPCR検査の結果を解析したように、まったく同じことをやってみる。
mikuhatsune.hatenadiary.com

PCRの感度を30%〜70%、特異度を>99% とすると、感染者割合は0.405%(95%信用区間 0.015%〜2.35%)、370人に検査したので1人(95%信用区間 0〜10人)は本当は陽性の患者がいると思われる。
PCR検査で安心できるね()
f:id:MikuHatsune:20200509174534p:plain

新型肺炎COVID-19の集中治療を要する患者の推移をSIRモデルを使ってrstanで推定する

こんなことをした。
mikuhatsune.hatenadiary.com
集中治療学会が、人工呼吸器を要している患者、ECMOをしている患者、ECMOで死亡した患者、ECMOから回復した患者、など日ベースで公開してる。
これに、毎日の感染者や死亡者のデータをくっつけて、SIRを使って人工呼吸器を要する患者やECMOを要する患者を推定する。
Japan Coronavirus: 9,231 Cases and 190 Deaths - Worldometer

SIRモデルとしてはこんなのを想定する。
SIRは普通にSIRモデルだが、Iから人工呼吸器管理になる人V、VからECMO管理になる人E、ECMOから死亡する人DEとして
Iから死亡する人はDIとすると、推移にパラメータpは9つ必要で
\frac{dS}{dt}=-\frac{p_1SI}{N}
\frac{dI}{dt}=\frac{p_1SI}{N}-(p_2I+p_6)I
\frac{dR}{dt}=p_2I+p_8V+p_9E
\frac{dV}{dt}=p_3I - (p_7+p_8)V
\frac{dE}{dt}=p_4V-p_9E
\frac{dD_I}{dt}=p_6I+p_7V
\frac{dD_E}{dt}=p_5E
となる。ただしN=S+I+R とする。
最初の人口は1.2億人で始める。
f:id:MikuHatsune:20200506171908p:plain

5月5日までのデータで既にかなりずれているようにも思えるが、SIRモデルの簡潔さからいうと仕方ないように思う。
f:id:MikuHatsune:20200506172907p:plain

pはこんな感じになった。例えばp_4 は人工呼吸器からECMOにいくパラメータで、p_5 はECMOから死亡にいくパラメータだが、これは1日あたりなので、適当に積分したらいいのだろうか。
f:id:MikuHatsune:20200506172957p:plain

365日後まで推定すると、感染者(軽症も含む)のピークは9月2日にきて、人工呼吸器を要する患者のピークは10月6日、ECMOを行なっている患者のピークは11月22日にくる。集中治療としては10-11月中が最大の山場、となるのかもしれない。
f:id:MikuHatsune:20200506173109p:plain

mat <- read.csv("japan.csv", stringsAsFactors=FALSE)       # worldmeter のデータ
ecmo <- merge(l222[,c("date", "全国")], merge(l555[,c("date", "全国")], l333, by.x="date", by.y="date"), by.x="date", by.y="date")
colnames(ecmo) <- c("date", "ECMOall", "ventilation", "ECMOrecovery", "ECMOdeath", "ECMOdoing")

mat$Date <- gsub("Apr ", "04-", gsub("Mar ", "03-", sprintf("2020-%s", gsub("Feb ", "02-", mat$Date))))
m <- as.data.frame(merge(ecmo, mat, by.x="date", by.y="Date"))

Pop <- 120000000

DE <- m$ECMOdeath
DI <- m$Total.Deaths - DE
E <- m$ECMOdoing
I <- m$Total.Cases
V <- m$ventilation - E
R <- I - m$Active.Cases
S <- Pop - (I + R + E + V + DI + DE)
N <- cbind.data.frame(S=S, I=I, R=R, V=V, E=E, DI=DI, DE=DE)


library(igraph)
G <- graph_from_edgelist(rbind(c(1,2),c(2,3),c(2,4),c(2,6),c(4,5),c(4,3),c(4,6),c(5,7),c(5,3)))
G <- graph_from_edgelist(rbind(c(1,2),c(2,3),c(2,4),c(4,5),c(5,7),c(2,6),c(4,6),c(4,3),c(5,3)))
lmat <- layout.norm(rbind(c(-1,0),c(0,0),c(1,0),c(0,-1),c(0,-2),c(1,-1),c(0,-3)))

pmat <- mapply(function(z) colMeans(lmat[get.edgelist(G)[z,],]), seq(nrow(get.edgelist(G))))
ppos <- c(3, 3, 2, 2, 2, 3, 3, 1, 1)


V(G)$label <- c("S", "I", "R", "V", "E", "DI", "DE")
V(G)$label.cex <- 2
V(G)$label.font <- 2
V(G)$label.color <- "black"
V(G)$size <- 20
V(G)$color <- cols
V(G)$color[6] <- "grey"
E(G)$color <- "black"
E(G)$width <- 3
# svg("fig02.svg", 48/9, 12/9)
par(mar=c(0, 0, 0, 0))
plot(G, layout=lmat)
legend("bottomleft", legend=sprintf("%s:\t%s", c("S", "I", "R", "V", "E", "DI", "DE"), names(cols)), pch=15, col=cols, bty="n", cex=1.5)
for(i in 1:9){
  #txt <- as.expression(substitute(italic(g[x%->%~y]), list(x=i, y=i+1)))
  txt <- as.expression(substitute(italic(p[x]), list(x=i)))
  text(pmat[1, i], pmat[2, i], txt, xpd=TRUE, pos=ppos[i], cex=1.5)
}


library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

code <- "
data{
  int<lower=1> T;
  int<lower=1> Tf;
  int<lower=1> N; // population
  int<lower=0> S0[7]; //S, I, R, V, E, DI, DE
  int<lower=0> S[T, 7]; //S, I, R, V, E, DI, DE
}
parameters{
  vector<lower=0, upper=0.5>[9] p;
}
transformed parameters{
  vector<lower=0>[7] y[T+Tf];
  y[1] = to_vector(S0);
  for(i in 1:(T+Tf-1)){
    y[i+1][1] = y[i][1] - p[1]*y[i][1]*y[i][2]/sum(y[i]);
    y[i+1][2] = y[i][2] + p[1]*y[i][1]*y[i][2]/sum(y[i]) - (p[2]+p[3]+p[6])*y[i][2];
    y[i+1][3] = y[i][3] + p[2]*y[i][2]+p[8]*y[i][4]+p[9]*y[i][5];
    y[i+1][4] = y[i][4] + p[3]*y[i][2] - (p[4]+p[7]+p[8])*y[i][4];
    y[i+1][5] = y[i][5] + p[4]*y[i][4] - (p[5]+p[9])*y[i][5];
    y[i+1][6] = y[i][6] + p[6]*y[i][2] + p[7]*y[i][4];
    y[i+1][7] = y[i][7] + p[5]*y[i][5];
  }
}
model{
  for(i in 1:T){
    S[i,] ~ multinomial(y[i]/sum(y[i]));
  }
}
"

m0 <- stan_model(model_code=code)
standata <- list(T=nrow(N), S=N, S0=unlist(N[1,]), N=Pop, Tf=365)
fit <- sampling(m0, standata, iter=7000, warmup=5000, chain=4)
ex <- extract(fit, pars=head(fit@model_pars, -1))

alpha <- 0.05
cia <- c(alpha/2, 0.5, 1-alpha/2)
m <- abind::abind(mapply(function(z) apply(ex$y[,,z], 2, quantile, cia), seq(dim(ex$y)[3]), SIMPLIFY=FALSE), along=3)

p <- apply(ex$p, 2, quantile, cia)

cols <- c("Susceptible"="blue", "Infected"="yellow", "Recovery"="skyblue", Ventilation="darkgreen", ECMO="orange", "Death infection"="black", "Death ECMO"="red")

library(vioplot)
par(mar=c(4.5, 5, 2, 2), cex.lab=1.5)
vioplot(as.data.frame(ex$p), ylab=as.expression(substitute(italic(g))), xlab="Parameter", horizontal=TRUE, side="right", colMed=NA, las=1)
text(p[2,], seq(ncol(p)), sprintf("%.3f", p[2,]), pos=1)
abline(h=seq(ncol(p)), lty=3)

xl <- c(0, 82)
yl <- c(0, 18000)
par(mar=c(4.5, 5, 4, 3))
matplot(m[2,,-1], type="l", col=cols[-1], xlim=xl, ylim=yl, lwd=4, lty=1, xlab="Date", ylab="人数")
for(i in 2:ncol(N)){
  points(N[,i], pch=16, col=cols[i])
}
legend("topleft", legend=names(cols), col=cols, pch=15, xpd=TRUE, ncol=1, cex=1.3)

xl <- c(0, dim(m)[2])
yl <- c(0, max(m[2,,-1]))
matplot(m[2,,-1], type="l", col=cols[-1], xlim=xl, ylim=yl, lwd=3, lty=1, xlab="", ylab="人数")
for(i in 2:ncol(N)){
  points(N[,i], pch=16, col=cols[i])
}

Mday <- apply(m[2,,], 2, which.max)
xl <- c(0, dim(m)[2])
yl <- c(0, 1)
xd <- 10
yd <- 0.03 + 0.05*(0:2)
d <- seq(as.Date("2020-02-15"), by="day", length.out=dim(m)[2])
r <- rle(format(d, "%h"))
par(mar=c(4.5, 5, 4, 3))
matplot(m[2,,]/Pop, type="l", col=cols, xlim=xl, ylim=yl, lwd=3, lty=1, xlab="Date", ylab="Proportion", xaxt="n", las=1, cex.lab=1.5)
axis(1, at=cumsum(r$lengths), labels=r$values, cex.axis=0.75)
legend("left", legend=names(cols), col=cols, pch=15, xpd=TRUE, ncol=1, cex=1.3)
j <- 0
pa <- par()$usr
for(i in c(5, 4, 2)){
  j <- j + 1
  # txt <- sprintf("Max %s on %d (%.2f)", names(cols)[i], Mday[i], m[2, Mday[i], i]/Pop)
  txt <- sprintf("Max %s on %s (%.2f)", names(cols)[i], format(d[Mday[i]], "%m/%d"), m[2, Mday[i], i]/Pop)
  segments(Mday[i], pa[4], y1=pa[4]+yd[j], xpd=TRUE, col=cols[i], lwd=3)
  segments(Mday[i], pa[4]+yd[j], x1=Mday[i]+xd, xpd=TRUE, col=cols[i], lwd=3)
  text(Mday[i]+xd, pa[4]+yd[j], txt, xpd=TRUE, pos=4, cex=1.2)
  abline(v=Mday[i], lty=3)
}
# worldmeter のデータ取り
library(stringr)
country <- "japan"
html <- sprintf("https://www.worldometers.info/coronavirus/country/%s", country)
h <- paste(readLines(html), collapse="")

s0 <- str_extract_all(h, '<script type="text/javascript">.*?</script>')[[1]]
s1 <- mapply(function(z) str_extract_all(z, " categories: \\[.*?\\]"), s0)
s2 <- lapply(mapply(function(z) str_replace(str_replace_all(z, "categories|[\":\\[\\]]", ""), "^ +", ""), s1), unique)
s3 <- unname(sapply(s2, strsplit, ","))


d1 <- mapply(function(z) str_extract_all(z, "data: \\[.*?\\]"), s0)
d2 <- lapply(mapply(function(z) str_replace(str_replace_all(z, "data|[\":\\[\\]]", ""), "^ +", ""), d1), unique)
d3 <- unname(lapply(lapply(d2, strsplit, ","), lapply, as.numeric))

n1 <- mapply(str_extract_all, s0, " title: \\{\\s+text: '.+?'")
n2 <- unname(rapply(n1, function(z) str_remove_all(str_extract(z, "text: '.+?'"), "text: |'"), how="replace"))

mat <- cbind.data.frame(s3[[1]], d3[[2]][[1]], d3[[3]][[1]], d3[[4]][[1]], d3[[5]][[1]], d3[[6]][[1]], d3[[7]][[1]], d3[[8]][[1]])
colnames(mat) <- c("Date", n2[[2]][[1]], n2[[3]][1], n2[[4]][1], n2[[5]][1], n2[[6]][1], n2[[7]][1], "Recoveries")

write.table(mat, sprintf("%s.csv", country), sep=",", row.names=FALSE)