新型肺炎COVID-19が流行してから全身麻酔件数がどれほど減ったかをrstanで推定する

こんなツイーヨを観測した。

コロナウィルス騒動で全国のオペ室今どうなってるのかアンケート7回目

全国の麻酔科の先生におたずねします。
今（4/24）現在、オペ室は、
— DAJ (@dajhiroki) 2020年4月24日

コロナウィルス騒動で全国のオペ室今どうなってるのかアンケート6回目

全国の麻酔科の先生におたずねします。
今（4/17）現在、オペ室は、
— DAJ (@dajhiroki) 2020年4月17日

コロナウィルス騒動で全国のオペ室今どうなってるのかアンケート5回目

全国の麻酔科の先生におたずねします。
今（4/10）現在、オペ室は、
— DAJ (@dajhiroki) 2020年4月10日

コロナウィルス騒動で全国のオペ室今どうなってるのかアンケート4回目

全国の麻酔科の先生におたずねします。
今（4/6）現在、オペ室は、
— DAJ (@dajhiroki) 2020年4月5日

コロナウィルス騒動で全国のオペ室今どうなってるのかアンケート３回目

全国の麻酔科の先生におたずねします。
今（3/27）現在、オペ室は、
— DAJ (@dajhiroki) 2020年3月27日

コロナウィルス騒動、全国のオペ室はどうなっているのかアンケート第二弾

麻酔科の先生にお尋ねします。
手術室は、
— DAJ (@dajhiroki) 2020年3月19日

コロナウィルス騒動でオペ室っていまみなさんのところではどんな感じになってるんでしょうか…

麻酔科の先生にお尋ねします。
今（3/13）現在、
— DAJ (@dajhiroki) 2020年3月13日

全身麻酔の導入において、気管挿管は飛沫感染のリスクが高いので、不要不急の手術は控えるように麻酔科学会、各種外科系学会から提言が出ている。
とはいっても、不要不急の手術ってなんぞや、ということで、手術をしないと死ぬような緊急疾患は麻酔をするとして、癌手術なんかもいずれは手術できなくなりそうである。

上記では、ツイッターで適当に麻酔科から「手術室の運営方針として、手術をどれくらい減らしているか」を聞いていて、(1)制限なし（100%運用）、(2)やや制限（半分以上として70%運用）、(3)制限（半分以下として40%運用）、(4)予定手術はしない、緊急手術のみ対応（12%）、と勝手にした。というのも、制限具合をベータ分布からサンプリングする上でいい感じに設定する必要があるからである。
平均 $m$ 、分散 $v=0.01$ として、各設定についてベータ分布のパラメータ $\alpha$ 、 $\beta$ は、平均 $\frac{\alpha}{\alpha + \beta}$ 、分散 $\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$ となるから、
$\alpha=\frac{m^2(1-m)}{v-m}$ 、 $\beta=(1-m)(\frac{m(1-m)}{v}-1)$
ベータ分布を使って、平均的な手術室稼働率の低下具合はこんな感じとする。
f:id:MikuHatsune:20200426225806p:plain

上記のアンケート結果を使って、手術室稼働率の低下具合の週ごとの推移はこんな感じになる。通常通り運営していた施設はどんどん減って、少し制限（半分以上）しながら運営している施設が増えている。大きく制限（半分以下）している施設も遅れて増えてきているが、予定手術はしない、という施設はまだそこまで増加していない。しかし、手術しない、という施設が増えるのも時間の問題のようである。
f:id:MikuHatsune:20200426225909p:plain
推定の問題として、アンケート回答は4つなので、どの回答をするかというパラメータ $\theta$ はsimplex であるが、通常通り、という回答をする施設、というか割合は確実に減るはずなので、あるアンケート回答週 $t$ において、通常通り、という回答が最初の項だとすると $\theta_{t,1}<\theta_{t+1,1}$ というようにしてorderedでパラメータを作りたかったが、うまくいかなかった。
$\theta_t$ については時系列サンプリングにしている。あまりよくなかった。

さて、こういう事態になって、実際の全身麻酔件数はどれくらい減っているのかを推定する。全国での全身麻酔の件数は統計があるようなのでここからexcelデータを拾ってくる。
すべて | 統計データを探す | 政府統計の総合窓口

週ごとにアンケートを取ってくれているようなので、少なくとも週に1件以上は全身麻酔をしているデータを対象としようと思ったが、適当に年間100件以上している施設だけ解析に使った。1989件あった。
さて、各施設 $H_j$ が、ある週 $t$ においてアンケート回答結果が $d=\{1,2,3,4\}$ だったとき、 $p$ を手術室運営%とすると
$d=1$ ならば通常通り100%の手術室運営なので、 $p \sim U(0.95, 1.05)$
$d=2$ ならば、やや制限（半分以上として70%運用）として、 $p\sim \beta(0.7)$
$d=3$ ならば、制限（半分以下として40%運用）として、 $p\sim \beta(0.4)$
$d=4$ ならば、予定手術はしない、緊急手術のみ対応として、 $p\sim \beta(0.12)$
とサンプリングして、前回実績数に $p$ をかけて手術件数を減ずることにする。
結果、こんな感じになった。
最初の頃に比べて、手術件数は2/3になっているようである。
f:id:MikuHatsune:20200426231152p:plain
これの推定の問題として、最初に例えばやや制限（半分以上として70%運用）と答えたとして、次の週以降はやや制限、制限、予定はしない、のいずれかしか答えないだろうが、通常通り、も答えうる、というのがよろしくないので、もう少しサンプリングの制限を加えた推定というかコードにしたいが、とりあえずここまで。

qdate <- c("03-13", "03-20", "03-27", "04-06", "04-10", "04-17", "04-27")
dat <- rbind(
c(89.3, 5.4, 3.6, 1.8),
c(76.6, 17, 2.1, 4.3),
c(82.2, 6.8, 4.1, 6.8),
c(64.4, 23.1, 6.9, 5.6),
c(45.7, 38.5, 13.1, 2.7),
c(26.7, 37.5, 26.7, 9.1),
c(24.7, 37.8, 29.5, 8)
)/100
dat <- dat/rowSums(dat)
colnames(dat) <- c("通常通り", "少し制限（半分以上）", "大きく制限（半分以下）", "原則予定手術なし")

N <- c(113, 47, 73, 160, 221, 232, 288)

datN <- round(dat * N)

# GA のデータは、エクセルの全身麻酔、というところだけあればよい
# 列名をGAにしている
# 下から5行は DPC対象病院群ということになっているので削除
GA <- read.csv("generalanesth.csv")
GAN <- na.omit(head(GA$GA, -5))
GAN12 <- round(GAN[GAN > 100]/12)

m <- c(0.7, 0.4, 0.12)
v <- 0.01
b <- cbind(m^2*(1-m)/v-m, (1-m)*(m*(1-m)/v-1))
x <- seq(0, 1, length=300)
d <- mapply(function(z) dbeta(x, b[z, 1], b[z, 2]), seq(m))
cols <- c("green", "orange", "red")
par(mar=c(5, 5, 2, 2), las=1, cex.lab=1.5)
matplot(x*100, d, type="l", lwd=3, lty=1, xlab="平時と比べて手術件数の減少割合[%]", ylab="確率密度", col=cols)
legend("topright", legend=sprintf("平均%d%sの減少", m*100, "%"), col=cols, pch=15, cex=1.5)


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

code <- "
data {
  int<lower=1> T; // 時間の数
  matrix<lower=0, upper=1>[T, 4] M;
  int<lower=0> N_hospital; // 全身麻酔をしている病院の数
  int<lower=0> Ope[N_hospital]; // 各病院の全身麻酔件数
  int<lower=0> N[T, 4]; // アンケート回答数
  matrix<lower=0>[3, 2] b; // beta distribution
}
parameters {
  simplex[4] y[T];
  real<lower=0> s[4];
}
model {
  for(i in 3:T){
    for(j in 1:4){
      y[i][j] ~ normal(2*y[i-1][j] - y[i-2][j], s[j]);
    }
  }
  for(i in 1:T){
    N[i,] ~ multinomial(y[i]);
  }
}
generated quantities{
  matrix[T, N_hospital] R;
  int<lower=0> idx;
  for(i in 1:T){
    for(j in 1:N_hospital){
      idx = categorical_rng(y[i]);
      // multinomial_rng で整数をサンプリングしたいが、
      // matrix がrealしかもたないのでこうした
      if(idx == 1){
        R[i, j] = Ope[j] * uniform_rng(0.95, 1.05);
      } else if(idx == 2){
        R[i, j] = Ope[j] * beta_rng(b[1, 1], b[1, 2]);
      } else if(idx == 3){
        R[i, j] = Ope[j] * beta_rng(b[2, 1], b[2, 2]);
      } else {
        R[i, j] = Ope[j] * beta_rng(b[3, 1], b[3, 2]);
      }
    }
  }
}
"

m0 <- stan_model(model_code=code)
standata <- list(T=nrow(dat), M=dat, N_hospital=length(GAN12), Ope=GAN12, N=datN, b=b)
fit <- sampling(m0, standata, iter=2000, warmup=1000, 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(4), SIMPLIFY=FALSE), along=3)

# 透明にしながら実線も描けるようにする
cols256 <- list("green", "yellow", "orange", "red")
cols256 <- lapply(lapply(cols256, col2rgb), c)
cols <- mapply(function(z) rgb(z[1], z[2], z[3], maxColorValue=256), cols256)
colsa <- mapply(function(z) rgb(z[1], z[2], z[3], 20, maxColorValue=256), cols256)

yl <- c(0, 1)
par(mar=c(5, 5, 2, 2), cex.lab=1.5, las=1)
matplot(m[2,,], type="n", col=colsa, lty=1, lwd=3, ylim=yl, xaxt="n", xlab="週", ylab="アンケート回答割合")
axis(1, at=seq(qdate), labels=gsub("-", "/", qdate))
legend("topright", legend=colnames(dat), ncol=1, col=cols, pch=15, yjust=-0.1, xpd=TRUE, bty="n", cex=1.4)
for(i in seq(dim(m)[3])){
  lines(c(seq(N)), m[2,,i], col=cols[i], lwd=4)
  polygon(c(c(seq(N)), rev(c(seq(N)))), c(m[1,,i], rev(m[3,,i])), col=colsa[i], border=NA)
}


library(vioplot)
Z <- as.data.frame(t(mapply(function(z) rowSums(ex$R[z,,]), seq(dim(ex$R)[1]))))
Z <- Z/4
par(mar=c(5, 5, 2, 2), cex.lab=1.5, las=0)
plot(c(0.5, length(qdate)+0.5), c(0, 1), type="n", xaxt="n", ylim=c(0, max(Z)), xlab="週", ylab="全施設での全身麻酔件数 [週あたり]")
axis(1, at=seq(qdate), labels=gsub("-", "/", qdate))
vioplot(Z, add=TRUE)
text(seq(qdate), apply(Z, 2, min), sprintf("%.1f%s", colMeans(Z)/colMeans(Z)[1]*100, "%"), pos=1)