typo and add multilevel to README

MangiolaLaboratory · Sep 4, 2023 · 1611e48 · 1611e48
1 parent 76d1a19
commit 1611e48
Show file tree

Hide file tree

Showing 8 changed files with 65 additions and 32 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -8,17 +8,17 @@ S3method(sccomp_estimate,data.frame)
 S3method(sccomp_predict,data.frame)
 S3method(sccomp_remove_outliers,sccomp_tbl)
 S3method(sccomp_replicate,data.frame)
+S3method(sccomp_test,data.frame)
 S3method(simulate_data,data.frame)
-S3method(test_contrasts,data.frame)
 export(plot_summary)
 export(remove_unwanted_variation)
 export(sccomp_boxplot)
 export(sccomp_estimate)
 export(sccomp_predict)
 export(sccomp_remove_outliers)
 export(sccomp_replicate)
+export(sccomp_test)
 export(simulate_data)
-export(test_contrasts)
 import(Rcpp)
 import(SeuratObject)
 import(dplyr)
@@ -66,6 +66,7 @@ importFrom(purrr,map_int)
 importFrom(purrr,map_lgl)
 importFrom(purrr,pmap)
 importFrom(purrr,prepend)
+importFrom(purrr,reduce)
 importFrom(purrr,when)
 importFrom(readr,read_file)
 importFrom(rlang,":=")

diff --git a/R/functions_multi_beta_binomial.R b/R/functions_multi_beta_binomial.R
@@ -300,7 +300,7 @@ sccomp_glm_data_frame_counts = function(.data,
     add_attr(parse_formula(formula_composition), "factors" ) |> 
 
     # Print estimates
-    test_contrasts() |>
+    sccomp_test() |>
 
     # drop hypothesis testing as the estimation exists without probabilities.
     # For hypothesis testing use sccomp_test
@@ -424,7 +424,7 @@ multi_beta_binomial_glm = function(.data,
     add_attr(formula_composition, "formula_composition") |>
     add_attr(formula_variability, "formula_variability") |>
 
-    test_contrasts(
+    sccomp_test(
       contrasts = contrasts,
       percent_false_positive = percent_false_positive,
       test_composition_above_logit_fold_change = test_composition_above_logit_fold_change

diff --git a/R/methods.R b/R/methods.R
@@ -733,7 +733,7 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate,
 
 
     # Print estimates
-    test_contrasts() |>
+    sccomp_test() |>
 
     # drop hypothesis testing as the estimation exists without probabilities.
     # For hypothesis testing use sccomp_test
@@ -748,7 +748,7 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate,
 
 }
 
-#' test_contrasts
+#' sccomp_test
 #'
 #' @description This function test contrasts from a sccomp result.
 #'
@@ -775,14 +775,14 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate,
 #'     cores = 1
 #'   ) |>
 #'
-#'   test_contrasts("typecancer - typebenign")
+#'   sccomp_test("typecancer - typebenign")
 #'
-test_contrasts <- function(.data,
+sccomp_test <- function(.data,
                            contrasts = NULL,
                            percent_false_positive = 5,
                            test_composition_above_logit_fold_change = 0.2,
                            pass_fit = TRUE) {
-  UseMethod("test_contrasts", .data)
+  UseMethod("sccomp_test", .data)
 }
 
 
@@ -792,7 +792,7 @@ test_contrasts <- function(.data,
 #' @export
 #' 
 #'
-test_contrasts.data.frame = function(.data,
+sccomp_test.data.frame = function(.data,
                                      contrasts = NULL,
                                      percent_false_positive = 5,
                                      test_composition_above_logit_fold_change = 0.2,

diff --git a/R/utilities.R b/R/utilities.R
@@ -668,7 +668,7 @@ get_random_intercept_design2 = function(.data_, .sample, formula_composition ){
             mutate_if(is.integer, ~1) |> 
             pivot_longer(-!!.sample, names_to = "factor"),
 
-          by = join_by(sample_, factor)
+          by = join_by(!!.sample, factor)
         ) |> 
 
          # Create unique name

diff --git a/inst/stan/glm_multi_beta_binomial.stan b/inst/stan/glm_multi_beta_binomial.stan
@@ -133,8 +133,9 @@ transformed data{
   Q_ast = qr_thin_Q(X) * sqrt(N - 1);
   R_ast_inverse = inverse(qr_thin_R(X) / sqrt(N - 1));
   // If I get crazy diagonal matrix omit it
-  if(max(R_ast_inverse)>1000 || N_random_intercepts>0){
-    print("sccomp says: The QR deconposition resulted in extreme values, probably for the correlation structure of your design matrix. Omitting QR decomposition.");
+  if(N_random_intercepts>0) { 
+    if(max(R_ast_inverse)>1000 )
+      print("sccomp says: The QR deconposition resulted in extreme values, probably for the correlation structure of your design matrix. Omitting QR decomposition.");
     Q_ast = X;
     R_ast_inverse = diag_matrix(rep_vector(1.0, C));
   }

diff --git a/man/fragments/intro.Rmd b/man/fragments/intro.Rmd
@@ -63,7 +63,7 @@ single_cell_object |>
     bimodal_mean_variability_association = TRUE,
     cores = 1 
   ) |> 
-	sccomp_test()
+	sccomp_test(test_composition_above_logit_fold_change = 0.2)
 
 ```
 
@@ -72,14 +72,15 @@ single_cell_object |>
 ```{r, message=FALSE, warning=FALSE}
 
 counts_obj |>
-  sccomp_glm( 
+  sccomp_estimate( 
     formula_composition = ~ type, 
     .sample = sample,
     .cell_group = cell_group,
     .count = count, 
     bimodal_mean_variability_association = TRUE,
     cores = 1 
-  )
+  ) |> 
+	sccomp_test(test_composition_above_logit_fold_change = 0.2)
 
 ```
 
@@ -90,14 +91,17 @@ Of the output table, the estimate columns start with the prefix `c_` indicate `c
 ```{r, message=FALSE, warning=FALSE}
 
 seurat_obj |>
-  sccomp_glm( 
+  sccomp_estimate( 
     formula_composition = ~ 0 + type, 
-    contrasts =  c("typecancer - typehealthy", "typehealthy - typecancer"),
     .sample = sample,
     .cell_group = cell_group, 
     bimodal_mean_variability_association = TRUE,
     cores = 1 
-  )
+  ) |> 
+	sccomp_test(
+	  contrasts =  c("typecancer - typehealthy", "typehealthy - typecancer"),
+	  test_composition_above_logit_fold_change = 0.2
+	)
 
 ```
 
@@ -115,11 +119,10 @@ library(loo)
 # Fit first model
 model_with_factor_association = 
   seurat_obj |>
-  sccomp_glm( 
+  sccomp_estimate( 
     formula_composition = ~ type, 
     .sample =  sample, 
     .cell_group = cell_group, 
-    check_outliers = FALSE, 
     bimodal_mean_variability_association = TRUE,
     cores = 1, 
     enable_loo = TRUE
@@ -128,11 +131,10 @@ model_with_factor_association =
 # Fit second model
 model_without_association = 
   seurat_obj |>
-  sccomp_glm( 
+  sccomp_estimate( 
     formula_composition = ~ 1, 
     .sample =  sample, 
     .cell_group = cell_group, 
-    check_outliers = FALSE, 
     bimodal_mean_variability_association = TRUE,
     cores = 1 , 
     enable_loo = TRUE
@@ -154,14 +156,17 @@ We can model the cell-group variability also dependent on the type, and so test
 
 res = 
   seurat_obj |>
-  sccomp_glm( 
+  sccomp_estimate( 
     formula_composition = ~ type, 
     formula_variability = ~ type,
     .sample = sample,
     .cell_group = cell_group,
     bimodal_mean_variability_association = TRUE,
     cores = 1 
-  )
+  ) |> 
+	sccomp_test(
+	  test_composition_above_logit_fold_change = 0.2
+	)
 
 res
 
@@ -227,3 +232,29 @@ This plot is provided only if differential variability has been tested. The diff
 ```{r}
 plots$credible_intervals_2D
 ```
+
+# Multilevel modelling
+
+`sccomp` is cabable of estimating population (i.e. fixed) and group (i.e. random) effects. The formulation is analogous to the `lme4` package and `brms`.
+
+!! For now, only one grouping is allowed (e.g. group2__).
+
+```{r, message=FALSE, warning=FALSE}
+
+res = 
+  seurat_obj |>
+  sccomp_estimate( 
+    formula_composition = ~ type + (type | group2__), 
+    formula_variability = ~ type,
+    .sample = sample,
+    .cell_group = cell_group,
+    bimodal_mean_variability_association = TRUE,
+    cores = 1 
+  ) |> 
+	sccomp_test(
+	  test_composition_above_logit_fold_change = 0.2
+	)
+
+res
+
+```
diff --git a/man/test_contrasts.Rd → man/sccomp_test.Rd b/man/test_contrasts.Rd → man/sccomp_test.Rd
diff --git a/tests/testthat/test-sccomp_.R b/tests/testthat/test-sccomp_.R
@@ -465,7 +465,7 @@ test_that("test constrasts",{
 
   new_test =
     estimate |>
-    test_contrasts() |>
-    test_contrasts()
+    sccomp_test() |>
+    sccomp_test()
 
 })