## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  eval=FALSE
)

## ----warning=FALSE, message=FALSE---------------------------------------------
# library(RaCE.NMA)
# library(ggplot2)
# library(dplyr)
# library(cowplot)
# library(gridExtra)
# library(reshape2)

## ----fig.asp=0.35-------------------------------------------------------------
# s <- 0.1
# g1a<-ggplot()+scale_x_continuous(limits=c(-3*s,1+3*s))+
#   stat_function(fun = dnorm, args = list(mean = 0, sd =s))+
#   stat_function(fun = dnorm, args = list(mean = 1, sd =s))+
#   theme_bw()+labs(x=" ",y="Density",subtitle=expression(hat(sigma)~"=0.1"))+
#   theme(panel.grid.minor = element_blank(),panel.grid.major.y = element_blank())
# s <- 0.3
# g1b<-ggplot()+scale_x_continuous(limits=c(-3*s,1+3*s))+
#   stat_function(fun = dnorm, args = list(mean = 0, sd =s))+
#   stat_function(fun = dnorm, args = list(mean = 1, sd =s))+
#   theme_bw()+labs(x="Posterior Intervention Effects",y=NULL,subtitle=expression(hat(sigma)~"=0.3"))+
#   theme(panel.grid.minor = element_blank(),panel.grid.major.y = element_blank())
# s <- 0.5
# g1c<-ggplot()+scale_x_continuous(limits=c(-3*s,1+3*s))+
#   stat_function(fun = dnorm, args = list(mean = 0, sd =s))+
#   stat_function(fun = dnorm, args = list(mean = 1, sd =s))+
#   theme_bw()+labs(x=" ",y=NULL,subtitle=expression(hat(sigma)~"=0.5"))+
#   theme(panel.grid.minor = element_blank(),panel.grid.major.y = element_blank())
# grid.arrange(g1a,g1b,g1c,nrow=1)

## ----cache=TRUE---------------------------------------------------------------
# ## Perform simulation study
# set.seed(1)
# results <- matrix(NA,nrow=0,ncol=6)
# for(iter in 1:20){
#   for(J in c(6,12,18)){
#     for(K in c(J/3,2*J/3,J)){
#       for(s in c(0.1,0.3,0.5)){
#         if(K==J){
#           mu_hat <- 1:K
#         }else{
#           mu_hat <- sample(1:K,J,replace=T)
#           while(length(unique(mu_hat))<K){mu_hat <- sample(1:K,J,replace=T)}
#         }
# 
#         mcmc <- mcmc_raceNMA(mu_hat = mu_hat,s=rep(s,J),mu=mean(mu_hat),
#                              sigma0=max(1,4*sd(mu_hat)),nu0=mu_hat,
#                              tau=1,iter=3000,nu_iter=2,chains=4,
#                              burn_prop=0.5,thin=3,verbose = FALSE)
# 
#         equal <- posterior_equal <- matrix(NA,nrow=J,ncol=J)
#         for(i in 1:(J-1)){for(j in (i+1):J){
#           equal[i,j] <- ifelse(mu_hat[i]==mu_hat[j],1,0) # test if treatments are truly equal in mean
#           posterior_equal[i,j] <- mean( # assess if treatments are rank-clustered
#             mcmc[,paste0("G",i)] == mcmc[,paste0("G",j)]
#             )
#         }}
# 
#         results <- rbind(
#           results,
#           c(iter,J,K,s,
#             mean(posterior_equal[equal==1],na.rm=T),
#             mean(posterior_equal[equal==0],na.rm=T))
#         )
#       }
#     }
#   }
# }

## -----------------------------------------------------------------------------
# ## Plotting results from simulation study
# results <- as.data.frame(results)
# names(results) <- c("Iteration","J","K","s","Prob_Clustered","Prob_Distinct")
# results$s<-factor(results$s, levels=c(.1,.3,.5),
#                   labels=c(expression(hat(sigma)~"=0.1"),
#                            expression(hat(sigma)~"=0.3"),
#                            expression(hat(sigma)~"=0.5")))
# results$J <- factor(results$J,levels=c(6,12,18),
#                     labels=c(expression(J~"= 6"),
#                              expression(J~"= 12"),
#                              expression(J~"= 18")))
# results_melt <- melt(results,id.vars=1:4)
# results_melt <- results_melt[!is.nan(results_melt$value),] # drop "cluster" cases when K=J
# ggplot(results_melt,aes(x=factor(K),y=value,
#                         color=factor(variable,labels=c("Rank-Clustered","Distinct"))))+
#   facet_grid(s~J,scales="free_x",labeller=label_parsed)+
#   geom_boxplot(outlier.alpha=0.5,position="identity")+theme_bw()+
#   scale_color_manual(values=c("skyblue","darkblue"))+
#   labs(x="Number of Rank-Clusters, K",
#        y="Posterior Rank-Clustering Probability",
#        color=NULL)+
#   theme(legend.position = "bottom",
#         panel.grid.minor = element_blank(),
#         panel.grid.major.x = element_blank())

## ----fig.asp = 0.4------------------------------------------------------------
# J <- 4
# mu_hat <- c(-1,0,1,0)
# s <- c(.1,.1,.1,1)
# 
# set.seed(2)
# data <- matrix(data = rnorm(10000*4,mean=mu_hat,sd=s), ncol = 4, byrow = TRUE)
# 
# forestplot_muhat(data = data, order_by_average = FALSE)

## ----cache=TRUE---------------------------------------------------------------
# mcmc <- mcmc_raceNMA(mu_hat = mu_hat, s = s,
#                    iter = 50000, nu_iter = 2, chains = 4,
#                    burn_prop = 0.5,thin = 1,verbose = FALSE)
# 
# posterior_equal <- matrix(NA,nrow=J,ncol=J)
# for(i in 1:(J-1)){for(j in (i+1):J){
#   posterior_equal[i,j] <- mean(mcmc[,paste0("G",i)] == mcmc[,paste0("G",j)])
# }}
# round(posterior_equal,2)

## ----fig.asp-0.6--------------------------------------------------------------
# wang_ranks <- t(apply(data,1,function(mu){rank(mu)}))
# wang_ranks_probs <- apply(wang_ranks,2,function(rank){
#   sapply(1:J,function(j){mean(rank==j)})
# })
# race_ranks <- t(apply( mcmc[,paste0("mu",1:J)], 1, function(mu){
#   rank(mu,ties.method="min")
# }))
# race_ranks_probs <- apply(race_ranks,2,function(rank){
#   sapply(1:J,function(j){mean(rank==j)})
# })
# 
# g4a<-ggplot(melt(wang_ranks_probs),
#             aes(x=Var1,y=factor(Var2),fill=value)) +
#   geom_tile() + theme_minimal() +
#   scale_y_discrete(limits=rev) +
#   scale_x_continuous(breaks=1:4,limits=c(.5,4.5)) +
#   scale_fill_gradient(low="white",high="black",limits=c(0,1))+
#   labs(x="Rank",y="Treatment",fill="Probability")+
#   theme(panel.grid = element_blank(),legend.position = "bottom")+
#   geom_text(aes(x=Var1,y=factor(Var2),label=round(value,2)),
#         color=ifelse(melt(wang_ranks_probs)$value>0.4,"white","black"))
# 
# g4b<-ggplot(melt(race_ranks_probs),
#             aes(x=Var1,y=factor(Var2,levels=paste0("mu",1:J),
#                                 labels=paste0(1:J)),fill=value))+
#   geom_tile() + theme_minimal() +
#   scale_y_discrete(limits=rev) +
#   scale_x_continuous(breaks=1:4,limits=c(.5,4.5)) +
#   scale_fill_gradient(low="white",high="black",limits=c(0,1)) +
#   labs(x="Rank",y="Treatment",fill="Probability")+
#   theme(panel.grid = element_blank(),legend.position = "bottom")+
#   geom_text(aes(x=Var1,y=factor(Var2,levels=paste0("mu",1:J),labels=paste0(1:J)),
#                 label=round(value,2)),
#             color=ifelse(melt(race_ranks_probs)$value>0.4,"white","black"))
# 
# plot_grid(plot_grid(g4a+theme(legend.position = "none"),
#                     g4b+theme(legend.position = "none"),
#                     labels = c('A', 'B'), label_size = 12),
#           get_plot_component(g4a, 'guide-box-bottom', return_all = TRUE),
#           nrow=2,rel_heights = c(.9,.1) )

## -----------------------------------------------------------------------------
# data("wang_posterior") # loads posterior of non-baseline treatments
# 
# # define assumed mean and variance for baseline treatment (R-CHOP)
# mu_hat_baseline <- 0
# var_baseline <- min(apply(wang_posterior,2,var))/10
# 
# # calculate summary statistics, mu_hat and cov, for all treatments
# mu_hat <- c( mu_hat_baseline, apply(wang_posterior,2,mean) )
# cov <- cbind( c(var_baseline, rep(0, 10) ),
#               rbind(0, cov(wang_posterior)) )
# 
# # store treatment names
# treatments <- c("R-CHOP", names(wang_posterior))

## ----fig.asp = 0.4, warning=FALSE---------------------------------------------
# forestplot_muhat(data=wang_posterior,names=treatments[-1])

## -----------------------------------------------------------------------------
# mcmc <- mcmc_raceNMA(mu_hat = mu_hat, cov = cov,
#                    mu0 = mean(mu_hat), sigma0 = sqrt(10*var(mu_hat)),
#                    iter = 50000, nu_iter = 5, chains = 4,
#                    seed = 1, verbose = FALSE)

## ----fig.asp=0.45-------------------------------------------------------------
# # Figure 6
# g9a <- clusterplot_ranks(data=cbind(0,wang_posterior),
#                             names=treatments,label_ranks=1)+
#   theme(legend.position = "bottom")
# g9b <- clusterplot_ranks(mcmc=mcmc,names=treatments,label_ranks = 1)
# plot_grid(plot_grid(g9a+theme(legend.position = "none"),
#                     g9b+theme(legend.position = "none"),
#                     labels = c('A', 'B'), label_size = 12),
#           get_plot_component(g9a, 'guide-box-bottom', return_all = TRUE),
#           nrow=2,rel_heights = c(.9,.1))
# # Figure 7
# g10a <- cumulativeprobplot_ranks(data=cbind(0,wang_posterior),names=treatments)+
#   theme(legend.position = "bottom") +
#   guides(color = guide_legend(nrow = 2))
# g10b <- cumulativeprobplot_ranks(mcmc=mcmc,names=treatments)
# plot_grid(plot_grid(g10a+theme(legend.position = "none"), g10b+theme(legend.position = "none"),
#                         labels = c('A', 'B'), label_size = 12),
#               get_plot_component(g10a, 'guide-box-bottom', return_all = TRUE),
#           nrow=2,rel_heights = c(.85,.15) )
# # Table 2
# res_WANG <- calculate_SUCRA_MNBT(data = cbind(0,wang_posterior),
#                                  level = 0.50, names = treatments)
# res_RaCENMA <- calculate_SUCRA_MNBT(mcmc = mcmc, level = 0.50, names = treatments)
# table2 <- left_join(res_WANG,res_RaCENMA,by="Treatment")[,c(1,2,4,3,5)]
# names(table2) <- c("Treatment","SUCRA (Wang)","SUCRA (RaCE)","MNBT (Wang)","MNBT (RaCE)")
# print(table2)

## ----fig.asp=0.4--------------------------------------------------------------
# traceplot_K(mcmc)+
#   scale_x_continuous(breaks=seq(125000,250000,by=25000),
#                      labels=paste0(seq(125,250,by=25),"k"))
# traceplot_mu(mcmc,names=treatments)+
#   scale_x_continuous(breaks=seq(125000,250000,by=25000),
#                      labels=paste0(seq(125,250,by=25),"k"))
# forestplot_muhat(mcmc=mcmc,names=treatments)