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

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  warning=FALSE,
  comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)
```

# Overview

In this vignette, we provide a detailed tutorial of the key functionalities of the `RaCE.NMA` R package. The package has functions to complete three primary tasks in the analysis pipeline when utilizing the RaCE-NMA methodology:

1. Model estimation
2. Model assessment
3. Inference and summary statistics

After loading the necessary packages to run our analyses, we demonstrate each of the primary tasks in the analysis pipeline on two datasets: A toy dataset (as seen in Simulation Study 2 of the paper) and real data (as seen in the Case Study of the paper).

```{r warning=FALSE,message=FALSE}
library(RaCE.NMA)
library(reshape2) # data reformatting
library(ggplot2)  # creating nice figures
```

# Toy Data

First, we use a "toy" dataset that can be seen in Simulation Study 2 of the paper. The data is included in the `RaCE.NMA` package and can be loaded using the following code:

```{r}
data("toy_data")
head(toy_data)
```

The data contains $10{,}000$ rows and 4 columns. Each column represented a treatment, and each row a posterior draw of relative treatment effects for those treatments. We assume these posterior draws are based on a standard NMA model (i.e., no rank-clustering). 

To visualize these posterior draws, we create a forest plot (with 95\% CIs) via the aptly-named `create_forestplot` function:

```{r fig.asp=0.45, fig.cap="Figure 1: Forest plot of estimated relative treatment effects in a toy example dataset, based on a standard NMA model."}
forestplot_muhat(data = toy_data, level = 0.95, 
                  order_by_average = FALSE)
```

We observe that the first three treatments have narrow and well-separated credible intervals of relative efficacy. The fourth treatment, however, has a very wide interval to indicate high uncertainty in its relative efficacy.

Now we are ready to start the analysis pipeline.

## 1. Model Estimation

The `RaCE.NMA` packages contains a primary model-fitting function, `mcmc_RCMVN`. The function requires many inputs, which may be broken down into three categories:

  + Previous results from an NMA analysis: The user must provide results on the relative treatment effects from a previous NMA study. Various formats are allowed; see below.
  + Model hyperparameters: Parameters which correspond to "priors" in the Bayesian model. Vague priors are available by default.
  + MCMC settings: Specifications of an MCMC estimation procedure, including the number of independent chains, length of each chain, thinning and burn-in, etc. By default, the function runs 2 chains, each of length 15000 and removes the first half of iterations as burn-in.
  
### Results of a standard NMA model

If available, we recommend inputting the full posterior of the relative treatment effect from a previous NMA analysis into the `mcmc_RCMVN` function. In this case, the `posterior` argument should be provided an $N\times J$ matrix, where the (i,j) entry indicates the relative treatment effect of treatment $j$ in MCMC draw $i$. An example of this functionality is provided below:

```{r eval=FALSE}
mcmc_raceNMA(posterior = toy_data)
```

In the above code, default settings are used for all other inputs. (We recommend against this practice.)

Alternatively, one may provide summary statistics of the model posterior (instead of a full posterior itself). This is useful when applying the RaCE-NMA methodology to the results of a published paper post-hoc, when only summary statistics are available. In this setting, the user must provide two inputs:

  + `mu_hat`: A vector of length $J$ (the number of treatments), where the $j$th entry displays the average relative treatment effect of treatment $j$
  + `cov` or `s`: A $J\times J$ covariance matrix of relative treatment effects *or* a vector of length $J$ of treatment-specific standard deviations of relative treatment effects.
  
An example of this functionality is provided below:
```{r eval=TRUE}
mu_hat_toy <- apply(toy_data, 2, mean)
cov_toy <- cov(toy_data)
s_toy <- apply(toy_data, 2, sd)
```
```{r eval=FALSE}
mcmc_raceNMA(mu_hat = mu_hat_toy, cov = cov_toy) # example with mu_hat and cov
mcmc_raceNMA(mu_hat = mu_hat_toy, s   = s_toy  ) # example with mu_hat and s
```
  
If duplicative inputs are provided, the function will always use those provided the "richest" information. Specifically,

  + If `posterior` is provided, `mu_hat`, `cov`, and/or `s` will be ignored. 
  + If `cov` is provided, `s` is ignored. 


### Model hyperparameters

Next, we may specify the following hyperparameters:

  + `mu0`: The prior mean on the rank-clustered relative treatment effects. Defaults to `mean(mu_hat)`, which aims to be vague.
  + `sigma0`: The prior standard deviation on the rank-clustered relative treatment effects. Defaults to `sqrt(10*var(mu_hat))`, which aims to be vague.
  
The following code demonstrates how one could specify these hyperparameters directly in the estimation function:

```{r eval=FALSE}
mcmc_raceNMA(
  posterior = toy_data, # results of a standard NMA model
  mu0 = mean(mu_hat_toy), sigma0 = sqrt(10*var(mu_hat_toy)), # model hyperparameters
)
```

### MCMC settings

Last, we may specify parameters that control how the MCMC chains are run:

  + `tau`: The standard deviation of the Metropolis Hastings proposal distribution. Defaults to `min(|mu_hat_i-mu_hat_j|)`. 
  + `nu0`: How cluster-specific mean parameters are initialized. Defaults to `NULL`, which randomly samples from the prior distribution.
  + `iter`: The number of times the rank-clustering structure is sampled in each MCMC chain. Defaults to $4{,}000$.
  + `nu_iter`: The number of times each cluster-specific mean is sampled per sampling of the rank-clustering structure. Defaults to $5$. In total, there will be `iter`$\times$`nu_iter` samples from the posterior.
  + `chains`: The number of independent MCMC chains. Defaults to 2.
  + `burn_prop`: The proportion of MCMC samples in each chain to be removed as burn-in. Defaults to 0.5.
  + `thin`: A number indicating how often to thin samples. Defaults to 1, indicating no thinning.
  + `seed`: The random seed to run chains. Defaults to `NULL`, meaning the environment seed is inherited.
  + `cores`: How many CPU cores to use when drawing chains. The default, `cores=1`, draws chains sequentially. Values above 1 draw chains in parallel; `cores` is silently capped at the number of  `chains`. Results do not depend on `cores`.

The following code demonstrates how one could specify some of these settings directly in the estimation function:

```{r}
mcmc_results_toy <- mcmc_raceNMA(
  posterior = toy_data, # results of a standard NMA model
  mu0 = mean(mu_hat_toy), sigma0 = sqrt(10*var(mu_hat_toy)), # model hyperparameters
  iter = 20000, nu_iter = 3, # MCMC settings
  chains = 2, burn_prop = 0.5, thin = 2, seed = 1, cores=2
)
```

## 2. Model Assessment

The package contains three functions to assess if the MCMC chains have converged and mixed. The first two produce trace plots of `mu` and `K`, the primary model parameters. The third calculates the $\hat{R}$ statistic for each `mu` parameter. We display and interpret each plot below.

```{r, fig.cap="Figure 2: Trace plots of relative treatment effects in a toy example dataset, based on the RaCE-NMA model."}
traceplot_mu(mcmc_results_toy)
```

We observe acceptable mixing and convergence in the trace plots for $\mu$. Ideally, the two chains would like similar to each other, with each exhibiting a "white noise pattern". This is certainly true for $\mu_1$, $\mu_2$, and $\mu_3$. For $\mu_4$, we notice some "jumps" in value across iterations. Based on our toy data, in which treatment 4 has extremely high uncertainty, we expect this type of behavior. Furthermore, since the chains both exhibit similar patterns, we are satisfied.

```{r, fig.cap="Figure 3: Trace plot of the number of rank-clusters, K, in a toy example dataset, based on the RaCE-NMA model."}
traceplot_K(mcmc_results_toy)
```

We again observe acceptable mixing and convergence in the trace plot for $K$. That is, the two chains are indistinguishable and exhibit a seemingly random pattern between the values $3$ and $4$.

```{r, fig.cap="Table 1: R-hat statistics for the relative treatment effects in a toy example dataset, based on the RaCE-NMA model."}
calculate_Rhat(mcmc_results_toy, level = 0.9)
```
The $\hat{R}$ statistic is a standard convergence diagnostic statistic used to assess MCMC chains. A standard rule of thumb is to ensure that $\hat{R}<1.1$. We see that the point estimate and 90% confidence intervals for all $\mu$ parameters are $\leq1.1$, indicating acceptable convergence.

## 3. Inference and summary statistics

The package contains three functions to visualize model results. We demonstrate each in turn:

Posterior rank probabilities for each intervention. Ranks are displayed on the ordinal scale (1 = best). Interventions are grouped into rank clusters representing non-exclusive superiority sets; multiple interventions within the same cluster may share identical posterior ranks.

The first is to create a clustering probability matrix.  In the figure, each row represents a treatment and each column a rank level. The color of the corresponding cell is the estimated posterior probability that a treatment belongs to a specific rank. Note that the `label_ranks` input allows the user to specify which rank levels in the graphic are labeled with a *numerical* probability (corresponding to the color shown). See the figure and corresponding interpretation below.

```{r fig.asp=0.6, fig.cap="Figure 4: Posterior rank probabilities for each intervention in a toy example dataset, based on the RaCE-NMA model. Ranks are displayed on the ordinal scale (1 = best). Interventions are grouped into rank clusters representing non-exclusive superiority sets; multiple interventions within the same cluster may share identical posterior ranks."}
clusterplot_ranks(mcmc=mcmc_results_toy, label_ranks = 1:4)
```

The model estimates a 0.97 probability that Treatment 1 is in first place. Simultaneously, the model estimates a 0.26 probability that Treatment 4 is in first place. Thus, there is some probability, based on the model, that Treatments 1 and 4 are simultaneously the "best" treatments. That being said, Treatment 4 has substantial uncertainty in its rank, which aligns with the high uncertainty in original relative treatment effect data. Treatment 2 has a 0.66 probability of being in second place and Treatment 3 has a 0.80 probability of being in fourth (last) place.


The second function creates a cumulative ranking probability curve for each treatment. In words, the plot displays the probability that each treatment is *at least* in a specific rank place. For example, we see that Treatment 4 has a 0.80 probability of being in first or second place.

```{r fig.asp=0.4, fig.cap="Figure 5: Cumulative ranking probability curves for each intervention in a toy example dataset, based on the RaCE-NMA model."}
cumulativeprobplot_ranks(mcmc=mcmc_results_toy)
```

The third function creates a table displaying the SUCRA value and the Median Number of Better Treatments (MNBT) for each treatment (with a corresponding credible interval). The credible level may be specified using the `credible` argument.

```{r, fig.cap="Table 2: SUCRA and MNBT statistics for each intervention in a toy example dataset, based on the RaCE-NMA model."}
calculate_SUCRA_MNBT(mcmc=mcmc_results_toy, level = 0.5)
```

# Real Data (Wang et. al [2022])

Next, we demonstrate the package on a real NMA study provided by Wang et. al (2022) and analyzed as a Case Study in our paper. The data is included in the `RaCE.NMA` package and can be loaded using the following code:

```{r}
data("wang_posterior")
round(head(wang_posterior), 3)
```

It is important to note that `wang_posterior` only includes relative treatment effects for non-baseline treatments. An 11th treatment, "R-CHOP", was the baseline with a relative treatment effect of 0 by construction. Thus, for technical reasons described in the paper, we will need to define summary statistics that can be used when fitting the RaCE-NMA model.

```{r}
# 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
mu_hat_wang <- c( mu_hat_baseline, apply(wang_posterior,2,mean) )
cov_wang <- cbind( c(var_baseline, rep(0, 10) ),
              rbind(0, cov(wang_posterior)) )

# store treatment names
treatments <- c("R-CHOP", names(wang_posterior))
```

The following forest plot displays the relative treatment effects of the 11 treatments in the study, as estimated by Wang et. al (2022).

```{r fig.asp=0.5, fig.cap="Figure 6: Forest plot of estimated relative treatment effects in a network meta-analysis by Wang et. al (2022), based on a standard NMA model."}
forestplot_muhat(data = cbind(0,wang_posterior), 
                 names = treatments)
```

Now, we proceed through the analysis pipeline.

## 1. Model Estimation

We fit the model below using 2 chains, each with $20{,}000\times 4$ MCMC iterations. We parallelize to use 2 CPU cores, one per chain.

```{r cache=TRUE}
mcmc_results_casestudy <- mcmc_raceNMA(
  mu_hat = mu_hat_wang, cov = cov_wang,
  iter = 20000, nu_iter = 4, chains = 2, seed = 1, cores = 2
)
```

## 2. Model Assessment

Next, we view trace plots and convergence diagnostics.

```{r, fig.cap="Figure 7: Trace plots of relative treatment effects in a network meta-analysis by Wang et. al (2022), based on the RaCE-NMA model."}
traceplot_mu(mcmc_results_casestudy, names=treatments)
```
```{r, fig.cap="Figure 8: Trace plot of the number of rank-clusters, K, in a network meta-analysis by Wang et. al (2022), based on the RaCE-NMA model."}
traceplot_K(mcmc_results_casestudy)
```
```{r, fig.cap="Table 3: R-hat statistics for the relative treatment effects in a network meta-analysis by Wang et. al (2022), based on the RaCE-NMA model."}
calculate_Rhat(mcmc_results_casestudy, names=treatments, level = 0.8)
```

We observe that chains have probably not run long enough (refer to the Reproducibility tab for a longer run of chains). However, for demonstration purposes we continue with the next stage of the analysis pipeline below.

## 3. Inference and summary statistics

Disclaimer: The results below are based on chains which have not converged, but are shown for demonstration purposes. Please see the Reproducibility tab for accurate results based on this data.

```{r fig.asp=0.6, fig.cap="Figure 9: Posterior rank probabilities for each intervention in a network meta-analysis by Wang et. al (2022), based on the RaCE-NMA model. Ranks are displayed on the ordinal scale (1 = best). Interventions are grouped into rank clusters representing non-exclusive superiority sets; multiple interventions within the same cluster may share identical posterior ranks."}
clusterplot_ranks(mcmc=mcmc_results_casestudy, names = treatments, label_ranks = 1)
```
```{r fig.asp=0.4, fig.cap="Figure 10: Cumulative ranking probability curves for each intervention in a network meta-analysis by Wang et. al (2022), based on the RaCE-NMA model."}
cumulativeprobplot_ranks(mcmc=mcmc_results_casestudy, names=treatments)
```
```{r, fig.cap="Table 4: SUCRA and MNBT statistics for each intervention in a network meta-analysis by Wang et. al (2022), based on the RaCE-NMA model."}
calculate_SUCRA_MNBT(mcmc=mcmc_results_casestudy,names=treatments,level=0.50)
```
