---
title: "An introduction to BASiCStan"
subtitle: "A Stan implementation of BASiCS"
output: rmarkdown::html_vignette
vignette: >
    %\VignetteIndexEntry{An introduction to BASiCStan}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
---

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

# Outline

This package implements alternative inference methods for
[BASiCS](https://bioconductor.org/packages/release/bioc/html/BASiCS.html).
The original package uses adaptive Metropolis within Gibbs sampling, while
BASiCStan uses Stan to enable the use of maximum *a posteriori* estimation,
variational inference, and Hamiltonian Monte Carlo. These each
have advantages for different use cases.

# Use

To use `BASiCStan`, we can first use `BASICS_MockSCE()` to generate an
toy example dataset. We will also set a seed for reproducibility.

```{r setup}
suppressPackageStartupMessages(library("BASiCStan"))
set.seed(42)
sce <- BASiCS_MockSCE()
```

The interface for running MCMC using BASiCS and using alternative inference
methods using Stan is very similar. It is worth noting that the joint prior
between mean and over-dispersion parameters, corresponding to
`Regression = TRUE`, is the only supported mode in `BASiCStan()`.

```{r basics, results = "hide"}
amwg_fit <- BASiCS_MCMC(
    sce,
    N = 200,
    Thin = 10,
    Burn = 100,
    Regression = TRUE
)
stan_fit <- BASiCStan(sce, Method = "sampling", iter = 10)
```

The output of `BASiCStan()` is an object of class `BASiCS_Chain`, similar to
the output of `BASiCS_MCMC()`. Therefore, you could use these as you would an
object generated using Metropolis within Gibbs sampling.
For example, we can plot the
relationship between mean and over-dispersion estimated using the joint
regression prior:

```{r}
BASiCS_ShowFit(amwg_fit)
BASiCS_ShowFit(stan_fit)
```

# Stan MCMC diagnostics

Using Stan has advantages outside of the variety of inference engines available.
By returning a Stan object that we can later convert to a `BASiCS_Chain` object,
we can leverage an even broader set of functionality.
For example, Stan has the ability to easily generate MCMC diagnostics
using simple functions. For example, `summary()` outputs a number of useful
per-chain and across-chain diagnostics:

```{r}
stan_obj <- BASiCStan(sce, ReturnBASiCS = FALSE, Method = "sampling", iter = 10)
head(summary(stan_obj)$summary)
```

We can then convert this object to a `BASiCS_Chain` and carry on a workflow as
before:

```{r}
Stan2BASiCS(stan_obj)
```

# Session info

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