---
title: "Introduction to BASSLINE"
date: '`r Sys.Date()`'
author:
  - name: Nathan Constantine-Cooke
    email: nathan.constantine-cooke@ed.ac.uk
  - name: Catalina Vallejos
    email: catalina.vallejos@igmm.ed.ac.uk

output:
  rmarkdown::html_vignette:
    toc: true
    number_sections: false
    toc_depth: 2
    fig_width: 6
    fig_height: 4

bibliography: library.bib

vignette: >
  %\VignetteIndexEntry{Introduction to BASSLINE}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev="png"
)
```


## Introduction

```{r setup}
library(BASSLINE)
``` 


BASSLINE (**BA**ye**S**ian **S**urvival ana**L**ys**I**s usi**N**g shap**E**
mixtures of log-normal distributions) uses shape mixtures of log-normal 
distributions to fit data with fat tails and has been adapted from code produced
for *Objective Bayesian Survival Analysis Using Shape Mixtures of Log-Normal 
*Distributions [@Vallejos2015]. Some of the functions have been rewritten 
in C++ for increased performance.

5 distributions from the log-normal family are supported by `BASSLINE`:

* The log-normal distribution
* The log student's T distribution
* The log-logistic distribution
* The log-Laplace distribution
* The log-exponential power distribution

As well as MCMC (Markov chain Monte Carlo) algorithms for the 5
distributions, additional functions which allow log-marginal likelihood
estimators and deviance information  criteria to be calculated are provided.
Case deletion analysis and outlier detection are also supported.

This vignette demonstrates how to use the `BASSLINE` package to carry out
survival analysis using the included `cancer` data-set from the veterans
administration lung cancer trial. 


## Quick Start

Essential parameters for running the MCMC are : 

* `N`: total number of iterations
* `thin`: length of the thinning period (i.e. only every `thin` 
iterations will be stored in the output)
* `burn`: length of burn-in period (i.e. the initial `burn` 
iterations that will be discarded from the output)
* `Time`: Vector of survival times
* `Cens`: Vector indicating if observations are censored
* `X`: Matrix of covariates for each observation

Starting values for the $\beta$s and $\sigma^2$ are randomly sampled
from an appropriate distribution if not given by the user as arguments.
Additional arguments allow the type of prior to be specified, and the
observations to be specified as set or point observations. See the documentation
for any of the MCMC functions included with `BASSLINE` for more information on
these additional arguments.

```{r MCMC Help, eval = FALSE}
?MCMC_LN()
```

Note that BASSLINE does not support factors/ levels. Factors should be converted
to separate binary variables for each level which can be easily done via the
provided `BASSLINE_convert` function. For example:

```{r Original Table, echo = F}
df <- data.frame(Time = c(10,15,24,21), Cens = c(1,1,0,1),
                 treatment = as.factor(c("A", "B", "C", "A")))
```

```{r, echo = FALSE}
knitr::kable(df)
```

can be converted by simply passing the dataframe object to the function.  

```{r Converted Table}
converted <- BASSLINE_convert(df)
```

```{r echo = FALSE}
knitr::kable(converted)
```


## The Cancer Data Set

Included with `BASSLINE` is an example data set, `cancer`. This data has been
obtained from a study conducted by the US Veterans Administration where male
patients with advanced inoperable lung cancer were given either standard or
experimental chemotherapy treatment [@VACancerTrial]. 137 patients took part in
the trial, 9 of whom left the study before their death and are thus right
censored. Various covariates were also documented for each patient.

Viewing the first 5 observations shows the data set's format:

```{r Cancer View,  eval = F}
data(cancer)
head(cancer, 5)
```

```{r Cancer Table, echo = F}
knitr::kable(cancer[1:5, ])
```

The first column of `cancer` denotes the survival time for the observation. 
The second column denotes the censored status for the observation (0 for right
censored; 1 for not censored ). All remaining columns are covariates. 
Additional information can be found in the documentation for `cancer`. 

```{r Cancer Help, eval = F}
?cancer
```

```{r Time and Cens}
Time <- cancer[,"Time"]
Cens <- cancer[,"Cens"]
```

## MCMC chains 

MCMC chains can be easily generated by providing the aforementioned essential
parameters. As previous discussed, starting values are randomly sampled if not
provided by the user. The user will possibly obtain improved results from 
experimenting with these starting values. 

Please note N = 1000, as used in these examples, is not enough to reach
convergence and is only used as a demonstration. The user is advised to run
longer chains with a longer burn-in period for more accurate estimations
(especially for the log-exponential power model).

```{r MCMC}
# Log-normal
LN <- MCMC_LN(N = 1000, thin = 20, burn = 40, Time = Time,
              Cens = Cens, X = cancer[,3:11])

# Log-student's T
LST <- MCMC_LST(N = 1000, thin = 20, burn = 40 , Time = Time, Cens = Cens,
                X = cancer[,3:11])

# Log-Laplace
LLAP <- MCMC_LLAP(N = 1000, thin = 20, burn = 40,  Time = Time, Cens = Cens,
                  X = cancer[,3:11])

#Log-exponential power
LEP <- MCMC_LEP(N = 1000, thin = 20, burn = 40, Time = Time, Cens = Cens,
                X = cancer[,3:11])

#Log-logistic
LLOG <- MCMC_LLOG(N = 1000, thin = 20, burn = 40, Time = Time, Cens = Cens,
                  X = cancer[,3:11])
```

### Diagnostics

After generating MCMC chains, their suitability should be assessed. This can be 
done, in part, via the `Trace_Plot` function included with `BASSLINE` which will
plot a chain for a variable across (non-discarded) iterations. We will 
investigate the chain for $\beta_1$ from the log-normal model. 

```{r Trace Plot, fig.retina=2}
Trace_plot(1, LN)
```

For additional analysis of chains, the `coda` package is recommended which
offers many different functions:  

```{r ACF plot, message = F, fig.retina= 2, fig.height=9}
library(coda)
plot(as.mcmc(LN[,1:10])) 
```

ACF plots are also available via `coda`:

```{r}
acfplot(as.mcmc(LN[,1:10]))
```

## Deviance Information Criterion

The deviance information criterion (DIC), a hierarchical
modeling generalization of the Akaike information criterion [@DIC], can be
easily calculated for the 5 models.

If the deviance is defined as 

 $$D\left(\theta, y\right) = -2 \log\left(f\left(y|\theta\right)\right)$$
 
 then
 
 $$ DIC = D\left(\bar{\theta}\right) + 2p_D $$
 
 Where $p_D$ is the number of effective parameters. 
 
Each member of the log-normal family has a function to calculate DIC. We will 
present an example using the log-normal distribution and compare two models with
differing covariates (all of the available covariates vs. only the last 4
covariates).  
 
```{r DIC}
LN <- MCMC_LN(N = 1000, thin = 20, burn = 40, Time = Time,
              Cens = Cens, X = cancer[,3:11])
DIC_LN(Time = Time, Cens = Cens, X = cancer[,3:11], chain = LN)

# Reduced model
LN.2 <- MCMC_LN(N = 1000, thin = 20, burn = 40, Time = Time,
              Cens = Cens, X = cancer[,8:11])
DIC_LN(Time = Time, Cens = Cens, X = cancer[,8:11], chain = LN.2)
```

## Log-Marginal Likelihood

The log-marginal likelihood is given by:

$$ \log\left(m\left(t\right)\right) = 
\log\left(  
\int_{-\infty}^{\infty}
\int_{0}^{\infty}
\int_\Theta 
f\left(t \mid \beta, \sigma^2, \theta \right)
\pi\left(\beta, \sigma^2, \theta\right)
d\beta \: d\sigma^2 \: d\theta \right)$$

And can be easily estimated using the supplied function for each distribution
which is based on the work of Chib [@chib] and Chib and Jeliaskov [@chib2].
The function will return a list which includes the log-marginal likelihood
estimation, the log-likelihood ordinate, the log-prior ordinate, and the
log-posterior ordinates. Messages detailing the progress of the algorithm are
provided to the user.

```{r LML}
LML_LN(thin = 20, Time, Cens = Cens, X = cancer[,3:11], chain = LN) 
```


## Leave-One-Out Cross-Validation Analysis 

Leave-one-out cross-validation analysis is also available for all 5 of the 
supported distributions. The functions returns matrices with $n$ rows.

Its first column contains the logarithm of the conditional predictive ordinate
(CPO) [@cpo]. Larger values of the CPO indicate better predictive accuracy.

The second and third columns contain the Kullback–Leibler (KL) divergence
between $$ \pi\left(\beta, \sigma^2, \theta \mid t_{-i}\right)$$ and 
$$ \pi\left(\beta, \sigma^2, \theta \mid t\right)$$ and its calibration index
$p_i$ [@CI] respectively. The later is used in order to evaluate the existence
of influential observations. If $p_i$ is substantially larger than
0.5, observation $i$ is declared influential. Suggested cut-off values are
0.8 and 0.9.

```{r Case Deletion}
LN.CD <- CaseDeletion_LN(Time, Cens = Cens, X = cancer[,3:11], chain = LN)
```

```{r}
knitr::kable(LN.CD[1:5,])
```

It is sensible to report the number of observations which are deemed
influential for a given cutoff which can be easily found by the user. 

```{r}
sum(LN.CD[,3] > 0.8)
```


## Outlier Detection

Outlier detection is available for a specified observation (obs) for the 
log-student's T, log-Laplace, log-logistic, log-exponential power models.
This returns a unique number corresponding to the Bayes Factor associated with
the test $M_0 : \Lambda_\text{obs} = \lambda_\text{ref}$ versus
$M_0 : \Lambda_\text{obs} \neq \lambda_\text{ref}$ (with all other
$\Lambda_j,  j \neq \text{obs}$ free). The value of $\lambda_\text{ref}$ is 
required as input.

The recommended value of $\lambda_\text{ref}$ is 1 with the exception of the
log-logistic model where we recommend 0.4 instead. The calculations which 
support these recommendations can be found in the original paper
[@Vallejos2015]. 


```{r}
OD.LST <- rep(0, 5)
for(i in 1 : 5)
{
	OD.LST[i] <- BF_lambda_obs_LST(N = 100, thin = 20 , burn = 1, ref = 1,
	                               obs = i, Time = Time, Cens = Cens,
	                               X = cancer[,3:11], chain = LST)
}
```

## References
