From b6f979298f75def0d0f890de517eac8790579318 Mon Sep 17 00:00:00 2001
From: Stefano Mangiola <mangiolastefano@gmail.com>
Date: Wed, 27 May 2026 14:13:09 +0930
Subject: [PATCH 1/4] Update version to 2.1.33 and enhance Stan models with
 non-negative constraints for random effect parameters. This change ensures
 backward compatibility while improving the handling of random effects in the
 glm_multi_beta_binomial models.

---
 DESCRIPTION                                          | 2 +-
 inst/stan/glm_multi_beta_binomial.stan               | 9 ++++++++-
 inst/stan/glm_multi_beta_binomial_generate_data.stan | 6 ++++++
 inst/stan/glm_multi_beta_binomial_simulate_data.stan | 7 ++++++-
 4 files changed, 21 insertions(+), 3 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index da8a8cb7..0c7d39c4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: sccomp
 Type: Package
 Title: Differential Composition and Variability Analysis for Single-Cell Data
-Version: 2.1.32
+Version: 2.1.33
 Date: 2026-05-11
 Authors@R: c(person("Stefano", "Mangiola", email = "stefano.mangiola@unimelb.edu.au", role = c("aut", "cre")), person("Alexandra J.", "Roth-Schulze", role = "aut"), person("Marie", "Trussart", role = "aut"), person("Enrique", "Zozaya-Valdés", role = "aut"), person("Mengyao", "Ma", role = "aut"), person("Zijie", "Gao", role = "aut"), person("Alan F.", "Rubin", role = "aut"), person("Terence P.", "Speed", role = "aut"), person("Heejung", "Shim", role = "aut"), person("Anthony T.", "Papenfuss", role = "aut"))
 Description: Comprehensive R package for differential composition and variability analysis in single-cell RNA sequencing, CyTOF, and microbiome data. Provides robust Bayesian modeling with outlier detection, random effects, and advanced statistical methods for cell type proportion analysis. Features include probabilistic outlier identification, mixed-effect modeling, differential variability testing, and comprehensive visualization tools. Perfect for cancer research, immunology, developmental biology, and single-cell genomics applications.
diff --git a/inst/stan/glm_multi_beta_binomial.stan b/inst/stan/glm_multi_beta_binomial.stan
index e74600de..b09a4cdf 100755
--- a/inst/stan/glm_multi_beta_binomial.stan
+++ b/inst/stan/glm_multi_beta_binomial.stan
@@ -430,7 +430,14 @@ parameters{
 
   // Shared hyperprior scalars (one mu, one sigma per slot)
   array[4 * (is_random_effect>0)] real random_effect_sigma_mu;
-  array[4 * (is_random_effect>0)] real random_effect_sigma_sigma;
+  // Constrained to be non-negative to break the sign-flip redundancy in
+  // sigma_vec[m] = exp((sigma_mu + sigma_sigma * sigma_raw[m]) / 3.0):
+  // (+sigma_sigma, +sigma_raw) and (-sigma_sigma, -sigma_raw) describe the
+  // same model, producing a bimodal posterior that HMC has to traverse
+  // through the degenerate region at 0 to mix. This is a reparameterisation,
+  // not a prior change: sigma_sigma = 0 (all cell types share one RE-SD) is
+  // still the mode of the half-normal prior below.
+  array[4 * (is_random_effect>0)] real<lower=0> random_effect_sigma_sigma;
 
   // For models with a single group (kept from the original design)
   array[is_random_effect>0] real zero_random_effect;
diff --git a/inst/stan/glm_multi_beta_binomial_generate_data.stan b/inst/stan/glm_multi_beta_binomial_generate_data.stan
index 4427d86c..8b3e3f23 100755
--- a/inst/stan/glm_multi_beta_binomial_generate_data.stan
+++ b/inst/stan/glm_multi_beta_binomial_generate_data.stan
@@ -117,6 +117,12 @@ parameters {
   array[M * (ncol_X_random_eff[4]>0)] cholesky_factor_corr[how_many_factors_in_random_design[4] * (ncol_X_random_eff[4]>0)] sigma_correlation_factor_4;
 
   array[4 * (is_random_effect>0)] real random_effect_sigma_mu;
+  // NOTE: deliberately *not* declared <lower=0> here even though the fitting
+  // model in glm_multi_beta_binomial.stan now constrains it. This program
+  // binds parameter values from previously-saved fits via generate_quantities;
+  // adding the constraint would reject fits produced before the sign-flip
+  // fix was introduced. The constraint only matters during sampling, so
+  // leaving this permissive is sound and preserves backward compatibility.
   array[4 * (is_random_effect>0)] real random_effect_sigma_sigma;
   array[is_random_effect>0] real zero_random_effect;
 }
diff --git a/inst/stan/glm_multi_beta_binomial_simulate_data.stan b/inst/stan/glm_multi_beta_binomial_simulate_data.stan
index a719643e..66f09ad8 100755
--- a/inst/stan/glm_multi_beta_binomial_simulate_data.stan
+++ b/inst/stan/glm_multi_beta_binomial_simulate_data.stan
@@ -34,7 +34,12 @@ parameters{
   
   // sd of random intercept
   array[is_random_effect>0] real random_effect_sigma_mu;
-  array[is_random_effect>0] real random_effect_sigma_sigma;
+  // <lower=0> mirrors the fitting model in glm_multi_beta_binomial.stan: it
+  // breaks the (+sigma_sigma, +sigma_raw) / (-sigma_sigma, -sigma_raw) sign
+  // redundancy so prior-predictive draws match the constrained parameterisation
+  // the fit actually uses. It is not a minimum-variance assumption: sigma_sigma
+  // = 0 (homogeneous RE-SDs across categories) remains the prior mode.
+  array[is_random_effect>0] real<lower=0> random_effect_sigma_sigma;
 
 	// Covariance
   array[M-1 * (is_random_effect>0)] vector[how_many_factors_in_random_design[1]]  random_effect_sigma_raw;

From 3d2e95a3b78042571b824899d6774d393a7b0c8c Mon Sep 17 00:00:00 2001
From: Stefano Mangiola <mangiolastefano@gmail.com>
Date: Mon, 1 Jun 2026 17:02:29 +1000
Subject: [PATCH 2/4] Fix sigma_sigma to 1 to avoid poor fit

Refactor Stan models to ensure backward compatibility with random effect parameters. Update comments for clarity and maintain consistency in handling sigma_sigma across models, fixing its value to 1.0 during sampling while retaining the original structure for saved fits.
---
 inst/stan/common_functions.stan               |  2 +-
 inst/stan/glm_multi_beta_binomial.stan        | 26 ++++++++-----------
 ...glm_multi_beta_binomial_generate_data.stan | 19 ++++++--------
 ...glm_multi_beta_binomial_simulate_data.stan |  6 +----
 4 files changed, 21 insertions(+), 32 deletions(-)

diff --git a/inst/stan/common_functions.stan b/inst/stan/common_functions.stan
index dc2a35df..2e28d89f 100644
--- a/inst/stan/common_functions.stan
+++ b/inst/stan/common_functions.stan
@@ -110,7 +110,7 @@
       array[,] int group_factor_indexes,                  // factor x group -> col index
       array[] vector raw_vec,                             // length ncol, each length M
       real sigma_mu,                                      // hierarchical hyperprior mean
-      real sigma_sigma,                                   // hierarchical hyperprior SD
+      real sigma_sigma,                                   // hierarchical hyperprior SD (callers pass 1.0; see fitting model)
       array[] vector sigma_raw,                           // length M, each length n_factors
       array[] matrix corr_chol                            // length M, each n_factors x n_factors
     ) {
diff --git a/inst/stan/glm_multi_beta_binomial.stan b/inst/stan/glm_multi_beta_binomial.stan
index b09a4cdf..2562d763 100755
--- a/inst/stan/glm_multi_beta_binomial.stan
+++ b/inst/stan/glm_multi_beta_binomial.stan
@@ -405,7 +405,8 @@ parameters{
   //   * sigma_correlation_factor_k: per-category Cholesky of correlation matrix
   //                                  (n_factors[k] x n_factors[k]; 1x1 = no LKJ work)
   //
-  // Hyperprior scalars sigma_mu / sigma_sigma are shared in length-4 arrays.
+  // Hyperprior scalar sigma_mu is shared (length-4 array). Between-category
+  // log-SD spread is fixed at sigma_sigma = 1.0 in build_re_block (see below).
   // ----------------------------------------------------------------------
 
   // Slot 1
@@ -428,16 +429,12 @@ parameters{
   array[M * (ncol_X_random_eff[4]>0)] vector[how_many_factors_in_random_design[4]]  random_effect_sigma_raw_4;
   array[M * (ncol_X_random_eff[4]>0)] cholesky_factor_corr[how_many_factors_in_random_design[4] * (ncol_X_random_eff[4]>0)] sigma_correlation_factor_4;
 
-  // Shared hyperprior scalars (one mu, one sigma per slot)
+  // Shared hyperprior scalars (one mu per slot)
   array[4 * (is_random_effect>0)] real random_effect_sigma_mu;
-  // Constrained to be non-negative to break the sign-flip redundancy in
-  // sigma_vec[m] = exp((sigma_mu + sigma_sigma * sigma_raw[m]) / 3.0):
-  // (+sigma_sigma, +sigma_raw) and (-sigma_sigma, -sigma_raw) describe the
-  // same model, producing a bimodal posterior that HMC has to traverse
-  // through the degenerate region at 0 to mix. This is a reparameterisation,
-  // not a prior change: sigma_sigma = 0 (all cell types share one RE-SD) is
-  // still the mode of the half-normal prior below.
-  array[4 * (is_random_effect>0)] real<lower=0> random_effect_sigma_sigma;
+  // Retained for backward compatibility with saved fits / R interfaces; not
+  // used in the likelihood (build_re_block is called with sigma_sigma = 1.0).
+  // Pinned so draws still expose the name at 1.0 without a prior.
+  array[4 * (is_random_effect>0)] real<lower=1, upper=1> random_effect_sigma_sigma;
 
   // For models with a single group (kept from the original design)
   array[is_random_effect>0] real zero_random_effect;
@@ -486,7 +483,7 @@ transformed parameters{
     random_effect_1 = build_re_block(
       M, n_groups[1], how_many_factors_in_random_design[1], ncol_X_random_eff[1],
       group_factor_indexes_for_covariance_1, raw_vec,
-      random_effect_sigma_mu[1], random_effect_sigma_sigma[1],
+      random_effect_sigma_mu[1], 1.0,
       random_effect_sigma_raw_1, sigma_correlation_factor_1
     );
   }
@@ -497,7 +494,7 @@ transformed parameters{
     random_effect_2 = build_re_block(
       M, n_groups[2], how_many_factors_in_random_design[2], ncol_X_random_eff[2],
       group_factor_indexes_for_covariance_2, raw_vec,
-      random_effect_sigma_mu[2], random_effect_sigma_sigma[2],
+      random_effect_sigma_mu[2], 1.0,
       random_effect_sigma_raw_2, sigma_correlation_factor_2
     );
   }
@@ -508,7 +505,7 @@ transformed parameters{
     random_effect_3 = build_re_block(
       M, n_groups[3], how_many_factors_in_random_design[3], ncol_X_random_eff[3],
       group_factor_indexes_for_covariance_3, raw_vec,
-      random_effect_sigma_mu[3], random_effect_sigma_sigma[3],
+      random_effect_sigma_mu[3], 1.0,
       random_effect_sigma_raw_3, sigma_correlation_factor_3
     );
   }
@@ -519,7 +516,7 @@ transformed parameters{
     random_effect_4 = build_re_block(
       M, n_groups[4], how_many_factors_in_random_design[4], ncol_X_random_eff[4],
       group_factor_indexes_for_covariance_4, raw_vec,
-      random_effect_sigma_mu[4], random_effect_sigma_sigma[4],
+      random_effect_sigma_mu[4], 1.0,
       random_effect_sigma_raw_4, sigma_correlation_factor_4
     );
   }
@@ -623,7 +620,6 @@ model{
   // ----------------------------------------------------------------------
   if (is_random_effect > 0) {
     random_effect_sigma_mu ~ std_normal();
-    random_effect_sigma_sigma ~ std_normal();
     zero_random_effect ~ std_normal();
   }
 
diff --git a/inst/stan/glm_multi_beta_binomial_generate_data.stan b/inst/stan/glm_multi_beta_binomial_generate_data.stan
index 8b3e3f23..50e4400a 100755
--- a/inst/stan/glm_multi_beta_binomial_generate_data.stan
+++ b/inst/stan/glm_multi_beta_binomial_generate_data.stan
@@ -117,13 +117,10 @@ parameters {
   array[M * (ncol_X_random_eff[4]>0)] cholesky_factor_corr[how_many_factors_in_random_design[4] * (ncol_X_random_eff[4]>0)] sigma_correlation_factor_4;
 
   array[4 * (is_random_effect>0)] real random_effect_sigma_mu;
-  // NOTE: deliberately *not* declared <lower=0> here even though the fitting
-  // model in glm_multi_beta_binomial.stan now constrains it. This program
-  // binds parameter values from previously-saved fits via generate_quantities;
-  // adding the constraint would reject fits produced before the sign-flip
-  // fix was introduced. The constraint only matters during sampling, so
-  // leaving this permissive is sound and preserves backward compatibility.
-  array[4 * (is_random_effect>0)] real random_effect_sigma_sigma;
+  // Retained for backward compatibility when binding saved draws; build_re_block
+  // uses sigma_sigma = 1.0 (same as the fitting model). Unconstrained so old fits
+  // with arbitrary stored values still load.
+  array[4 * (is_random_effect>0)] real<lower=0> random_effect_sigma_sigma;
   array[is_random_effect>0] real zero_random_effect;
 }
 transformed parameters {
@@ -155,7 +152,7 @@ transformed parameters {
     random_effect_1 = build_re_block(
       M, n_groups[1], how_many_factors_in_random_design[1], ncol_X_random_eff[1],
       group_factor_indexes_for_covariance_1, raw_vec,
-      random_effect_sigma_mu[1], random_effect_sigma_sigma[1],
+      random_effect_sigma_mu[1], 1.0,
       random_effect_sigma_raw_1, sigma_correlation_factor_1
     );
   }
@@ -165,7 +162,7 @@ transformed parameters {
     random_effect_2 = build_re_block(
       M, n_groups[2], how_many_factors_in_random_design[2], ncol_X_random_eff[2],
       group_factor_indexes_for_covariance_2, raw_vec,
-      random_effect_sigma_mu[2], random_effect_sigma_sigma[2],
+      random_effect_sigma_mu[2], 1.0,
       random_effect_sigma_raw_2, sigma_correlation_factor_2
     );
   }
@@ -175,7 +172,7 @@ transformed parameters {
     random_effect_3 = build_re_block(
       M, n_groups[3], how_many_factors_in_random_design[3], ncol_X_random_eff[3],
       group_factor_indexes_for_covariance_3, raw_vec,
-      random_effect_sigma_mu[3], random_effect_sigma_sigma[3],
+      random_effect_sigma_mu[3], 1.0,
       random_effect_sigma_raw_3, sigma_correlation_factor_3
     );
   }
@@ -185,7 +182,7 @@ transformed parameters {
     random_effect_4 = build_re_block(
       M, n_groups[4], how_many_factors_in_random_design[4], ncol_X_random_eff[4],
       group_factor_indexes_for_covariance_4, raw_vec,
-      random_effect_sigma_mu[4], random_effect_sigma_sigma[4],
+      random_effect_sigma_mu[4], 1.0,
       random_effect_sigma_raw_4, sigma_correlation_factor_4
     );
   }
diff --git a/inst/stan/glm_multi_beta_binomial_simulate_data.stan b/inst/stan/glm_multi_beta_binomial_simulate_data.stan
index 66f09ad8..c515bf1c 100755
--- a/inst/stan/glm_multi_beta_binomial_simulate_data.stan
+++ b/inst/stan/glm_multi_beta_binomial_simulate_data.stan
@@ -34,11 +34,7 @@ parameters{
   
   // sd of random intercept
   array[is_random_effect>0] real random_effect_sigma_mu;
-  // <lower=0> mirrors the fitting model in glm_multi_beta_binomial.stan: it
-  // breaks the (+sigma_sigma, +sigma_raw) / (-sigma_sigma, -sigma_raw) sign
-  // redundancy so prior-predictive draws match the constrained parameterisation
-  // the fit actually uses. It is not a minimum-variance assumption: sigma_sigma
-  // = 0 (homogeneous RE-SDs across categories) remains the prior mode.
+  // Retained for backward compatibility; fitting model fixes sigma_sigma at 1.0.
   array[is_random_effect>0] real<lower=0> random_effect_sigma_sigma;
 
 	// Covariance

From f11ca53de0711e32ea6ef1e50134231a21ab416e Mon Sep 17 00:00:00 2001
From: Stefano Mangiola <mangiolastefano@gmail.com>
Date: Mon, 1 Jun 2026 17:32:39 +1000
Subject: [PATCH 3/4] fix the prior drop

---
 inst/stan/glm_multi_beta_binomial.stan               | 5 +++--
 inst/stan/glm_multi_beta_binomial_generate_data.stan | 4 ++--
 inst/stan/glm_multi_beta_binomial_simulate_data.stan | 3 ++-
 3 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/inst/stan/glm_multi_beta_binomial.stan b/inst/stan/glm_multi_beta_binomial.stan
index 2562d763..bdb73a12 100755
--- a/inst/stan/glm_multi_beta_binomial.stan
+++ b/inst/stan/glm_multi_beta_binomial.stan
@@ -433,8 +433,8 @@ parameters{
   array[4 * (is_random_effect>0)] real random_effect_sigma_mu;
   // Retained for backward compatibility with saved fits / R interfaces; not
   // used in the likelihood (build_re_block is called with sigma_sigma = 1.0).
-  // Pinned so draws still expose the name at 1.0 without a prior.
-  array[4 * (is_random_effect>0)] real<lower=1, upper=1> random_effect_sigma_sigma;
+  // Still sampled via half-normal prior below so the name stays in draws.
+  array[4 * (is_random_effect>0)] real<lower=0> random_effect_sigma_sigma;
 
   // For models with a single group (kept from the original design)
   array[is_random_effect>0] real zero_random_effect;
@@ -620,6 +620,7 @@ model{
   // ----------------------------------------------------------------------
   if (is_random_effect > 0) {
     random_effect_sigma_mu ~ std_normal();
+    random_effect_sigma_sigma ~ std_normal();  // half-normal via <lower=0>
     zero_random_effect ~ std_normal();
   }
 
diff --git a/inst/stan/glm_multi_beta_binomial_generate_data.stan b/inst/stan/glm_multi_beta_binomial_generate_data.stan
index 50e4400a..89a215cd 100755
--- a/inst/stan/glm_multi_beta_binomial_generate_data.stan
+++ b/inst/stan/glm_multi_beta_binomial_generate_data.stan
@@ -118,8 +118,8 @@ parameters {
 
   array[4 * (is_random_effect>0)] real random_effect_sigma_mu;
   // Retained for backward compatibility when binding saved draws; build_re_block
-  // uses sigma_sigma = 1.0 (same as the fitting model). Unconstrained so old fits
-  // with arbitrary stored values still load.
+  // uses sigma_sigma = 1.0 (same as the fitting model). <lower=0> matches the
+  // fitting model; values from saved draws are bound but not used in the build.
   array[4 * (is_random_effect>0)] real<lower=0> random_effect_sigma_sigma;
   array[is_random_effect>0] real zero_random_effect;
 }
diff --git a/inst/stan/glm_multi_beta_binomial_simulate_data.stan b/inst/stan/glm_multi_beta_binomial_simulate_data.stan
index c515bf1c..c98bc401 100755
--- a/inst/stan/glm_multi_beta_binomial_simulate_data.stan
+++ b/inst/stan/glm_multi_beta_binomial_simulate_data.stan
@@ -34,7 +34,8 @@ parameters{
   
   // sd of random intercept
   array[is_random_effect>0] real random_effect_sigma_mu;
-  // Retained for backward compatibility; fitting model fixes sigma_sigma at 1.0.
+  // Retained for backward compatibility; fitting model fixes sigma_sigma at 1.0
+  // in build_re_block but keeps this parameter with <lower=0> for draw structure.
   array[is_random_effect>0] real<lower=0> random_effect_sigma_sigma;
 
 	// Covariance

From dc59ffb3d04f1df1faaf56858836d1196ea3d90a Mon Sep 17 00:00:00 2001
From: Stefano Mangiola <mangiolastefano@gmail.com>
Date: Thu, 4 Jun 2026 17:35:21 +0930
Subject: [PATCH 4/4] Update Stan models to clarify hyperprior handling

This commit modifies the handling of hyperprior parameters in the Stan models, specifically setting sigma_mu to 0.0 and sigma_sigma to 1.0 for backward compatibility. Comments have been updated for clarity, ensuring that the changes are well-documented and consistent across the glm_multi_beta_binomial models.
---
 inst/stan/common_functions.stan               |  2 +-
 inst/stan/glm_multi_beta_binomial.stan        | 20 +++++++++----------
 ...glm_multi_beta_binomial_generate_data.stan | 11 +++++-----
 ...glm_multi_beta_binomial_simulate_data.stan |  4 ++--
 4 files changed, 18 insertions(+), 19 deletions(-)

diff --git a/inst/stan/common_functions.stan b/inst/stan/common_functions.stan
index 2e28d89f..5e8ac484 100644
--- a/inst/stan/common_functions.stan
+++ b/inst/stan/common_functions.stan
@@ -109,7 +109,7 @@
       int ncol,                                           // # cols in this block's design matrix
       array[,] int group_factor_indexes,                  // factor x group -> col index
       array[] vector raw_vec,                             // length ncol, each length M
-      real sigma_mu,                                      // hierarchical hyperprior mean
+      real sigma_mu,                                      // hierarchical hyperprior mean (callers pass 0.0; see fitting model)
       real sigma_sigma,                                   // hierarchical hyperprior SD (callers pass 1.0; see fitting model)
       array[] vector sigma_raw,                           // length M, each length n_factors
       array[] matrix corr_chol                            // length M, each n_factors x n_factors
diff --git a/inst/stan/glm_multi_beta_binomial.stan b/inst/stan/glm_multi_beta_binomial.stan
index bdb73a12..b7243074 100755
--- a/inst/stan/glm_multi_beta_binomial.stan
+++ b/inst/stan/glm_multi_beta_binomial.stan
@@ -405,8 +405,8 @@ parameters{
   //   * sigma_correlation_factor_k: per-category Cholesky of correlation matrix
   //                                  (n_factors[k] x n_factors[k]; 1x1 = no LKJ work)
   //
-  // Hyperprior scalar sigma_mu is shared (length-4 array). Between-category
-  // log-SD spread is fixed at sigma_sigma = 1.0 in build_re_block (see below).
+  // Hyperprior scalars retained for backward compatibility. build_re_block uses
+  // fixed sigma_mu = 0.0 and sigma_sigma = 1.0 (see transformed parameters).
   // ----------------------------------------------------------------------
 
   // Slot 1
@@ -429,11 +429,11 @@ parameters{
   array[M * (ncol_X_random_eff[4]>0)] vector[how_many_factors_in_random_design[4]]  random_effect_sigma_raw_4;
   array[M * (ncol_X_random_eff[4]>0)] cholesky_factor_corr[how_many_factors_in_random_design[4] * (ncol_X_random_eff[4]>0)] sigma_correlation_factor_4;
 
-  // Shared hyperprior scalars (one mu per slot)
+  // Shared hyperprior scalars (one mu, one sigma per slot)
+  // Retained for backward compatibility; not used in the likelihood
+  // (build_re_block is called with sigma_mu = 0.0, sigma_sigma = 1.0).
+  // Priors below keep both names in draws.
   array[4 * (is_random_effect>0)] real random_effect_sigma_mu;
-  // Retained for backward compatibility with saved fits / R interfaces; not
-  // used in the likelihood (build_re_block is called with sigma_sigma = 1.0).
-  // Still sampled via half-normal prior below so the name stays in draws.
   array[4 * (is_random_effect>0)] real<lower=0> random_effect_sigma_sigma;
 
   // For models with a single group (kept from the original design)
@@ -483,7 +483,7 @@ transformed parameters{
     random_effect_1 = build_re_block(
       M, n_groups[1], how_many_factors_in_random_design[1], ncol_X_random_eff[1],
       group_factor_indexes_for_covariance_1, raw_vec,
-      random_effect_sigma_mu[1], 1.0,
+      0.0, 1.0,
       random_effect_sigma_raw_1, sigma_correlation_factor_1
     );
   }
@@ -494,7 +494,7 @@ transformed parameters{
     random_effect_2 = build_re_block(
       M, n_groups[2], how_many_factors_in_random_design[2], ncol_X_random_eff[2],
       group_factor_indexes_for_covariance_2, raw_vec,
-      random_effect_sigma_mu[2], 1.0,
+      0.0, 1.0,
       random_effect_sigma_raw_2, sigma_correlation_factor_2
     );
   }
@@ -505,7 +505,7 @@ transformed parameters{
     random_effect_3 = build_re_block(
       M, n_groups[3], how_many_factors_in_random_design[3], ncol_X_random_eff[3],
       group_factor_indexes_for_covariance_3, raw_vec,
-      random_effect_sigma_mu[3], 1.0,
+      0.0, 1.0,
       random_effect_sigma_raw_3, sigma_correlation_factor_3
     );
   }
@@ -516,7 +516,7 @@ transformed parameters{
     random_effect_4 = build_re_block(
       M, n_groups[4], how_many_factors_in_random_design[4], ncol_X_random_eff[4],
       group_factor_indexes_for_covariance_4, raw_vec,
-      random_effect_sigma_mu[4], 1.0,
+      0.0, 1.0,
       random_effect_sigma_raw_4, sigma_correlation_factor_4
     );
   }
diff --git a/inst/stan/glm_multi_beta_binomial_generate_data.stan b/inst/stan/glm_multi_beta_binomial_generate_data.stan
index 89a215cd..061fc8cd 100755
--- a/inst/stan/glm_multi_beta_binomial_generate_data.stan
+++ b/inst/stan/glm_multi_beta_binomial_generate_data.stan
@@ -118,8 +118,7 @@ parameters {
 
   array[4 * (is_random_effect>0)] real random_effect_sigma_mu;
   // Retained for backward compatibility when binding saved draws; build_re_block
-  // uses sigma_sigma = 1.0 (same as the fitting model). <lower=0> matches the
-  // fitting model; values from saved draws are bound but not used in the build.
+  // uses sigma_mu = 0.0 and sigma_sigma = 1.0 (same as the fitting model).
   array[4 * (is_random_effect>0)] real<lower=0> random_effect_sigma_sigma;
   array[is_random_effect>0] real zero_random_effect;
 }
@@ -152,7 +151,7 @@ transformed parameters {
     random_effect_1 = build_re_block(
       M, n_groups[1], how_many_factors_in_random_design[1], ncol_X_random_eff[1],
       group_factor_indexes_for_covariance_1, raw_vec,
-      random_effect_sigma_mu[1], 1.0,
+      0.0, 1.0,
       random_effect_sigma_raw_1, sigma_correlation_factor_1
     );
   }
@@ -162,7 +161,7 @@ transformed parameters {
     random_effect_2 = build_re_block(
       M, n_groups[2], how_many_factors_in_random_design[2], ncol_X_random_eff[2],
       group_factor_indexes_for_covariance_2, raw_vec,
-      random_effect_sigma_mu[2], 1.0,
+      0.0, 1.0,
       random_effect_sigma_raw_2, sigma_correlation_factor_2
     );
   }
@@ -172,7 +171,7 @@ transformed parameters {
     random_effect_3 = build_re_block(
       M, n_groups[3], how_many_factors_in_random_design[3], ncol_X_random_eff[3],
       group_factor_indexes_for_covariance_3, raw_vec,
-      random_effect_sigma_mu[3], 1.0,
+      0.0, 1.0,
       random_effect_sigma_raw_3, sigma_correlation_factor_3
     );
   }
@@ -182,7 +181,7 @@ transformed parameters {
     random_effect_4 = build_re_block(
       M, n_groups[4], how_many_factors_in_random_design[4], ncol_X_random_eff[4],
       group_factor_indexes_for_covariance_4, raw_vec,
-      random_effect_sigma_mu[4], 1.0,
+      0.0, 1.0,
       random_effect_sigma_raw_4, sigma_correlation_factor_4
     );
   }
diff --git a/inst/stan/glm_multi_beta_binomial_simulate_data.stan b/inst/stan/glm_multi_beta_binomial_simulate_data.stan
index c98bc401..27ddd8da 100755
--- a/inst/stan/glm_multi_beta_binomial_simulate_data.stan
+++ b/inst/stan/glm_multi_beta_binomial_simulate_data.stan
@@ -34,8 +34,8 @@ parameters{
   
   // sd of random intercept
   array[is_random_effect>0] real random_effect_sigma_mu;
-  // Retained for backward compatibility; fitting model fixes sigma_sigma at 1.0
-  // in build_re_block but keeps this parameter with <lower=0> for draw structure.
+  // Retained for backward compatibility; fitting model fixes sigma_mu = 0.0 and
+  // sigma_sigma = 1.0 in build_re_block but keeps these parameters for draw structure.
   array[is_random_effect>0] real<lower=0> random_effect_sigma_sigma;
 
 	// Covariance