## ----eval=FALSE---------------------------------------------------------------
#  if (!requireNamespace("BiocManager")) install.packages("BiocManager")
#  BiocManager::install("sccomp")

## ----eval=FALSE---------------------------------------------------------------
#  devtools::install_github("stemangiola/sccomp")

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

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

## ----message=FALSE, warning=FALSE---------------------------------------------

sccomp_result = 
  counts_obj |>
  sccomp_estimate( 
    formula_composition = ~ type, 
    .sample = sample,
    .cell_group = cell_group,
    .count = count, 
    bimodal_mean_variability_association = TRUE,
    cores = 1 
  ) |> 
  sccomp_remove_outliers(cores = 1) |> # Optional
  sccomp_test()


## ----out.height="200%"--------------------------------------------------------
plots = sccomp_result |> plot() 

## -----------------------------------------------------------------------------
plots$boxplot

## -----------------------------------------------------------------------------
plots$credible_intervals_1D

## -----------------------------------------------------------------------------
plots$credible_intervals_2D

## ----message=FALSE, warning=FALSE---------------------------------------------

seurat_obj |>
  sccomp_estimate( 
    formula_composition = ~ 0 + type, 
    .sample = sample,
    .cell_group = cell_group, 
    bimodal_mean_variability_association = TRUE,
    cores = 1 
  ) |> 
  sccomp_test( contrasts =  c("typecancer - typehealthy", "typehealthy - typecancer"))



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

## ----message=FALSE, warning=FALSE---------------------------------------------

res = 
  seurat_obj |>
  sccomp_estimate( 
    formula_composition = ~ type, 
    formula_variability = ~ type,
    .sample = sample,
    .cell_group = cell_group,
    bimodal_mean_variability_association = TRUE,
    cores = 1 
  )

res


## -----------------------------------------------------------------------------
res %>% attr("fit") %>% rstan::traceplot("beta[2,1]")

## -----------------------------------------------------------------------------
plots = res |> sccomp_test() |> plot()

plots$credible_intervals_1D

## -----------------------------------------------------------------------------
plots$credible_intervals_2D

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