#############################
## Forecasting fit analysis
#############################
library(zikaInfer)
library(ggplot2)
library(extrafont)
library(cowplot)
library(hexbin)
library(RColorBrewer)
chainWD <- "~/net/home/zika/outputs/"
# Heatmaps ----------------------------------------------------------------
heatmap_theme <- theme(axis.text.x=element_text(size=8,family="Arial",color="black"),
plot.title=element_text(size=8,family="Arial",color="black"),
axis.text.y=element_text(size=8,family="Arial",color="black"),
axis.title.x=element_text(size=10,family="Arial",color="black"),
axis.title.y=element_text(size=10,family="Arial",color="black"),
legend.text=element_text(size=8,family="Arial",color="black"),
legend.title=element_text(size=8,family="Arial",color="black"),
strip.text=element_text(size=8,family="Arial")) + theme_bw()
# Reading in data ---------------------------------------------------------
bahia_chain <- zikaInfer::load_mcmc_chains(location = paste0(chainWD,"bahia_forecast"),
asList = FALSE,convertMCMC = FALSE,unfixed = FALSE,thin = 1,burnin = 50000)
nejm_chain <- zikaInfer::load_mcmc_chains(location = paste0(chainWD,"northeast_forecast"),
asList = FALSE,convertMCMC = FALSE,unfixed = FALSE,thin = 1,burnin = 50000)
bahia_chain$propn_increase <- bahia_chain$incPropn2/bahia_chain$incPropn
nejm_chain$propn_increase <- nejm_chain$incPropn2/nejm_chain$incPropn
bahia_dat<- make_micro_dat(bahia_chain,10000)
bahia_dat$Model <- "Bahia, Brazil (reports)"
nejm_dat<- make_micro_dat(nejm_chain,10000)
nejm_dat$Model <- "Northeast Brazil"
dat <- rbind(bahia_dat, nejm_dat)
dat$Model <- as.factor(dat$Model)
# Microceph risk curves ---------------------------------------------------
p <- ggplot(data=dat) +
geom_vline(xintercept=c(14*7, 28*7),col="grey40",lty="dashed",size=1) +
geom_ribbon(aes(x=weeks,ymax=upper,ymin=lower,fill=Model,group=Model),alpha=0.6) +
geom_line(aes(x=weeks,y=means,group=Model),col="black",lty="dashed") +
facet_wrap(~Model,ncol=1) +
ylab("Relative probability of \nmicrocephaly given ZIKV infection")+
xlab("Week of gestation at time of infection") +
theme_bw() +
theme(axis.text.x=element_text(size=8,family="Arial",color="black"),
panel.grid.minor=element_blank(),
plot.title=element_text(size=8,family="Arial",color="black"),
axis.text.y=element_text(size=8,family="Arial",color="black"),
axis.title.x=element_text(size=10,family="Arial",color="black"),
axis.title.y=element_text(size=10,family="Arial",color="black"),
legend.position="bottom",
legend.title=element_blank(),
legend.text=element_text(size=8,family="Arial",color="black"),
legend.justification=c(0,0)) +
scale_y_continuous(expand=c(0,0),limits=c(0,0.05)) +
scale_x_continuous(expand=c(0,0),limits=c(0,40*7),breaks=seq(0,280,by=7*4), labels=seq(0,40,by=4))
p
# Densities
bahia_chain$version <- "Bahia, Brazil (reports)"
nejm_chain$version <- "Northeast Brazil"
b_chain <- reshape2::melt(bahia_chain, id.vars="version")
n_chain <- reshape2::melt(nejm_chain, id.vars="version")
chains <- rbind(b_chain,n_chain)
chains$version <- as.factor(chains$version)
use_vars <- c("birth_reduction","propn_increase","incPropn","incPropn2","abortion_rate")
tmp_chain <- chains[chains$variable %in% use_vars,]
fake_dat <- data.frame(version=c(rep("Bahia - reports",12),
rep("Northeast Brazil NEJM", 12)),
variable=rep(c("propn2","propn2","birth_reduction","birth_reduction","propn_increase",
"propn_increase","incPropn","incPropn","incPropn2","incPropn2",
"abortion_rate","abortion_rate"),2),
value=rep(c(0,2,0,1,0,500,0,0.01,0,1,0,1),2)
)
fake_dat <- fake_dat[fake_dat$variable %in% use_vars,]
convert_names <- c("propn2"="Microcephaly reporting rate 2016",
"birth_reduction"="Proportion reduction in affected births",
"propn_increase"="Increase in ZIKV reporting rate",
"incPropn"="ZIKV incidence reporting rate 2015",
"incPropn2"="ZIKV incidence reporting rate 2016",
"abortion_rate"="Proportion of microcephaly affected births aborted")
tmp_chain$variable <- convert_names[as.character(tmp_chain$variable)]
fake_dat$variable <- convert_names[as.character(fake_dat$variable)]
use_vars2 <- c("birth_reduction","propn_increase","abortion_rate","propn")
convert_names <- c("propn"="Microcephaly reporting\nrate 2015/2016",
"birth_reduction"="Proportion reduction\nin affected births",
"propn_increase"="Increase in ZIKV reporting rate",
"incPropn"="ZIKV incidence reporting rate 2015",
"incPropn2"="ZIKV incidence reporting rate 2016",
"abortion_rate"="Proportion of microcephaly\naffected births aborted")
tmp_chain <- b_chain[b_chain$variable %in% use_vars2,]
tmp_chain$variable <- convert_names[as.character(tmp_chain$variable)]
quantiles <- reshape2::melt(plyr::ddply(tmp_chain, "variable",function(x) quantile(x$value,c(0.025,0.5,0.975))))
colnames(quantiles) <- c("variable","quantile","value")
densities_bahia <- ggplot(tmp_chain) +
geom_density(aes(x=value,y=..density..,fill=variable)) +
geom_vline(data=quantiles[quantiles$quantile %in% c("2.5%","97.5%"),],aes(xintercept=value),col="black",lty="dashed")+
geom_vline(data=quantiles[quantiles$quantile %in% c("50%"),],aes(xintercept=value),col="black")+
facet_wrap(~variable,scale="free") +
heatmap_theme +
ylab("Density") +
xlab("Value") +
theme_bw() +
theme(legend.position="none")
densities_bahia
tmp_chain <- n_chain[n_chain$variable %in% use_vars2,]
tmp_chain$variable <- convert_names[as.character(tmp_chain$variable)]
quantiles <- reshape2::melt(plyr::ddply(tmp_chain, "variable",function(x) quantile(x$value,c(0.025,0.5,0.975))))
colnames(quantiles) <- c("variable","quantile","value")
densities_nejm <- ggplot(tmp_chain) +
geom_density(aes(x=value,y=..density..,fill=variable)) +
geom_vline(data=quantiles[quantiles$quantile %in% c("2.5%","97.5%"),],aes(xintercept=value),col="black",lty="dashed")+
geom_vline(data=quantiles[quantiles$quantile %in% c("50%"),],aes(xintercept=value),col="black")+
facet_wrap(~variable,scale="free") +
heatmap_theme +
ylab("Density") +
xlab("Value") +
theme_bw() +
theme(legend.position="none")
densities_nejm
# Corrplots Bahia ---------------------------------------------------------------
heatmap_theme <- heatmap_theme + theme(
axis.text.x=element_text(size=8,family="Arial",color="black"),
axis.text.y=element_text(size=8,family="Arial",color="black"),
legend.text=element_text(size=6,family="Arial",color="black")
)
p1 <- ggplot() + geom_hex(data=bahia_chain,aes(x=birth_reduction,y=propn_increase,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.005) +
ylab("Relative increase in\n ZIKV reporting rate") +
xlab("Proportion decrease in\nZIKV affected births") +
labs(fill="Density") + heatmap_theme
p2 <- ggplot() + geom_hex(data=bahia_chain,aes(x=birth_reduction,y=propn,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.005) +
ylab("Microcephaly reporting\nrate in first wave") +
xlab("Proportion decrease in\nZIKV affected births") +
labs(fill="Density") + heatmap_theme
p3 <- ggplot() + geom_hex(data=bahia_chain,aes(x=birth_reduction,y=abortion_rate,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.005) +
ylab("Proportion of microcephaly\naffected births aborted") +
xlab("Proportion decrease in\nZIKV affected births") +
labs(fill="Density") + heatmap_theme
p4 <- ggplot() + geom_hex(data=bahia_chain,aes(x=propn_increase,y=propn,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.007) +
ylab("Microcephaly reporting\nrate in first wave") +
xlab("Relative increase in\n ZIKV reporting rate") +
labs(fill="Density") + heatmap_theme
p5 <- ggplot() + geom_hex(data=bahia_chain,aes(x=propn_increase,y=abortion_rate,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.006) +
ylab("Proportion of microcephaly\naffected births aborted") +
xlab("Relative increase in\nZIKV reporting rate") +
labs(fill="Density") + heatmap_theme
p6 <- ggplot() + geom_hex(data=bahia_chain,aes(x=propn,y=abortion_rate,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.006) +
ylab("Proportion of microcephaly\naffected births aborted") +
xlab("Microcephaly reporting\nrate in first wave") +
labs(fill="Density") + heatmap_theme
cor_bahia <- plot_grid(p1,p2,p3,p4,p5,p6,ncol=3)
cor_bahia
# Corrplots NEJM ----------------------------------------------------------
p1 <- ggplot() + geom_hex(data=nejm_chain,aes(x=birth_reduction,y=propn_increase,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.005) +
ylab("Relative increase in\n ZIKV reporting rate") +
xlab("Proportion decrease in\nZIKV affected births") +
labs(fill="Density") + heatmap_theme
p2 <- ggplot() + geom_hex(data=nejm_chain,aes(x=birth_reduction,y=propn,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.0035) +
ylab("Microcephaly reporting\nrate in first wave") +
xlab("Proportion decrease in\nZIKV affected births") +
labs(fill="Density") + heatmap_theme
p3 <- ggplot() + geom_hex(data=nejm_chain,aes(x=birth_reduction,y=abortion_rate,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.002) +
ylab("Proportion of microcephaly\naffected births aborted") +
xlab("Proportion decrease in\nZIKV affected births") +
labs(fill="Density") + heatmap_theme
p4 <- ggplot() + geom_hex(data=nejm_chain,aes(x=propn_increase,y=propn,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.0075) +
ylab("Microcephaly reporting\nrate in first wave") +
xlab("Relative increase in\n ZIKV reporting rate") +
labs(fill="Density") + heatmap_theme
p5 <- ggplot() + geom_hex(data=nejm_chain,aes(x=propn_increase,y=abortion_rate,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.003) +
ylab("Proportion of microcephaly\naffected births aborted") +
xlab("Relative increase in\nZIKV reporting rate") +
labs(fill="Density") + heatmap_theme
p6 <- ggplot() + geom_hex(data=nejm_chain,aes(x=propn,y=abortion_rate,fill=..density..))+
scale_fill_gradient2(low="#5E4FA2",mid="#FAFDB8",high="#9E0142",midpoint= 0.0025) +
ylab("Proportion of microcephaly\naffected births aborted") +
xlab("Microcephaly reporting\nrate in first wave") +
labs(fill="Density") + heatmap_theme
cor_nejm <- plot_grid(p1,p2,p3,p4,p5,p6,ncol=3)
cor_nejm
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.