## ----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---------------------------------------------------------------
#  
#  single_cell_object |>
#    sccomp_glm(
#      formula_composition = ~ type,
#      .sample =  sample,
#      .cell_group = cell_group,
#      bimodal_mean_variability_association = TRUE,
#      cores = 1
#    )
#  

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

counts_obj |>
  sccomp_glm( 
    formula_composition = ~ type, 
    .sample = sample,
    .cell_group = cell_group,
    .count = count, 
    bimodal_mean_variability_association = TRUE,
    cores = 1 
  )


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

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


## ---- eval=FALSE, message=FALSE, warning=FALSE--------------------------------
#  library(loo)
#  
#  # Fit first model
#  model_with_factor_association =
#    seurat_obj |>
#    sccomp_glm(
#      formula_composition = ~ type,
#      .sample =  sample,
#      .cell_group = cell_group,
#      check_outliers = FALSE,
#      bimodal_mean_variability_association = TRUE,
#      cores = 1,
#      enable_loo = TRUE
#    )
#  
#  # Fit second model
#  model_without_association =
#    seurat_obj |>
#    sccomp_glm(
#      formula_composition = ~ 1,
#      .sample =  sample,
#      .cell_group = cell_group,
#      check_outliers = FALSE,
#      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_glm( 
    formula_composition = ~ type, 
    formula_variability = ~ type,
    .sample = sample,
    .cell_group = cell_group,
    bimodal_mean_variability_association = TRUE,
    cores = 1 
  )

res


## ---- out.height="200%"-------------------------------------------------------
plots = plot_summary(res) 

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

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

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

## -----------------------------------------------------------------------------
plots = plot_summary(res)

plots$credible_intervals_1D

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

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