## ----eval=FALSE---------------------------------------------------------------
#  if (!requireNamespace("BiocManager")) install.packages("BiocManager")
#  
#  # Step 1
#  BiocManager::install("sccomp")
#  
#  # Step 2
#  install.packages("cmdstanr", repos = c("https://stan-dev.r-universe.dev/", getOption("repos")))
#  
#  # Step 3
#  cmdstanr::check_cmdstan_toolchain(fix = TRUE) # Just checking system setting
#  cmdstanr::install_cmdstan()

## ----eval=FALSE---------------------------------------------------------------
#  
#  # Step 1
#  devtools::install_github("MangiolaLaboratory/sccomp")
#  
#  # Step 2
#  install.packages("cmdstanr", repos = c("https://stan-dev.r-universe.dev/", getOption("repos")))
#  
#  # Step 3
#  cmdstanr::check_cmdstan_toolchain(fix = TRUE) # Just checking system setting
#  cmdstanr::install_cmdstan()

## ----message=FALSE, warning=FALSE---------------------------------------------
library(dplyr)
library(sccomp)
library(ggplot2)
library(forcats)
library(tidyr)
data("seurat_obj")
data("sce_obj")
data("counts_obj")

## ----eval=FALSE---------------------------------------------------------------
#  
#  sccomp_result =
#    sce_obj |>
#    sccomp_estimate(
#      formula_composition = ~ type,
#      .sample =  sample,
#      .cell_group = cell_group,
#      cores = 1
#    ) |>
#    sccomp_remove_outliers(cores = 1) |> # Optional
#    sccomp_test()
#  

## ----message=FALSE, warning=FALSE, eval = instantiate::stan_cmdstan_exists()----
#  
#  sccomp_result =
#    counts_obj |>
#    sccomp_estimate(
#      formula_composition = ~ type,
#      .sample = sample,
#      .cell_group = cell_group,
#      .count = count,
#      cores = 1, verbose = FALSE
#    ) |>
#    sccomp_remove_outliers(cores = 1, verbose = FALSE) |> # Optional
#    sccomp_test()

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  sccomp_result

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  sccomp_result |>
#    sccomp_boxplot(factor = "type")

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  sccomp_result |>
#    plot_1D_intervals()

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  sccomp_result |>
#    plot_2D_intervals()

## ----out.height="200%", eval=FALSE--------------------------------------------
#  sccomp_result |> plot()

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  
#  res =
#      seurat_obj |>
#      sccomp_estimate(
#        formula_composition = ~ type + continuous_covariate,
#        .sample = sample, .cell_group = cell_group,
#        cores = 1, verbose=FALSE
#      )
#  
#  res

## -----------------------------------------------------------------------------
seurat_obj[[]] |> as_tibble()

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  res =
#    seurat_obj |>
#    sccomp_estimate(
#      formula_composition = ~ type + (1 | group__),
#      .sample = sample,
#      .cell_group = cell_group,
#      bimodal_mean_variability_association = TRUE,
#      cores = 1, verbose = FALSE
#    )
#  
#  res

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  res =
#    seurat_obj |>
#    sccomp_estimate(
#        formula_composition = ~ type + (type | group__),
#        .sample = sample,
#        .cell_group = cell_group,
#        bimodal_mean_variability_association = TRUE,
#        cores = 1, verbose = FALSE
#      )
#  
#  res

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  res =
#    seurat_obj |>
#    sccomp_estimate(
#        formula_composition = ~ type + (type | group__) + (1 | group2__),
#        .sample = sample,
#        .cell_group = cell_group,
#        bimodal_mean_variability_association = TRUE,
#        cores = 1, verbose = FALSE
#      )
#  
#  res

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  sccomp_result |>
#     sccomp_proportional_fold_change(
#       formula_composition = ~  type,
#       from =  "healthy",
#       to = "cancer"
#      ) |>
#    select(cell_group, statement)

## ----message=FALSE, warning=FALSE, eval = instantiate::stan_cmdstan_exists()----
#  
#  seurat_obj |>
#    sccomp_estimate(
#      formula_composition = ~ 0 + type,
#      .sample = sample,
#      .cell_group = cell_group,
#      cores = 1, verbose = FALSE
#    ) |>
#    sccomp_test( contrasts =  c("typecancer - typehealthy", "typehealthy - typecancer"))
#  
#  

## ----message=FALSE, warning=FALSE, eval = instantiate::stan_cmdstan_exists()----
#  library(loo)
#  
#  # Fit first model
#  model_with_factor_association =
#    seurat_obj |>
#    sccomp_estimate(
#      formula_composition = ~ type,
#      .sample =  sample,
#      .cell_group = cell_group,
#      inference_method = "hmc",
#      enable_loo = TRUE
#    )
#  
#  # Fit second model
#  model_without_association =
#    seurat_obj |>
#    sccomp_estimate(
#      formula_composition = ~ 1,
#      .sample =  sample,
#      .cell_group = cell_group,
#      inference_method = "hmc",
#      enable_loo = TRUE
#    )
#  
#  # Compare models
#  loo_compare(
#     attr(model_with_factor_association, "fit")$loo(),
#     attr(model_without_association, "fit")$loo()
#  )
#  

## ----message=FALSE, warning=FALSE, eval = instantiate::stan_cmdstan_exists()----
#  
#  res =
#    seurat_obj |>
#    sccomp_estimate(
#      formula_composition = ~ type,
#      formula_variability = ~ type,
#      .sample = sample,
#      .cell_group = cell_group,
#      cores = 1, verbose = FALSE
#    )
#  
#  res
#  

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  plots = res |> sccomp_test() |> plot()
#  
#  plots$credible_intervals_1D

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  plots$credible_intervals_2D

## ----eval = instantiate::stan_cmdstan_exists()--------------------------------
#  
#  library(cmdstanr)
#  library(posterior)
#  library(bayesplot)
#  
#  # Assuming res contains the fit object from cmdstanr
#  fit <- res |> attr("fit")
#  
#  # Extract draws for 'beta[2,1]'
#  draws <- as_draws_array(fit$draws("beta[2,1]"))
#  
#  # Create a traceplot for 'beta[2,1]'
#  mcmc_trace(draws, pars = "beta[2,1]")
#  

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