Skip to content
This repository was archived by the owner on Nov 20, 2025. It is now read-only.

susie.ash fixes, citation fixes #35

Merged
alexmccreight merged 1 commit into from
Nov 3, 2025
Merged

susie.ash fixes, citation fixes #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 29 additions & 10 deletions R/sufficient_stats_methods.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -332,30 +332,49 @@ update_variance_components.ss <- function(data, params, model, ...) {
omega <- matrix(rep(omega_res$diagXtOmegaX, L), nrow = L, ncol = data$p, byrow = TRUE) +
matrix(rep(1 / model$V, data$p), nrow = L, ncol = data$p, byrow = FALSE)

# Update the sparse effect variance
sparse_var <- mean(colSums(model$alpha * model$V))

# Update sigma2 and tau2 via MoM
mom_result <- mom_unmappable(data, params, model, omega, model$tau2)

# If tau2 = 0 skip ash entirely
if (mom_result$tau2 == 0) {
return(list(
sigma2 = mom_result$sigma2,
tau2 = 0,
theta = rep(0, data$p)
))
}

# Remove the sparse effects
b <- colSums(model$alpha * model$mu)
residuals <- data$y - data$X %*% b

# Build ASH grid based on sparse variance and scaling factor
if (sparse_var * params$ash_scaling_factor > 1e-8) {
grid_factors <- exp(seq(log(0.01), log(10), length.out = 20 - 1))
est_sa2 <- c(0, sparse_var * params$ash_scaling_factor * grid_factors)
# Get minimum of active sparse effect variances as upper bound
active_V <- model$V[model$V > 0.01]
if (length(active_V) > 0) {
min_sparse_var <- min(active_V)
} else {
est_sa2 <- c(0, (2^(0.05*(1:20-1)) - 1)^4)
est_sa2 <- est_sa2 * (0.01 / max(est_sa2))
min_sparse_var <- 0.01 # Default when no active sparse effects
}

# Region 1: Very dense near zero (tau2/100 to tau2/10)
small_grid <- exp(seq(log(mom_result$tau2/100), log(mom_result$tau2/10), length.out = 10))

# Region 2: Moderate density around tau2 (tau2/10 to tau2*3)
medium_grid <- exp(seq(log(mom_result$tau2/10), log(mom_result$tau2*3), length.out = 6))

# Region 3: Sparse coverage towards minimum sparse effect
large_grid <- exp(seq(log(mom_result$tau2*3), log(min_sparse_var), length.out = 3))

est_sa2 <- c(0, unique(c(small_grid, medium_grid, large_grid)))

# Call mr.ash with residuals
mrash_output <- mr.ash.alpha::mr.ash(
X = data$X,
y = residuals,
sa2 = est_sa2,
intercept = FALSE,
standardize = FALSE,
sigma2 = model$sigma2,
sigma2 = mom_result$sigma2,
update.sigma2 = params$estimate_residual_variance,
max.iter = 3000
)
Expand Down
3 changes: 1 addition & 2 deletions R/susie.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -279,7 +279,6 @@ susie <- function(X, y, L = min(10, ncol(X)),
estimate_prior_variance = TRUE,
estimate_prior_method = c("optim", "EM", "simple"),
unmappable_effects = c("none", "inf", "ash"),
ash_scaling_factor = 0.15,
check_null_threshold = 0,
prior_tol = 1e-9,
residual_variance_upperbound = Inf,
Expand Down Expand Up @@ -311,7 +310,7 @@ susie <- function(X, y, L = min(10, ncol(X)),
prior_weights, null_weight, standardize, intercept,
estimate_residual_variance, estimate_residual_method,
estimate_prior_variance, estimate_prior_method,
unmappable_effects, ash_scaling_factor, check_null_threshold, prior_tol,
unmappable_effects, check_null_threshold, prior_tol,
residual_variance_upperbound, model_init, coverage,
min_abs_corr, compute_univariate_zscore, na.rm,
max_iter, tol, convergence_method, verbose, track_fit,
Expand Down
2 changes: 0 additions & 2 deletions R/susie_constructors.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@ individual_data_constructor <- function(X, y, L = min(10, ncol(X)),
estimate_prior_variance = TRUE,
estimate_prior_method = "optim",
unmappable_effects = "none",
ash_scaling_factor = 0.15,
check_null_threshold = 0,
prior_tol = 1e-9,
residual_variance_upperbound = Inf,
Expand Down Expand Up @@ -141,7 +140,6 @@ individual_data_constructor <- function(X, y, L = min(10, ncol(X)),
estimate_prior_variance = estimate_prior_variance,
estimate_prior_method = estimate_prior_method,
unmappable_effects = unmappable_effects,
ash_scaling_factor = ash_scaling_factor,
check_null_threshold = check_null_threshold,
prior_tol = prior_tol,
residual_variance_upperbound = residual_variance_upperbound,
Expand Down
10 changes: 5 additions & 5 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,11 +85,11 @@ If you use any of the summary data methods such as `susie_ss` or
If you use the Servin-Stephens prior on residual variance estimates
(`estimate_residual_method = "Servin_Stephens"`), please also cite:

> Denault, W.R.P., Carbonetto, P., Li, R., Consortium, A.D.F.G., Wang,
> G. & Stephens, M. (2025). Accounting for uncertainty in residual
> variances improves calibration of the "Sum of Single Effects" model
> for small sample sizes. *bioRxiv*, 2025-05. Under review for *Nature
> Methods*.
> Denault, W.R.P., Carbonetto, P., Li, R., Alzheimer's Disease Functional
> Genomics Consortium, Wang, G. & Stephens, M. (2025). Accounting for
> uncertainty in residual variances improves calibration of the "Sum of
> Single Effects" model for small sample sizes. *bioRxiv*, 2025-05.
> Under review for *Nature Methods*.

If you use infinitesimal effects modeling (`unmappable_effects = "inf"`),
please also cite:
Expand Down
12 changes: 6 additions & 6 deletions inst/CITATION
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,18 +68,18 @@ bibentry(header = paste("If estimate_residual_method = 'Servin_Stephens' is used
author = c(person("William", "R.P. Denault"),
person("Peter", "Carbonetto"),
person("Ruixi", "Li"),
person("A.D.F.G.", "Consortium"),
person("Alzheimer's Disease Functional Genomics Consortium"),
person("Gao", "Wang"),
person("Matthew", "Stephens")),
journal = "bioRxiv",
year = "2025",
note = "Under review for Nature Methods",
textVersion =
paste("Denault, W.R.P., Carbonetto, P., Li, R., Consortium, A.D.F.G.,",
"Wang, G. & Stephens, M. (2025). Accounting for uncertainty in",
"residual variances improves calibration of the 'Sum of Single",
"Effects' model for small sample sizes. bioRxiv, 2025-05.",
"Under review for Nature Methods."))
paste("Denault, W.R.P., Carbonetto, P., Li, R., Alzheimer's Disease",
"Functional Genomics Consortium, Wang, G. & Stephens, M. (2025).",
"Accounting for uncertainty in residual variances improves",
"calibration of the 'Sum of Single Effects' model for small",
"sample sizes. bioRxiv, 2025-05. Under review for Nature Methods."))

bibentry(header = paste("If unmappable_effects = 'inf' is used,",
"please also cite:"),
Expand Down