---
title: "Reproducibility"
output:
  rmarkdown::html_document:
    toc: true
    toc_depth: 3
    toc_float: true  
number_sections: true  
vignette: >
  %\VignetteIndexEntry{reproducibility}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  eval=FALSE
)
```

```{r warning=FALSE, message=FALSE}
library(RaCE.NMA)
library(ggplot2)
library(dplyr)
library(cowplot)
library(gridExtra)
library(reshape2)
```

In this vignette, we provide code to reproduce the analyses presented in the paper associated with this R package. Specifically, we reproduce two simulation studies (Sections 3.1 and 3.2) and a case study (Section 4). We assume knowledge of the associated R package functions; please refer to the "Tutorial" and "Reference" pages for more information. No figures are produced; see paper.

## Reproduction of Simulation Study 1 (Section 3.1)

The following code reproduces Figure 1.

```{r 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)
```

The following code reproduces Figure 2: First, we perform a simulation study. Second, we create the figure.

```{r 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))
        )
      }
    }
  }
}
```
```{r}
## 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())
```

## Reproduction of Simulation Study 2 (Section 3.2)

The following code reproduces Figure 3.

```{r 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)
```

The following code reproduces Table 1: First, we fit the RaCE model to the simulated data. Second, we calculate and display the values in Table 1.

```{r 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)
```
The following code reproduces Figure 4: First, we calculate the posterior rank-clustering probability under a traditional model (i.e., assumption of distinct ranks). Second, we calculate the posterior rank-clustering probability based on the RaCE model. Third, we produce the figure.

```{r 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) )
```

## Reproduction of Case Study (Section 4)

We first load the posteriors from Wang et al. (2022) and calculate the mean and covariances of the relative treatment effects.

```{r}
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))
```

The following code chunk reproduces Figure 5.

```{r fig.asp = 0.4, warning=FALSE}
forestplot_muhat(data=wang_posterior,names=treatments[-1])
```

Next, we fit the RaCE model to estimated relative treatment effects from Wang et al. (2022).

```{r}
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)
```

The following code chunk reproduces Figure 6, Figure 7, and Table 2: First, we calculate the values underpinning each graphic. Then, we create each in turn.

```{r 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)
```

The following code chunk reproduces the Figures 8--10 (Appendix B).

```{r 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)
```