In Alanocallaghan/BASiCStan: Stan implementation of BASiCS

knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

Outline

This package implements alternative inference methods for BASiCS. 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.

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().

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:

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:

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:

Stan2BASiCS(stan_obj)

Session info

sessionInfo()

Alanocallaghan/BASiCStan documentation built on Feb. 19, 2024, 10:40 p.m.