## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "HP-"
)

## ----setup, message=FALSE, warning=FALSE--------------------------------------
library(cophescan)
library(dplyr)
library(ggplot2)

## -----------------------------------------------------------------------------
data("cophe_multi_trait_data")
trait_dat = cophe_multi_trait_data$summ_stat$Trait_1
str(trait_dat)
querysnpid <- cophe_multi_trait_data$querysnpid
LD <- cophe_multi_trait_data$LD

## ----message=FALSE, warning=FALSE, results='hide'-----------------------------
## Hide print messages from coloc
res.multi.lbf <- list()
for (trait_idx in seq_along(cophe_multi_trait_data$summ_stat)){
  querytrait_ss <- cophe_multi_trait_data$summ_stat[[trait_idx]]
  # Here LD is the same
  querytrait_ss$LD <- LD
  trait_variant_pair <- paste0('Trait', trait_idx, '_', querysnpid)
  res.multi.lbf[[trait_variant_pair]] <- cophe.susie.lbf(querytrait_ss, querysnpid = querysnpid, querytrait = paste0('Trait_', trait_idx))
}

res.multi.lbf.df = bind_rows(res.multi.lbf)


## -----------------------------------------------------------------------------
head(res.multi.lbf.df)

## ----warning=F, message=F, fig.width = 8, fig.height=8, out.width = "75%",  fig.ncol = 3, fig.align = "center"----
covar=FALSE
covar_vec=rep(1, nrow(res.multi.lbf.df))
## Set covar to TRUE to include covariates, uncomment the following 2 lines 
# covar=TRUE
# covar_vec = cophe_multi_trait_data$covar_vec
cophe.hier.res <- run_metrop_priors(res.multi.lbf.df, avg_posterior=TRUE, avg_pik = TRUE, covar_vec = covar_vec, covar = covar, nits = 50000, thin = 5)
names(cophe.hier.res)

## ----message=FALSE------------------------------------------------------------
cophe.hier.res.chain.list <- lapply(1:4, function(x) 
  run_metrop_priors(res.multi.lbf.df, avg_posterior=TRUE, avg_pik = TRUE, 
                    covar_vec = covar_vec, covar = covar, nits = 50000, thin = 5))

## ----mcmcDiagnostics, warning=FALSE, message=FALSE, fig.width=5, fig.height=5, out.width="80%", fig.align="center"----
# Store user parameters
old_par <- par(no.readonly = TRUE)

# chain_colors <- c("#e63946c4", "#f1faee", "#a8dadc", "#457b9d" )
chain_colors <- c("#f4f1de", "#e07a5fb2", "#3d405bb2", "#81b29aa6")

layout(matrix(c(1, 2, 3, 4, 5, 5), ncol=2, byrow = TRUE), respect = TRUE, 
       heights = c(0.9, 0.9, 0.1))


matplot(sapply(cophe.hier.res.chain.list, function(x) x$ll), type = "l", 
        col = chain_colors, 
     main ="loglik", ylab = "ll", xlab = "Iteration", lty = 1)

y_ax <- c("alpha", "beta", "gamma")
num_pars <- ifelse(covar, 3, 2) 
for (idx in 1:num_pars) {
    matplot(sapply(cophe.hier.res.chain.list, function(x) x$parameters[idx, ]),
        type = "l", col = chain_colors,
        main = paste(y_ax[idx]), ylab = y_ax[idx], xlab = "Iteration", lty = 1
    )
}

if (!covar) {
    plot(1, type = "n", axes = FALSE, xlab = "", ylab = "")
}

par(mar=c(0, 0, 0, 0))
plot(1, type = "n", axes = FALSE, xlab = "", ylab = "")
legend("top", legend = paste("Chain", 1:4), col = chain_colors, lty = 1, lwd = 2, 
       horiz = TRUE, bty = "n")

# Reset user parameters
par(old_par)


## -----------------------------------------------------------------------------
res.post.prob = cbind(cophe.hier.res.chain.list[[1]]$avg.posterior, cophe.hier.res$data)

## -----------------------------------------------------------------------------
res.hier.predict <- cophe.hyp.predict(as.data.frame(res.post.prob ))
col_disp <- c( "PP.Hn", "PP.Ha", "PP.Hc", "nsnps", "querysnp", "querytrait",
               "typeBF",  "grp", "cophe.hyp.call")
knitr::kable(res.hier.predict[, col_disp], row.names = FALSE, digits=3)

## ----cophePlots, message=FALSE------------------------------------------------
res.plots = cophe_plot(res.hier.predict, traits.dat = cophe_multi_trait_data$summ_stat, querysnpid = querysnpid, query_trait_names = paste0('Trait_', 1:24))

# if (!require(ggpubr)) {
  # install.packages("ggpubr") 
# }
# ggpubr::ggarrange(res.plots$pval, res.plots$ppHa, res.plots$ppHc, ncol = 2, 
#                 nrow = 2)

## ----warning=FALSE, message=FALSE, echo=FALSE, out.width = "40%"--------------
res.plots$pval + theme(legend.position="bottom")
res.plots$ppHa + theme(legend.position="bottom")
res.plots$ppHc + theme(legend.position="bottom")