diff --git a/NAMESPACE b/NAMESPACE
index 436eb956..79782749 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -8,8 +8,8 @@ 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)
@@ -17,8 +17,8 @@ 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
index 7e9a8275..2c026484 100644
--- 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
index e943a0ab..de3aae48 100644
--- 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
index 38b05144..c96fd07d 100644
--- 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
index cbb6797b..aaf6eba4 100644
--- 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
index ae06fe9f..feaab835 100644
--- 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
+
+```
\ No newline at end of file
diff --git a/man/test_contrasts.Rd b/man/sccomp_test.Rd
similarity index 92%
rename from man/test_contrasts.Rd
rename to man/sccomp_test.Rd
index f6fef968..d3998042 100644
--- a/man/test_contrasts.Rd
+++ b/man/sccomp_test.Rd
@@ -1,10 +1,10 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/methods.R
-\name{test_contrasts}
-\alias{test_contrasts}
-\title{test_contrasts}
+\name{sccomp_test}
+\alias{sccomp_test}
+\title{sccomp_test}
 \usage{
-test_contrasts(
+sccomp_test(
   .data,
   contrasts = NULL,
   percent_false_positive = 5,
@@ -41,6 +41,6 @@ data("counts_obj")
     cores = 1
   ) |>
 
-  test_contrasts("typecancer - typebenign")
+  sccomp_test("typecancer - typebenign")
 
 }
diff --git a/tests/testthat/test-sccomp_.R b/tests/testthat/test-sccomp_.R
index d368dd0b..79a3effa 100644
--- 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()
 
 })