diff --git a/DESCRIPTION b/DESCRIPTION
index baa32511..f87a088c 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.27
+Version: 2.1.28
 Date: 2024-01-15
 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/NAMESPACE b/NAMESPACE
index 206c10fe..6193be61 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,8 +12,8 @@ S3method(sccomp_proportional_fold_change,sccomp_tbl)
 S3method(sccomp_remove_outliers,sccomp_tbl)
 S3method(sccomp_remove_unwanted_effects,sccomp_tbl)
 S3method(sccomp_replicate,sccomp_tbl)
+S3method(sccomp_simulate,sccomp_tbl)
 S3method(sccomp_test,sccomp_tbl)
-S3method(simulate_data,tbl)
 export(clear_draw_files)
 export(clear_stan_model_cache)
 export(plot_1D_intervals)
@@ -27,6 +27,7 @@ export(sccomp_remove_outliers)
 export(sccomp_remove_unwanted_effects)
 export(sccomp_remove_unwanted_variation)
 export(sccomp_replicate)
+export(sccomp_simulate)
 export(sccomp_stan_models_cache_dir)
 export(sccomp_test)
 export(sccomp_theme)
@@ -119,6 +120,7 @@ importFrom(rlang,quo_is_symbolic)
 importFrom(rlang,quo_name)
 importFrom(rlang,quo_squash)
 importFrom(rlang,set_names)
+importFrom(rlang,sym)
 importFrom(scales,trans_new)
 importFrom(stats,C)
 importFrom(stats,as.formula)
diff --git a/R/simulate_data.R b/R/simulate_data.R
index c30c54e0..8128084a 100644
--- a/R/simulate_data.R
+++ b/R/simulate_data.R
@@ -1,4 +1,4 @@
-#' simulate_data
+#' sccomp_simulate
 #'
 #' @description This function simulates data from a fitted model.
 #'
@@ -7,17 +7,21 @@
 #' @importFrom magrittr divide_by
 #' @importFrom magrittr multiply_by
 #' @importFrom magrittr equals
-#' @importFrom rlang quo_is_null
+#' @importFrom rlang quo_is_null sym
 #' @importFrom SingleCellExperiment colData
 #' @importFrom parallel detectCores
+#' @importFrom lifecycle deprecate_warn
+#' @importFrom tidyr crossing expand_grid cross_join
 #'
-#' @param .data A tibble including a cell_group name column | sample name column | read counts column | factor columns | Pvalue column | a significance column
-#' @param .estimate_object The result of sccomp_estimate execution. This is used for sampling from real-data properties.
+#' @param fit The result of sccomp_estimate execution. This is used for sampling from real-data properties.
 #' @param formula_composition A formula. The formula describing the model for differential abundance, for example ~treatment
+#' @param new_data A tibble including sample-specific columns (sample identifier and factor columns from formula). If coefficients are provided separately, cell_group column is optional. Otherwise, should include cell_group column.
 #' @param formula_variability A formula. The formula describing the model for differential variability, for example ~treatment
+#' @param coefficients A data frame/tibble with cell-type specific coefficients. Must contain a column matching the cell_group column name, and columns matching the design matrix column names (e.g., "(Intercept)", "typeB"). If NULL, posterior beta_raw will be used.
+#' @param mean_dispersion_slope Optional numeric value for the slope parameter of the mean-dispersion association. If NULL, the slope from the fitted model (prec_coeff[2]) will be used.
+#' @param mean_dispersion_intercept Optional numeric value for the intercept parameter of the mean-dispersion association. If NULL, the intercept from the fitted model (prec_coeff[1]) will be used. This ensures that when only the slope is changed, the intercept remains the same.
 #' @param .sample A column name as symbol. The sample identifier
 #' @param .cell_group A column name as symbol. The cell_group identifier
-#' @param .coefficients The column names for coefficients, for example, c(b_0, b_1)
 #' @param variability_multiplier A real scalar. This can be used for artificially increasing the variability of the simulation for benchmarking purposes.
 #' @param number_of_draws An integer. How may copies of the data you want to draw from the model joint posterior distribution.
 #' @param mcmc_seed An integer. Used for Markov-chain Monte Carlo reproducibility. By default a random number is sampled from 1 to 999999. This itself can be controlled by set.seed()
@@ -26,6 +30,7 @@
 #' @param cache_stan_model A character string specifying the cache directory for compiled Stan models. 
 #'                        The sccomp version will be automatically appended to ensure version isolation.
 #'                        Default is `sccomp_stan_models_cache_dir` which points to `~/.sccomp_models`.
+#' @param .coefficients (Deprecated) The column names for coefficients in new_data, for example, c(b_0, b_1). Use the 'coefficients' parameter instead. This parameter is placed last for backward compatibility.
 #' 
 #' @return A tibble (`tbl`) with the following columns:
 #' \itemize{
@@ -67,27 +72,30 @@
 #' #     counts_obj = counts_obj |> mutate(b_0 = 0, b_1 = 0)
 #'
 #' #     # Simulate data
-#' #     simulate_data(counts_obj, estimate, ~type, ~1, sample, cell_group, c(b_0, b_1))
+#' #     sccomp_simulate(estimate, ~type, ~1, counts_obj, sample, cell_group, c(b_0, b_1))
 #' #   }
 #' # }
-simulate_data <- function(.data,
-                          .estimate_object,
+sccomp_simulate <- function(fit,
                           formula_composition,
                           formula_variability = NULL,
+                          new_data = NULL,
+                          coefficients = NULL,
+                          mean_dispersion_slope = NULL,
+                          mean_dispersion_intercept = NULL,
                           .sample = NULL,
                           .cell_group = NULL,
-                          .coefficients = NULL,
-                          variability_multiplier = 5,
+                          variability_multiplier = 1,
                           number_of_draws = 1,
                           mcmc_seed = sample_seed(),
                           cores = detectCores(),
                           sig_figs = 9,
-                          cache_stan_model = sccomp_stan_models_cache_dir) {
+                          cache_stan_model = sccomp_stan_models_cache_dir,
+                          .coefficients = NULL) {
   
   # Run the function
   check_and_install_cmdstanr()
   
-  UseMethod("simulate_data", .data)
+  UseMethod("sccomp_simulate", fit)
 }
 
 #' @export
@@ -99,87 +107,568 @@ simulate_data <- function(.data,
 #' @importFrom readr read_file
 #' @importFrom tibble column_to_rownames
 #'
-simulate_data.tbl = function(.data,
-                             .estimate_object,
+sccomp_simulate.sccomp_tbl = function(fit,
                              formula_composition,
                              formula_variability = NULL,
+                                     new_data = NULL,
+                                     coefficients = NULL,
+                                     mean_dispersion_slope = NULL,
+                                     mean_dispersion_intercept = NULL,
                              .sample = NULL,
                              .cell_group = NULL,
-                             .coefficients = NULL,
-                             variability_multiplier = 5,
+                             variability_multiplier = 1,
                              number_of_draws = 1,
                              mcmc_seed = sample_seed(),
                              cores = detectCores(),
                              sig_figs = 9,
-                             cache_stan_model = sccomp_stan_models_cache_dir) {
+                                     cache_stan_model = sccomp_stan_models_cache_dir,
+                                     .coefficients = NULL) {
+  
   
+  model_data = attr(fit, "model_input")
+  original_data = attr(fit, "model_input")
   
+  # Get sample and cell_group from fit attributes if not provided
+  if(is.null(.sample)) {
+    .sample_attr = attr(fit, ".sample")
+    if(!is.null(.sample_attr)) {
+      .sample = .sample_attr
+    } else {
+      stop("sccomp says: .sample must be provided or available in fit attributes")
+    }
+  } else {
   .sample = enquo(.sample)
+  }
+  
+  if(is.null(.cell_group)) {
+    .cell_group_attr = attr(fit, ".cell_group")
+    if(!is.null(.cell_group_attr)) {
+      .cell_group = .cell_group_attr
+    } else {
+      stop("sccomp says: .cell_group must be provided or available in fit attributes")
+    }
+  } else {
   .cell_group = enquo(.cell_group)
-  .coefficients = enquo(.coefficients)
+  }
+  
+  # Handle coefficients parameter
+  # If coefficients table is provided, use it; otherwise check for .coefficients column names (backward compatibility)
+  .coefficients_quo = enquo(.coefficients)
+  
+  # Deprecate .coefficients parameter in favor of coefficients table
+  if(!rlang::quo_is_null(.coefficients_quo)) {
+    lifecycle::deprecate_warn(
+      "2.2.0",
+      "sccomp_simulate(.coefficients)",
+      details = "sccomp says: .coefficients parameter is deprecated. Please use the 'coefficients' parameter with a data frame/tibble instead. See ?sccomp_simulate for details."
+    )
+  }
   
-  #Check column class
+  # Use new_data if provided, otherwise use count_data from fit
+  if(is.null(new_data)) {
+    .data = attr(fit, "count_data")
+  } else {
+    .data = new_data
+  }
+  
+  # Check column class
+  # If coefficients table is provided separately, new_data might not have cell_group column
+  # In that case, skip cell_group validation for new_data
+  if(is.null(coefficients) || quo_name(.cell_group) %in% colnames(.data)) {
   check_if_columns_right_class(.data, !!.sample, !!.cell_group)
+  }
+  # If coefficients are provided separately, we'll get cell_group info from coefficients table
+  
+  # Get formulas if not provided
+  if(is.null(formula_variability)) {
+    formula_variability = attr(fit, "formula_variability")
+  }
+  original_formula_composition = attr(fit, "formula_composition")
   
-  model_data = attr(.estimate_object, "model_input")
+  # Validate dimensions before proceeding
+  # Use prepare_replicate_data machinery to get X_which for proper beta subsetting
+  # This avoids padding beta with zeros and matches the approach used in sccomp_predict
+  original_X = original_data$X
+  original_Xa = original_data$Xa
   
-  # # Select model based on noise model
-  # if(attr(.estimate_object, "noise_model") == "multi_beta_binomial") my_model = stanmodels$glm_multi_beta_binomial_simulate_data
-  # else if(attr(.estimate_object, "noise_model") == "dirichlet_multinomial") my_model = get_model_from_data("model_glm_dirichlet_multinomial_generate_quantities.rds", glm_dirichlet_multinomial_generate_quantities)
-  # else if(attr(.estimate_object, "noise_model") == "logit_normal_multinomial") my_model = get_model_from_data("glm_multinomial_logit_linear_simulate_data.stan", read_file("~/PostDoc/sccomp/dev/stan_models/glm_multinomial_logit_linear_simulate_data.stan"))
+  # Prepare data using the same machinery as replicate_data
+  prepared_data = prepare_replicate_data(
+    X = original_X,
+    Xa = original_Xa,
+    N = nrow(.data |> distinct(!!.sample)),
+    intercept_in_design = original_data$intercept_in_design,
+    X_random_effect = original_data$X_random_effect,
+    X_random_effect_2 = original_data$X_random_effect_2,
+    .sample = !!.sample,
+    .cell_group = !!.cell_group,
+    .count = quo(count),  # Dummy - not used for simulation
+    formula_composition = formula_composition,
+    original_formula_composition = original_formula_composition,
+    formula_variability = formula_variability,
+    new_data = .data |>
+      nest(data___ = -!!.sample) |>
+      mutate(.exposure = sample(model_data$exposure, size = n(), replace = TRUE )) |>
+      unnest(data___) |>
+      select(!!.sample, any_of(c(parse_formula(formula_composition), parse_formula(formula_variability)))) |>
+      distinct(),
+    original_count_data = attr(fit, "count_data") |>
+      distinct(!!.sample) |>
+      mutate(dummy = 1)
+  )
   
+  # Get X_which for beta subsetting
+  X_which = prepared_data$X_which
+  XA_which = prepared_data$XA_which
+  
+  # Create data_for_model with simulated dimensions
+  # Use prepared_data$X and prepared_data$Xa (which match X_which) instead of creating new ones
+  # If coefficients are provided separately, expand .data to include all cell_groups
+  if(!is.null(coefficients) && !quo_name(.cell_group) %in% colnames(.data)) {
+    cell_group_colname = quo_name(.cell_group)
+    all_cell_groups = coefficients[[cell_group_colname]]
+    # Create a tibble with all cell_groups and cross join with .data
+    cell_group_tibble = tibble(!!.cell_group := all_cell_groups)
+    .data_expanded = 
+      .data |>
+      cross_join(cell_group_tibble)
+  } else {
+    .data_expanded = .data
+  }
+  
+  # If coefficients table is provided, don't pass coefficients to data_to_simulation_covariates
+  # (we'll extract them separately later)
+  # Use an empty quosure if coefficients table is provided
+  if(!is.null(coefficients)) {
+    coefficients_quo_for_data_prep = quo(NULL)
+  } else {
+    coefficients_quo_for_data_prep = .coefficients_quo
+  }
   
   data_for_model =
-    .data |>
+    .data_expanded |>
     nest(data___ = -!!.sample) |>
     mutate(.exposure = sample(model_data$exposure, size = n(), replace = TRUE )) |>
     unnest(data___) |>
-    data_simulation_to_model_input(
+    data_to_simulation_covariates(
       formula_composition,
-      #formula_variability,
-      !!.sample, !!.cell_group, .exposure, !!.coefficients
+      formula_variability,
+      !!.sample, !!.cell_group, .exposure, !!coefficients_quo_for_data_prep
     )
   names(data_for_model)  = names(data_for_model) |> stringr::str_c("_simulated")
   
-  # Drop data from old input
-  original_data = .estimate_object |> attr("model_input")
-  original_data = original_data[(names(original_data) %in% c("C", "M", "A", "ncol_X_random_eff", "is_random_effect", "how_many_factors_in_random_design"  ))]
+  # Update X_simulated and Xa_simulated to use the prepared matrices (which match X_which)
+  data_for_model$X_simulated = prepared_data$X
+  data_for_model$Xa_simulated = prepared_data$Xa
+  data_for_model$C_simulated = ncol(prepared_data$X)
+  data_for_model$A_simulated = ncol(prepared_data$Xa)
+  
+  # Extract coefficients based on prepared_data$X columns (which match X_which)
+  # Coefficients can be provided as:
+  # 1. A separate table (coefficients parameter) - cell-type specific
+  # 2. Column names in new_data (.coefficients parameter) - backward compatibility
+  # 3. NULL - use posterior beta_raw
+  
+  # Get design matrix column names from prepared_data$X (which matches X_which)
+  X_simulated_colnames = colnames(prepared_data$X)
+  
+  # Get cell group names
+  # If coefficients table is provided, get cell_group names from there
+  # Otherwise, get from .data_expanded (which includes cell_group)
+  if(!is.null(coefficients)) {
+    cell_group_colname = quo_name(.cell_group)
+    
+    # Validate coefficients table has cell_group column
+    if(!cell_group_colname %in% colnames(coefficients)) {
+      stop("sccomp says: coefficients table must contain a column matching the cell_group column name (", cell_group_colname, ")")
+    }
+    
+    cell_group_names = 
+      coefficients |>
+      pull(!!sym(cell_group_colname)) |>
+      unique() |>
+      sort()
+  } else {
+    cell_group_names = 
+      .data_expanded |>
+      distinct(!!.cell_group) |>
+      arrange(!!.cell_group) |>
+      pull(!!.cell_group)
+  }
   
-  # [1]  5.6260004 -0.6940178
-  # prec_sd  = 0.816423129
+  # Initialize beta_simulated_aligned
+  beta_simulated_aligned = NULL
+  
+  # Case 1: coefficients table is provided
+  if(!is.null(coefficients)) {
+    # Create beta matrix: rows = cell_types, cols = design columns (matching X_simulated_colnames)
+    beta_simulated_aligned = matrix(0, nrow = length(cell_group_names), ncol = length(X_simulated_colnames))
+    rownames(beta_simulated_aligned) = cell_group_names
+    colnames(beta_simulated_aligned) = X_simulated_colnames
+    
+    # Extract coefficients from coefficients table
+    # Use column name directly instead of quosure
+    coefficients_data = 
+      coefficients |>
+      select(all_of(cell_group_colname), any_of(X_simulated_colnames)) |>
+      arrange(!!sym(cell_group_colname))
+    
+    # Map coefficients to design matrix columns
+    # Coefficient column names should match design matrix column names
+    for(col_name in X_simulated_colnames) {
+      if(col_name %in% colnames(coefficients_data)) {
+        # Match by cell_group and extract coefficient values
+        for(i in 1:length(cell_group_names)) {
+          cell_name = cell_group_names[i]
+          matching_row = coefficients_data[[cell_group_colname]] == cell_name
+          if(any(matching_row)) {
+            beta_simulated_aligned[i, col_name] = coefficients_data[matching_row, col_name][[1]]
+          }
+        }
+      }
+    }
+    
+    # Note: Normalization to sum-to-zero happens after transpose (see below)
+    # beta_simulated_aligned has [cell_types, design_columns], we need to normalize columns
+    # but we'll normalize rows after transpose to [design_columns, cell_types]
+    data_for_model$beta_simulated = beta_simulated_aligned
+    
+  } else if(!rlang::quo_is_null(.coefficients_quo)) {
+    # Case 2: .coefficients column names provided in new_data (backward compatibility)
+    # Get coefficient column names
+    coeff_names = quo_names(.coefficients_quo)
+    
+    # Create beta matrix: rows = cell_types, cols = design columns (matching X_simulated_colnames)
+    beta_simulated_aligned = matrix(0, nrow = length(cell_group_names), ncol = length(X_simulated_colnames))
+    rownames(beta_simulated_aligned) = cell_group_names
+    colnames(beta_simulated_aligned) = X_simulated_colnames
+    
+    # Extract coefficients from data
+    coefficients_data = 
+      .data |>
+      select(!!.cell_group, all_of(coeff_names)) |>
+      distinct() |>
+      arrange(!!.cell_group)
+    
+    # Map coefficients to design matrix columns
+    # If coefficient names match design column names, use them directly
+    # Otherwise, map by position (assuming same order)
+    for(i in 1:length(X_simulated_colnames)) {
+      col_name = X_simulated_colnames[i]
+      # Try to find matching coefficient column
+      if(col_name %in% coeff_names) {
+        beta_simulated_aligned[, col_name] = coefficients_data[[col_name]]
+      } else if(i <= length(coeff_names)) {
+        # Map by position if no name match
+        beta_simulated_aligned[, i] = coefficients_data[[coeff_names[i]]]
+      }
+    }
+    
+    data_for_model$beta_simulated = beta_simulated_aligned
+  }
+  # Case 3: coefficients is NULL and .coefficients is NULL - use posterior beta_raw (no action needed)
+  
+  # Issue 3: Validate dimensions
+  if(data_for_model$C_simulated > original_data$C) {
+    stop("sccomp says: C_simulated (", data_for_model$C_simulated, ") cannot be larger than C (", original_data$C, ") from the fitted model. The simulated design matrix has more columns than the original model.")
+  }
+  if(data_for_model$M_simulated > original_data$M) {
+    stop("sccomp says: M_simulated (", data_for_model$M_simulated, ") cannot be larger than M (", original_data$M, ") from the fitted model. The simulated data has more cell groups than the original model.")
+  }
+  if(data_for_model$A_simulated > original_data$A) {
+    stop("sccomp says: A_simulated (", data_for_model$A_simulated, ") cannot be larger than A (", original_data$A, ") from the fitted model. The simulated variability design has more columns than the original model.")
+  }
+  
+  # Extract necessary data from original model - reuse pattern from replicate_data
+  # These fields are always present in model_input from sccomp_estimate (set in data_spread_to_model_input)
+  original_data_subset = original_data[(names(original_data) %in% c("C", "M", "A", "A_intercept_columns", "intercept_in_design", "bimodal_mean_variability_association", "ncol_X_random_eff", "is_random_effect", "how_many_factors_in_random_design", "n_groups", "group_factor_indexes_for_covariance", "group_factor_indexes_for_covariance_2"))]
+  
+  # Create proper random effect design matrices
+  # Stan requires matrices to have at least 1 column, so we use max(1, ncol) to ensure valid dimensions
+  ncol_re1 = max(1, original_data_subset$ncol_X_random_eff[1])
+  ncol_re2 = max(1, original_data_subset$ncol_X_random_eff[2])
+  
+  # Initialize with empty matrices
+  X_random_effect_simulated = matrix(0, nrow = nrow(data_for_model$X_simulated), ncol = ncol_re1)
+  X_random_effect_2_simulated = matrix(0, nrow = nrow(data_for_model$X_simulated), ncol = ncol_re2)
+  
+  # Create proper random effect design matrices if they exist in the original model
+  if(original_data_subset$is_random_effect > 0 && 
+     !is.null(original_data$X_random_effect) && 
+     original_data_subset$ncol_X_random_eff[1] > 0) {
+    
+    # Get random effect elements from formula
+    random_effect_elements = parse_formula_random_effect(formula_composition)
+    original_grouping_names = original_formula_composition |> formula_to_random_effect_formulae() |> pull(grouping)
+    
+    if(length(original_grouping_names) > 0 && 
+       (random_effect_elements$grouping %in% original_grouping_names[1]) |> any()) {
+      
+      # Create random effect design for simulated data
+      random_effect_grouping = 
+        formula_composition |>
+        formula_to_random_effect_formulae() |>
+        mutate(design = map2(
+          formula, grouping,
+          ~  get_random_effect_design3(.data, .x, .y, !!.sample,
+                                       accept_NA_as_average_effect = TRUE )
+        ))
+      
+      if((random_effect_grouping$grouping %in% original_grouping_names[1]) |> any()) {
+        X_random_effect_simulated =
+          random_effect_grouping |>
+          filter(grouping==original_grouping_names[1]) |> 
+          mutate(design_matrix = map(
+            design,
+            ~ ..1 |>
+              select(!!.sample, group___label, value) |>
+              filter(group___label %in% colnames(original_data$X_random_effect)) |> 
+              pivot_wider(names_from = group___label, values_from = value) |>
+              mutate(across(everything(), ~ .x |> replace_na(0)))
+          )) |>
+          pull(design_matrix) |> 
+          _[[1]] |> 
+          column_to_rownames(quo_name(.sample))
+        
+        # Ensure matrix has correct dimensions (pad with zeros if needed)
+        if(ncol(X_random_effect_simulated) < ncol_re1) {
+          padding = matrix(0, nrow = nrow(X_random_effect_simulated), 
+                          ncol = ncol_re1 - ncol(X_random_effect_simulated))
+          X_random_effect_simulated = cbind(X_random_effect_simulated, padding)
+        } else if(ncol(X_random_effect_simulated) > ncol_re1) {
+          X_random_effect_simulated = X_random_effect_simulated[, 1:ncol_re1, drop = FALSE]
+        }
+      }
+    }
+  }
+  
+  # Create second random effect design matrix if it exists
+  if(original_data_subset$is_random_effect > 0 && 
+     !is.null(original_data$X_random_effect_2) && 
+     original_data_subset$ncol_X_random_eff[2] > 0) {
+    
+    random_effect_elements = parse_formula_random_effect(formula_composition)
+    original_grouping_names = original_formula_composition |> formula_to_random_effect_formulae() |> pull(grouping)
+    
+    if(length(original_grouping_names) > 1 && 
+       (random_effect_elements$grouping %in% original_grouping_names[2]) |> any()) {
+      
+      random_effect_grouping = 
+        formula_composition |>
+        formula_to_random_effect_formulae() |>
+        mutate(design = map2(
+          formula, grouping,
+          ~  get_random_effect_design3(.data, .x, .y, !!.sample,
+                                       accept_NA_as_average_effect = TRUE )
+        ))
+      
+      if((random_effect_grouping$grouping %in% original_grouping_names[2]) |> any()) {
+        X_random_effect_2_simulated =
+          random_effect_grouping |>
+          filter(grouping==original_grouping_names[2]) |> 
+          mutate(design_matrix = map(
+            design,
+            ~ ..1 |>
+              select(!!.sample, group___label, value) |>
+              filter(group___label %in% colnames(original_data$X_random_effect_2)) |> 
+              pivot_wider(names_from = group___label, values_from = value) |>
+              mutate(across(everything(), ~ .x |> replace_na(0)))
+          )) |>
+          pull(design_matrix) |> 
+          _[[1]] |> 
+          column_to_rownames(quo_name(.sample))
+        
+        # Ensure matrix has correct dimensions
+        if(ncol(X_random_effect_2_simulated) < ncol_re2) {
+          padding = matrix(0, nrow = nrow(X_random_effect_2_simulated), 
+                          ncol = ncol_re2 - ncol(X_random_effect_2_simulated))
+          X_random_effect_2_simulated = cbind(X_random_effect_2_simulated, padding)
+        } else if(ncol(X_random_effect_2_simulated) > ncol_re2) {
+          X_random_effect_2_simulated = X_random_effect_2_simulated[, 1:ncol_re2, drop = FALSE]
+        }
+      }
+    }
+  }
   
   mod_rng = load_model("glm_multi_beta_binomial_simulate_data", threads = cores, cache_dir = cache_stan_model)
   
+  # Check if coefficients are provided and validate/normalize them
+  user_provided_beta = 0
+  # When coefficients are not provided, pass empty array with correct dimensions for Stan
+  # Stan expects array[C_simulated] vector[M_simulated], so we create a list of length C_simulated
+  # Each element is an empty vector of length M_simulated
+  beta_simulated_provided = lapply(1:data_for_model$C_simulated, function(i) {
+    rep(0.0, data_for_model$M_simulated)
+  })
+  
+  if(!is.null(data_for_model$beta_simulated) && nrow(data_for_model$beta_simulated) > 0 && ncol(data_for_model$beta_simulated) > 0) {
+    # Coefficients are provided - validate dimensions
+    # Note: beta_simulated has been aligned with prepared_data$X columns above
+    # beta_simulated has rows=cell_types, cols=coefficients (matching X_simulated_colnames)
+    # Stan expects rows=design_columns (C), cols=cell_types (M), so we need to transpose
+    beta_provided = t(data_for_model$beta_simulated)
+    
+    # Validate dimensions after transpose
+    # beta_provided should match length_X_which (number of design columns in X_which)
+    expected_rows = length(X_which)
+    if(nrow(beta_provided) != expected_rows) {
+      stop("sccomp says: beta_simulated must have ", expected_rows, " design columns (matching X_which), but has ", nrow(beta_provided), ". This may happen if the design matrix columns don't match between the formula and the original model.")
+    }
+    if(ncol(beta_provided) != data_for_model$M_simulated) {
+      stop("sccomp says: beta_simulated must have ", data_for_model$M_simulated, " cell groups (after transpose), but has ", ncol(beta_provided))
+    }
+    
+    # Validate that coefficients sum to zero for each design column
+    # User-provided coefficients must sum to exactly zero (within floating-point tolerance)
+    # This is required for compositional models where log-ratios must sum to zero
+    max_abs_val = max(abs(beta_provided))
+    tolerance = .Machine$double.eps * 100 * max(max_abs_val, 1)
+    
+    # Check row sums using tidyverse approach
+    row_sums = beta_provided |>
+      as_tibble(.name_repair = "minimal") |>
+      rowwise() |>
+      mutate(row_sum = sum(c_across(everything()))) |>
+      pull(row_sum)
+    
+    # Find columns that don't sum to zero
+    invalid_cols = which(abs(row_sums) > tolerance)
+    if(length(invalid_cols) > 0) {
+      col_name = X_simulated_colnames[invalid_cols[1]]
+      row_sum = row_sums[invalid_cols[1]]
+      stop(
+        "sccomp says: User-provided coefficients for design column '", col_name, 
+        "' do not sum to zero (sum = ", round(row_sum, 10), "). ",
+        "Coefficients in compositional models must sum to zero. ",
+        "Please normalize your coefficients (e.g., set the last element to -sum(others)) ",
+        "or subtract the mean from each coefficient vector."
+      )
+    }
+    
+    # Convert to list format for Stan (array[length_X_which] vector[M_simulated])
+    # NOTE: We use vector[M_simulated] instead of sum_to_zero_vector[M_simulated] in Stan
+    # to avoid floating-point precision issues. Small precision errors in sum-to-zero constraint
+    # are acceptable since Stan doesn't enforce strict constraints with vector types.
+    # However, user-provided coefficients must sum to zero within reasonable tolerance (checked above).
+    # Stan expects a list/array where each element is a vector of length M_simulated
+    # In R, we pass this as a list of vectors
+    # The length matches X_which, not C_simulated
+    beta_simulated_provided = lapply(1:nrow(beta_provided), function(c) {
+      as.numeric(beta_provided[c, ])
+    })
+    
+    user_provided_beta = 1
+  }
+  
+  # Handle optional mean-dispersion slope and intercept parameters
+  # Similar to user_provided_beta, allow user to override prec_coeff[2] (slope) and prec_coeff[1] (intercept)
+  user_provided_prec_coeff_slope = 0
+  prec_coeff_slope_provided = 0.0  # Default value (will be ignored if user_provided_prec_coeff_slope == 0)
+  user_provided_prec_coeff_intercept = 0
+  prec_coeff_intercept_provided = 0.0  # Default value (will be ignored if user_provided_prec_coeff_intercept == 0)
+  
+  # If slope is provided but intercept is not, automatically preserve the mean intercept from the fit
+  if(!is.null(mean_dispersion_slope) && is.null(mean_dispersion_intercept)) {
+    fit_obj = attr(fit, "fit")
+    prec_coeff_summary = fit_obj$summary("prec_coeff")
+    mean_dispersion_intercept = prec_coeff_summary$mean[1]
+  }
+  
+  if(!is.null(mean_dispersion_slope)) {
+    if(!is.numeric(mean_dispersion_slope) || length(mean_dispersion_slope) != 1) {
+      stop("sccomp says: mean_dispersion_slope must be a single numeric value.")
+    }
+    user_provided_prec_coeff_slope = 1
+    prec_coeff_slope_provided = as.numeric(mean_dispersion_slope)
+  }
+  
+  if(!is.null(mean_dispersion_intercept)) {
+    if(!is.numeric(mean_dispersion_intercept) || length(mean_dispersion_intercept) != 1) {
+      stop("sccomp says: mean_dispersion_intercept must be a single numeric value.")
+    }
+    user_provided_prec_coeff_intercept = 1
+    prec_coeff_intercept_provided = as.numeric(mean_dispersion_intercept)
+  }
+  
+  # Issue 2: Xa_simulated is already in data_for_model, so we don't need to add it again
+  # Combine all data for Stan model
+  stan_data = data_for_model |> 
+    c(original_data_subset) |> 
+    c(list(
+      variability_multiplier = variability_multiplier,
+      X_random_effect_simulated = X_random_effect_simulated,
+      X_random_effect_2_simulated = X_random_effect_2_simulated,
+      user_provided_beta = user_provided_beta,
+      beta_simulated_provided = beta_simulated_provided,
+      user_provided_prec_coeff_slope = user_provided_prec_coeff_slope,
+      prec_coeff_slope_provided = prec_coeff_slope_provided,
+      user_provided_prec_coeff_intercept = user_provided_prec_coeff_intercept,
+      prec_coeff_intercept_provided = prec_coeff_intercept_provided,
+      length_X_which = length(X_which),
+      length_XA_which = length(XA_which),
+      X_which = X_which,
+      XA_which = XA_which
+    ))
+  
+  # Get posterior draws - reuse pattern from replicate_data
+  number_of_draws_in_the_fit = attr(fit, "fit") |> get_output_samples()
+  number_of_draws = min(number_of_draws, number_of_draws_in_the_fit)
+  
+  draws_matrix = attr(fit, "fit")$draws(format = "matrix")
+  
+  if(number_of_draws > nrow(draws_matrix)) {
+    number_of_draws = nrow(draws_matrix)
+  }
+  
+  # Sample draws if needed (reuse pattern from replicate_data)
+  if(number_of_draws < nrow(draws_matrix)) {
+    draws_matrix = draws_matrix[sample(seq_len(nrow(draws_matrix)), size = number_of_draws),, drop = FALSE]
+  }
+  
+  # Note: We no longer normalize sum-to-zero parameters since we use vector[M] instead of sum_to_zero_vector[M]
+  # Small precision errors in sum-to-zero constraint are acceptable since Stan doesn't enforce strict constraints
+  
+  # Generate quantities - reuse pattern from replicate_data
   fit = mod_rng |> sample_safe(
     generate_quantities_fx,
-    attr(.estimate_object , "fit")$draws(format = "matrix"),
-    
-    # This is for the new data generation with selected factors to do adjustment
-    data = data_for_model |> c(original_data) |> c(list(variability_multiplier = variability_multiplier)),
+    draws_matrix,
+    data = stan_data,
     seed = mcmc_seed,
-    parallel_chains = attr(.estimate_object , "fit")$metadata()$threads_per_chain, 
+    parallel_chains = attr(fit, "fit")$metadata()$threads_per_chain, 
     threads_per_chain = cores,
     sig_figs = sig_figs
   )
   
+  # Parse generated quantities - reuse parse_generated_quantities from replicate_data
+  # Get cell group names from the simulated data (same as used in data_to_simulation_covariates)
+  # Use .data_expanded which includes cell_group (either from original data or expanded from coefficients)
+  cell_group_names = 
+    .data_expanded |>
+    distinct(!!.cell_group) |>
+    arrange(!!.cell_group) |>
+    pull(!!.cell_group)
+  
+  sample_names = rownames(data_for_model$X_simulated)
+  
   parsed_fit =
     fit |>
     parse_generated_quantities(number_of_draws = number_of_draws) |>
     
-    # Get sample name
+    # Get sample name - reuse pattern from sccomp_replicate
     nest(data = -N) |>
     arrange(N) |>
-    mutate(!!.sample := rownames(data_for_model$X_simulated)) |>
+    mutate(!!.sample := sample_names) |>
     unnest(data) |>
     
-    # get cell type name
+    # get cell type name - reuse pattern from sccomp_replicate
     nest(data = -M) |>
-    mutate(!!.cell_group := colnames(data_for_model$beta_simulated)) |>
+    mutate(!!.cell_group := cell_group_names) |>
     unnest(data) |>
     
     select(-N, -M)
   
-  .data |>
+  # Join with original data - reuse pattern from sccomp_predict
+  # Use .data_expanded which includes cell_group (either from original data or expanded from coefficients)
+  .data_expanded |>
     left_join(
       parsed_fit,
       by = c(quo_name(.sample), quo_name(.cell_group))
@@ -188,14 +677,75 @@ simulate_data.tbl = function(.data,
   
 }
 
+#' DEPRECATED: simulate_data
+#'
+#' @description This function is DEPRECATED. Please use \code{\link{sccomp_simulate}} instead.
+#'
+#' @param .data A tibble including a cell_group name column | sample name column | read counts column | factor columns | Pvalue column | a significance column
+#' @param .estimate_object The result of sccomp_estimate execution. This is used for sampling from real-data properties.
+#' @param formula_composition A formula. The formula describing the model for differential abundance, for example ~treatment
+#' @param formula_variability A formula. The formula describing the model for differential variability, for example ~treatment
+#' @param .sample A column name as symbol. The sample identifier
+#' @param .cell_group A column name as symbol. The cell_group identifier
+#' @param .coefficients (Deprecated) The column names for coefficients in new_data, for example, c(b_0, b_1). Use the 'coefficients' parameter in sccomp_simulate() instead.
+#' @param variability_multiplier A real scalar. This can be used for artificially increasing the variability of the simulation for benchmarking purposes.
+#' @param number_of_draws An integer. How may copies of the data you want to draw from the model joint posterior distribution.
+#' @param mcmc_seed An integer. Used for Markov-chain Monte Carlo reproducibility. By default a random number is sampled from 1 to 999999. This itself can be controlled by set.seed()
+#' @param cores Integer, the number of cores to be used for parallel calculations.
+#' @param sig_figs Number of significant figures to use for Stan model output. Default is 9.
+#' @param cache_stan_model A character string specifying the cache directory for compiled Stan models. 
+#'                        The sccomp version will be automatically appended to ensure version isolation.
+#'                        Default is `sccomp_stan_models_cache_dir` which points to `~/.sccomp_models`.
+#'
+#' @export
+#'
+simulate_data <- function(.data,
+                          .estimate_object,
+                          formula_composition,
+                          formula_variability = NULL,
+                          .sample = NULL,
+                          .cell_group = NULL,
+                          .coefficients = NULL,
+                          variability_multiplier = 1,
+                          number_of_draws = 1,
+                          mcmc_seed = sample_seed(),
+                          cores = detectCores(),
+                          sig_figs = 9,
+                          cache_stan_model = sccomp_stan_models_cache_dir) {
+  
+  lifecycle::deprecate_warn(
+    "2.1.26",
+    "sccomp::simulate_data()",
+    details = "sccomp says: simulate_data is deprecated. Please use sccomp_simulate() instead."
+  )
+  
+  sccomp_simulate(
+    fit = .estimate_object,
+    formula_composition = formula_composition,
+    formula_variability = formula_variability,
+    new_data = .data,
+    .sample = .sample,
+    .cell_group = .cell_group,
+    .coefficients = .coefficients,
+    mean_dispersion_slope = mean_dispersion_slope,
+    mean_dispersion_intercept = mean_dispersion_intercept,
+    variability_multiplier = variability_multiplier,
+    number_of_draws = number_of_draws,
+    mcmc_seed = mcmc_seed,
+    cores = cores,
+    sig_figs = sig_figs,
+    cache_stan_model = cache_stan_model
+  )
+}
+
 #' @importFrom purrr when
 #' @importFrom stats model.matrix
 #'
 #' @keywords internal
 #' @noRd
 #'
-data_simulation_to_model_input =
-  function(.data, formula, .sample, .cell_type, .exposure, .coefficients, truncation_ajustment = 1, approximate_posterior_inference ){
+data_to_simulation_covariates =
+  function(.data, formula, formula_variability = ~ 1, .sample, .cell_type, .exposure, .coefficients, truncation_ajustment = 1, approximate_posterior_inference ){
     
     # Define the variables as NULL to avoid CRAN NOTES
     sd <- NULL
@@ -209,12 +759,15 @@ data_simulation_to_model_input =
     .coefficients = enquo(.coefficients)
     
     factor_names = parse_formula(formula)
+    factor_names_variability = parse_formula(formula_variability)
     
     sample_data =
       .data %>%
-      select(!!.sample, any_of(factor_names)) %>%
+      select(!!.sample, any_of(c(factor_names, factor_names_variability))) %>%
       distinct() %>%
       arrange(!!.sample)
+    
+    # Create composition design matrix
     X =
       sample_data %>%
       model.matrix(formula, data=.) %>%
@@ -230,11 +783,25 @@ data_simulation_to_model_input =
         .x
       }
     
-    if(factor_names == "1") XA = X[,1, drop=FALSE]
-    else XA = X[,c(1,2), drop=FALSE]
+    # Create variability design matrix (Issue 2: Use formula_variability)
+    Xa =
+      sample_data %>%
+      model.matrix(formula_variability, data=.) %>%
+      apply(2, function(x) {
+        
+        if(sd(x)==0 ) x
+        else x |> scale(scale=FALSE)
+        
+      } ) %>%
+      {
+        .x = (.)
+        rownames(.x) = sample_data %>% pull(!!.sample)
+        .x
+      }
     
-    XA = XA |> 
-      as_tibble()  |> 
+    # Unique variability design (for compatibility)
+    XA = Xa %>%
+      as_tibble() %>%
       distinct()
     
     cell_cluster_names =
@@ -243,21 +810,44 @@ data_simulation_to_model_input =
       arrange(!!.cell_type) %>%
       pull(!!.cell_type)
     
-    coefficients =
+    # Extract coefficients
+    # .coefficients is a quosure pointing to column names like c(b_0, b_1)
+    # Use quo_names to extract the actual column names from the quosure
+    # If .coefficients is NULL quosure, skip coefficient extraction
+    if(rlang::quo_is_null(.coefficients)) {
+      coeff_names = character(0)
+      coefficients = matrix(nrow = 0, ncol = 0)
+    } else {
+      coeff_names = quo_names(.coefficients)
+      
+      if(length(coeff_names) > 0) {
+        # Pivot to long format: cell_type | coefficient_name | value
+        coefficients_long =
       .data %>%
-      select(!!.cell_type, !!.coefficients) %>%
-      unnest(!!.coefficients) %>%
+          select(!!.cell_type, all_of(coeff_names)) %>%
       distinct() %>%
       arrange(!!.cell_type) %>%
-      pivot_wider(names_from = quo_name(.cell_type), values_from = !!.coefficients) %>%
-      column_to_rownames(quo_name(.cell_type)) %>%
-      t()
+          pivot_longer(cols = all_of(coeff_names), names_to = "coefficient_name", values_to = "value")
+        
+        # Pivot to wide format: coefficient_name | (cell_type_1) | (cell_type_2) | ...
+        coefficients_wide = coefficients_long %>%
+          pivot_wider(names_from = quo_name(.cell_type), values_from = value) %>%
+          column_to_rownames("coefficient_name") %>%
+          as.matrix()
+        
+        # Transpose to get: rows = coefficients, columns = cell_types
+        coefficients = t(coefficients_wide)
+      } else {
+        coefficients = matrix(nrow = 0, ncol = 0)
+      }
+    }
     
     list(
       N = .data %>% distinct(!!.sample) %>% nrow(),
       M = .data %>% distinct(!!.cell_type) %>% nrow(),
       exposure = .data %>% distinct(!!.sample, !!.exposure) %>% arrange(!!.sample) %>% pull(!!.exposure),
       X = X,
+      Xa = Xa,  # Issue 2: Add Xa for variability design
       XA = XA,
       C = ncol(X),
       A =  ncol(XA),
diff --git a/R/utilities.R b/R/utilities.R
index 3c6ca7b4..1fed492b 100755
--- a/R/utilities.R
+++ b/R/utilities.R
@@ -664,7 +664,62 @@ calculate_na_fraction_contribution = function(my_design_matrix, na_cols, design_
     bind_rows()
 }
 
-
+#' Normalize sum_to_zero_vector parameters in posterior draws
+#'
+#' @description
+#' This function normalizes sum_to_zero_vector parameters in Stan posterior draws to ensure
+#' they sum to exactly zero. This fixes floating-point precision issues when using
+#' generate_quantities with sum_to_zero_vector types.
+#'
+#' @details
+#' Stan's sum_to_zero_vector type has a strict constraint that the sum must be exactly zero.
+#' When posterior draws are saved and reloaded, floating-point precision can cause the sum
+#' to deviate slightly from zero, causing generate_quantities to fail. This function fixes
+#' this by setting the last element of each vector to exactly the negative sum of all others,
+#' guaranteeing a zero sum.
+#'
+#' @param draws A matrix of posterior draws from a Stan model
+#' @param param_pattern A regular expression pattern matching the parameter names to normalize
+#'   (e.g., "^beta_raw\\[" for beta_raw parameters)
+#'
+#' @return A matrix with normalized parameters that sum to exactly zero
+#'
+#' @keywords internal
+#' @noRd
+#'
+normalize_sum_to_zero_params = function(draws, param_pattern) {
+  param_cols = grep(param_pattern, colnames(draws))
+  if(length(param_cols) == 0) return(draws)
+  
+  # Parse parameter names to extract indices
+  param_names = colnames(draws)[param_cols]
+  # Extract first index (grouping dimension) and second index (M dimension)
+  first_indices = gsub(paste0("^", param_pattern, "\\[(\\d+),.*"), "\\1", param_names)
+  second_indices = gsub(paste0("^", param_pattern, "\\[\\d+,(\\d+)\\]"), "\\1", param_names)
+  
+  # Group by first index and normalize each group
+  unique_groups = unique(first_indices)
+  for(group_idx in unique_groups) {
+    group_cols = param_cols[first_indices == group_idx]
+    if(length(group_cols) > 1) {
+      # Get the M indices for this group to find the last element
+      m_indices = as.integer(second_indices[first_indices == group_idx])
+      last_m_idx = max(m_indices)
+      last_col_idx = group_cols[m_indices == last_m_idx]
+      
+      # Normalize each row's vector to sum to exactly zero
+      for(i in seq_len(nrow(draws))) {
+        vec = as.numeric(draws[i, group_cols])
+        # Calculate sum of all but the last element
+        sum_others = sum(vec[-which(group_cols == last_col_idx)])
+        # Set last element to exactly negative sum of others (guarantees zero sum)
+        vec[which(group_cols == last_col_idx)] = -sum_others
+        draws[i, group_cols] = vec
+      }
+    }
+  }
+  return(draws)
+}
 
 #' @importFrom purrr when
 #' @importFrom stats model.matrix
diff --git a/dev/sccomp_simulate.qmd b/dev/sccomp_simulate.qmd
new file mode 100644
index 00000000..aeeb9182
--- /dev/null
+++ b/dev/sccomp_simulate.qmd
@@ -0,0 +1,859 @@
+---
+title: "Data Simulation with sccomp_simulate"
+author: "sccomp"
+date: today
+format:
+  html:
+    toc: true
+    code-fold: true
+    code-tools: true
+---
+
+## Introduction
+
+The `sccomp_simulate` function allows you to simulate compositional count data from a fitted sccomp model. This vignette demonstrates three main use cases:
+
+1. **Posterior Predictive Checks**: Simulate data using the full fitted model (all parameters from posterior distribution)
+2. **New Data Simulation**: Simulate data for new sample conditions using the fitted model
+3. **Custom Coefficients**: Simulate data with user-specified coefficients for specific scenarios
+
+## Setup
+
+```{r setup, message=FALSE, warning=FALSE}
+# Load development version to ensure latest features
+if(requireNamespace("devtools", quietly = TRUE)) {
+  devtools::load_all(quiet = TRUE)
+}
+library(sccomp)
+library(dplyr)
+library(ggplot2)
+library(tidyr)
+```
+
+Load example data:
+
+```{r load-data, message=FALSE}
+data("counts_obj")
+```
+
+## 1. Posterior Predictive Check (Using Everything from Fit)
+
+The simplest use case is to simulate data using all parameters from the fitted model. This is useful for posterior predictive checks to validate model fit.
+
+### Fit the Model
+
+```{r fit-model, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+estimate = sccomp_estimate(
+  counts_obj,
+  formula_composition = ~ type, 
+  formula_variability = ~ 1,
+  sample = "sample", 
+  cell_group = "cell_group", 
+  abundance = "count",
+  
+  verbose = FALSE,
+  max_sampling_iterations = 2000
+)
+```
+
+### Simulate from Posterior Distribution
+
+When `new_data = NULL` and `coefficients = NULL`, the function uses the original data and samples from the posterior distribution of all model parameters:
+
+```{r posterior-predictive, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Simulate using all parameters from the fitted model
+simulated_data = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = NULL,        # Use original data
+  coefficients = NULL,    # Use posterior beta_raw
+  
+  number_of_draws = 10     # Generate 10 replicates
+)
+
+# View the simulated data
+head(simulated_data, 20)
+```
+
+The simulated data includes:
+- Original columns from the input data
+- `generated_counts`: Simulated cell counts
+- `replicate`: Which draw from the posterior (1, 2, 3, etc.)
+
+### Visualize Posterior Predictive Check
+
+```{r visualize-ppc, eval = instantiate::stan_cmdstan_exists(), message=FALSE, fig.height=6}
+# Compare observed vs simulated
+comparison = 
+  bind_rows(
+    counts_obj |> 
+      mutate(type = "observed") |>
+      select(sample, cell_group, count, type),
+    simulated_data |>
+      filter(replicate == 1) |>
+      mutate(type = "simulated") |>
+      select(sample, cell_group, generated_counts, type) |>
+      rename(count = generated_counts)
+  )
+
+comparison |>
+  ggplot(aes(x = cell_group, y = count, fill = type)) +
+  geom_boxplot(position = "dodge", alpha = 0.7, outlier.alpha = 0.3) +
+  scale_y_continuous(trans = "log1p") +
+  theme_minimal() +
+  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
+  labs(title = "Posterior Predictive Check: Observed vs Simulated",
+       subtitle = "All samples combined",
+       x = "Cell Group", y = "Count (log1p scale)")
+```
+
+## 2. Simulate for New Data (New Sample Conditions)
+
+You can simulate data for new sample conditions (e.g., new treatment groups, new time points) while using the fitted model parameters.
+
+### Create New Data
+
+When `coefficients = NULL`, `new_data` must include the `cell_group` column. We can expand the sample data to include all cell groups:
+
+```{r new-data-example, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Create new sample conditions
+sample_conditions = tibble(
+  sample = c("new_sample_1", "new_sample_2", "new_sample_3"),
+  type = factor(c("normal", "cancer", "cancer"))
+)
+
+# Expand to include all cell groups from the fitted model
+# Get unique cell groups from original data
+all_cell_groups = counts_obj |>
+  distinct(cell_group) |>
+  pull(cell_group)
+
+# Create full new_data with all combinations
+new_samples = sample_conditions |>
+  tidyr::crossing(cell_group = all_cell_groups)
+
+# Simulate for new samples using posterior parameters
+simulated_new = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = new_samples,
+  coefficients = NULL,    # Still use posterior beta_raw
+  
+  number_of_draws = 10
+)
+
+head(simulated_new, 20)
+```
+
+Note: When `coefficients = NULL`, the function automatically expands `new_data` to include all cell groups from the fitted model. The design matrix for the new samples is created based on the `formula_composition` and `formula_variability`.
+
+### Simulate for Subset of Samples
+
+You can also simulate for a subset of the original samples. When `coefficients = NULL`, the data must include `cell_group`:
+
+```{r subset-samples, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Select a subset of original samples (with cell_group)
+subset_samples = counts_obj |>
+  filter(sample %in% unique(counts_obj$sample)[1:3]) |>  # First 3 samples
+  select(sample, type, cell_group)
+
+# Simulate for subset
+simulated_subset = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = subset_samples,
+  coefficients = NULL,
+  cores = 1
+)
+
+head(simulated_subset, 15)
+```
+
+## 3. Simulate with Custom Coefficients
+
+You can specify custom coefficients to simulate specific scenarios (e.g., hypothetical treatment effects, what-if analyses).
+
+### Create Coefficients Table
+
+The coefficients table must contain:
+- A column matching the `cell_group` column name
+- Columns matching the design matrix column names (e.g., `(Intercept)`, `typecancer`)
+
+```{r create-coefficients, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Get design matrix column names from the fitted model
+# These typically include "(Intercept)" and factor levels
+design_cols = c("(Intercept)", "typecancer")
+
+# Create coefficients table
+# Example: Set all coefficients to zero (null scenario)
+coeffs_null = counts_obj |>
+  distinct(cell_group) |>
+  mutate(
+    `(Intercept)` = 0,
+    `typecancer` = 0
+  )
+
+coeffs_null
+```
+
+### Verify Zero Coefficients Give Equal Mean Proportions
+
+When all coefficients are zero (both intercept and treatment effects), all cell types should have equal mean proportions (1/M, where M is the number of cell types). This is because with zero log-odds, `softmax(0, 0, ..., 0) = (1/M, 1/M, ..., 1/M)`.
+
+Let's verify this by simulating 30 samples:
+
+```{r verify-zero-coefficients, eval = instantiate::stan_cmdstan_exists(), message=FALSE, fig.height=6}
+# Simulate with zero coefficients across 30 samples
+sample_data_many = tibble(
+  sample = paste0("sample_", 1:30),
+  type = factor(rep(c("normal", "cancer"), each = 15))
+)
+
+sim_zero = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = sample_data_many,
+  coefficients = coeffs_null,
+  
+  number_of_draws = 10
+)
+
+# Expected proportion (1/M where M is number of cell types)
+M = nrow(coeffs_null)
+expected_prop = 1 / M
+
+cat("Expected proportion per cell type (1/M):", round(expected_prop, 4), "\n")
+cat("Number of cell types (M):", M, "\n\n")
+
+# Visualize with boxplots - all cell types should have similar mean proportions
+sim_zero |>
+  ggplot(aes(x = cell_group, y = generated_proportions)) +
+  geom_boxplot(fill = "steelblue", alpha = 0.7, outlier.alpha = 0.3) +
+  geom_hline(
+    yintercept = expected_prop, 
+    linetype = "dashed", 
+    color = "red", 
+    linewidth = 1
+  ) +
+  theme_minimal() +
+  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
+  labs(
+    title = "Proportions with Zero Coefficients (30 samples)",
+    subtitle = paste("Expected mean:", round(expected_prop, 4), "(red dashed line) - All cell types should align"),
+    x = "Cell Group",
+    y = "Proportion",
+    caption = "With all coefficients = 0, all cell types should have equal mean proportions (1/M)"
+  )
+```
+
+**Note**: Individual samples may show variation due to:
+- Random effects (if present in the model)
+- Variability parameters (alpha) affecting the beta-binomial distribution
+- Sampling variability
+
+However, the **mean** across many samples should converge to 1/M for each cell type.
+
+### Simulate with Custom Coefficients
+
+```{r simulate-custom-coefficients, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Create new data for samples
+new_samples_coef = tibble(
+  sample = c("scenario_1", "scenario_2"),
+  type = factor(c("normal", "cancer"))
+)
+
+# Simulate with custom coefficients
+simulated_custom = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = new_samples_coef,
+  coefficients = coeffs_null,
+  
+  number_of_draws = 10
+)
+
+head(simulated_custom, 20)
+
+# Note: With zero coefficients, mean proportions should be equal across cell types
+# Individual samples may vary due to random effects and variability parameters,
+# but the mean across many samples should converge to 1/M for each cell type
+```
+
+### Simulate Different Scenarios
+
+You can create multiple coefficient tables to simulate different scenarios. When using `new_data`, it's useful to simulate from at least 30 samples to see clear patterns:
+
+```{r multiple-scenarios, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Create sample data with 30 samples for better visualization
+new_samples_scenarios = tibble(
+  sample = paste0("scenario_sample_", 1:30),
+  type = factor(rep(c("normal", "cancer"), each = 15))
+)
+
+# Scenario 1: Strong cancer effect for specific cell type
+coeffs_strong_effect = counts_obj |>
+  distinct(cell_group) |>
+  mutate(
+    `(Intercept)` = 0,
+    `typecancer` = if_else(cell_group == "B1", 2.0, 0.0)  # Strong effect for B1
+  )
+
+# Scenario 2: Moderate effect across all cell types
+coeffs_moderate = counts_obj |>
+  distinct(cell_group) |>
+  mutate(
+    `(Intercept)` = 0,
+    `typecancer` = 0.5  # Moderate effect for all
+  )
+
+# Simulate both scenarios
+scenario1 = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = new_samples_scenarios,
+  coefficients = coeffs_strong_effect,
+  cores = 1
+) |> mutate(scenario = "Strong B1 effect")
+
+scenario2 = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = new_samples_scenarios,
+  coefficients = coeffs_moderate,
+  cores = 1
+) |> mutate(scenario = "Moderate effect")
+
+# Compare scenarios - focus on cancer samples
+bind_rows(scenario1, scenario2) |>
+  filter(type == "cancer") |>
+  ggplot(aes(x = cell_group, y = generated_counts, fill = scenario)) +
+  geom_boxplot(position = "dodge", alpha = 0.7) +
+  scale_y_continuous(trans = "log1p") +
+  theme_minimal() +
+  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
+  labs(title = "Comparison of Different Coefficient Scenarios (30 cancer samples)",
+       x = "Cell Group", y = "Simulated Counts (log1p scale)")
+```
+
+## 4. Combining New Data with Custom Coefficients
+
+When both `new_data` and `coefficients` are provided, `new_data` can omit the `cell_group` column. The function automatically expands `new_data` to include all cell groups from the `coefficients` table using a cross join:
+
+```{r combine-newdata-coefficients, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Sample-specific data (no cell_group needed)
+sample_data = tibble(
+  sample = c("experiment_1", "experiment_2", "experiment_3"),
+  type = factor(c("normal", "cancer", "cancer"))
+)
+
+# Cell-type-specific coefficients
+coeffs_table = counts_obj |>
+  distinct(cell_group) |>
+  mutate(
+    `(Intercept)` = 0,
+    `typecancer` = 0
+  )
+
+# Simulate: new_data specifies samples, coefficients specifies cell types
+simulated_combined = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = sample_data,      # Sample-specific
+  coefficients = coeffs_table,  # Cell-type-specific
+  cores = 1
+)
+
+# The result has all combinations of samples x cell_groups
+simulated_combined |>
+  select(sample, type, cell_group, generated_counts) |>
+  head(20)
+```
+
+## 5. Advanced: Multiple Draws and Variability
+
+### Custom Mean-Dispersion Association Slope
+
+The `mean_dispersion_slope` parameter allows you to override the slope of the mean-dispersion association from the fitted model. This controls how variability (dispersion) relates to mean abundance. If `NULL`, the slope from the fitted model (`prec_coeff[2]`) is used.
+
+We can visualize the effect by comparing 2D plots from `plot_2D_intervals()` showing the mean-variance association:
+
+```{r mean-dispersion-slope, eval = instantiate::stan_cmdstan_exists(), message=FALSE, fig.height=8}
+# First, get the 2D plot from the original fit (shows the association)
+# The original estimate was fitted with formula_variability = ~ 1, but we need ~ type for 2D plots
+# So we'll refit with ~ type for variability
+estimate_with_variability = sccomp_estimate(
+  counts_obj,
+  formula_composition = ~ type,
+  formula_variability = ~ type,  # Include type in variability to get 2D plots for factor
+  sample = "sample",
+  cell_group = "cell_group",
+  abundance = "count",
+  verbose = FALSE
+)
+
+estimate_tested = estimate_with_variability |>
+  sccomp_test()
+
+plot_original = plot_2D_intervals(estimate_tested) +
+  labs(title = "Original Fit: Mean-Dispersion Association")
+
+# Simulate data with zero slope (no association)
+# When mean_dispersion_slope is provided, the intercept is automatically preserved from the original fit
+# Use formula_variability = ~ type to get 2D plots for the factor as well
+# Use multiple replicates and treat them as separate samples for better fit
+sim_no_association = sccomp_simulate(
+  estimate_with_variability,  # Use the refitted estimate with ~ type variability
+  formula_composition = ~ type,
+  formula_variability = ~ type,  # Include type in variability to get 2D plots for factor
+  new_data = NULL,
+  coefficients = NULL,
+  mean_dispersion_slope = 0.0,  # No association (intercept is automatically preserved)
+  
+  number_of_draws = 10  # Use multiple replicates
+)
+
+# Get the original fit parameters for comparison
+fit_original = attr(estimate_tested, "fit")
+prec_coeff_original = fit_original$summary("prec_coeff")
+cat("Original fit: intercept =", round(prec_coeff_original$mean[1], 2), 
+    ", slope =", round(prec_coeff_original$mean[2], 2), "\n")
+cat("Simulation used: intercept =", round(prec_coeff_original$mean[1], 2), 
+    ", slope = 0.0\n")
+
+# Fit the simulated data (with zero association)
+# Treat each replicate as a separate sample to get more data
+sim_data_for_fit = sim_no_association |>
+  mutate(sample = paste0(sample, "_rep", replicate)) |>
+  select(sample, cell_group, type, count = generated_counts)
+
+estimate_simulated = sccomp_estimate(
+  sim_data_for_fit,
+  formula_composition = ~ type,
+  formula_variability = ~ type,  # Match the variability formula used in simulation
+  sample = "sample",
+  cell_group = "cell_group",
+  abundance = "count",
+  
+  verbose = FALSE
+) |>
+  sccomp_test()
+
+# Get the fitted parameters from simulated data
+fit_simulated = attr(estimate_simulated, "fit")
+prec_coeff_simulated = fit_simulated$summary("prec_coeff")
+cat("Fitted from simulated data: intercept =", round(prec_coeff_simulated$mean[1], 2), 
+    ", slope =", round(prec_coeff_simulated$mean[2], 2), "\n")
+
+plot_simulated = plot_2D_intervals(estimate_simulated) +
+  labs(title = "Simulated with Slope = 0: No Mean-Dispersion Association")
+
+# Display plots side by side
+library(patchwork)
+plot_original / plot_simulated
+
+# Note: The regression line in the plot shows the ESTIMATED parameters from fitting the simulated data,
+# not the simulation parameters themselves. The simulation correctly used:
+# - intercept = original fit intercept (preserved automatically)
+# - slope = 0.0 (as specified)
+# However, when we fit the simulated data, the model estimates new parameters from the data structure.
+# The fitted intercept may differ from the simulation intercept because:
+# - The simulated data has a different structure (slope = 0 changes the alpha-beta relationship)
+# - The model estimates parameters that best fit the observed data
+# The "(Intercept, adjusted)" facet shows the variability effect after adjusting for
+# the estimated association, which should be closer to horizontal when slope = 0.
+```
+
+### Simulate with Non-Zero Coefficients and No Association
+
+We can also simulate with specific coefficients (e.g., a treatment effect) while removing the mean-dispersion association:
+
+```{r simulate-with-coefficients-no-association, eval = instantiate::stan_cmdstan_exists(), message=FALSE, fig.height=10}
+# Create coefficients table with a non-zero effect for typecancer
+# We'll use a moderate effect size
+# Coefficients must sum to zero for compositional models
+coeffs_no_assoc = estimate_tested |>
+  distinct(cell_group) |>
+  mutate(
+    `(Intercept)` = 0,
+    typecancer_raw = if_else(cell_group == "B1", 1.5, 0)  # B1 has positive effect
+  ) |>
+  mutate(
+    typecancer = typecancer_raw - mean(typecancer_raw)  # Normalize to sum to zero
+  ) |>
+  select(-typecancer_raw)
+
+# Simulate with these coefficients and slope = 0 (no association)
+# Use formula_variability = ~ type to get 2D plots for the factor as well
+sim_coeffs_no_assoc = sccomp_simulate(
+  estimate_with_variability,  # Use the refitted estimate with ~ type variability
+  formula_composition = ~ type,
+  formula_variability = ~ type,  # Include type in variability to get 2D plots for factor
+  new_data = NULL,
+  coefficients = coeffs_no_assoc,
+  mean_dispersion_slope = 0.0,  # No association
+  
+  number_of_draws = 10
+)
+
+cat("Simulation with coefficients and no association completed\n")
+cat("Coefficients used: B1 typecancer =", coeffs_no_assoc$typecancer[coeffs_no_assoc$cell_group == "B1"], "\n")
+
+# Fit the simulated data
+sim_data_coeffs = sim_coeffs_no_assoc |>
+  mutate(sample = paste0(sample, "_rep", replicate)) |>
+  select(sample, cell_group, type, count = generated_counts)
+
+estimate_simulated_coeffs = sccomp_estimate(
+  sim_data_coeffs,
+  formula_composition = ~ type,
+  formula_variability = ~ type,  # Match the variability formula used in simulation
+  sample = "sample",
+  cell_group = "cell_group",
+  abundance = "count",
+  
+  verbose = FALSE
+) |>
+  sccomp_test()
+
+# Get fitted parameters
+fit_simulated_coeffs = attr(estimate_simulated_coeffs, "fit")
+prec_coeff_simulated_coeffs = fit_simulated_coeffs$summary("prec_coeff")
+cat("Fitted from simulated data (with coefficients): intercept =", round(prec_coeff_simulated_coeffs$mean[1], 2), 
+    ", slope =", round(prec_coeff_simulated_coeffs$mean[2], 2), "\n")
+
+plot_simulated_coeffs = plot_2D_intervals(estimate_simulated_coeffs) +
+  labs(title = "Simulated with Coefficients (b1≠0) and Slope = 0: No Mean-Dispersion Association")
+
+# Display all three plots: original, simulated without coefficients, simulated with coefficients
+plot_original / plot_simulated / plot_simulated_coeffs
+```
+
+The three plots show:
+1. **Top**: Original fit with mean-dispersion association (sloped regression line)
+2. **Middle**: Simulated with zero coefficients and slope = 0 (no association, horizontal line)
+3. **Bottom**: Simulated with non-zero coefficients (B1 has positive effect) and slope = 0 (no association, horizontal line)
+
+Even with non-zero coefficients, when `mean_dispersion_slope = 0`, there should be no association between mean and dispersion, resulting in a horizontal regression line.
+```
+
+### Multiple Posterior Draws
+
+You can simulate multiple draws from the posterior distribution:
+
+```{r multiple-draws, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Simulate 5 draws from posterior
+simulated_multi = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = NULL,
+  coefficients = NULL,
+  
+  number_of_draws = 10
+)
+
+# Check replicate column
+simulated_multi |>
+  distinct(replicate) |>
+  count()
+```
+
+### Variability Multiplier
+
+You can artificially increase variability for benchmarking or sensitivity analysis:
+
+```{r variability-multiplier, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Standard simulation
+sim_normal = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = NULL,
+  coefficients = NULL,
+  
+  variability_multiplier = 1  # Default
+) |> mutate(variability = "Normal")
+
+# Increased variability
+sim_high = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = NULL,
+  coefficients = NULL,
+  
+  variability_multiplier = 10  # 10x variability
+) |> mutate(variability = "High")
+
+# Compare
+bind_rows(sim_normal, sim_high) |>
+  filter(replicate == 1) |>
+  ggplot(aes(x = cell_group, y = generated_counts, fill = variability)) +
+  geom_boxplot(position = "dodge", alpha = 0.7) +
+  scale_y_continuous(trans = "log1p") +
+  theme_minimal() +
+  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
+  labs(title = "Effect of Variability Multiplier",
+       x = "Cell Group", y = "Simulated Counts (log1p scale)")
+```
+
+## 6. Comparing Simulated Scenarios: Significant Cell Type Analysis
+
+A powerful use case is comparing simulations with different coefficient scenarios, especially for cell types that show significant differences. This helps visualize the impact of treatment effects vs. null scenarios.
+
+### Find a Significantly Different Cell Type
+
+First, let's identify which cell types are significantly different between conditions:
+
+```{r find-significant, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Test for significance
+test_results = estimate |>
+  sccomp_test()
+
+# Find significantly different cell types (e.g., FDR < 0.05)
+significant_cells = test_results |>
+  filter(c_FDR < 0.05 | is.na(c_FDR)) |>
+  arrange(c_FDR) |>
+  head(3)
+
+# Display significant cell types
+significant_cells |>
+  select(cell_group, parameter, c_effect, c_lower, c_upper, c_FDR)
+
+# Pick the most significant one for demonstration
+if(nrow(significant_cells) > 0) {
+  selected_cell = significant_cells$cell_group[1]
+  cat("\nSelected cell type for comparison:", selected_cell, "\n")
+} else {
+  # If no significant, pick one with largest effect
+  selected_cell = test_results |>
+    arrange(desc(abs(c_effect))) |>
+    head(1) |>
+    pull(cell_group)
+  cat("\nNo significant cell types found. Using cell type with largest effect:", selected_cell, "\n")
+}
+```
+
+### Simulate with Posterior Parameters (Real Effect)
+
+```{r simulate-posterior, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Simulate using all parameters from fit (includes real effects)
+sim_posterior = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = NULL,
+  coefficients = NULL,  # Use posterior beta_raw (real effects)
+  
+  number_of_draws = 10
+)
+
+# Calculate proportions for the selected cell type
+# Use generated_proportions directly from the simulation output
+proportions_posterior = sim_posterior |>
+  filter(cell_group == selected_cell) |>
+  select(sample, type, replicate, generated_proportions) |>
+  rename(proportion = generated_proportions) |>
+  mutate(scenario = "Posterior (real effects)")
+```
+
+### Simulate with Different Coefficient Scenarios
+
+```{r simulate-scenarios, eval = instantiate::stan_cmdstan_exists(), message=FALSE}
+# Create coefficients tables for different scenarios
+coeffs_null = counts_obj |>
+  distinct(cell_group) |>
+  mutate(
+    `(Intercept)` = 0,
+    `typecancer` = 0  # No effect
+  )
+
+coeffs_positive = counts_obj |>
+  distinct(cell_group) |>
+  mutate(
+    `(Intercept)` = 0,
+    `typecancer` = if_else(cell_group == selected_cell, 1.5, 0.0)  # Positive effect for selected cell type
+  ) |>
+  mutate(`typecancer` = `typecancer` - mean(`typecancer`))  # Normalize to sum to zero
+
+coeffs_negative = counts_obj |>
+  distinct(cell_group) |>
+  mutate(
+    `(Intercept)` = 0,
+    `typecancer` = if_else(cell_group == selected_cell, -1.5, 0.0)  # Negative effect for selected cell type
+  ) |>
+  mutate(`typecancer` = `typecancer` - mean(`typecancer`))  # Normalize to sum to zero
+
+# Get sample data
+sample_data = counts_obj |>
+  distinct(sample, type)
+
+# Simulate with different coefficients
+sim_null = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = sample_data,
+  coefficients = coeffs_null,  # Zero coefficients (no effect)
+  
+  number_of_draws = 10
+)
+
+sim_positive = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = sample_data,
+  coefficients = coeffs_positive,  # Positive coefficients
+  
+  number_of_draws = 10
+)
+
+sim_negative = sccomp_simulate(
+  estimate,
+  formula_composition = ~ type,
+  formula_variability = ~ 1,
+  new_data = sample_data,
+  coefficients = coeffs_negative,  # Negative coefficients
+  
+  number_of_draws = 10
+)
+
+# Calculate proportions for the selected cell type
+# Use generated_proportions directly from the simulation output
+proportions_null = sim_null |>
+  filter(cell_group == selected_cell) |>
+  select(sample, type, replicate, generated_proportions) |>
+  rename(proportion = generated_proportions) |>
+  mutate(scenario = "Zero (null)")
+
+proportions_positive = sim_positive |>
+  filter(cell_group == selected_cell) |>
+  select(sample, type, replicate, generated_proportions) |>
+  rename(proportion = generated_proportions) |>
+  mutate(scenario = "Positive (+1.5)")
+
+proportions_negative = sim_negative |>
+  filter(cell_group == selected_cell) |>
+  select(sample, type, replicate, generated_proportions) |>
+  rename(proportion = generated_proportions) |>
+  mutate(scenario = "Negative (-1.5)")
+```
+
+### Compare Proportions Between Scenarios
+
+```{r compare-scenarios, eval = instantiate::stan_cmdstan_exists(), message=FALSE, fig.height=6}
+# Combine all scenarios
+comparison = bind_rows(
+  proportions_posterior,
+  proportions_null,
+  proportions_positive,
+  proportions_negative
+)
+
+# Calculate mean proportions by type and scenario
+summary_comparison = comparison |>
+  group_by(type, scenario) |>
+  summarise(
+    mean_proportion = mean(proportion, na.rm = TRUE),
+    sd_proportion = sd(proportion, na.rm = TRUE),
+    lower = quantile(proportion, 0.025, na.rm = TRUE),
+    upper = quantile(proportion, 0.975, na.rm = TRUE),
+    .groups = "drop"
+  )
+
+# Visualize comparison
+p1 = comparison |>
+  ggplot(aes(x = type, y = proportion, fill = scenario)) +
+  geom_boxplot(position = "dodge", alpha = 0.7, outlier.alpha = 0.3) +
+  facet_wrap(~scenario, ncol = 2) +
+  theme_minimal() +
+  labs(
+    title = paste("Proportion Comparison:", selected_cell),
+    subtitle = "Posterior (real) vs Zero vs Positive (+1.5) vs Negative (-1.5) coefficients",
+    x = "Type", 
+    y = "Proportion",
+    fill = "Scenario"
+  ) +
+  theme(legend.position = "none")
+
+# Show mean differences
+p2 = summary_comparison |>
+  ggplot(aes(x = type, y = mean_proportion, fill = scenario)) +
+  geom_col(position = "dodge", alpha = 0.7) +
+  geom_errorbar(
+    aes(ymin = lower, ymax = upper),
+    position = position_dodge(width = 0.9),
+    width = 0.2
+  ) +
+  theme_minimal() +
+  labs(
+    title = "Mean Proportions with 95% CI",
+    x = "Type",
+    y = "Mean Proportion",
+    fill = "Scenario"
+  )
+
+# Display plots
+p1
+p2
+
+# Calculate difference between types for each scenario
+# First, get mean proportion per type, scenario, and replicate
+type_means = comparison |>
+  group_by(type, scenario, replicate) |>
+  summarise(mean_prop = mean(proportion, na.rm = TRUE), .groups = "drop")
+
+# Get type levels to use correct column names after pivot
+type_levels = sort(unique(as.character(type_means$type)))
+type1 = type_levels[1]
+type2 = type_levels[2]
+
+# Calculate difference between types for each replicate
+differences = type_means |>
+  pivot_wider(names_from = type, values_from = mean_prop) |>
+  mutate(
+    difference = .data[[type2]] - .data[[type1]],
+    abs_difference = abs(difference)
+  )
+
+# Show difference summary
+cat("\nDifference between", type2, "and", type1, ":\n")
+differences |>
+  group_by(scenario) |>
+  summarise(
+    mean_diff = mean(difference, na.rm = TRUE),
+    sd_diff = sd(difference, na.rm = TRUE),
+    lower_diff = quantile(difference, 0.025, na.rm = TRUE),
+    upper_diff = quantile(difference, 0.975, na.rm = TRUE),
+    .groups = "drop"
+  ) |>
+  print()
+```
+
+This comparison clearly shows:
+- **Posterior simulation**: Reflects the actual treatment effects estimated from the data
+- **Zero coefficients simulation**: Shows what would happen if there were no treatment effects (null scenario)
+
+The difference between these scenarios demonstrates the magnitude of the treatment effect for the selected cell type.
+
+## 7. Summary
+
+The `sccomp_simulate` function provides flexible data simulation capabilities:
+
+- **Posterior Predictive Checks**: Use `new_data = NULL, coefficients = NULL` to simulate from the full posterior
+- **New Conditions**: Provide `new_data` with new sample conditions while using posterior parameters
+- **Custom Scenarios**: Provide `coefficients` table to specify custom effect sizes
+- **Combined**: Use both `new_data` and `coefficients` for maximum flexibility
+- **Effect Comparison**: Compare different coefficient scenarios to visualize treatment effects
+
+All simulations respect the compositional nature of the data and include proper uncertainty quantification through the posterior distribution.
+
diff --git a/inst/NEWS.rd b/inst/NEWS.rd
index aaba2b49..8faf7173 100644
--- a/inst/NEWS.rd
+++ b/inst/NEWS.rd
@@ -1,12 +1,11 @@
 \name{NEWS}
 \title{News for Package \pkg{sccomp}}
 
-\section{News in version 2.1.25}{
+\section{News in version 2.1.28}{
 \itemize{
-    \item Improved plotting significance colouring controls. Added \code{significance_statistic} to \code{sccomp_boxplot()} with default \code{c("pH0", "FDR")}, so colouring now defaults to posterior probability while still supporting FDR-based colouring.
-    \item Bayesian FDR messaging is now shown only when FDR is selected, avoiding FDR-specific text when probability-based significance is used.
-    \item Fixed the package version dot-numbering (\url{https://github.com/MangiolaLaboratory/sccomp/issues/256}).
-}}
+    \item **Simulation improvements:** Changed default \code{variability_multiplier} parameter from 5 to 1 in \code{sccomp_simulate()} and \code{simulate_data()}. The new default preserves the fitted model's variability by default, rather than artificially increasing it. Users can still specify custom values greater than 1 to increase variability for benchmarking or sensitivity analysis.
+    \item Enhanced simulation documentation. Improved clarity on how \code{mean_dispersion_slope} parameter affects both intercept and factor-level variability parameters in the variability formula. Updated vignette examples to use consistent \code{number_of_draws} values.
+    \item Clarified behavior of mean-dispersion association parameters in simulation. The \code{mean_dispersion_slope} affects the slope for both intercept columns (alpha1) and factor columns (alpha2, e.g., type), with intercept only added to intercept columns as per the model specification.
 
 \section{News in version 2.1.27}{
 \itemize{
@@ -16,6 +15,14 @@
     \item **Technical implementation:** Modified the Stan model \code{glm_multi_beta_binomial_generate_data.stan} to save the linear predictors as \code{mu_unconstrained} before applying the softmax transformation. Updated \code{sccomp_predict()} to extract these values using both summary and draw modes, and updated \code{sccomp_calculate_residuals()} to include them in the output. The implementation includes backward compatibility - if the unconstrained predictors are not available (e.g., with older cached models), the functions gracefully handle their absence.
     \item **Documentation and testing:** Updated function documentation to describe the new output columns. Added comprehensive unit tests in \code{test-unconstrained-predictors.R} covering all output modes (summary with robust/non-robust statistics, draws mode) and verifying that unconstrained predictors are distinct from proportions and properly formatted. All 43 unit tests passing.
     \item **Why this feature is useful:** The unconstrained predictors provide insight into the model's linear scale, which can be valuable for: (1) understanding effect sizes on the logit scale before transformation, (2) performing downstream analyses that require linear-scale values, (3) debugging model convergence issues, (4) comparing effect magnitudes across different cell types on a common scale, and (5) advanced statistical analyses that benefit from access to both transformed and untransformed predictions.
+
+\section{News in version 2.1.25}{
+\itemize{
+    \item Improved plotting significance colouring controls. Added \code{significance_statistic} to \code{sccomp_boxplot()} with default \code{c("pH0", "FDR")}, so colouring now defaults to posterior probability while still supporting FDR-based colouring.
+    \item Bayesian FDR messaging is now shown only when FDR is selected, avoiding FDR-specific text when probability-based significance is used.
+    \item Fixed the package version dot-numbering (\url{https://github.com/MangiolaLaboratory/sccomp/issues/256}).
+}}
+
 }}
 
 \section{News in version 2.1.22}{
diff --git a/inst/stan/glm_multi_beta_binomial_simulate_data.stan b/inst/stan/glm_multi_beta_binomial_simulate_data.stan
index a719643e..15590a23 100755
--- a/inst/stan/glm_multi_beta_binomial_simulate_data.stan
+++ b/inst/stan/glm_multi_beta_binomial_simulate_data.stan
@@ -1,3 +1,7 @@
+functions{
+  #include common_functions.stan
+}
+
 data{
   int N_simulated; // Number of subjects
   int M_simulated; // Number of categories
@@ -5,80 +9,291 @@ data{
   int A_simulated;
   array[N_simulated] int exposure_simulated;
   matrix[N_simulated, C_simulated] X_simulated;
-  matrix[A_simulated, A_simulated] XA_simulated;
-  matrix[C_simulated,M_simulated] beta_simulated;
+  matrix[N_simulated, A_simulated] Xa_simulated; // Variability design matrix for simulated data
+  matrix[A_simulated, A_simulated] XA_simulated; // Unique variability design (for compatibility with old code)
+
   real<lower=0> variability_multiplier;
+  
+  // Optional: provided coefficients to override posterior beta_raw
+  int<lower=0, upper=1> user_provided_beta; // 1 if beta_simulated is provided, 0 to use posterior beta_raw
+  // NOTE: Using vector[M_simulated] instead of sum_to_zero_vector[M_simulated] to avoid floating-point precision issues
+  // Floating-point precision can cause sum_to_zero_vector to fail the strict sum-to-zero constraint when values are passed from R
+  // The sum-to-zero constraint is a modeling constraint (compositional data must sum to 1, so log-ratios sum to 0)
+  // but does not need to be enforced at the Stan type level for generate_quantities
+  // Small precision errors are acceptable since Stan doesn't enforce strict constraints with vector types
+  array[C_simulated] vector[M_simulated] beta_simulated_provided; // Provided coefficients
+  
+  // Optional: provided slope and intercept for mean-dispersion association
+  int<lower=0, upper=1> user_provided_prec_coeff_slope; // 1 if prec_coeff_slope is provided, 0 to use posterior prec_coeff[2]
+  real prec_coeff_slope_provided; // Provided slope for mean-dispersion association (alpha = beta * slope + intercept)
+  int<lower=0, upper=1> user_provided_prec_coeff_intercept; // 1 if prec_coeff_intercept is provided, 0 to use posterior prec_coeff[1]
+  real prec_coeff_intercept_provided; // Provided intercept for mean-dispersion association (alpha = beta * slope + intercept)
 
-	int M;
+  int M;
 	int C;
 	int A;
+	int A_intercept_columns; // How many intercept columns in variability design
+	int intercept_in_design; // Whether intercept is in design
+	int bimodal_mean_variability_association; // Whether to use bimodal association
+	
+	// Indices for subsetting beta and alpha (matching replicate_data approach)
+	int length_X_which;
+	int length_XA_which;
+	array[length_X_which] int X_which;
+	array[length_XA_which] int XA_which;
 	array[2] int ncol_X_random_eff;
 	int is_random_effect;
   array[2] int how_many_factors_in_random_design;
+  
+  // Random effects design matrices (if random effects exist, otherwise empty matrices)
+  // Stan requires at least 1 column, so we use max(1, ncol) - the conditional logic will handle when ncol is 0
+  matrix[N_simulated, max(1, ncol_X_random_eff[1])] X_random_effect_simulated;
+  matrix[N_simulated, max(1, ncol_X_random_eff[2])] X_random_effect_2_simulated;
+  
+  // Covariance setup for random effects (if random effects exist)
+  array[2] int n_groups;
+  array[how_many_factors_in_random_design[1], n_groups[1]] int group_factor_indexes_for_covariance;
+  array[how_many_factors_in_random_design[2], n_groups[2]] int group_factor_indexes_for_covariance_2;
+}
 
-
+transformed data{
+  // Issue 5: Ensure we have valid dimensions for random effect matrices
+  // These are used in generated quantities to subset matrices correctly
+  // ncol_re1 and ncol_re2 represent the actual matrix dimensions (always >= 1)
+  // ncol_X_random_eff[1] and ncol_X_random_eff[2] represent the number of random effect parameters (can be 0)
+  int ncol_re1 = max(1, ncol_X_random_eff[1]);
+  int ncol_re2 = max(1, ncol_X_random_eff[2]);
 }
 
 parameters{
-
-  matrix[C, M-1] beta_raw_raw; // matrix with C rows and number of cells (-1) columns
-  matrix[A, M] alpha; // Variability
+  // NOTE: Using vector[M] instead of sum_to_zero_vector[M] to avoid floating-point precision issues
+  // Floating-point precision can cause sum_to_zero_vector to fail the strict sum-to-zero constraint
+  // when posterior draws are reloaded from disk. The sum-to-zero constraint is a modeling constraint
+  // (compositional data must sum to 1, so log-ratios sum to 0) but does not need to be enforced
+  // at the Stan type level for generate_quantities. Small precision errors are acceptable.
+  array[C] vector[M] beta_raw; // Each row is a vector of length M
+  matrix[A, M] alpha; // Variability - kept in parameters for generate_quantities compatibility
+  // Note: alpha from posterior is NOT used - alpha_simulated is computed from beta in generated quantities
   // To exclude
   array[2] real prec_coeff;
   real<lower=0> prec_sd;
   real<lower=0, upper=1> mix_p;
   
-  // Random intercept // matrix with N_groupings rows and number of cells (-1) columns
-  matrix[ncol_X_random_eff[1] * (is_random_effect>0), M-1] random_effect_raw;
-  matrix[ncol_X_random_eff[2] * (is_random_effect>0), M-1] random_effect_raw_2;
+  // Random intercept - using vector instead of sum_to_zero_vector to avoid precision issues
+  array[ncol_X_random_eff[1] * (is_random_effect>0)] vector[M] random_effect_raw;
+  array[ncol_X_random_eff[2] * (ncol_X_random_eff[2]>0)] vector[M] random_effect_raw_2;
   
   // sd of random intercept
-  array[is_random_effect>0] real random_effect_sigma_mu;
-  array[is_random_effect>0] real random_effect_sigma_sigma;
+  array[2 * (is_random_effect>0)] real random_effect_sigma_mu;
+  array[2 * (is_random_effect>0)] real random_effect_sigma_sigma;
 
 	// Covariance
-  array[M-1 * (is_random_effect>0)] vector[how_many_factors_in_random_design[1]]  random_effect_sigma_raw;
-	array[M-1 * (is_random_effect>0)] cholesky_factor_corr[how_many_factors_in_random_design[1] * (is_random_effect>0)] sigma_correlation_factor;
+  array[M * (is_random_effect>0)] vector[how_many_factors_in_random_design[1]]  random_effect_sigma_raw;
+	array[M * (is_random_effect>0)] cholesky_factor_corr[how_many_factors_in_random_design[1] * (is_random_effect>0)] sigma_correlation_factor;
 
 	// Covariance
-  array[M-1 * (is_random_effect>0)] vector[how_many_factors_in_random_design[2]]  random_effect_sigma_raw_2;
-	array[M-1 * (is_random_effect>0)] cholesky_factor_corr[how_many_factors_in_random_design[2] * (is_random_effect>0)] sigma_correlation_factor_2;
+  array[M * (is_random_effect>0)] vector[how_many_factors_in_random_design[2]]  random_effect_sigma_raw_2;
+	array[M * (is_random_effect>0)] cholesky_factor_corr[how_many_factors_in_random_design[2] * (is_random_effect>0)] sigma_correlation_factor_2;
 
   // If I have just one group
   array[is_random_effect>0] real zero_random_effect;
   
 }
 
+transformed parameters{
+  // Convert vector to matrix
+  // If user_provided_beta is 1, use provided coefficients
+  // Otherwise, use beta_raw from posterior draws
+  // Note: We use vector[M] instead of sum_to_zero_vector[M] to avoid floating-point precision issues
+  // Small precision errors in sum-to-zero constraint are acceptable since Stan doesn't enforce strict constraints
+  matrix[C,M] beta;
+  if(user_provided_beta == 1) {
+    // Use provided coefficients - map them to the full beta matrix using X_which
+    // First, create beta matrix from beta_raw (needed for rows not in X_which)
+    for(c in 1:C) {
+      beta[c,] = to_row_vector(beta_raw[c]);
+    }
+    // Then override the rows specified by X_which with provided coefficients
+    for(i in 1:length_X_which) {
+      int c = X_which[i];
+      // Provided coefficients are vector[M_simulated]
+      // Map them to full M dimensions (no padding needed - X_which ensures correct mapping)
+      for(m in 1:M_simulated) {
+        beta[c, m] = beta_simulated_provided[i][m];
+      }
+      // If M_simulated < M, pad remaining with zeros
+      if(M_simulated < M) {
+        for(m in (M_simulated + 1):M) {
+          beta[c, m] = 0.0;
+        }
+      }
+    }
+  } else {
+    // Use beta_raw from posterior draws
+    for(c in 1:C) {
+      beta[c,] = to_row_vector(beta_raw[c]);
+    }
+  }
+  
+  // Non centered parameterisation SD of random effects (matching main model)
+  array[M * (ncol_X_random_eff[1]> 0)] vector[how_many_factors_in_random_design[1]] random_effect_sigma;
+  if(ncol_X_random_eff[1]> 0) for(m in 1:(M)) random_effect_sigma[m] = random_effect_sigma_mu[1] + random_effect_sigma_sigma[1] * random_effect_sigma_raw[m];
+  if(ncol_X_random_eff[1]> 0) for(m in 1:(M)) random_effect_sigma[m] = exp(random_effect_sigma[m]/3.0);
+  
+  // Non centered parameterisation SD of random effects 2
+  array[M * (ncol_X_random_eff[2]> 0)] vector[how_many_factors_in_random_design[2]] random_effect_sigma_2;
+  if(ncol_X_random_eff[2]> 0) for(m in 1:(M)) random_effect_sigma_2[m] = random_effect_sigma_mu[2] + random_effect_sigma_sigma[2] * random_effect_sigma_raw_2[m];
+  if(ncol_X_random_eff[2]> 0) for(m in 1:(M)) random_effect_sigma_2[m] = exp(random_effect_sigma_2[m]/3.0);
+    
+  matrix[ncol_X_random_eff[1] * (is_random_effect>0), M] random_effect;
+  matrix[ncol_X_random_eff[2] * (is_random_effect>0), M] random_effect_2;
+  
+  // random intercept
+  if(ncol_X_random_eff[1]> 0){
+    
+    // Convert vector array
+    array[ncol_X_random_eff[1]] vector[M] random_effect_raw_vec;
+    for(i in 1:ncol_X_random_eff[1]) {
+      random_effect_raw_vec[i] = random_effect_raw[i];
+    }
+    
+    // Covariate setup
+    random_effect =
+    get_random_effect_matrix(
+      M,
+      n_groups[1],
+      how_many_factors_in_random_design[1],
+      is_random_effect,
+      ncol_X_random_eff[1],
+      group_factor_indexes_for_covariance,
+      random_effect_raw_vec,
+      random_effect_sigma,
+      sigma_correlation_factor
+      );
+      
+  }
+  
+  // random intercept
+  if(ncol_X_random_eff[2]>0 ){
+
+    // Convert vector array
+    array[ncol_X_random_eff[2]] vector[M] random_effect_raw_2_vec;
+    for(i in 1:ncol_X_random_eff[2]) {
+      random_effect_raw_2_vec[i] = random_effect_raw_2[i];
+    }
+
+    // Covariate setup
+    random_effect_2 =
+    get_random_effect_matrix(
+      M,
+      n_groups[2],
+      how_many_factors_in_random_design[2],
+      is_random_effect,
+      ncol_X_random_eff[2],
+      group_factor_indexes_for_covariance_2,
+      random_effect_raw_2_vec,
+      random_effect_sigma_2,
+      sigma_correlation_factor_2
+      );
+
+  }
+}
+
 generated quantities{
   array[N_simulated, M_simulated] int counts_uncorrected;
   matrix[N_simulated, M_simulated] counts;
-  matrix[A_simulated,M_simulated] alpha_simulated;
-  matrix[M_simulated,N_simulated] mu = (X_simulated * beta_simulated)';
+  matrix[M_simulated,N_simulated] mu;
   matrix[M_simulated,N_simulated] precision;
 
-  matrix[A_simulated,M_simulated] beta_intercept_slope;
   // Vector of the generated exposure_simulated
   array[N_simulated] real generated_exposure;
 
-  // matrix[A_simulated,M_simulated] alpha_intercept_slope;
-
-  // All this because if A_simulated ==1 we have ocnversion problems
-  // This works only with two discrete groups
-  if(A_simulated == 1) beta_intercept_slope = to_matrix(beta_simulated[A_simulated,], A_simulated, M_simulated, 0);
-  else beta_intercept_slope = (XA_simulated * beta_simulated[1:A_simulated,]);
-  // if(A_simulated == 1)  alpha_intercept_slope = alpha;
-  // else alpha_intercept_slope = (XA_simulated * alpha);
+  // Use the transformed beta from transformed parameters
+  // Subset beta using X_which (matching replicate_data approach) - no padding needed
+  matrix[length_X_which, M_simulated] my_beta = beta[X_which, 1:M_simulated];
+  
+  // Compute alpha from beta using the association relationship (matching root model)
+  // First compute full alpha_simulated from beta, then subset using XA_which
+  matrix[A_simulated, M_simulated] alpha_simulated;
+  
+  // Use user-provided slope if available, otherwise use posterior prec_coeff[2]
+  real prec_coeff_slope = user_provided_prec_coeff_slope ? prec_coeff_slope_provided : prec_coeff[2];
+  // Use user-provided intercept if available, otherwise use posterior prec_coeff[1]
+  real prec_coeff_intercept = user_provided_prec_coeff_intercept ? prec_coeff_intercept_provided : prec_coeff[1];
+  
+  // Build alpha from beta following the root model's association
+  // Map A_simulated columns to corresponding beta columns via XA_which -> X_which
+  if(A_simulated == 1) {
+    // Simple case: single intercept column (variability ~ 1)
+    // Find intercept column in X_which
+    int beta_col = intercept_in_design && length_X_which > 0 ? X_which[1] : 1;
+    for(m in 1:M_simulated) {
+      alpha_simulated[1, m] = beta[beta_col, m] * prec_coeff_slope + prec_coeff_intercept;
+    }
+  } else {
+    // Multiple columns: handle intercept and non-intercept columns separately
+    int A_intercept_columns_sim = min(A_intercept_columns, A_simulated);
+    
+    // Intercept columns: alpha = beta * prec_coeff_slope + prec_coeff_intercept
+    for(a in 1:A_intercept_columns_sim) {
+      int alpha_col_idx = XA_which[a];
+      // Find corresponding beta column - intercept columns map to first columns in X_which
+      int beta_col = intercept_in_design && length_X_which > 0 ? X_which[min(a, length_X_which)] : 1;
+      for(m in 1:M_simulated) {
+        alpha_simulated[a, m] = beta[beta_col, m] * prec_coeff_slope + prec_coeff_intercept;
+      }
+    }
+    
+    // Non-intercept columns: alpha = beta * prec_coeff_slope
+    if(A_simulated > A_intercept_columns_sim) {
+      for(a in (A_intercept_columns_sim + 1):A_simulated) {
+        int alpha_col_idx = XA_which[a];
+        // Find corresponding beta column - map via XA_which to X_which
+        int beta_col_idx = a;
+        int beta_col = beta_col_idx <= length_X_which ? X_which[beta_col_idx] : X_which[length_X_which];
+        for(m in 1:M_simulated) {
+          alpha_simulated[a, m] = beta[beta_col, m] * prec_coeff_slope;
+        }
+      }
+    }
+  }
+  
+  // Subset alpha_simulated using XA_which to get my_alpha (matching replicate_data approach)
+  matrix[length_XA_which, M_simulated] my_alpha = alpha_simulated[XA_which,];
 
-  // PRECISION REGRESSION
-  for(a in 1:A_simulated) for(m in 1:M_simulated) alpha_simulated[a,m] = normal_rng( beta_intercept_slope[a,m] * prec_coeff[2] + prec_coeff[1], prec_sd);
+  // Calculate mu using the subsetted beta (matching replicate_data approach)
+  mu = (X_simulated * my_beta)';
+  
+  // Add random effects if present
+  // Issue 5: Clarify dimension handling - ncol_X_random_eff[1] is the actual number of random effect parameters
+  // The matrix X_random_effect_simulated has ncol_re1 columns (which is max(1, ncol_X_random_eff[1]))
+  // We only use the first ncol_X_random_eff[1] columns if ncol_X_random_eff[1] > 0
+  if(ncol_X_random_eff[1]> 0 && is_random_effect > 0) {
+    // Extract the relevant subset of random effects (ncol_X_random_eff[1] rows, M_simulated columns)
+    matrix[ncol_X_random_eff[1], M_simulated] random_effect_subset = random_effect[1:ncol_X_random_eff[1], 1:M_simulated];
+    // Use only the first ncol_X_random_eff[1] columns of the design matrix
+    mu = mu + (X_random_effect_simulated[,1:ncol_X_random_eff[1]] * random_effect_subset)';
+  }
+  
+  if(ncol_X_random_eff[2]>0 && is_random_effect > 0) {
+    // Extract the relevant subset of random effects (ncol_X_random_eff[2] rows, M_simulated columns)
+    matrix[ncol_X_random_eff[2], M_simulated] random_effect_2_subset = random_effect_2[1:ncol_X_random_eff[2], 1:M_simulated];
+    // Use only the first ncol_X_random_eff[2] columns of the design matrix
+    mu = mu + (X_random_effect_2_simulated[,1:ncol_X_random_eff[2]] * random_effect_2_subset)';
+  }
 
-  precision = (X_simulated[,1:A_simulated] * alpha_simulated)';
+  // Calculate precision using the actual alpha from posterior
+  // Issue 2: Use Xa_simulated (variability design matrix) with my_alpha (matching replicate_data approach)
+  precision = (Xa_simulated * my_alpha)';
 
-  // Precision adjustment
+  // Precision adjustment (variability multiplier)
   precision = precision - log(variability_multiplier);
 
+  // Convert to proportions
   for(i in 1:N_simulated) mu[,i] = softmax(mu[,i]);
-  for(i in 1:cols(mu)) {
+  
+  // Generate counts
+  for(i in 1:N_simulated) {
     counts_uncorrected[i,] = beta_binomial_rng(
       exposure_simulated[i],
       mu[,i] .* exp(precision[,i]),
diff --git a/man/sccomp_simulate.Rd b/man/sccomp_simulate.Rd
new file mode 100644
index 00000000..26b48504
--- /dev/null
+++ b/man/sccomp_simulate.Rd
@@ -0,0 +1,114 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/simulate_data.R
+\name{sccomp_simulate}
+\alias{sccomp_simulate}
+\title{sccomp_simulate}
+\usage{
+sccomp_simulate(
+  fit,
+  formula_composition,
+  formula_variability = NULL,
+  new_data = NULL,
+  coefficients = NULL,
+  mean_dispersion_slope = NULL,
+  mean_dispersion_intercept = NULL,
+  .sample = NULL,
+  .cell_group = NULL,
+  variability_multiplier = 1,
+  number_of_draws = 1,
+  mcmc_seed = sample_seed(),
+  cores = detectCores(),
+  sig_figs = 9,
+  cache_stan_model = sccomp_stan_models_cache_dir,
+  .coefficients = NULL
+)
+}
+\arguments{
+\item{fit}{The result of sccomp_estimate execution. This is used for sampling from real-data properties.}
+
+\item{formula_composition}{A formula. The formula describing the model for differential abundance, for example ~treatment}
+
+\item{formula_variability}{A formula. The formula describing the model for differential variability, for example ~treatment}
+
+\item{new_data}{A tibble including sample-specific columns (sample identifier and factor columns from formula). If coefficients are provided separately, cell_group column is optional. Otherwise, should include cell_group column.}
+
+\item{coefficients}{A data frame/tibble with cell-type specific coefficients. Must contain a column matching the cell_group column name, and columns matching the design matrix column names (e.g., "(Intercept)", "typeB"). If NULL, posterior beta_raw will be used.}
+
+\item{mean_dispersion_slope}{Optional numeric value for the slope parameter of the mean-dispersion association. If NULL, the slope from the fitted model (prec_coeff\link{2}) will be used.}
+
+\item{mean_dispersion_intercept}{Optional numeric value for the intercept parameter of the mean-dispersion association. If NULL, the intercept from the fitted model (prec_coeff\link{1}) will be used. This ensures that when only the slope is changed, the intercept remains the same.}
+
+\item{.sample}{A column name as symbol. The sample identifier}
+
+\item{.cell_group}{A column name as symbol. The cell_group identifier}
+
+\item{variability_multiplier}{A real scalar. This can be used for artificially increasing the variability of the simulation for benchmarking purposes.}
+
+\item{number_of_draws}{An integer. How may copies of the data you want to draw from the model joint posterior distribution.}
+
+\item{mcmc_seed}{An integer. Used for Markov-chain Monte Carlo reproducibility. By default a random number is sampled from 1 to 999999. This itself can be controlled by set.seed()}
+
+\item{cores}{Integer, the number of cores to be used for parallel calculations.}
+
+\item{sig_figs}{Number of significant figures to use for Stan model output. Default is 9.}
+
+\item{cache_stan_model}{A character string specifying the cache directory for compiled Stan models.
+The sccomp version will be automatically appended to ensure version isolation.
+Default is \code{sccomp_stan_models_cache_dir} which points to \verb{~/.sccomp_models}.}
+
+\item{.coefficients}{(Deprecated) The column names for coefficients in new_data, for example, c(b_0, b_1). Use the 'coefficients' parameter instead. This parameter is placed last for backward compatibility.}
+}
+\value{
+A tibble (\code{tbl}) with the following columns:
+\itemize{
+\item \strong{sample} - A character column representing the sample name.
+\item \strong{type} - A factor column representing the type of the sample.
+\item \strong{phenotype} - A factor column representing the phenotype in the data.
+\item \strong{count} - An integer column representing the original cell counts.
+\item \strong{cell_group} - A character column representing the cell group identifier.
+\item \strong{b_0} - A numeric column representing the first coefficient used for simulation.
+\item \strong{b_1} - A numeric column representing the second coefficient used for simulation.
+\item \strong{generated_proportions} - A numeric column representing the generated proportions from the simulation.
+\item \strong{generated_counts} - An integer column representing the generated cell counts from the simulation.
+\item \strong{replicate} - An integer column representing the replicate number for each draw from the posterior distribution.
+}
+}
+\description{
+This function simulates data from a fitted model.
+}
+\examples{
+
+# print("cmdstanr is needed to run this example.")
+# Note: Before running the example, ensure that the 'cmdstanr' package is installed:
+# install.packages("cmdstanr", repos = c("https://stan-dev.r-universe.dev/", getOption("repos")))
+
+# \donttest{
+#   if (instantiate::stan_cmdstan_exists()) {
+#     data("counts_obj")
+#     library(dplyr)
+
+#     estimate = sccomp_estimate(
+#       counts_obj,
+#       ~ type, ~1, "sample", "cell_group", "count",
+#       cores = 1
+#     )
+
+#     # Set coefficients for cell_groups. In this case all coefficients are 0 for simplicity.
+#     counts_obj = counts_obj |> mutate(b_0 = 0, b_1 = 0)
+
+#     # Simulate data
+#     sccomp_simulate(
+#       estimate,
+#       formula_composition = ~type,
+#       formula_variability = ~1,
+#       new_data = counts_obj,
+#       .sample = sample,
+#       .cell_group = cell_group,
+#       .coefficients = c("b_0", "b_1")
+#     )
+#   }
+# }
+}
+\references{
+S. Mangiola, A.J. Roth-Schulze, M. Trussart, E. Zozaya-Valdés, M. Ma, Z. Gao, A.F. Rubin, T.P. Speed, H. Shim, & A.T. Papenfuss, sccomp: Robust differential composition and variability analysis for single-cell data, Proc. Natl. Acad. Sci. U.S.A. 120 (33) e2203828120, https://doi.org/10.1073/pnas.2203828120 (2023).
+}
diff --git a/man/simulate_data.Rd b/man/simulate_data.Rd
index 237b6a65..feb1c6c2 100644
--- a/man/simulate_data.Rd
+++ b/man/simulate_data.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/simulate_data.R
 \name{simulate_data}
 \alias{simulate_data}
-\title{simulate_data}
+\title{DEPRECATED: simulate_data}
 \usage{
 simulate_data(
   .data,
@@ -12,7 +12,7 @@ simulate_data(
   .sample = NULL,
   .cell_group = NULL,
   .coefficients = NULL,
-  variability_multiplier = 5,
+  variability_multiplier = 1,
   number_of_draws = 1,
   mcmc_seed = sample_seed(),
   cores = detectCores(),
@@ -33,7 +33,7 @@ simulate_data(
 
 \item{.cell_group}{A column name as symbol. The cell_group identifier}
 
-\item{.coefficients}{The column names for coefficients, for example, c(b_0, b_1)}
+\item{.coefficients}{(Deprecated) The column names for coefficients in new_data, for example, c(b_0, b_1). Use the 'coefficients' parameter in sccomp_simulate() instead.}
 
 \item{variability_multiplier}{A real scalar. This can be used for artificially increasing the variability of the simulation for benchmarking purposes.}
 
@@ -49,49 +49,6 @@ simulate_data(
 The sccomp version will be automatically appended to ensure version isolation.
 Default is \code{sccomp_stan_models_cache_dir} which points to \verb{~/.sccomp_models}.}
 }
-\value{
-A tibble (\code{tbl}) with the following columns:
-\itemize{
-\item \strong{sample} - A character column representing the sample name.
-\item \strong{type} - A factor column representing the type of the sample.
-\item \strong{phenotype} - A factor column representing the phenotype in the data.
-\item \strong{count} - An integer column representing the original cell counts.
-\item \strong{cell_group} - A character column representing the cell group identifier.
-\item \strong{b_0} - A numeric column representing the first coefficient used for simulation.
-\item \strong{b_1} - A numeric column representing the second coefficient used for simulation.
-\item \strong{generated_proportions} - A numeric column representing the generated proportions from the simulation.
-\item \strong{generated_counts} - An integer column representing the generated cell counts from the simulation.
-\item \strong{replicate} - An integer column representing the replicate number for each draw from the posterior distribution.
-}
-}
 \description{
-This function simulates data from a fitted model.
-}
-\examples{
-
-# print("cmdstanr is needed to run this example.")
-# Note: Before running the example, ensure that the 'cmdstanr' package is installed:
-# install.packages("cmdstanr", repos = c("https://stan-dev.r-universe.dev/", getOption("repos")))
-
-# \donttest{
-#   if (instantiate::stan_cmdstan_exists()) {
-#     data("counts_obj")
-#     library(dplyr)
-
-#     estimate = sccomp_estimate(
-#       counts_obj,
-#       ~ type, ~1, "sample", "cell_group", "count",
-#       cores = 1
-#     )
-
-#     # Set coefficients for cell_groups. In this case all coefficients are 0 for simplicity.
-#     counts_obj = counts_obj |> mutate(b_0 = 0, b_1 = 0)
-
-#     # Simulate data
-#     simulate_data(counts_obj, estimate, ~type, ~1, sample, cell_group, c(b_0, b_1))
-#   }
-# }
-}
-\references{
-S. Mangiola, A.J. Roth-Schulze, M. Trussart, E. Zozaya-Valdés, M. Ma, Z. Gao, A.F. Rubin, T.P. Speed, H. Shim, & A.T. Papenfuss, sccomp: Robust differential composition and variability analysis for single-cell data, Proc. Natl. Acad. Sci. U.S.A. 120 (33) e2203828120, https://doi.org/10.1073/pnas.2203828120 (2023).
+This function is DEPRECATED. Please use \code{\link{sccomp_simulate}} instead.
 }
diff --git a/tests/testthat/test-cache_stan_model.R b/tests/testthat/test-cache_stan_model.R
index 5b846c53..af0ef9e8 100644
--- a/tests/testthat/test-cache_stan_model.R
+++ b/tests/testthat/test-cache_stan_model.R
@@ -5,13 +5,15 @@ test_that("cache_stan_model parameter works correctly", {
   # Test that cache_stan_model parameter is available in all relevant functions
   expect_true("cache_stan_model" %in% names(formals(sccomp_estimate)))
   expect_true("cache_stan_model" %in% names(formals(sccomp_replicate)))
-  expect_true("cache_stan_model" %in% names(formals(simulate_data)))
+  expect_true("cache_stan_model" %in% names(formals(sccomp_simulate)))
+  expect_true("cache_stan_model" %in% names(formals(simulate_data)))  # Deprecated but still has parameter
   expect_true("cache_stan_model" %in% names(formals(sccomp_remove_outliers)))
   
   # Test that cache_stan_model defaults to sccomp_stan_models_cache_dir
   expect_equal(formals(sccomp_estimate)$cache_stan_model, quote(sccomp_stan_models_cache_dir))
   expect_equal(formals(sccomp_replicate)$cache_stan_model, quote(sccomp_stan_models_cache_dir))
-  expect_equal(formals(simulate_data)$cache_stan_model, quote(sccomp_stan_models_cache_dir))
+  expect_equal(formals(sccomp_simulate)$cache_stan_model, quote(sccomp_stan_models_cache_dir))
+  expect_equal(formals(simulate_data)$cache_stan_model, quote(sccomp_stan_models_cache_dir))  # Deprecated but still has parameter
   expect_equal(formals(sccomp_remove_outliers)$cache_stan_model, quote(sccomp_stan_models_cache_dir))
   
   # Test that sccomp:::load_model function handles cache_stan_model correctly
diff --git a/tests/testthat/test-simulate_data.R b/tests/testthat/test-simulate_data.R
new file mode 100644
index 00000000..aad54aca
--- /dev/null
+++ b/tests/testthat/test-simulate_data.R
@@ -0,0 +1,407 @@
+# Unit tests for sccomp_simulate function
+# Testing simulation functionality with various model configurations
+
+library(testthat)
+library(sccomp)
+library(dplyr)
+
+# Helper function to skip tests if cmdstan is not available
+skip_cmdstan <- function() {
+  if (!instantiate::stan_cmdstan_exists()) {
+    skip("CmdStan not available")
+  }
+}
+
+test_that("sccomp_simulate works with simple model", {
+  skip_cmdstan()
+  
+  # Load test data
+  data("counts_obj")
+  
+  # Fit a simple model
+  estimate = sccomp_estimate(
+    counts_obj,
+    ~ type, ~1, "sample", "cell_group", "count",
+    cores = 1
+  )
+  
+  # Create coefficients table (all zeros for simplicity)
+  coeffs_table = counts_obj |>
+    distinct(cell_group) |>
+    mutate(`(Intercept)` = 0, `typecancer` = 0)
+  
+  # Create sample-specific new_data
+  new_data_samples = counts_obj |>
+    distinct(sample, type)
+  
+  # Simulate data
+  result = sccomp_simulate(
+    estimate,
+    ~type, 
+    ~1, 
+    new_data = new_data_samples,
+    coefficients = coeffs_table,
+    cores = 1
+  )
+  
+  # Check that result is a tibble
+  expect_s3_class(result, "tbl")
+  
+  # Check that result has expected columns
+  expect_true("sample" %in% colnames(result))
+  expect_true("cell_group" %in% colnames(result))
+  
+  # Check that result has data
+  expect_true(nrow(result) > 0)
+})
+
+test_that("sccomp_simulate works with variability formula", {
+  skip_cmdstan()
+  
+  # Load test data
+  data("counts_obj")
+  
+  # Fit a model with variability formula
+  estimate = sccomp_estimate(
+    counts_obj,
+    ~ type, ~type, "sample", "cell_group", "count",
+    cores = 1
+  )
+  
+  # Create coefficients table
+  coeffs_table = counts_obj |>
+    distinct(cell_group) |>
+    mutate(`(Intercept)` = 0, `typecancer` = 0)
+  
+  # Create sample-specific new_data
+  new_data_samples = counts_obj |>
+    distinct(sample, type)
+  
+  # Simulate data with variability formula
+  result = sccomp_simulate(
+    estimate,
+    ~type, 
+    ~type,  # Use variability formula
+    new_data = new_data_samples,
+    coefficients = coeffs_table,
+    cores = 1
+  )
+  
+  # Check that result is a tibble
+  expect_s3_class(result, "tbl")
+  
+  # Check that result has expected columns
+  expect_true("sample" %in% colnames(result))
+  expect_true("cell_group" %in% colnames(result))
+})
+
+test_that("sccomp_simulate validates dimensions correctly", {
+  skip_cmdstan()
+  
+  # Load test data
+  data("counts_obj")
+  
+  # Fit a simple model
+  estimate = sccomp_estimate(
+    counts_obj,
+    ~ type, ~1, "sample", "cell_group", "count",
+    cores = 1
+  )
+  
+  # Create coefficients table with more cell groups than original (should fail)
+  coeffs_table_expanded = counts_obj |>
+    bind_rows(
+      counts_obj |>
+        mutate(cell_group = paste0(cell_group, "_new"))
+    ) |>
+    distinct(cell_group) |>
+    mutate(`(Intercept)` = 0, `typecancer` = 0)
+  
+  # Create sample-specific new_data
+  new_data_samples = counts_obj |>
+    distinct(sample, type)
+  
+  # This should fail with dimension validation error
+  expect_error(
+    sccomp_simulate(
+      estimate,
+      ~type, 
+      ~1, 
+      new_data = new_data_samples,
+      coefficients = coeffs_table_expanded,
+      cores = 1
+    ),
+    "M_simulated.*cannot be larger than M"
+  )
+})
+
+test_that("sccomp_simulate works with number_of_draws parameter", {
+  skip_cmdstan()
+  
+  # Load test data
+  data("counts_obj")
+  
+  # Fit a simple model
+  estimate = sccomp_estimate(
+    counts_obj,
+    ~ type, ~1, "sample", "cell_group", "count",
+    cores = 1
+  )
+  
+  # Create coefficients table
+  coeffs_table = counts_obj |>
+    distinct(cell_group) |>
+    mutate(`(Intercept)` = 0, `typecancer` = 0)
+  
+  # Create sample-specific new_data
+  new_data_samples = counts_obj |>
+    distinct(sample, type)
+  
+  # Simulate with multiple draws
+  result = sccomp_simulate(
+    estimate,
+    ~type, 
+    ~1, 
+    new_data = new_data_samples,
+    coefficients = coeffs_table,
+    number_of_draws = 3,
+    cores = 1
+  )
+  
+  # Check that result is a tibble
+  expect_s3_class(result, "tbl")
+  
+  # Check that we have data (may have replicate column if multiple draws)
+  expect_true(nrow(result) > 0)
+})
+
+test_that("sccomp_simulate works with variability_multiplier", {
+  skip_cmdstan()
+  
+  # Load test data
+  data("counts_obj")
+  
+  # Fit a simple model
+  estimate = sccomp_estimate(
+    counts_obj,
+    ~ type, ~1, "sample", "cell_group", "count",
+    cores = 1
+  )
+  
+  # Create coefficients table
+  coeffs_table = counts_obj |>
+    distinct(cell_group) |>
+    mutate(`(Intercept)` = 0, `typecancer` = 0)
+  
+  # Create sample-specific new_data
+  new_data_samples = counts_obj |>
+    distinct(sample, type)
+  
+  # Simulate with different variability multipliers
+  result1 = sccomp_simulate(
+    estimate,
+    ~type, 
+    ~1, 
+    new_data = new_data_samples,
+    coefficients = coeffs_table,
+    variability_multiplier = 1,
+    cores = 1
+  )
+  
+  result2 = sccomp_simulate(
+    estimate,
+    ~type, 
+    ~1, 
+    new_data = new_data_samples,
+    coefficients = coeffs_table,
+    variability_multiplier = 10,
+    cores = 1
+  )
+  
+  # Both should work
+  expect_s3_class(result1, "tbl")
+  expect_s3_class(result2, "tbl")
+})
+
+test_that("sccomp_simulate handles missing formula_variability gracefully", {
+  skip_cmdstan()
+  
+  # Load test data
+  data("counts_obj")
+  
+  # Fit a model with variability formula
+  estimate = sccomp_estimate(
+    counts_obj,
+    ~ type, ~type, "sample", "cell_group", "count",
+    cores = 1
+  )
+  
+  # Create coefficients table
+  coeffs_table = counts_obj |>
+    distinct(cell_group) |>
+    mutate(`(Intercept)` = 0, `typecancer` = 0)
+  
+  # Create sample-specific new_data
+  new_data_samples = counts_obj |>
+    distinct(sample, type)
+  
+  # Simulate without specifying formula_variability (should use from estimate)
+  result = sccomp_simulate(
+    estimate,
+    ~type, 
+    NULL,  # Let it use the formula from estimate
+    new_data = new_data_samples,
+    coefficients = coeffs_table,
+    cores = 1
+  )
+  
+  # Should work
+  expect_s3_class(result, "tbl")
+})
+
+test_that("sccomp_simulate works with subset of original data", {
+  skip_cmdstan()
+  
+  # Load test data
+  data("counts_obj")
+  
+  # Fit a model
+  estimate = sccomp_estimate(
+    counts_obj,
+    ~ type, ~1, "sample", "cell_group", "count",
+    cores = 1
+  )
+  
+  # Create coefficients table
+  coeffs_table = counts_obj |>
+    distinct(cell_group) |>
+    mutate(`(Intercept)` = 0, `typecancer` = 0)
+  
+  # Use subset of data (fewer samples)
+  new_data_samples_subset = counts_obj |>
+    filter(sample %in% unique(sample)[1:3]) |>
+    distinct(sample, type)
+  
+  # Simulate data
+  result = sccomp_simulate(
+    estimate,
+    ~type, 
+    ~1, 
+    new_data = new_data_samples_subset,
+    coefficients = coeffs_table,
+    cores = 1
+  )
+  
+  # Should work with subset
+  expect_s3_class(result, "tbl")
+  expect_true(nrow(result) >= nrow(new_data_samples_subset))
+})
+
+test_that("sccomp_simulate works with separate coefficients table", {
+  skip_cmdstan()
+  
+  # Load test data
+  data("counts_obj")
+  
+  # Fit a model
+  estimate = sccomp_estimate(
+    counts_obj,
+    ~ type, ~1, "sample", "cell_group", "count",
+    cores = 1
+  )
+  
+  # Create separate coefficients table (cell-type specific)
+  coeffs_table = counts_obj |>
+    distinct(cell_group) |>
+    mutate(`(Intercept)` = 0, `typecancer` = 0)
+  
+  # Create sample-specific new_data (no cell_group needed)
+  new_data_samples = counts_obj |>
+    distinct(sample, type)
+  
+  # Simulate data with separate coefficients table
+  result = sccomp_simulate(
+    estimate,
+    ~type, 
+    ~1, 
+    new_data = new_data_samples,
+    coefficients = coeffs_table,
+    cores = 1
+  )
+  
+  # Should work
+  expect_s3_class(result, "tbl")
+  expect_true("sample" %in% colnames(result))
+  expect_true("cell_group" %in% colnames(result))
+  expect_true(nrow(result) > 0)
+  
+  # Check that all cell_groups from coefficients table are present
+  expect_true(all(coeffs_table$cell_group %in% result$cell_group))
+})
+
+test_that("sccomp_simulate works with separate coefficients table and subset data", {
+  skip_cmdstan()
+  
+  # Load test data
+  data("counts_obj")
+  
+  # Fit a model
+  estimate = sccomp_estimate(
+    counts_obj,
+    ~ type, ~1, "sample", "cell_group", "count",
+    cores = 1
+  )
+  
+  # Create separate coefficients table
+  coeffs_table = counts_obj |>
+    distinct(cell_group) |>
+    mutate(`(Intercept)` = 0, `typecancer` = 0)
+  
+  # Create subset of samples
+  new_data_samples = counts_obj |>
+    filter(sample %in% unique(sample)[1:3]) |>
+    distinct(sample, type)
+  
+  # Simulate data
+  result = sccomp_simulate(
+    estimate,
+    ~type, 
+    ~1, 
+    new_data = new_data_samples,
+    coefficients = coeffs_table,
+    cores = 1
+  )
+  
+  # Should work
+  expect_s3_class(result, "tbl")
+  expect_true(nrow(result) > 0)
+})
+
+test_that("sccomp_simulate uses posterior beta when coefficients are NULL", {
+  skip_cmdstan()
+  
+  # Load test data
+  data("counts_obj")
+  
+  # Fit a model
+  estimate = sccomp_estimate(
+    counts_obj,
+    ~ type, ~1, "sample", "cell_group", "count",
+    cores = 1
+  )
+  
+  # Simulate without providing coefficients (should use posterior beta_raw)
+  result = sccomp_simulate(
+    estimate,
+    ~type, 
+    ~1, 
+    new_data = NULL,  # Use original data
+    coefficients = NULL,  # No coefficients provided
+    cores = 1
+  )
+  
+  # Should work
+  expect_s3_class(result, "tbl")
+  expect_true(nrow(result) > 0)
+})
+