inst/doc/banocc-vignette.R

In banocc: Bayesian ANalysis Of Compositional Covariance

## ---- eval=FALSE, echo=FALSE--------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly=TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("BAnOCC")

## ----load, eval=TRUE----------------------------------------------------------
library(banocc)

## ----run-banocc-help, eval=FALSE----------------------------------------------
#  ?run_banocc
#  ?get_banocc_output

## ----rerun, cache=TRUE, echo=FALSE--------------------------------------------
rerun <- 0

## ----basic-run-banocc, eval=TRUE, cache=TRUE, dependson=c('rerun'), results="hide"----
data(compositions_null)
compiled_banocc_model <- rstan::stan_model(model_code = banocc::banocc_model) 
b_fit     <- banocc::run_banocc(C = compositions_null, compiled_banocc_model=compiled_banocc_model)
b_output <- banocc::get_banocc_output(banoccfit=b_fit)

## ----input-hyperparameters, eval=TRUE, cache=TRUE, dependson=c('rerun'), results="hide"----
p <- ncol(compositions_null)
b_fit_hp <- banocc::run_banocc(C = compositions_null, 
                               compiled_banocc_model = compiled_banocc_model,
                               n = rep(0, p),
                               L = 10 * diag(p), 
                               a = 0.5,
                               b = 0.01)

## ----sampling-sampling, eval=TRUE, cache=TRUE, dependson=c('rerun'), results="hide"----
b_fit_sampling <- banocc::run_banocc(C = compositions_null,
                                     compiled_banocc_model = compiled_banocc_model,
                                     chains = 2,
                                     iter = 11,
                                     warmup = 5,
                                     thin = 2)

## ----sampling-cores, eval=FALSE-----------------------------------------------
#  # This code is not run
#  b_fit_cores <- banocc::run_banocc(C = compositions_null,
#                                    compiled_banocc_model = compiled_banocc_model,
#                                    chains = 2,
#                                    cores = 2)

## ----sampling-init, eval=TRUE, cache=TRUE, dependson=c('rerun'), results="hide"----
init <- list(list(m = rep(0, p),
                  O = diag(p),
                  lambda = 0.02),
             list(m = runif(p),
                  O = 10 * diag(p),
                  lambda = runif(1, 0.1, 2)))
b_fit_init <- banocc::run_banocc(C = compositions_null,
                                 compiled_banocc_model = compiled_banocc_model,
                                 chains = 2,
                                 init = init)

## ---- eval=FALSE--------------------------------------------------------------
#  ?stan

## ----output-ci, eval=TRUE, cache=TRUE, dependson=c('rerun'), results="hide"----
# Get 90% credible intervals
b_out_90 <- banocc::get_banocc_output(banoccfit=b_fit,
                                      conf_alpha = 0.1)
# Get 99% credible intervals
b_out_99 <- banocc::get_banocc_output(banoccfit=b_fit,
                                      conf_alpha = 0.01)

## ----eval-convergence, eval=TRUE, cache=TRUE, dependson=c('rerun'), results="hide"----
# Default is to evaluate convergence
b_out_ec <- banocc::get_banocc_output(banoccfit=b_fit)

# This can be turned off using `eval_convergence`
b_out_nec <- banocc::get_banocc_output(banoccfit=b_fit,
                                       eval_convergence = FALSE)

## ----show-eval-convergence, eval=TRUE-----------------------------------------
# Iterations are too few, so estimates are missing
b_out_ec$Estimates.median

# Convergence was not evaluated, so estimates are not missing
b_out_nec$Estimates.median

## ----output-extra, eval=TRUE, cache=TRUE, dependson=c('rerun'), results="hide"----
# Get the smallest credible interval width that includes zero
b_out_min_width <- banocc::get_banocc_output(banoccfit=b_fit,
                                             get_min_width = TRUE)

# Get the scaled neighborhood criterion
b_out_snc <- banocc::get_banocc_output(banoccfit=b_fit,
                                       calc_snc = TRUE)

## ----traceplot, eval=TRUE, cache=TRUE-----------------------------------------
# The inverse covariances of feature 1 with all other features
rstan::traceplot(b_fit$Fit, pars=paste0("O[1,", 2:9, "]"))

## ----traceplot-warmup, eval=TRUE, cache=TRUE----------------------------------
# The inverse covariances of feature 1 with all other features, including warmup
rstan::traceplot(b_fit$Fit, pars=paste0("O[1,", 2:9, "]"),
                 inc_warmup=TRUE)

## ----Rhat, eval=TRUE, cache=TRUE----------------------------------------------
# This returns a named vector with the Rhat values for all parameters
rhat_all <- rstan::summary(b_fit$Fit)$summary[, "Rhat"]

# To see the Rhat values for the inverse covariances of feature 1
rhat_all[paste0("O[1,", 2:9, "]")]

## ---- eval=FALSE, echo=FALSE--------------------------------------------------
#  # The above plot is from Dropbox/hutlab/Emma/paper1/writeup/figures/supplemental/lambbda_behavior.png

## ----model, eval=FALSE--------------------------------------------------------
#  # This code is not run
#  cat(banocc::banocc_model)