---
title: "Bayesian FE & RE NMA for HR data (via gemtc package) - Result Generation"
author: "Sandro Gsteiger, Nick Howlett"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Bayesian FE & RE NMA for HR data (via gemtc package) - Result Generation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## Introduction
This vignette provides a short example of a Bayesian NMA for HR data. The model fit relies on the `gemtc` package, pre- and post-processing is done with `gemtcPlus`.

## Prepare the environment
```{r setup, warning=FALSE, message=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(dplyr)
library(gemtc)
library(gemtcPlus)
library(ggmcmc)

```


## Load in the data
```{r}
# load example data
data("hr_data", package = "gemtcPlus")
```

## Plan the model 
```{r}
#Plan model
model_plan <- plan_hr(bth.model = "FE", 
                          n.chain = 3, 
                          n.iter = 6000, 
                          thin = 1,
                          n.adapt = 1000, 
                          link = "identity",
                          linearModel = "fixed")
```

## Ready the data
```{r, echo=FALSE}
model_input <-  nma_pre_proc(data = hr_data, plan = model_plan)

```

## Fit the model 
```{r}
model  <- nma_fit(model_input = model_input)
```

## Post processing
```{r}
ggs_traceplot(ggs(model$samples))
ggs_density(ggs(model$samples))
summary(model)
plot(model_input$fitting_data)
get_mtc_sum(model)
```


# Random Effects example
Update model plan and re-run fit.
```{r}
#Plan model
model_plan <- plan_hr(bth.model = "RE", 
                      n.chain = 3,
                      n.iter = 6000,
                      thin = 1,
                      n.adapt = 1000,
                      link = "identity",
                      linearModel = "random", 
                      bth.prior = mtc.hy.prior(type = "var", distr = "dlnorm",-4.18, 1 / 1.41 ^ 2)
)


```

## Fit the model 
```{r}
model  <- nma_fit(model_input = model_input)
```


# Post processing
## Inspect convergence
The `ggmcmc` package provides ggplot2 versions of all major convergence plots and diagnostics. 

__Figure__ Traceplot
```{r, results='asis', fig.height=9, fig.width=6}
ggs_traceplot(ggs(model$samples))
```

__Figure__ Densityplot
```{r, results='asis', fig.height=9, fig.width=6}
ggs_density(ggs(model$samples))
```

__Figure__ Brooks-Gelman-Rubin convergence diagnostic (Rhat)
```{r, results='asis', fig.height=6, fig.width=6}
ggs_Rhat(ggs(model$samples))
```

__Figure__ Auto-correlation plot
```{r, results='asis', fig.height=9, fig.width=6}
ggs_autocorrelation(ggs(model$samples))
```

__Figure__ Running means
```{r, results='asis', fig.height=9, fig.width=6}
ggs_running(ggs(model$samples))
```


# Produce outputs of interest

## Posterior summaries (log-scale)
The contrasts in this model are log-hazard ratios (which correspond to differences in log-hazard rates).

Unfortunately, `gemtc` does not provide an estimate of the effective sample size (`n.eff`). Instead, a time-series SE is given. As a rule of thumb, the length of the MCMC is sufficient if the time-series SE is smaller than 2%(-5%) of the posterior SD.

```{r}
summary(model)
```

In the example here, the chain length seems borderline (sufficient for posterior means and medians, but rather a bit too small for stable 95% credible intervals).

To judge overall model fit, the residual deviance should be compared to the number of independent data points (which can be done via a small utility function in `gemtcPlus`).

```{r}
get_mtc_sum(model)
```


## Hazard ratio estimates

Assume new treatment is "A" and is to be compared vs all other treatments.

__Table__ Hazard ratios A vs other treatments
```{r}
HR <- get_mtc_newVsAll(model, new.lab = "A", transform = "exp", digits = 2)
HR
```

__Table__ Probability A better than other treatments (better meaning smaller HR)
```{r}
get_mtc_probBetter(model, new.lab = "A", smaller.is.better = TRUE, sort.by = "effect")
```

__Figure__ Forest plot A vs other treatments
```{r, results='asis', fig.height=3, fig.width=6, warning=FALSE}
plot_mtc_forest(x = HR, lab = "Hazard ratio A vs others", sort.by = "effect")  
```

__Table__ Cross-tabulation of HRs
```{r, results='asis'}
ctab <- round(exp(relative.effect.table(model)), 2)
pander::pandoc.table(as.data.frame(ctab), split.tables = Inf)
```


## Treatment rankings
```{r}
rk  <- rank.probability(model, preferredDirection = -1)
mrk <- reshape2::melt(rk[,], varnames = c("Treatment", "Rank"), value.name = "Probability")

fig <- ggplot(data = mrk) +
  geom_line(aes(Rank, Probability, color = Treatment, linetype = Treatment), size = 1.5) +
  theme_bw()
```

__Figure__ Rankogram
```{r, results='asis', fig.height=3, fig.width=6}
plot(fig)
```

__Table__ Rank probabilities
```{r}
colnames(rk) <- paste("Rank", 1:ncol(rk))
rk
```


# Extract model code (e.g. for Appendix)
```{r}
cat(model$model$code)
```


# Session info
BEE repository: `r getwd()`

```{r}
date()
sessionInfo()
```
