From 2051dfac0bc0f7d02854e6d34d23a9e608728bbb Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Fri, 9 Jan 2026 13:46:59 +1030 Subject: [PATCH 01/12] Enhance simulate_data function with variability handling and dimension validation - Added support for variability formulas in the simulate_data function, allowing users to specify a formula for variability alongside the main formula. - Implemented dimension validation to ensure that simulated data does not exceed the dimensions of the original model, preventing potential errors during simulation. - Introduced a new utility function to normalize sum-to-zero vector parameters in posterior draws, addressing floating-point precision issues. - Updated Stan model to accommodate new design matrices for random effects and variability, ensuring compatibility with existing functionality. - Added comprehensive unit tests to validate the new features and ensure robustness against various model configurations. --- DESCRIPTION | 2 +- R/simulate_data.R | 254 ++++++++++++++-- R/utilities.R | 57 +++- ...glm_multi_beta_binomial_simulate_data.stan | 220 ++++++++++++-- tests/testthat/test-simulate_data.R | 280 ++++++++++++++++++ 5 files changed, 747 insertions(+), 66 deletions(-) create mode 100644 tests/testthat/test-simulate_data.R diff --git a/DESCRIPTION b/DESCRIPTION index 49d52056..fe358fe4 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.23 +Version: 2.1.24 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/R/simulate_data.R b/R/simulate_data.R index c30c54e0..0a1d1195 100644 --- a/R/simulate_data.R +++ b/R/simulate_data.R @@ -122,13 +122,16 @@ simulate_data.tbl = function(.data, check_if_columns_right_class(.data, !!.sample, !!.cell_group) model_data = attr(.estimate_object, "model_input") + original_data = attr(.estimate_object, "model_input") - # # 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")) - + # Get formulas if not provided + if(is.null(formula_variability)) { + formula_variability = attr(.estimate_object, "formula_variability") + } + original_formula_composition = attr(.estimate_object, "formula_composition") + # Validate dimensions before proceeding + # Issue 3: Add dimension validation data_for_model = .data |> nest(data___ = -!!.sample) |> @@ -136,49 +139,206 @@ simulate_data.tbl = function(.data, unnest(data___) |> data_simulation_to_model_input( formula_composition, - #formula_variability, + formula_variability, # Issue 2: Use formula_variability !!.sample, !!.cell_group, .exposure, !!.coefficients ) 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" ))] + # 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] + } + } + } + } - # [1] 5.6260004 -0.6940178 - # prec_sd = 0.816423129 + # 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) + # 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 + )) + + # Get posterior draws - reuse pattern from replicate_data + number_of_draws_in_the_fit = attr(.estimate_object, "fit") |> get_output_samples() + number_of_draws = min(number_of_draws, number_of_draws_in_the_fit) + + draws_matrix = attr(.estimate_object, "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] + } + + # Normalize all sum_to_zero_vector parameters to ensure they sum to exactly zero + # This fixes floating-point precision issues when using generate_quantities + # Uses the utility function from utilities.R + draws_matrix = normalize_sum_to_zero_params(draws_matrix, "^beta_raw\\[") + draws_matrix = normalize_sum_to_zero_params(draws_matrix, "^random_effect_raw\\[") + draws_matrix = normalize_sum_to_zero_params(draws_matrix, "^random_effect_raw_2\\[") + + # 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(.estimate_object, "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_simulation_to_model_input) + cell_group_names = + .data |> + 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) + # Join with original data - reuse pattern from sccomp_predict .data |> left_join( parsed_fit, @@ -195,7 +355,7 @@ simulate_data.tbl = function(.data, #' @noRd #' data_simulation_to_model_input = - function(.data, formula, .sample, .cell_type, .exposure, .coefficients, truncation_ajustment = 1, approximate_posterior_inference ){ + 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 +369,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 +393,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 +420,34 @@ 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 + coeff_names = quo_names(.coefficients) + + # 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) 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 cd66d98c..28f89f0f 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/inst/stan/glm_multi_beta_binomial_simulate_data.stan b/inst/stan/glm_multi_beta_binomial_simulate_data.stan index a719643e..57d6554b 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,232 @@ 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 variability_multiplier; - 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 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 + // Must use sum_to_zero_vector to match the main model's parameter types for generate_quantities + // The precision issue is handled by normalizing in transformed parameters + array[C] sum_to_zero_vector[M] beta_raw; // Each row is a sum_to_zero_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 prec_sd; real 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 - must match main model type + array[ncol_X_random_eff[1] * (is_random_effect>0)] sum_to_zero_vector[M] random_effect_raw; + array[ncol_X_random_eff[2] * (ncol_X_random_eff[2]>0)] sum_to_zero_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 sum_to_zero_vector to regular matrix (matching main model) + matrix[C,M] beta; + 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 sum_to_zero_vector array to vector array for function call + 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] = to_vector(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 sum_to_zero_vector array to vector array for function call + 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] = to_vector(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 to match simulated data dimensions + // Assumes C_simulated <= C and M_simulated <= M (simulation typically uses same or subset of factors) + matrix[C_simulated, M_simulated] beta_simulated = beta[1:C_simulated, 1:M_simulated]; + + // Compute alpha from beta using the association relationship (matching root model) + // This matches how alpha is related to beta in the main model's priors + // The root model uses: alpha = beta * prec_coeff[2] + prec_coeff[1] for intercept columns + // and alpha = beta * prec_coeff[2] for non-intercept columns + matrix[A_simulated, M_simulated] alpha_simulated; + + // Build alpha from beta following the root model's association + // Note: A_simulated corresponds to variability design columns, which should match composition design columns + // We use beta_simulated indices that correspond to the variability design columns + if(A_simulated == 1) { + // Simple case: single intercept column (variability ~ 1) + // Use beta[1] which corresponds to the intercept in composition design + for(m in 1:M_simulated) { + alpha_simulated[1, m] = beta_simulated[1, m] * prec_coeff[2] + prec_coeff[1]; + } + } 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[2] + prec_coeff[1] + // Note: beta indices should match alpha indices (both use same design structure) + for(a in 1:A_intercept_columns_sim) { + // Ensure we don't go out of bounds for beta_simulated + int beta_idx = min(a, C_simulated); + for(m in 1:M_simulated) { + alpha_simulated[a, m] = beta_simulated[beta_idx, m] * prec_coeff[2] + prec_coeff[1]; + } + } + + // Non-intercept columns: alpha = beta * prec_coeff[2] + if(A_simulated > A_intercept_columns_sim) { + for(a in (A_intercept_columns_sim + 1):A_simulated) { + // Ensure we don't go out of bounds for beta_simulated + int beta_idx = min(a, C_simulated); + for(m in 1:M_simulated) { + alpha_simulated[a, m] = beta_simulated[beta_idx, m] * prec_coeff[2]; + } + } + } + } - // 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 posterior beta + mu = (X_simulated * beta_simulated)'; + + // 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) instead of subsetting X_simulated + precision = (Xa_simulated * alpha_simulated)'; - // 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/tests/testthat/test-simulate_data.R b/tests/testthat/test-simulate_data.R new file mode 100644 index 00000000..8cfd4a60 --- /dev/null +++ b/tests/testthat/test-simulate_data.R @@ -0,0 +1,280 @@ +# Unit tests for simulate_data 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("simulate_data 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 + ) + + # Set coefficients for cell_groups (all zeros for simplicity) + counts_obj_with_coefs = counts_obj |> + mutate(b_0 = 0, b_1 = 0) + + # Simulate data + result = simulate_data( + counts_obj_with_coefs, + estimate, + ~type, + ~1, + sample, + cell_group, + c(b_0, b_1), + 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 same number of rows as input (or more if multiple draws) + expect_true(nrow(result) >= nrow(counts_obj_with_coefs)) +}) + +test_that("simulate_data 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 + ) + + # Set coefficients + counts_obj_with_coefs = counts_obj |> + mutate(b_0 = 0, b_1 = 0) + + # Simulate data with variability formula + result = simulate_data( + counts_obj_with_coefs, + estimate, + ~type, + ~type, # Use variability formula + sample, + cell_group, + c(b_0, b_1), + 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("simulate_data 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 data with more cell groups than original (should fail) + counts_obj_expanded = counts_obj |> + bind_rows( + counts_obj |> + mutate(cell_group = paste0(cell_group, "_new")) + ) |> + mutate(b_0 = 0, b_1 = 0) + + # This should fail with dimension validation error + expect_error( + simulate_data( + counts_obj_expanded, + estimate, + ~type, + ~1, + sample, + cell_group, + c(b_0, b_1), + cores = 1 + ), + "M_simulated.*cannot be larger than M" + ) +}) + +test_that("simulate_data 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 + ) + + # Set coefficients + counts_obj_with_coefs = counts_obj |> + mutate(b_0 = 0, b_1 = 0) + + # Simulate with multiple draws + result = simulate_data( + counts_obj_with_coefs, + estimate, + ~type, + ~1, + sample, + cell_group, + c(b_0, b_1), + 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("simulate_data 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 + ) + + # Set coefficients + counts_obj_with_coefs = counts_obj |> + mutate(b_0 = 0, b_1 = 0) + + # Simulate with different variability multipliers + result1 = simulate_data( + counts_obj_with_coefs, + estimate, + ~type, + ~1, + sample, + cell_group, + c(b_0, b_1), + variability_multiplier = 1, + cores = 1 + ) + + result2 = simulate_data( + counts_obj_with_coefs, + estimate, + ~type, + ~1, + sample, + cell_group, + c(b_0, b_1), + variability_multiplier = 10, + cores = 1 + ) + + # Both should work + expect_s3_class(result1, "tbl") + expect_s3_class(result2, "tbl") +}) + +test_that("simulate_data 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 + ) + + # Set coefficients + counts_obj_with_coefs = counts_obj |> + mutate(b_0 = 0, b_1 = 0) + + # Simulate without specifying formula_variability (should use from estimate) + result = simulate_data( + counts_obj_with_coefs, + estimate, + ~type, + NULL, # Let it use the formula from estimate + sample, + cell_group, + c(b_0, b_1), + cores = 1 + ) + + # Should work + expect_s3_class(result, "tbl") +}) + +test_that("simulate_data 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 + ) + + # Use subset of data (fewer samples) + counts_obj_subset = counts_obj |> + filter(sample %in% unique(sample)[1:3]) |> + mutate(b_0 = 0, b_1 = 0) + + # Simulate data + result = simulate_data( + counts_obj_subset, + estimate, + ~type, + ~1, + sample, + cell_group, + c(b_0, b_1), + cores = 1 + ) + + # Should work with subset + expect_s3_class(result, "tbl") + expect_true(nrow(result) >= nrow(counts_obj_subset)) +}) + From d39569a584abb978f726d63956b7b5df7d780180 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Mon, 12 Jan 2026 11:38:01 +1030 Subject: [PATCH 02/12] New simulation function and deprecations - Introduced `sccomp_simulate` function for simulating data from fitted models, enhancing usability and flexibility. - Deprecated the `simulate_data` function, providing a warning to users to transition to `sccomp_simulate`. - Updated documentation to reflect changes and provide clear guidance on the new function's usage. - Enhanced unit tests to validate the new simulation functionality and ensure compatibility with existing workflows. - Adjusted Stan model to accommodate new parameters for user-provided coefficients, improving simulation accuracy and performance. --- DESCRIPTION | 2 +- NAMESPACE | 4 +- R/simulate_data.R | 442 +++++++++++++++--- ...glm_multi_beta_binomial_simulate_data.stan | 116 +++-- man/sccomp_simulate.Rd | 100 ++++ man/simulate_data.Rd | 49 +- tests/testthat/test-simulate_data.R | 285 +++++++---- 7 files changed, 781 insertions(+), 217 deletions(-) create mode 100644 man/sccomp_simulate.Rd diff --git a/DESCRIPTION b/DESCRIPTION index fe358fe4..aba87cc5 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.24 +Version: 2.1.25 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 0a1d1195..6db8beaa 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,19 @@ #' @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 .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 +28,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 +70,28 @@ #' # 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, .sample = NULL, .cell_group = NULL, - .coefficients = NULL, variability_multiplier = 5, 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,51 +103,264 @@ 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, .sample = NULL, .cell_group = NULL, - .coefficients = NULL, variability_multiplier = 5, 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) + } - #Check column class - check_if_columns_right_class(.data, !!.sample, !!.cell_group) + # Handle coefficients parameter + # If coefficients table is provided, use it; otherwise check for .coefficients column names (backward compatibility) + .coefficients_quo = enquo(.coefficients) - model_data = attr(.estimate_object, "model_input") - original_data = attr(.estimate_object, "model_input") + # 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." + ) + } + + # 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(.estimate_object, "formula_variability") + formula_variability = attr(fit, "formula_variability") } - original_formula_composition = attr(.estimate_object, "formula_composition") + original_formula_composition = attr(fit, "formula_composition") # Validate dimensions before proceeding - # Issue 3: Add dimension validation + # 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 + + # 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, # Issue 2: Use 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") + # 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) + } + + # 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.") @@ -268,6 +485,46 @@ simulate_data.tbl = function(.data, 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)) + } + + # 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. + # 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 + } + # 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 |> @@ -275,14 +532,20 @@ simulate_data.tbl = function(.data, c(list( variability_multiplier = variability_multiplier, X_random_effect_simulated = X_random_effect_simulated, - X_random_effect_2_simulated = X_random_effect_2_simulated + X_random_effect_2_simulated = X_random_effect_2_simulated, + user_provided_beta = user_provided_beta, + beta_simulated_provided = beta_simulated_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(.estimate_object, "fit") |> get_output_samples() + 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(.estimate_object, "fit")$draws(format = "matrix") + draws_matrix = attr(fit, "fit")$draws(format = "matrix") if(number_of_draws > nrow(draws_matrix)) { number_of_draws = nrow(draws_matrix) @@ -293,12 +556,8 @@ simulate_data.tbl = function(.data, draws_matrix = draws_matrix[sample(seq_len(nrow(draws_matrix)), size = number_of_draws),, drop = FALSE] } - # Normalize all sum_to_zero_vector parameters to ensure they sum to exactly zero - # This fixes floating-point precision issues when using generate_quantities - # Uses the utility function from utilities.R - draws_matrix = normalize_sum_to_zero_params(draws_matrix, "^beta_raw\\[") - draws_matrix = normalize_sum_to_zero_params(draws_matrix, "^random_effect_raw\\[") - draws_matrix = normalize_sum_to_zero_params(draws_matrix, "^random_effect_raw_2\\[") + # 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( @@ -306,15 +565,16 @@ simulate_data.tbl = function(.data, 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_simulation_to_model_input) + # 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 |> + .data_expanded |> distinct(!!.cell_group) |> arrange(!!.cell_group) |> pull(!!.cell_group) @@ -339,7 +599,8 @@ simulate_data.tbl = function(.data, select(-N, -M) # Join with original data - reuse pattern from sccomp_predict - .data |> + # 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)) @@ -348,13 +609,72 @@ 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 = 5, + 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.2.0", + "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, + 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 = +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 @@ -423,24 +743,34 @@ data_simulation_to_model_input = # 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 - coeff_names = quo_names(.coefficients) - - # Pivot to long format: cell_type | coefficient_name | value - coefficients_long = + # 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, all_of(coeff_names)) %>% + select(!!.cell_type, all_of(coeff_names)) %>% distinct() %>% arrange(!!.cell_type) %>% - 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) + 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(), diff --git a/inst/stan/glm_multi_beta_binomial_simulate_data.stan b/inst/stan/glm_multi_beta_binomial_simulate_data.stan index 57d6554b..a250e746 100755 --- a/inst/stan/glm_multi_beta_binomial_simulate_data.stan +++ b/inst/stan/glm_multi_beta_binomial_simulate_data.stan @@ -13,6 +13,15 @@ data{ matrix[A_simulated, A_simulated] XA_simulated; // Unique variability design (for compatibility with old code) real variability_multiplier; + + // Optional: provided coefficients to override posterior beta_raw + int 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 int M; int C; @@ -20,6 +29,12 @@ data{ 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; @@ -45,9 +60,12 @@ transformed data{ } parameters{ - // Must use sum_to_zero_vector to match the main model's parameter types for generate_quantities - // The precision issue is handled by normalizing in transformed parameters - array[C] sum_to_zero_vector[M] beta_raw; // Each row is a sum_to_zero_vector of length M + // 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 @@ -55,9 +73,9 @@ parameters{ real prec_sd; real mix_p; - // Random intercept - must match main model type - array[ncol_X_random_eff[1] * (is_random_effect>0)] sum_to_zero_vector[M] random_effect_raw; - array[ncol_X_random_eff[2] * (ncol_X_random_eff[2]>0)] sum_to_zero_vector[M] 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[2 * (is_random_effect>0)] real random_effect_sigma_mu; @@ -77,10 +95,38 @@ parameters{ } transformed parameters{ - // Convert sum_to_zero_vector to regular matrix (matching main model) + // 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; - for(c in 1:C) { - beta[c,] = to_row_vector(beta_raw[c]); + 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) @@ -99,10 +145,10 @@ transformed parameters{ // random intercept if(ncol_X_random_eff[1]> 0){ - // Convert sum_to_zero_vector array to vector array for function call + // 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] = to_vector(random_effect_raw[i]); + random_effect_raw_vec[i] = random_effect_raw[i]; } // Covariate setup @@ -124,10 +170,10 @@ transformed parameters{ // random intercept if(ncol_X_random_eff[2]>0 ){ - // Convert sum_to_zero_vector array to vector array for function call + // 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] = to_vector(random_effect_raw_2[i]); + random_effect_raw_2_vec[i] = random_effect_raw_2[i]; } // Covariate setup @@ -157,53 +203,55 @@ generated quantities{ array[N_simulated] real generated_exposure; // Use the transformed beta from transformed parameters - // Subset beta to match simulated data dimensions - // Assumes C_simulated <= C and M_simulated <= M (simulation typically uses same or subset of factors) - matrix[C_simulated, M_simulated] beta_simulated = beta[1:C_simulated, 1:M_simulated]; + // 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) - // This matches how alpha is related to beta in the main model's priors - // The root model uses: alpha = beta * prec_coeff[2] + prec_coeff[1] for intercept columns - // and alpha = beta * prec_coeff[2] for non-intercept columns + // First compute full alpha_simulated from beta, then subset using XA_which matrix[A_simulated, M_simulated] alpha_simulated; // Build alpha from beta following the root model's association - // Note: A_simulated corresponds to variability design columns, which should match composition design columns - // We use beta_simulated indices that correspond to the variability design columns + // Map A_simulated columns to corresponding beta columns via XA_which -> X_which if(A_simulated == 1) { // Simple case: single intercept column (variability ~ 1) - // Use beta[1] which corresponds to the intercept in composition design + // 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_simulated[1, m] * prec_coeff[2] + prec_coeff[1]; + alpha_simulated[1, m] = beta[beta_col, m] * prec_coeff[2] + prec_coeff[1]; } } 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[2] + prec_coeff[1] - // Note: beta indices should match alpha indices (both use same design structure) for(a in 1:A_intercept_columns_sim) { - // Ensure we don't go out of bounds for beta_simulated - int beta_idx = min(a, C_simulated); + 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_simulated[beta_idx, m] * prec_coeff[2] + prec_coeff[1]; + alpha_simulated[a, m] = beta[beta_col, m] * prec_coeff[2] + prec_coeff[1]; } } // Non-intercept columns: alpha = beta * prec_coeff[2] if(A_simulated > A_intercept_columns_sim) { for(a in (A_intercept_columns_sim + 1):A_simulated) { - // Ensure we don't go out of bounds for beta_simulated - int beta_idx = min(a, C_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_simulated[beta_idx, m] * prec_coeff[2]; + alpha_simulated[a, m] = beta[beta_col, m] * prec_coeff[2]; } } } } + + // 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,]; - // Calculate mu using the posterior beta - mu = (X_simulated * beta_simulated)'; + // 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 @@ -224,8 +272,8 @@ generated quantities{ } // Calculate precision using the actual alpha from posterior - // Issue 2: Use Xa_simulated (variability design matrix) instead of subsetting X_simulated - precision = (Xa_simulated * alpha_simulated)'; + // Issue 2: Use Xa_simulated (variability design matrix) with my_alpha (matching replicate_data approach) + precision = (Xa_simulated * my_alpha)'; // Precision adjustment (variability multiplier) precision = precision - log(variability_multiplier); diff --git a/man/sccomp_simulate.Rd b/man/sccomp_simulate.Rd new file mode 100644 index 00000000..dad29915 --- /dev/null +++ b/man/sccomp_simulate.Rd @@ -0,0 +1,100 @@ +% 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, + .sample = NULL, + .cell_group = NULL, + variability_multiplier = 5, + 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{.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, ~type, ~1, counts_obj, 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). +} diff --git a/man/simulate_data.Rd b/man/simulate_data.Rd index 237b6a65..2ad45a13 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, @@ -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-simulate_data.R b/tests/testthat/test-simulate_data.R index 8cfd4a60..aad54aca 100644 --- a/tests/testthat/test-simulate_data.R +++ b/tests/testthat/test-simulate_data.R @@ -1,4 +1,4 @@ -# Unit tests for simulate_data function +# Unit tests for sccomp_simulate function # Testing simulation functionality with various model configurations library(testthat) @@ -12,7 +12,7 @@ skip_cmdstan <- function() { } } -test_that("simulate_data works with simple model", { +test_that("sccomp_simulate works with simple model", { skip_cmdstan() # Load test data @@ -25,19 +25,22 @@ test_that("simulate_data works with simple model", { cores = 1 ) - # Set coefficients for cell_groups (all zeros for simplicity) - counts_obj_with_coefs = counts_obj |> - mutate(b_0 = 0, b_1 = 0) + # 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 = simulate_data( - counts_obj_with_coefs, - estimate, + result = sccomp_simulate( + estimate, ~type, ~1, - sample, - cell_group, - c(b_0, b_1), + new_data = new_data_samples, + coefficients = coeffs_table, cores = 1 ) @@ -48,11 +51,11 @@ test_that("simulate_data works with simple model", { expect_true("sample" %in% colnames(result)) expect_true("cell_group" %in% colnames(result)) - # Check that result has same number of rows as input (or more if multiple draws) - expect_true(nrow(result) >= nrow(counts_obj_with_coefs)) + # Check that result has data + expect_true(nrow(result) > 0) }) -test_that("simulate_data works with variability formula", { +test_that("sccomp_simulate works with variability formula", { skip_cmdstan() # Load test data @@ -65,19 +68,22 @@ test_that("simulate_data works with variability formula", { cores = 1 ) - # Set coefficients - counts_obj_with_coefs = counts_obj |> - mutate(b_0 = 0, b_1 = 0) + # 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 = simulate_data( - counts_obj_with_coefs, - estimate, + result = sccomp_simulate( + estimate, ~type, ~type, # Use variability formula - sample, - cell_group, - c(b_0, b_1), + new_data = new_data_samples, + coefficients = coeffs_table, cores = 1 ) @@ -89,7 +95,7 @@ test_that("simulate_data works with variability formula", { expect_true("cell_group" %in% colnames(result)) }) -test_that("simulate_data validates dimensions correctly", { +test_that("sccomp_simulate validates dimensions correctly", { skip_cmdstan() # Load test data @@ -102,31 +108,34 @@ test_that("simulate_data validates dimensions correctly", { cores = 1 ) - # Create data with more cell groups than original (should fail) - counts_obj_expanded = counts_obj |> + # 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")) ) |> - mutate(b_0 = 0, b_1 = 0) + 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( - simulate_data( - counts_obj_expanded, - estimate, + sccomp_simulate( + estimate, ~type, ~1, - sample, - cell_group, - c(b_0, b_1), + new_data = new_data_samples, + coefficients = coeffs_table_expanded, cores = 1 ), "M_simulated.*cannot be larger than M" ) }) -test_that("simulate_data works with number_of_draws parameter", { +test_that("sccomp_simulate works with number_of_draws parameter", { skip_cmdstan() # Load test data @@ -139,19 +148,22 @@ test_that("simulate_data works with number_of_draws parameter", { cores = 1 ) - # Set coefficients - counts_obj_with_coefs = counts_obj |> - mutate(b_0 = 0, b_1 = 0) + # 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 = simulate_data( - counts_obj_with_coefs, - estimate, + result = sccomp_simulate( + estimate, ~type, ~1, - sample, - cell_group, - c(b_0, b_1), + new_data = new_data_samples, + coefficients = coeffs_table, number_of_draws = 3, cores = 1 ) @@ -163,7 +175,7 @@ test_that("simulate_data works with number_of_draws parameter", { expect_true(nrow(result) > 0) }) -test_that("simulate_data works with variability_multiplier", { +test_that("sccomp_simulate works with variability_multiplier", { skip_cmdstan() # Load test data @@ -176,31 +188,32 @@ test_that("simulate_data works with variability_multiplier", { cores = 1 ) - # Set coefficients - counts_obj_with_coefs = counts_obj |> - mutate(b_0 = 0, b_1 = 0) + # 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 = simulate_data( - counts_obj_with_coefs, - estimate, + result1 = sccomp_simulate( + estimate, ~type, ~1, - sample, - cell_group, - c(b_0, b_1), + new_data = new_data_samples, + coefficients = coeffs_table, variability_multiplier = 1, cores = 1 ) - result2 = simulate_data( - counts_obj_with_coefs, - estimate, + result2 = sccomp_simulate( + estimate, ~type, ~1, - sample, - cell_group, - c(b_0, b_1), + new_data = new_data_samples, + coefficients = coeffs_table, variability_multiplier = 10, cores = 1 ) @@ -210,7 +223,7 @@ test_that("simulate_data works with variability_multiplier", { expect_s3_class(result2, "tbl") }) -test_that("simulate_data handles missing formula_variability gracefully", { +test_that("sccomp_simulate handles missing formula_variability gracefully", { skip_cmdstan() # Load test data @@ -223,19 +236,22 @@ test_that("simulate_data handles missing formula_variability gracefully", { cores = 1 ) - # Set coefficients - counts_obj_with_coefs = counts_obj |> - mutate(b_0 = 0, b_1 = 0) + # 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 = simulate_data( - counts_obj_with_coefs, - estimate, + result = sccomp_simulate( + estimate, ~type, NULL, # Let it use the formula from estimate - sample, - cell_group, - c(b_0, b_1), + new_data = new_data_samples, + coefficients = coeffs_table, cores = 1 ) @@ -243,7 +259,7 @@ test_that("simulate_data handles missing formula_variability gracefully", { expect_s3_class(result, "tbl") }) -test_that("simulate_data works with subset of original data", { +test_that("sccomp_simulate works with subset of original data", { skip_cmdstan() # Load test data @@ -256,25 +272,136 @@ test_that("simulate_data works with subset of original data", { cores = 1 ) + # Create coefficients table + coeffs_table = counts_obj |> + distinct(cell_group) |> + mutate(`(Intercept)` = 0, `typecancer` = 0) + # Use subset of data (fewer samples) - counts_obj_subset = counts_obj |> + new_data_samples_subset = counts_obj |> filter(sample %in% unique(sample)[1:3]) |> - mutate(b_0 = 0, b_1 = 0) + distinct(sample, type) # Simulate data - result = simulate_data( - counts_obj_subset, - estimate, + result = sccomp_simulate( + estimate, ~type, ~1, - sample, - cell_group, - c(b_0, b_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(counts_obj_subset)) + 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) }) From 769d220c0c0e676a49d3414f4211d80ffeaf9690 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Mon, 12 Jan 2026 17:54:45 +1030 Subject: [PATCH 03/12] Enhance coefficient validation in sccomp_simulate function - Added validation to ensure user-provided coefficients sum to zero within a specified tolerance, crucial for compositional models. - Implemented checks using tidyverse to identify and report invalid coefficients, improving error handling and user feedback. - Updated comments to clarify the importance of coefficient normalization in the context of simulation. --- R/simulate_data.R | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/R/simulate_data.R b/R/simulate_data.R index 6db8beaa..548f6048 100644 --- a/R/simulate_data.R +++ b/R/simulate_data.R @@ -511,10 +511,38 @@ sccomp_simulate.sccomp_tbl = function(fit, 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 From 255a2dd421658cd3479706a374e054107184ae27 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Tue, 13 Jan 2026 11:29:12 +1030 Subject: [PATCH 04/12] Enhance sccomp_simulate function with mean-dispersion parameters - Added optional parameters `mean_dispersion_slope` and `mean_dispersion_intercept` to allow users to specify the slope and intercept for the mean-dispersion association, enhancing flexibility in simulations. - Updated the `simulate_data` function to reflect changes in parameter defaults, setting `variability_multiplier` to 1 for consistency. - Enhanced documentation to include new parameters and their usage, ensuring clarity for users. - Adjusted the Stan model to accommodate the new mean-dispersion parameters, improving simulation accuracy. --- R/simulate_data.R | 48 +- dev/sccomp_simulate.qmd | 859 ++++++++++++++++++ ...glm_multi_beta_binomial_simulate_data.stan | 21 +- man/sccomp_simulate.Rd | 8 +- man/simulate_data.Rd | 2 +- 5 files changed, 928 insertions(+), 10 deletions(-) create mode 100644 dev/sccomp_simulate.qmd diff --git a/R/simulate_data.R b/R/simulate_data.R index 548f6048..1bf392f1 100644 --- a/R/simulate_data.R +++ b/R/simulate_data.R @@ -18,6 +18,8 @@ #' @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 variability_multiplier A real scalar. This can be used for artificially increasing the variability of the simulation for benchmarking purposes. @@ -78,9 +80,11 @@ sccomp_simulate <- function(fit, formula_variability = NULL, new_data = NULL, coefficients = NULL, + mean_dispersion_slope = NULL, + mean_dispersion_intercept = NULL, .sample = NULL, .cell_group = NULL, - variability_multiplier = 5, + variability_multiplier = 1, number_of_draws = 1, mcmc_seed = sample_seed(), cores = detectCores(), @@ -108,9 +112,11 @@ sccomp_simulate.sccomp_tbl = function(fit, formula_variability = NULL, new_data = NULL, coefficients = NULL, + mean_dispersion_slope = NULL, + mean_dispersion_intercept = NULL, .sample = NULL, .cell_group = NULL, - variability_multiplier = 5, + variability_multiplier = 1, number_of_draws = 1, mcmc_seed = sample_seed(), cores = detectCores(), @@ -553,6 +559,36 @@ sccomp_simulate.sccomp_tbl = function(fit, 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 |> @@ -563,6 +599,10 @@ sccomp_simulate.sccomp_tbl = function(fit, 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, @@ -666,7 +706,7 @@ simulate_data <- function(.data, .sample = NULL, .cell_group = NULL, .coefficients = NULL, - variability_multiplier = 5, + variability_multiplier = 1, number_of_draws = 1, mcmc_seed = sample_seed(), cores = detectCores(), @@ -687,6 +727,8 @@ simulate_data <- function(.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, diff --git a/dev/sccomp_simulate.qmd b/dev/sccomp_simulate.qmd new file mode 100644 index 00000000..4f7ce892 --- /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 3 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/stan/glm_multi_beta_binomial_simulate_data.stan b/inst/stan/glm_multi_beta_binomial_simulate_data.stan index a250e746..15590a23 100755 --- a/inst/stan/glm_multi_beta_binomial_simulate_data.stan +++ b/inst/stan/glm_multi_beta_binomial_simulate_data.stan @@ -22,6 +22,12 @@ data{ // 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 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 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 C; @@ -210,6 +216,11 @@ generated quantities{ // 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) { @@ -217,23 +228,23 @@ generated quantities{ // 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[2] + prec_coeff[1]; + 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[2] + prec_coeff[1] + // 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[2] + prec_coeff[1]; + alpha_simulated[a, m] = beta[beta_col, m] * prec_coeff_slope + prec_coeff_intercept; } } - // Non-intercept columns: alpha = beta * prec_coeff[2] + // 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]; @@ -241,7 +252,7 @@ generated quantities{ 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[2]; + alpha_simulated[a, m] = beta[beta_col, m] * prec_coeff_slope; } } } diff --git a/man/sccomp_simulate.Rd b/man/sccomp_simulate.Rd index dad29915..5f8f9a70 100644 --- a/man/sccomp_simulate.Rd +++ b/man/sccomp_simulate.Rd @@ -10,9 +10,11 @@ sccomp_simulate( formula_variability = NULL, new_data = NULL, coefficients = NULL, + mean_dispersion_slope = NULL, + mean_dispersion_intercept = NULL, .sample = NULL, .cell_group = NULL, - variability_multiplier = 5, + variability_multiplier = 1, number_of_draws = 1, mcmc_seed = sample_seed(), cores = detectCores(), @@ -32,6 +34,10 @@ sccomp_simulate( \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} diff --git a/man/simulate_data.Rd b/man/simulate_data.Rd index 2ad45a13..feb1c6c2 100644 --- a/man/simulate_data.Rd +++ b/man/simulate_data.Rd @@ -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(), From cecf2efdcf8a5cd9f5b238320d9223c13d65cecf Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Tue, 13 Jan 2026 11:30:56 +1030 Subject: [PATCH 05/12] Update NEWS for version 2.1.26 with simulation improvements - Changed default `variability_multiplier` from 5 to 1 in `sccomp_simulate()` and `simulate_data()`, preserving model variability by default. - Enhanced simulation documentation for clarity on `mean_dispersion_slope` effects and updated vignette examples for consistency. - Clarified behavior of mean-dispersion association parameters in simulation, detailing their impact on intercept and factor-level variability. --- inst/NEWS.rd | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/inst/NEWS.rd b/inst/NEWS.rd index 51992142..bcbec1b8 100644 --- a/inst/NEWS.rd +++ b/inst/NEWS.rd @@ -1,6 +1,13 @@ \name{NEWS} \title{News for Package \pkg{sccomp}} +\section{News in version 2.1.26}{ +\itemize{ + \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.22}{ \itemize{ \item Added automatic cleanup of Stan draw CSV files. New \code{cleanup_draw_files} parameter (default TRUE) in \code{sccomp_estimate()} and \code{sccomp_remove_outliers()} automatically removes large draw files after analysis completion, significantly reducing disk space usage. From 585977f875f245a5b4913cc2290e38e136e80047 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Tue, 13 Jan 2026 11:33:12 +1030 Subject: [PATCH 06/12] set version as 2.1.24 --- DESCRIPTION | 2 +- inst/NEWS.rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index aba87cc5..fe358fe4 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.25 +Version: 2.1.24 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/inst/NEWS.rd b/inst/NEWS.rd index bcbec1b8..37378241 100644 --- a/inst/NEWS.rd +++ b/inst/NEWS.rd @@ -1,7 +1,7 @@ \name{NEWS} \title{News for Package \pkg{sccomp}} -\section{News in version 2.1.26}{ +\section{News in version 2.1.24}{ \itemize{ \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. From 1491686c2f1342e7b03ce2e222f60bd8fc0fd8a2 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Tue, 13 Jan 2026 11:39:09 +1030 Subject: [PATCH 07/12] Update R/simulate_data.R Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- R/simulate_data.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/simulate_data.R b/R/simulate_data.R index 1bf392f1..8128084a 100644 --- a/R/simulate_data.R +++ b/R/simulate_data.R @@ -714,7 +714,7 @@ simulate_data <- function(.data, cache_stan_model = sccomp_stan_models_cache_dir) { lifecycle::deprecate_warn( - "2.2.0", + "2.1.26", "sccomp::simulate_data()", details = "sccomp says: simulate_data is deprecated. Please use sccomp_simulate() instead." ) From 14e4501cabd81447afc6d15bb1709ff916e3f8fb Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Tue, 13 Jan 2026 11:39:36 +1030 Subject: [PATCH 08/12] Update DESCRIPTION Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index fe358fe4..09b65fb3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: sccomp Type: Package Title: Differential Composition and Variability Analysis for Single-Cell Data Version: 2.1.24 -Date: 2024-01-15 +Date: 2026-01-13 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. License: GPL-3 From 6dc7b2e82b6ac1ed079bdbe452f18718824b3ead Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Tue, 13 Jan 2026 11:40:14 +1030 Subject: [PATCH 09/12] Update man/sccomp_simulate.Rd Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- man/sccomp_simulate.Rd | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/man/sccomp_simulate.Rd b/man/sccomp_simulate.Rd index 5f8f9a70..1488ca14 100644 --- a/man/sccomp_simulate.Rd +++ b/man/sccomp_simulate.Rd @@ -97,7 +97,15 @@ This function simulates data from a fitted model. # counts_obj = counts_obj |> mutate(b_0 = 0, b_1 = 0) # # Simulate data -# sccomp_simulate(estimate, ~type, ~1, counts_obj, sample, cell_group, c(b_0, b_1)) +# 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) +# ) # } # } } From a2660be62cdbf69e68505692fe08cb611be8f99a Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Tue, 13 Jan 2026 11:40:28 +1030 Subject: [PATCH 10/12] Update dev/sccomp_simulate.qmd Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- dev/sccomp_simulate.qmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dev/sccomp_simulate.qmd b/dev/sccomp_simulate.qmd index 4f7ce892..aeeb9182 100644 --- a/dev/sccomp_simulate.qmd +++ b/dev/sccomp_simulate.qmd @@ -69,7 +69,7 @@ simulated_data = sccomp_simulate( new_data = NULL, # Use original data coefficients = NULL, # Use posterior beta_raw - number_of_draws = 10 # Generate 3 replicates + number_of_draws = 10 # Generate 10 replicates ) # View the simulated data From 01e76261188fce696403fa75545a61a4bb0c4d08 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Tue, 20 Jan 2026 11:33:18 +1030 Subject: [PATCH 11/12] Update tests for cache_stan_model parameter in sccomp functions - Updated unit tests to check for the presence of the `cache_stan_model` parameter in `sccomp_simulate`, reflecting the recent function introduction. - Maintained checks for the deprecated `simulate_data` function to ensure backward compatibility. - Verified that the default value of `cache_stan_model` is set correctly across relevant functions. --- tests/testthat/test-cache_stan_model.R | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) 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 From ff3a83addf6151c8007f85acd9d2f0c4f4629786 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Wed, 18 Mar 2026 19:29:51 +1030 Subject: [PATCH 12/12] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- man/sccomp_simulate.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/sccomp_simulate.Rd b/man/sccomp_simulate.Rd index 1488ca14..26b48504 100644 --- a/man/sccomp_simulate.Rd +++ b/man/sccomp_simulate.Rd @@ -104,7 +104,7 @@ This function simulates data from a fitted model. # new_data = counts_obj, # .sample = sample, # .cell_group = cell_group, -# coefficients = c(b_0, b_1) +# .coefficients = c("b_0", "b_1") # ) # } # }