susieR: Sum of Single Effects (SuSiE) Model

The architecture is built around three primary object types:

Data Object

• Purpose: Holds processed, algorithm-ready input data.
• S3 Classes: individual, ss, or rss_lambda (this determines which backend methods are used).
• Properties: Immutable. Contains input matrices (e.g., X/y, XtX/Xty/yty, or z/R), metadata (n, p), and scaling attributes.

Params Object

• Purpose: Holds all algorithm parameters and user settings.
• Properties: Immutable. Contains settings for L (number of effects), max_iter, tol, convergence_method, estimate_prior_method, estimate_residual_method, and model options like unmappable_effects, refine, standardize, and intercept.
• Note: For rss_lambda data, certain parameters are fixed (e.g., unmappable_effects = "none", standardize = FALSE).

Model Object

• Purpose: Holds the fitted state, results, and algorithm outputs.
• Properties: Mutable. Updated during the fitting process. Contains model matrices (alpha, mu, mu2, V, sigma2), fitted values, algorithm outputs (ELBO, niter, converged), and final results (credible sets, PIPs, intercept, z-scores).

susieR 2.0 uses a unified architecture to support multiple data types (individual-level, sufficient statistics, and regression summary statistics) through a single algorithmic pipeline. It relies on S3 method dispatch to route different data representations to the appropriate backend computations.

The workflow follows this pattern:

1. Interface: You call high-level functions like susie(), susie_ss(), or susie_rss().
2. Constructor: These functions call constructors that validate your inputs and create immutable Data and Params objects.
3. Workhorse: The validated objects are passed to a workhorse function that orchestrates the main SuSiE algorithm.
4. IBSS Core: The core Iterative Bayesian Stepwise Selection (IBSS) algorithm executes the fitting process.
5. Backend Methods: S3 dispatch ensures the correct mathematical operations are performed based on whether the data is individual, ss, or rss_lambda type.

When individual-level genotype data is unavailable, susie_rss() can use a stochastic LD sketch to approximate the LD matrix. This is useful when you only have summary statistics or want to avoid the computational cost of the full LD matrix.

To use a stochastic sketch:

1. Generate a sketch matrix $U = X_{\text{std}}^ op W$, where $W \sim N(0, I_B / n)$ and $B$ is the sketch size.
2. Pass the transpose of this sketch, X = t(U), to the susie_rss() function. The function treats this $B \times p$ matrix as a genotype-like input that approximates the LD structure.

Important: Variance Inflation Correction Using a noisy sketch $\hat{R}$ can inflate the variance of the Single-Effect Regression (SER) residuals, leading to an elevated False Discovery Rate (FDR). To correct this, set the stochastic_ld_sample parameter to the sketch size $B$. This activates a dynamic variance inflation correction that accounts for the noise in the sketch at every iteration.

# Example: Using a stochastic sketch with correction
# Assuming U is your pre-computed sketch matrix (p x B)
# and B is the sketch size

susie_rss(
  z = z_scores,
  n = n_samples,
  L = L_signals,
  X = t(U),                # Pass the sketch as X
  stochastic_ld_sample = B  # Activate the dynamic correction
)

To execute the benchmark comparing standard SuSiE (Gaussian) against SuSiE-SS (Servin-Stephens), use jupyter nbconvert from the inst/notebooks/ directory. The simulation is designed to be incremental and uses a checkpoint system to avoid redundant computations.

Execution Command

jupyter nbconvert --to notebook --execute \
  --ExecutePreprocessor.timeout=0 \
  --output small_sample_benchmark_executed.ipynb \
  small_sample_benchmark.ipynb

Key Behaviors

• Timeout: The --ExecutePreprocessor.timeout=0 flag is required because the simulation can take several hours.
• Incremental Runs: If you increase n_rep in the notebook configuration, the system will only compute the new replicates and append them to the existing results in benchmark_results/.
• Starting Fresh: To invalidate existing results (e.g., after changing data or noise calibration), delete the benchmark_results/_meta.rds file. The system uses an MD5 hash of all simulation parameters to detect configuration changes and will automatically clear old results if a mismatch is detected.