COVID-19: estimate total number of positive persons in Japan

data_20200402.csv

DP.csv

fit-summary.csv

model.bugs

model {
  for (a in 1:A) {
    n.pos[a] ~ dbinom(p.pos, N[a])
    n.sym[a] ~ dbinom(q[a], n.pos[a])
    n.asy[a] <- n.pos[a] - n.sym[a]
    Y.sym[a] ~ dbinom(p.obs.sym, n.sym[a])
    Y.asy[a] ~ dbinom(p.obs.asy, n.asy[a])
    Y.sym.DP[a] ~ dbinom(q[a], N.pos.DP[a])
    logit(q[a]) <- q.x[a]
  }
  n.pos.sum <- sum(n.pos)
  p.pos     ~ dbeta(0.5, 1)
  p.obs.sym ~ dbeta(4, 2)
  p.obs.asy ~ dbeta(2, 4)
  q.x[1] ~ dnorm(0, 1e-4)
  q.x[2] ~ dnorm(0, 1e-4)
  for (a in 3:A) {
    q.x[a] ~ dnorm(2*q.x[a-1] - q.x[a-2], tau)
  }
  tau <- 1 / (sigma*sigma)
  sigma ~ dunif(0, 5)
}

pop.csv

run-model.R

library(dplyr)
library(rjags)
library(R2WinBUGS)

A <- 10
d_pop <- read.csv('pop.csv', stringsAsFactors = FALSE) %>%
  mutate(Age10 = rep(1:10, each=2)) %>%
  group_by(Age10) %>%
  summarise(Value = sum(Value)) %>% ungroup()
d_DP <- read.csv('DP.csv')
d <- read.csv('data_20200402.csv') %>% mutate(Y.sym.new = Y.sym + Y.unk)

data <- list(A = A,
             N = d_pop$Value*1000,
             Y.sym = d$Y.sym.new,
             Y.asy = d$Y.asy,
             N.pos.DP = d_DP$N,
             Y.sym.DP = d_DP$Y)
inits <- lapply(1:4, function(i) list(p.pos = 0.0001, q.x = rep(0,A), p.obs.sym = 0.5, p.obs.asy = 0.1,
                                      .RNG.name = 'base::Wichmann-Hill', .RNG.seed = i))
m <- jags.model('model.bugs', data, inits, n.chains=4, n.adapt=50000)
update(m, 50000)
post.list <- coda.samples(m, c('n.pos.sum', 'n.pos', 'p.pos', 'p.obs.sym', 'p.obs.asy', 'q'), n.iter=500000, thin=100)

mcmc.list2bugs <- function(mcmc.list){
  b1 <- mcmc.list[[1]]
  m1 <- as.matrix(b1)
  mall <- matrix(numeric(0), 0, ncol(m1))
  n.chains <- length(mcmc.list)
  for (i in 1:n.chains) {
    mall <- rbind(mall, as.matrix(mcmc.list[[i]]))
  }
  sims.array <- array(mall, dim = c(nrow(m1), n.chains, ncol(m1)))
  dimnames(sims.array) <- list(NULL, NULL, colnames(m1))
  mcpar <- attr(b1, "mcpar")
  as.bugs.array(
    sims.array = sims.array,
    model.file = NULL,
    program = NULL,
    DIC = FALSE,
    DICOutput = NULL,
    n.iter = mcpar[2],
    n.burnin = mcpar[1] - mcpar[3],
    n.thin = mcpar[3]
  )
}
post.bugs <- mcmc.list2bugs(post.list)
options(scipen = 10)
write.table(data.frame(post.bugs$summary, check.names = FALSE),
            file = 'fit-summary.csv', sep = ',', quote = TRUE, col.names = NA)
# save.image(file = 'result.RData')


pdf('traceplot.pdf')
plot(post.list)
dev.off()

library(ggplot2)
post.mcmc <- as.mcmc(post.bugs$sims.matrix)
d.plot <- data.frame(n.pos.sum = as.vector(post.mcmc[,'n.pos.sum']))

p <- ggplot(d.plot, aes(x=n.pos.sum)) +
  theme(text=element_text(size=18)) +
  geom_histogram(aes(y=..density..), binwidth=500, fill='white', color='black') +
  geom_density(fill='black', alpha=0.3) +
  xlim(NA, 20000) +
  labs(x='推定総陽性者数 in 日本 (4/2時点)', y='確率密度')
ggsave(p, file = 'fig-positive-persons-total.png', dpi=300, w=7, h=5)

age10 <- c('0-9', '10-19', '20-29', '30-39', '40-49', '50-59', '60-69', '70-79', '80-89', '90-99')
qua <- apply(post.mcmc[, sprintf('n.pos[%d]', 1:A)], 2, quantile, probs=c(0.025, 0.5, 0.975)) %>% t()
d.qua <- data.frame(Age10=factor(age10, levels=age10), qua, check.names = FALSE)

p <- ggplot() +
  theme(text=element_text(size=18), axis.text.x=element_text(angle=40, vjust=1, hjust=1)) +
  geom_line(data=d.qua, aes(x=Age10, y=`50%`), group=1) +
  geom_point(data=d.qua, aes(x=Age10, y=`50%`), size=2) +
  geom_ribbon(data=d.qua, aes(x=Age10, ymin=`2.5%`, ymax=`97.5%`), group=1, fill='black', alpha=0.3) +
  labs(x='年代', y='推定陽性者数 in 日本 (4/2時点)')
ggsave(p, file = 'fig-positive-persons-by-age.png', dpi=300, w=7, h=5)