library(ggplot2)
library(EpiGLM) ##For data(estData)

load("estBoot-post.RData")

hatnames <- gsub("theta_", "hat_", parnames)
hats <- list()
hats$tr <- hatnames[grep("^hat_ctr", hatnames)]
hats$gr <- hatnames[grep("^hat_cgr", hatnames)]
hats$ip <- hatnames[grep("^hat_cgr|^hat_ctr", hatnames, invert=TRUE)]

cnames <- c(resL[[1]]$coefnames,
            list(mmat_ip=rownames(resL[[1]]$ip_fit$coefs)))
## Check that all models are the one that is assumed
terms <- c("(Intercept)", "facility2", "is_smallTRUE", "seasonspring-summer", 
"CaseTypestaff")
stopifnot(all(sapply(cnames, function(x) all(x == terms))))

hats$int <- hatnames[grep("r.1$|\\.Intercept\\.$", hatnames)]
hats$facility2 <- hatnames[grep("r.2$|facility2$", hatnames)]
hats$is_small <- hatnames[grep("r.3$|is_smallTRUE$", hatnames)]
hats$spring_summer <- hatnames[grep("r.4$|seasonspring.summer$", hatnames)]
hats$staff <- hatnames[grep("r.5$|CaseTypestaff$", hatnames)]

hatdat <- resdf[, unlist(hatnames)]
hatdat <- stack(hatdat)
hatdat$linear_pred <- hatdat$term <- NA

predlabs <- c("transmission", "growth_rate", "symptomatic_period")

test <- hatdat$ind %in% hats$tr
hatdat$linear_pred[test] <- predlabs[1]
test <- hatdat$ind %in% hats$gr
hatdat$linear_pred[test] <- predlabs[2]
test <- hatdat$ind %in% hats$ip
hatdat$linear_pred[test] <- predlabs[3]

termlabs <- c("intercept", "facility2", "is_small", "spring_summer", "staff")
termlabs[-1] <- c("nursing home", "small", "spring-summer", "staff")
test <- hatdat$ind %in% hats$int
hatdat$term[test] <- termlabs[1]
test <- hatdat$ind %in% hats$facility2
hatdat$term[test] <- termlabs[2]
test <- hatdat$ind %in% hats$is_small
hatdat$term[test] <- termlabs[3]
test <- hatdat$ind %in% hats$spring_summer
hatdat$term[test] <- termlabs[4]
test <- hatdat$ind %in% hats$staff
hatdat$term[test] <- termlabs[5]

table(is.na(hatdat$values), hatdat$linear_pred)
hatdat <- hatdat[!is.na(hatdat$values), ]
stopifnot(!any(is.na(hatdat)))

melt_wrapper <- function(x, fun=identity, ...) {
  names(x) <- predlabs
  x <- lapply(x, fun, ...)
  x <- as.data.frame(x)
  x <- cbind(x, termlabs)
  x <- reshape::melt.data.frame(x, id="termlabs")
  names(x)[1:2] <- c("term", "linear_pred")
  x
}
data(estData)
mle <- c(estData$fit$admb$coeflist, list(ip=coef(estData$ip_fit$m)))
mle <- melt_wrapper(mle)
tmpf <- function(x, suf) {
  colname <- paste(x, suf, sep="_")
  ret <- res_tab[, colname]
  as.numeric(ret)
}  
foo <- hats[c("tr", "gr", "ip")]
# Because res_tab has dumb names, at time of writing
names(res_tab) <- gsub("^theta", "hat", names(res_tab))
bootLower <- melt_wrapper(foo, fun=tmpf, suf="bootLower")
bootUpper <- melt_wrapper(foo, fun=tmpf, suf="bootUpper")
rm(foo, tmpf, melt_wrapper)
bootLimit <- bootLower
bootLimit$value <- NULL
bootLimit$lower <- bootLower$value
bootLimit$upper <- bootUpper$value

spruce_labels <- function(x) {
  x$linear_pred <- factor(x$linear_pred,
    levels=c("transmission", "growth_rate", "symptomatic_period"),
    labels=c("transmission rate", "growth rate", "symptomatic period"))
  x
}
mle <- spruce_labels(mle)
hatdat <- spruce_labels(hatdat)
bootLimit <- spruce_labels(bootLimit)

test <- hatdat$term != "intercept"
hatnoint <- hatdat[test, ]
test <- mle$term != "intercept"
mlenoint <- mle[test, ]
test <- bootLimit$term != "intercept"
bootLimitnoint <- bootLimit[test, ]

the_breaks <- c(0.25, 0.5, 1, 2, 4, 8, 16)
g <- ggplot()
g <- g + geom_rect(data=bootLimitnoint,
                   aes(ymin=-Inf, ymax=Inf, xmin=lower, xmax=upper),
                   fill="#CCCCCC")
g <- g + geom_vline(data=mlenoint, aes(xintercept=value),
                    color="white")
for (i in unique(hatdat$linear_pred)){
  test <- hatnoint$linear_pred == i
  dat <- hatnoint[test, ]
  for (j in unique(hatdat$term)[-1]){
    test <- dat$term == j
    dat2 <- dat[test, ]
    binwidth <- diff(range(dat2$values)) / 30
    g <- g + geom_histogram(data=dat2, aes(x=values), binwidth=binwidth)
  }
}
g <- g + scale_x_continuous(breaks=log(the_breaks), labels=the_breaks)
g <- g + facet_grid(linear_pred~term)
g <- g + coord_flip()
g <- g + geom_vline(xintercept=0)
g <- g + theme_bw() 
g <- g + theme(panel.grid.minor=element_blank())
g <- g + theme(panel.grid.major=element_blank())
g <- g + theme(panel.background=element_blank())
g <- g + theme(strip.background=element_blank())
g <- g + theme(panel.border=element_blank())
g <- g + theme(axis.line=element_line())
g <- g + theme(strip.text.y=element_text())
g <- g + labs(x="Multiplicative effect", y="Count")
ggsave(g, filename="fig1.eps", width=190/25.4)
