---
title: 'Bayesian RE fractional polynomial NMA for grouped survival data'
author: "Sandro Gsteiger"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Bayesian RE NMA using FP model (example)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

_Tested on BEE R 3.6.1_

# Introduction
This vignette provides a short example of a Bayesian random effects fractional polynomial NMA model for grouped survival data. A random effect is put on the intercept term only (i.e. on the proportional hazards part of the polynomial) but not on the time-varying terms of the hazard ratio.



## Prepare the environment
```{r, warning = FALSE, results = "hide", message=FALSE}
library(dplyr)
library(gemtc)        
library(gemtcPlus)     
library(ggmcmc)       

```


## Load in the data
```{r, warning = FALSE}
data("grouped_TTE")

```

# First order random intercept model (log-normal distribution for RE variance)

## Plan the model 
```{r, warning = FALSE}
model_plan <-  plan_fp(model.pars = list(exponents = 0, t.eval = "midpoint"),
                       bth.model = "REINT", ref.std = "STUDY2", nma.ref.trt = "B",
                       bth.prior = list(type = "var", distr = "dlnorm", 
                                        meanlog = -2.71, 
                                        sdlog = 1.74),
                       n.chains = 2,
                       n.iter = 6000,
                       n.burnin = 1000,
                       n.thin = 1,
                       descr = "First order random intercept model (LN prior for RE Var)",
                       descr_s = "FP (1o, REint, LN)"
                       )
```

## Ready the data 
```{r, warning = FALSE}
# Returns list list contaiing a jags list ready for input to `nma_fit` and a network object
model_input <- nma_pre_proc(grouped_TTE, model_plan)
```

__Figure__ Network graph
```{r, fig.height=7, fig.width=7} 
plot(model_input$network, displaylabels = TRUE)
```

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


## Post processing
```{r, warning = FALSE}

# Prepare plot data
nodes <- colnames(as.mcmc(model)[[1]])
sel <- grep("d[2,", nodes, fixed = TRUE)
plot_data <- ggs(as.mcmc(model)[, sel])

```

Produce diagnostic plots to further assess convergence. Here: select the contrasts trt 2 vs trt 1 for visibility.


__Figure__ Traceplot
```{r, results='asis', fig.height=5, fig.width=6}
ggs_traceplot(plot_data)
```

__Figure__ Densityplot
```{r, results='asis', fig.height=5, fig.width=6}
ggs_density(plot_data)
```

__Figure__ Auto-correlation plot
```{r, results='asis', fig.height=5, fig.width=6}
ggs_autocorrelation(plot_data)
```

__Figure__ Running means
```{r, results='asis', fig.height=5, fig.width=6}
ggs_running(plot_data)
```


The diagnostic plots show that the chains are _much_ (!) too short in this example (for real use, run longer chains and consider thinning).


Save the results for later use.
```{r}
res_re1 <- model
rm(model)
```


# Second order random intercept model (log-normal distribution for RE variance)

## Plan the model 
```{r, warning = FALSE}
model_plan <-  plan_fp(model.pars = list(exponents = c(0, 1), 
                                         t.eval = "midpoint"),
                       bth.model = "REINT", ref.std = "STUDY2", nma.ref.trt = "B",
                       bth.prior = list(type = "var", distr = "dlnorm", 
                                        meanlog = -2.71, 
                                        sdlog = 1.74),
                       n.chains = 2,
                       n.iter = 6000,
                       n.burnin = 1000,
                       n.thin = 1,
                       descr = "Second order random intercept model (LN prior for RE Var)",
                       descr_s = "FP (2o, REint, LN)"
                       )
```

## Ready the data 
```{r, warning = FALSE}
# Returns list list contaiing a jags list ready for input to `nma_fit` and a network object
model_input <- nma_pre_proc(grouped_TTE, model_plan)
```
## Fit the model 
```{r, warning = FALSE}
model  <- nma_fit(model_input = model_input)
```




## Post processing
```{r, warning = FALSE}

# Prepare plot data
nodes <- colnames(as.mcmc(model)[[1]])
sel <- grep("d[2,", nodes, fixed = TRUE)
plot_data <- ggs(as.mcmc(model)[, sel])

```

__Figure__ Traceplot
```{r, results='asis', fig.height=5, fig.width=6}
ggs_traceplot(plot_data)
```

The chains are obviously MUCH too short!


Save the results for later use.
```{r}
res_re2 <- model
rm(model)
```


# First order random intercept model (uniform distribution for RE standard deviation)

## Plan the model 
```{r, warning = FALSE}
model_plan <-  plan_fp(model.pars = list(exponents = 0, t.eval = "midpoint"),
                       bth.model = "REINT", ref.std = "STUDY2", nma.ref.trt = "B",
                       bth.prior = list(type = "sd", distr = "unif", min = 0, max = 2),
                       n.chains = 2,
                       n.iter = 6000,
                       n.burnin = 1000,
                       n.thin = 1,
                       descr = "First order random intercept model (unif[0,2] prior)",
                       descr_s = "FP (1o, REint, U[0,2])"
                       )
```

## Ready the data 
```{r, warning = FALSE}
# Returns list list contaiing a jags list ready for input to `nma_fit` and a network object
model_input <- nma_pre_proc(grouped_TTE, model_plan)
```


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

## Post processing
Do convergence assessment as above (not shown).

The chains are obviously MUCH too short!

Save the results for later use.
```{r}
res_re3 <- model
rm(model)
```


# Second order random intercept model (uniform distribution for RE standard deviation)

## Plan the model 
```{r, warning = FALSE}
model_plan <-  plan_fp(model.pars = list(exponents = c(0, 1), t.eval = "midpoint"),
                       bth.model = "REINT", ref.std = "STUDY2", nma.ref.trt = "B",
                       bth.prior = list(type = "sd", distr = "unif", min = 0, max = 2),
                       n.chains = 2,
                       n.iter = 6000,
                       n.burnin = 1000,
                       n.thin = 1,
                       descr = "Second order random intercept model (unif[0,2] prior)",
                       descr_s = "FP (2o, REint, U[0,2])"
                       )
```

## Ready the data 
```{r, warning = FALSE}
# Returns list list contaiing a jags list ready for input to `nma_fit` and a network object
model_input <- nma_pre_proc(grouped_TTE, model_plan)
```


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

## Post processing
Do convergence assessment as above (not shown).

The chains are obviously MUCH too short!

Save the results for later use.
```{r}
res_re4 <- model
rm(model)
```



# Produce outputs of interest

Start with an object collecting all fits.
```{r}
all_res <- list(res_re1, res_re2, res_re3, res_re4)
```


## Model comparison
```{r, results='asis'}
dcompare <- get_fp_comparison(all_res)
cat("__Table__ Model comparison")
pander::pandoc.table(dcompare, row.names = FALSE, split.tables = Inf)
```


## Hazard ratio estimates
```{r, results='asis', warning=FALSE, fig.height=6, fig.width=6}
# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("### ", title, "  \n")

  ## Tables: Hazard ratio estimates for each segment

  HR_rev <- get_fp_HR(x = seq(0.1, 24, 0.1), 
                      fit = res_i, 
                      trt.nos = 1:res_i$data.jg$Ntrt, 
                      ref.no = 2, # use this treatment in the list of treatments as reference for the HRs 
                      revert = TRUE # revert to get HRs ref vs other treatments
                      ) 
  
  fig1 <- plot_fp_HR(HR_rev, xlab = "Month", breaks = c(0.25, 0.5, 1, 2))
  fig2 <- plot_fp_HR(HR_rev, xlab = "Month", breaks = c(0.25, 0.5, 1, 2), facet = TRUE) # !! HERE: ADD CIs !!
  
  
  dHRtab <- HR_rev %>%
    mutate(Month = round(x, 1)) %>%        # trick, otherwise equality testing fails to pick out all timepoints in filter step
    filter(Month %in% seq(3, 24, 3)) %>%
    mutate(Comparison = lab, Month, HR = round(median, 3), lCI = round(lCI, 3), uCI = round(uCI, 3)) %>%
    select(Comparison, Month, HR, lCI, uCI)
  
  
  cat("__Figure__ Hazard ratios treatment A vs other treatments\n")
  plot(fig1)
  cat("\n\n")
  cat("__Figure__ Hazard ratios treatment A vs other treatments (multi-panel)\n")
  plot(fig2)
  cat("\n\n")
  cat("__Table__ Hazard ratios treatment A vs other treatments\n")
  pander::pandoc.table(dHRtab, row.names = FALSE, split.tables = Inf)
  cat("\n\n")
  cat("\n\n")
  
  rm(HR_rev)
  rm(dHRtab)
  rm(fig1)
  rm(fig2)
  rm(res_i)
}
```


## Survivor function estimates
The NMA baseline estimate from the `ref_trt` arm from `ref_std` is used. These are combined with the time-varying hazard-ratio functions from the NMA to obtain the survivor functions for the other interventions.

```{r, results='asis', warning=FALSE, fig.height=4, fig.width=6}
ref_trt <- "B"
ref_std <- "STUDY2"
hor <- 60

# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("### ", title, "  \n")
  
  ## Plots of survivor functions over time ("NMA result"), ref study/arm and timehorizons specified in settings function
  sel_ref <- which(attr(res_i$data.jg, "d_arms")$study == ref_std & attr(res_i$data.jg, "d_arms")$treatment == ref_trt)
  id_ref_std <- attr(res_i$data.jg, "d_arms")$studyn[sel_ref]
  id_ref_arm <- attr(res_i$data.jg, "d_arms")$arm[sel_ref]
    
  S_extrap <- get_fp_S(fit = res_i, 
                       ref.std = id_ref_std, 
                       ref.arm = id_ref_arm, 
                       time = seq(0.1, hor, 0.1))

  fig <- ggplot(data = S_extrap) +        
    geom_line(aes(x = time, y = S, col = treatment, linetype = treatment)) +
    ylim(0, 1) +
    xlab("Month") + ylab("Survival probability") +
    theme_bw() +
    theme(legend.title = element_blank())
  cat("__Figure__ Survivor function estimates (time horizon:", hor, "months) \n")
  plot(fig)
  cat("\n\n")
  
  fig <- ggplot(data = S_extrap) + 
    facet_wrap(~treatment) +
    geom_ribbon(aes(x = time, ymin = lCrI, ymax = uCrI), fill = "lightblue", alpha = 0.8) +
    geom_line(aes(x = time, y = S)) +
    ylim(0, 1) +
    xlab("Month") + ylab("Survival probability") +
    theme_bw()
  cat("__Figure__ Survivor function estimates by treatment (time horizon:", hor, "months) \n")
  plot(fig)
  cat("\n\n")
  
  rm(list = c("S_extrap", "fig"))
  rm(res_i)
}
```


## Model fit: observed KM data vs estimated S(t)
For every arm in every study, the study baseline hazard estimate is combined with the corresponding contrast estimate (both from the NMA) to obtain the estimated survivor functions.

```{r, results='asis', warning=FALSE, fig.height=6, fig.width=6}
hor <- 36
# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("### ", title, "  \n")
  
  gof <- get_fp_GoF(fit = res_i, time = seq(0.1, hor, 0.1))

  fig <- ggplot() + 
    geom_line(data = gof %>% filter(type == "nma"), aes(x = time, y = S, col = treatment)) +
    geom_line(data = gof %>% filter(type == "obs"), aes(x = time, y = S, col = treatment), linetype = "dashed") +
    facet_wrap(~study, ncol = 2) +
    ylim(0, 1) + xlim(0, 36) +
    xlab("Month") + ylab("Survival probability") +
    theme_bw() +
    theme(legend.position = "top", legend.title = element_blank())

  cat("__Figure__ Goodness-of-fit: estimated (solid lines) and observed (dashed) survivor functions for each study\n")
  plot(fig)
  cat("\n\n")
  
  rm(list = c("gof", "fig"))
  rm(res_i)
}
```


## Parameter estimates (baseline and contrasts)

```{r, results='asis', warning=FALSE, fig.height=6, fig.width=6}
# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("### ", title, "  \n")

  cest <- get_fp_contrasts(res_i)

  cat("\n\n")
  cat("__Table__ Contrast estimates in fractional polynomial vs network reference\n")
  pander::pandoc.table(cest, row.names = FALSE, split.tables = Inf)
  cat("\n\n")

  rm(cest)
  rm(res_i)
}
```



## Posterior correlation plots

```{r, results='asis', warning=FALSE, fig.height=6, fig.width=6}
# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("### ", title, "  \n")

  corrs <- get_fp_corrs(res_i)

  for(j in 1:res_i$data.jg$Ntrt){
    cat("\n\n")
    cat("__Table__ Posterior correlations of (multivariate) contrasts for", dimnames(corrs)$treatment[j],"vs reference\n")
    pander::pandoc.table(corrs[j,,], row.names = dimnames(corrs[j,,])[[1]], split.tables = Inf)
    cat("\n\n")
  }
  
  rm(corrs)
  rm(res_i)
}
```


# Appendix

```{r, results='asis', warning=FALSE}
# loop through fits
for(i in seq_along(all_res)){
  res_i <- all_res[[i]]
  title <- res_i$descr
  cat("## ", title, "  \n\n")

  jginfo <- get_jags_info(res_i, include.comments = TRUE)
  cat("```\n", jginfo, "\n```\n\n")
  
  rm(jginfo)
  rm(out)
}  
```


# Session info
```{r}
date()
sessionInfo()
```