## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

# Checks whether CmdStan is available and the model-fitting chunks can be
# evaluated only when CmdStan is present. Otherwise the pre-computed fits in
# inst/extdata are used.
cmdstan_ok <- !identical(Sys.getenv("DIPPER_NO_CMDSTAN"), "true") &&
    requireNamespace("instantiate", quietly = TRUE) &&
    instantiate::stan_cmdstan_exists()

cache_file <- function(x) {
    system.file("extdata", x, package = "DiPPER")
}

have_cache <- nzchar(cache_file("fit_hintikka.rds")) &&
    nzchar(cache_file("fit_vatanen.rds"))

# Run the DiPPER models?
run_fit <- cmdstan_ok

# Can results be shown?
show_res <- cmdstan_ok || have_cache

## ----install_stan, eval=FALSE-------------------------------------------------
# # Install cmdstanr package
# install.packages("cmdstanr",
#                  repos = c("https://stan-dev.r-universe.dev",
#                            getOption("repos")))
# 
# # Check C++ toolchain and install the backend
# cmdstanr::check_cmdstan_toolchain(fix = TRUE)
# cmdstanr::install_cmdstan()

## ----install_bioc, eval=FALSE-------------------------------------------------
# # Install BiocManager if it is not already installed
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# 
# BiocManager::install("DiPPER")

## ----quick_start, eval=run_fit------------------------------------------------
# library(DiPPER)
# 
# # Load dataset (TreeSummarizedExperiment object)
# data("tse_hintikka")
# 
# # Run DiPPER
# fit <- dipper(
#     tse = tse_hintikka,
#     assay.type = "counts",
#     data.type = "counts",
#     formula = ~ Fat + XOS,
#     read.depth = TRUE
# )

## ----quick_start_cached, eval=!run_fit && have_cache, echo=FALSE--------------
# CmdStan is not available on this machine: use the pre-computed fit instead
library(DiPPER)
fit <- readRDS(cache_file("fit_hintikka.rds"))

## ----quick_print, eval=show_res-----------------------------------------------
print(fit)

## ----quick_results, eval=show_res---------------------------------------------
# Extract the results
res <- summary(fit)

# Print the first seven rows
res[1:7, ]

## ----quick_plot, eval=show_res, out.width="100%"------------------------------
plot(fit)

## ----load_vatanen-------------------------------------------------------------
library(DiPPER)

# Load dataset (TreeSummarizedExperiment object, agglomerated to genus level)
data("VatanenT_2016_subset")

## ----vatanen_metadata---------------------------------------------------------
VatanenT_2016_subset |>
    SummarizedExperiment::colData() |>
    head()

## ----vatanen_abundance--------------------------------------------------------
vatanen_assay <- SummarizedExperiment::assay(VatanenT_2016_subset,
                                             "relative_abundance")
vatanen_assay[c("Bifidobacterium", "Escherichia", "Veillonella"), 1:5]

## ----assay_vatanen, eval=run_fit----------------------------------------------
# fit_long <- dipper(
#     tse = VatanenT_2016_subset,
#     assay.type = "relative_abundance", # Name of the assay
#     data.type = "relabundance",
#     formula = ~ age_point + antibiotics + gender + (1 | subject_id),
#     read.depth = FALSE
# )

## ----assay_vatanen_cached, eval=!run_fit && have_cache, echo=FALSE------------
# CmdStan is not available on this machine: use the pre-computed fit instead
fit_long <- readRDS(cache_file("fit_vatanen.rds"))

## ----plot_vatanen, eval=show_res, out.width="100%"----------------------------
# Print the details of the run
print(fit_long)

# Extract the results for age_point
res_long <- summary(fit_long)
res_long[1:7, ]

# Create a forest plot for significant taxa
plot(fit_long)

## ----session_info-------------------------------------------------------------
sessionInfo()