# Reproducibility

    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.

    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.

    ## 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,cores=4)

            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())

## Reproduction of Simulation Study 2 (Section 3.2)

The following code reproduces Figure 3.

    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.

    mcmc <- mcmc_raceNMA(mu_hat = mu_hat, s = s,
                       iter = 50000, nu_iter = 2, chains = 4,
                       burn_prop = 0.5,thin = 1,cores = 4, 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.

    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.

    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.

    forestplot_muhat(data=wang_posterior,names=treatments[-1])

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

    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, cores = 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.

    # 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 B1–B3 (Appendix B).

    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)