新型肺炎の年齢階級別死亡率

新型コロナウイルス 国内感染の状況
こちらの2021年1月13日時点の各年代別感染者数と死亡者数のデータがあったので単純に推定してみる。

その1
drcパッケージにあるLL.4関数は、4パラメータモデルのロジスティック曲線に当てはめる。
これで簡単にやる

その2
rstan を使う。年齢を経るごとに死亡率が上がるとして、orderedで昇順の値のサンプリングをして、inv_logitにより0~1 の値になるようにする。すなわちこれが死亡率になり、やっていることはロジスティック曲線への変換である。

死亡率は以下のようになる(%)

      age             
 [1,]   5  0.007447185
 [2,]  15  0.007450361
 [3,]  25  0.007735825
 [4,]  35  0.013074665
 [5,]  45  0.059155769
 [6,]  55  0.310493872
 [7,]  65  1.314663597
 [8,]  75  4.443249524
 [9,]  85 12.025304864

rstan でやるよりかは、推定曲線へのフィッティングのほうが、predict関数でinterval="confidence"を指定すると区間推定値が出るので有用だと思う。
例えば100歳での死亡率は95%信頼区間で

     Lower Prediction      Upper 
  30.62692   33.85670   37.08648 

となる。

# 1/13
death <- cbind(c(0, 0, 2, 11, 35, 113, 328, 944, 2401),
               c(7160, 18377, 69311, 46455, 42939, 39148, 24823, 21243, 19967))
rownames(death) <- c("10代以下", sprintf("%d代", seq(10, 70, by=10)), "80代以上")

age <- seq(5, 85, by=10)
p <- death[,1]/death[,2]



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

code <- "
data{
  int N;
  int death[N];
  int positive[N];
}
parameters{
  ordered[N] theta;
}
model{
  for(i in 1:N){
    death[i] ~ binomial(positive[i], inv_logit(theta[i]));
  }
}
"

m0 <- stan_model(model_code=code)
standata <- list(N=nrow(death), death=death[,1], positive=death[,2])
fit <- sampling(m0, data=standata, iter=3000, warmup=1500, seed=1234)

ex <- extract(fit)
r <- 1/(1+exp(-ex$theta))

d <- drm(p ~ age, fct=LL.4())
x <- seq(0, 100, by=1)
y <- predict(d, as.data.frame(x), interval="confidence")

xl <- range(x)
yl <- range(0, r, y)

plot(d, log="", xlim=xl, ylim=yl, xlab="年齢", ylab="陽性者死亡率", bty="n", axes=FALSE)
lines(x, y[,2], col=2)
lines(x, y[,3], col=2)
vioplot(r, add=TRUE, at=age, colMed="yellow", wex=5)
axis(1, lwd=0, lwd.ticks=1, at=age, labels=age)
axis(2, lwd=0, lwd.ticks=1)