Last active
April 4, 2020 22:50
-
-
Save MatsuuraKentaro/b859b1036afe11d01fe1d001f3daf79f to your computer and use it in GitHub Desktop.
COVID-19: estimate total number of positive persons in Japan
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Age10 | Y.sym | Y.asy | Y.unk | |
---|---|---|---|---|
1 | 19 | 6 | 7 | |
2 | 20 | 12 | 18 | |
3 | 256 | 30 | 81 | |
4 | 212 | 23 | 104 | |
5 | 296 | 32 | 94 | |
6 | 334 | 20 | 58 | |
7 | 280 | 21 | 45 | |
8 | 212 | 19 | 55 | |
9 | 114 | 17 | 21 | |
10 | 20 | 6 | 4 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Age10 | N | Y | |
---|---|---|---|
1 | 1 | 0 | |
2 | 5 | 2 | |
3 | 28 | 25 | |
4 | 34 | 27 | |
5 | 27 | 19 | |
6 | 59 | 28 | |
7 | 177 | 76 | |
8 | 234 | 95 | |
9 | 52 | 27 | |
10 | 2 | 2 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
mean | sd | 2.5% | 25% | 50% | 75% | 97.5% | Rhat | n.eff | ||
---|---|---|---|---|---|---|---|---|---|---|
n.pos[1] | 571.3027 | 159.276253154546 | 389 | 463 | 528 | 633.24985203207 | 995 | 1.00273359464946 | 1400 | |
n.pos[2] | 646.54735 | 180.40521489152 | 442 | 524 | 597 | 718 | 1128 | 1.00280045249263 | 1300 | |
n.pos[3] | 788.81355 | 202.366841412625 | 563 | 651 | 732 | 868 | 1328 | 1.00266226760728 | 1400 | |
n.pos[4] | 847.64955 | 229.340902364302 | 592 | 691 | 783 | 938 | 1459 | 1.00267799682899 | 1400 | |
n.pos[5] | 1112.34675 | 296.987740993383 | 782 | 909 | 1029 | 1230 | 1905.02499360457 | 1.00270136266963 | 1400 | |
n.pos[6] | 988.9982 | 260.827368660906 | 699 | 811 | 915 | 1089 | 1690 | 1.00265209090134 | 1400 | |
n.pos[7] | 979.2296 | 259.509483916199 | 690 | 801 | 906 | 1081 | 1677 | 1.00278489415782 | 1300 | |
n.pos[8] | 943.2969 | 254.664178561851 | 660 | 769 | 871 | 1042 | 1618.02499247061 | 1.00267407148837 | 1400 | |
n.pos[9] | 528.11235 | 144.136180709102 | 364 | 430 | 489 | 585 | 911 | 1.00255621322743 | 1500 | |
n.pos[10] | 133.87285 | 37.5079554923187 | 88 | 108 | 125 | 149 | 234 | 1.00247803364773 | 1600 | |
n.pos.sum | 7540.1698 | 2010.14756242647 | 5345 | 6159 | 6963.99999999999 | 8319.99999999999 | 12933.1499660855 | 1.002724989488 | 1400 | |
p.obs.asy | 0.0462820189650233 | 0.0106371673964274 | 0.0249479044213861 | 0.0390468723561322 | 0.0467567974872957 | 0.0538146206566749 | 0.065973207967604 | 1.00256260563806 | 1500 | |
p.obs.sym | 0.731551421252528 | 0.149115496892973 | 0.406139527877314 | 0.631397971320101 | 0.753212794884425 | 0.849299316258468 | 0.958274541057583 | 1.00260743556801 | 1500 | |
p.pos | 0.0000597917406823182 | 0.0000159545997699604 | 0.0000423301153237783 | 0.0000488289060364238 | 0.0000552550203854312 | 0.0000660167601435308 | 0.000102447045958762 | 1.00275434384743 | 1400 | |
q[1] | 0.0638171116737871 | 0.0130368120212076 | 0.0409365927091878 | 0.0547055590843862 | 0.0629007434340131 | 0.0718747862452483 | 0.0922218566108755 | 1.00094340668855 | 20000 | |
q[2] | 0.0978514667313718 | 0.0157685607362956 | 0.0694697191165102 | 0.0868346238082514 | 0.0969061499115658 | 0.107940089403318 | 0.130931353227933 | 1.00127138890485 | 6700 | |
q[3] | 0.612663477703154 | 0.0360474418781044 | 0.542908700445883 | 0.588298893426953 | 0.612671574930056 | 0.636865874972209 | 0.683919508331088 | 1.00112138438452 | 11000 | |
q[4] | 0.556647738516761 | 0.0359573121573898 | 0.48782819876908 | 0.531921618190929 | 0.555981533417973 | 0.580743698485345 | 0.627679011526143 | 1.00094499280894 | 20000 | |
q[5] | 0.512964933900214 | 0.0322005449037579 | 0.451576030784536 | 0.491012656967577 | 0.512429210450138 | 0.534339023243265 | 0.578394497353506 | 1.00100000446474 | 20000 | |
q[6] | 0.559810425760402 | 0.0318546207652959 | 0.497800941613316 | 0.537972499282999 | 0.559783137745123 | 0.581225167956847 | 0.62299267042623 | 1.00100756756489 | 20000 | |
q[7] | 0.467753253141322 | 0.0253789012920255 | 0.418286653630966 | 0.450526451046241 | 0.467581382402763 | 0.484883620731577 | 0.517193106302658 | 1.00099492000193 | 20000 | |
q[8] | 0.409483681415604 | 0.0228188916634046 | 0.365658682907407 | 0.394176123092705 | 0.409003522101692 | 0.424840495650802 | 0.455343825123253 | 1.00118551930463 | 8700 | |
q[9] | 0.388747519361643 | 0.0321716946415977 | 0.328574619722074 | 0.366525406414358 | 0.387432313404832 | 0.409703518073575 | 0.454240612835382 | 1.00107984629458 | 14000 | |
q[10] | 0.279136724983478 | 0.0540495626286248 | 0.182315210175756 | 0.241238547979342 | 0.275819744497811 | 0.313424045532416 | 0.395513301469468 | 1.00093811682452 | 20000 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Age5 | Value | |
---|---|---|
0-4 | 4758 | |
5-9 | 5101 | |
10-14 | 5351 | |
15-19 | 5820 | |
20-24 | 6388 | |
25-29 | 6240 | |
30-34 | 6752 | |
35-39 | 7551 | |
40-44 | 8718 | |
45-49 | 9802 | |
50-54 | 8567 | |
55-59 | 7711 | |
60-64 | 7523 | |
65-69 | 8709 | |
70-74 | 8686 | |
75-79 | 7241 | |
80-84 | 5328 | |
85-89 | 3612 | |
90-94 | 1761 | |
95-99 | 479 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment