Repository navigation
Conversation
Tasks 2.1 (port remap_data, colScale) and 2.2 (port gen_wavelet_indx, DWT helper) landed in commit 10b7481. Marking complete in tasks.md. openspec list --changes now reports 9/71 tasks complete.
…ks 2.3-2.8)
Author pass for PR group 2 tasks 2.3-2.8. Implements the internal
mf_individual data class and the per-modality forward / inverse
DWT wrappers, with bit-identity to fsusieR / mvf.susie.alpha on the
DWT pipeline locked by the C2 and C3 fidelity tests.
## R/dwt.R (new)
- mf_dwt(Y_m, pos_m, max_padded_log2, filter_number, family,
verbose). Per-modality forward pipeline: remap_data ->
col_scale -> dwt_matrix -> pack as cbind(D, C) ->
gen_wavelet_indx. Univariate (T_m = 1) short-circuit:
centers/scales the single column, returns it as a length-1
"wavelet coefficient" with scale_index = 1L.
- mf_invert_dwt(D_packed, column_center, column_scale,
filter_number, family). Inverse via wavethresh::wr on a wd
skeleton (set $D and last entry of $C, then wr), reversing the
per-position centering and scaling. Univariate path is a single
multiply-add.
## R/data_class.R (new)
- create_mf_individual(X, Y, pos, max_padded_log2,
wavelet_filter_number, wavelet_family, standardize, intercept,
save_residuals, verbose). Internal constructor returning class
c("mf_individual", "individual"). Validates inputs, centers /
standardizes X, runs mf_dwt per modality, packages everything
the IBSS loop needs. Per design.md D2 / mf-data-class/spec.md.
- Univariate (T_m = 1), single-modality (M = 1), and ragged
multi-modality inputs all flow through one class with shape
branches, mirroring mvsusieR/refactor-s3's mv_individual.
## Tests
- tests/testthat/test_dwt.R: 30 expectations across 9 test_that
blocks. Forward shape / univariate short-circuit / DWT roundtrip
identity / Haar-wavelet bit-identity / C2 fidelity vs
fsusieR::susiF pipeline / C3 fidelity vs mvf.susie.alpha::DWT2
row-wise / verbose message forwarding. The Haar-wavelet
bit-identity test (filter_number = 1, tolerance 1e-12) pins the
absence of bug-in-our-port separately from wavethresh's ~2e-9
long-filter accumulation noise at filter_number = 10.
- tests/testthat/test_data_class.R: 52 expectations across 18
test_that blocks. Input validation / single-modality functional
/ ragged multi-modality / univariate short-circuit / mixed
univariate-and-functional / X centering / X scaling / zero-
variance X column / intercept = FALSE / save_residuals on / off
/ DWT cache contents agree with mf_dwt / default pos.
- All Phase 4 fidelity tests at tolerance 0 (bit-identical) or
1e-12 / 1e-8 per the binary tolerance philosophy. No graduated
tolerances.
## inst/notes/refactor-exceptions.md
PR group 2 entries appended for:
- mvf.susie.alpha/R/utils_wavelet_transform.R:24-50 (DWT2):
replaced-by-`R/utils_wavelet.R::dwt_matrix`.
- mvf.susie.alpha/R/utils_wavelet_transform.R:64-68 (pack_dwt):
replaced-by-inline (`cbind(D, C)` pack inlined into mf_dwt).
- mvf.susie.alpha/R/multfsusie.R:255-318 (per-modality DWT
pipeline): replaced-by-mf_dwt + create_mf_individual.
- fsusieR/R/operation_on_susiF_obj.R:1785-1796 (inverse-DWT
reconstruction skeleton): replaced-by-mf_invert_dwt.
## Manuscript-citation convention update
Reviewer pass on the author-pass diff caught that `@manuscript_ref`
is not a registered roxygen tag, so `roxygen2::roxygenise()` errors
on it. The fix switches the canonical citation form to the built-in
`@references` tag with a `Manuscript:` prefix:
```r
#' @references
#' Manuscript: methods/<file>.tex eq:<label>
```
This works with roxygen, renders cleanly in `?function`, stays
greppable for the reviewer-checklist scan, and avoids the cost of
writing a custom roxygen tag handler. Convention update propagated
to:
- CLAUDE.md (Phase 3 description)
- inst/notes/refactor-discipline.md (binding policy)
- inst/notes/test-fidelity-policy.md (binding policy)
- inst/notes/review-loop-methodology.md (reviewer item 3)
- inst/openspec/changes/.../proposal.md, design.md (D12),
specs/mf-ibss/spec.md
- tasks.md task 9.1 wording
The forbidden-token grep at task 9.1 now also forbids the literal
string `@manuscript_ref` to catch any reintroduction.
## inst/notes/review-loop-methodology.md "Reviewer model selection"
New subsection per Gao 2026-04-25: reviewer subagent on numerical
PRs SHALL use `model: "opus"`; documentation / lighter passes can
use sonnet. Records the rate-limit-fallback rule.
## Verification
devtools::test(): 162 / 162 expectations pass
roxygen2::roxygenise(): clean (no errors, no warnings)
R CMD check: 0 errors, 0 warnings, 1 NOTE
(declared Imports not yet used; resolves
in PR groups 3-6)
openspec validate --strict: valid
9/71 + 6 ticked = 15/71 tasks done.
## Review-loop status
- Author pass: tasks 2.3-2.8 + manuscript-citation convention
update.
- D13 audit substep: not run on R/dwt.R / R/data_class.R because
these are net-new files with no port-source code (the helpers
in R/utils_wavelet.R that they call were already audited).
- Reviewer pass 1 (sonnet): two findings; both fixed.
- Reviewer pass 2 (sonnet): clean.
- Reviewer pass 3 (opus): clean.
- Two consecutive clean reviewer passes achieved per
inst/notes/review-loop-methodology.md.
Tasks 2.3-2.8 ticked in tasks.md.
NAMESPACE re-generated by roxygen2 (only the autogenerated
header line; no exports to add yet).
Author pass for tasks 3.1-3.3. Implements the internal model class
constructor and the six per-effect getter / setter methods that the
IBSS loop will dispatch on. The constructor returns a list of class
c("mfsusie", "susie") so susieR generic accessors (susie_get_pip,
susie_get_cs) work where their semantics match the multi-modality
functional setting.
## R/model_class.R (new)
- initialize_susie_model.mf_individual(data, params, var_y, ...).
Pre-allocates the IBSS-iterated subset of the spec field set:
alpha (L x J, 1/J uniform), mu / mu2 (nested list[L][M] of
J x T_padded[m] zero matrices per design.md D2), V (length-L
per-effect prior scalar scaling, init 1), V_grid / pi_V /
null_prior_weight (lifted from params$prior; NULL until PR
group 4 builds mf_prior_scale_mixture), KL / lbf / lbf_variable
(NA at the right shape), sigma2 (params$residual_variance
with var_y fallback per design.md D2), pi (variable-selection
prior, length-J, defaults to uniform 1/J), fitted (list[M] of
zero matrices), intercept (length-M zero vector), L / J / M
storage helpers (all integer).
- Post-fit fields (elbo, niter, converged, pip, cs, residuals,
csd_X, wavelet_meta, call) are deferred to ibss_finalize.mf_individual
in PR group 7 per design.md D2.
- L stored as integer via as.integer(params$L). Verified that
susieR does NOT carry $L on its fit object (susieR uses dim(alpha)
to recover L), so no upstream coordination is needed for type
consistency in the C1 contract.
## R/model_methods.R (new)
Six accessors mirroring mvsusieR/R/model_methods.R signatures:
- get_alpha_l.mfsusie(model, l) -> J-vector (alpha[l, ])
- get_posterior_moments_l.mfsusie(model, l) -> list(mu, mu2) of
list[M] matrices
- get_posterior_mean_l.mfsusie(model, l) -> list[M] of
alpha[l, ] * mu[[l]][[m]]
- get_posterior_mean_sum.mfsusie(model) -> list[M], sum across L
- get_prior_variance_l.mfsusie(model, l) -> scalar V[l]
- set_prior_variance_l.mfsusie(model, l, V) -> updates V[l],
returns model
## Tests
tests/testthat/test_model_class.R: 71 expectations across 22
test_that blocks. Covers:
- class hierarchy c("mfsusie", "susie")
- alpha / mu / mu2 / fitted shape and zero-init
- KL / lbf / lbf_variable NA shape
- V (per-effect prior scalar) ones init
- sigma2 fallback to var_y when params$residual_variance is NULL
- sigma2 takes params$residual_variance when supplied
- pi defaults to uniform 1/J when prior_weights is NULL
- L / J / M integer storage
- degenerate cases: M = 1 functional, T_m = 1 univariate, ragged M = 3
- all six accessor contracts with deterministic injected mu / alpha
- set_prior_variance_l updates the right entry, leaves others alone
Test header carries the apple-to-apple comparison metadata per
design.md D12 (paradigm reference, not port source: mvsusieR/R/
mvsusie_constructors.R:341-427 and susieR/R/individual_data_methods.R).
## Review-loop status
This PR group is bespoke (paradigm-derived from mvsusieR but not a
direct line port), so the D13 port-quality audit substep does not
apply.
- Reviewer pass 1 (opus): 0 blockers, 2 notes — both forward-looking
concerns about consumers in PR groups 5-7, deferred:
1. R/model_class.R:111 intercept = rep(0, M) is length-M scalar
vector; if PR group 7 finalize wants per-(modality, position)
intercepts, the constructor will need to widen. Spec does not
pin the shape; not blocking.
2. R/model_methods.R:48-49 get_posterior_mean_sum reads model$mu[[1]]
to size the output, requiring L >= 1. Mfsusie() public API pins
L >= 1 default; no realistic caller hits L = 0. Not blocking.
- Reviewer pass 2 (opus): clean, both prior notes confirmed
reasonable to defer, no new findings.
- Two consecutive clean reviewer passes per
inst/notes/review-loop-methodology.md.
## Verification
devtools::test(): 233 / 233 expectations pass (was 162;
+71 in test_model_class.R)
roxygen2::roxygenise(): clean
R CMD check: 0 errors, 0 warnings, 1 NOTE (declared
Imports not yet used; resolves in PR groups
4-6)
openspec validate --strict: valid; 18 / 71 tasks done (was 15).
Tasks 3.1-3.3 ticked.
Per Gao 2026-04-25: splitting into _class.R + _methods.R is the mvsusieR / susieR convention because those packages have multiple constructors per class and ~15 methods. mfsusieR's scope (one constructor + 6 accessors per class) does not justify the split; a single file with a section divider is cleaner. Single-file pattern: - R/data_class.R holds create_mf_individual (one section, one constructor; future per-class methods on mf_individual will join this file with a section divider when PR group 5 lands). - R/model_class.R now holds initialize_susie_model.mf_individual AND the six accessors get_alpha_l / get_posterior_moments_l / get_posterior_mean_l / get_posterior_mean_sum / get_prior_variance_l / set_prior_variance_l, separated by a comment-block divider. CLAUDE.md "Naming rules" section gains a binding rule: class- related code lives in ONE file per class; do NOT split into _class.R + _methods.R unless / until the file genuinely outgrows single-file readability. Section divider format documented. Tests pass unchanged (233 / 233). roxygen clean.
The greedy outer loop landed on susieR master in this session (linear +L_greedy step, min(lbf) < lbf_min saturation criterion, warm-start across rounds via params$model_init). mfsusieR's mfsusie() puts L_greedy and lbf_min in the params list and lets susieR::susie_workhorse handle the rest. Updates: - design.md D5: drop the "Generalize susieR with L_greedy = NULL argument. Tracked outside this OpenSpec change but a hard prerequisite for Phase 3" language. Replace with: mvf.susie.alpha and fsusieR's per-iteration interleaved greedy is NOT ported; mfsusieR uses the susieR wrapper directly. - design.md D7 signature: add lbf_min = 0.1 to mfsusie() arg list. - design.md D7 renamed-args row for L_start -> L_greedy: note lbf_min is also exposed and passed through. - design.md D11d binary tolerance table: add two apple-to-orange rows for greedy-L vs legacy greedy. Add a paragraph clarifying that C2 / C3 fidelity tests run with greedy disabled on the reference side (fsusieR::susiF(greedy = FALSE, backfit = FALSE); mvf.susie.alpha::multfsusie(greedy = FALSE, backfit = FALSE)) AND mfsusieR L_greedy = NULL, so the algorithmic difference is taken out of the comparison. - design.md "External coordination: susieR L_greedy generalization" Migration Plan section: replaced with a short note that the upstream susieR change landed in this session. No graceful- degradation shim; no commit-hash pin (we ship against the local susieR fork until / unless upstream merges). - tasks.md group 0: dropped the two-task graceful-degradation prerequisite. Replaced with a single task 0.1 marked complete citing the susieR change landing. - specs/mf-ibss/spec.md "greedy R-search via L_greedy passthrough" requirement: replaced post-/pre-generalization scenario pair with a passthrough scenario, an L_greedy=NULL disable scenario, and an apple-to-orange smoke-test scenario for greedy vs legacy greedy. - specs/mf-public-api/spec.md formalArgs requirement: add lbf_min to the required argument list. EM-literature framing: The susieR wrapper maps onto Verbeek-Vlassis-Krose 2003 Greedy EM for Gaussian mixtures: full inference at fixed K, between-round dimension change with a Bayes-factor-based saturation criterion. The fsusieR/mvf.susie.alpha interleaved approach maps onto Celeux-et-al CEM2 in spirit but without CEM2's principled acceptance criterion. The susieR wrapper preserves variational ELBO monotonicity at each round; the interleaved approach loses it. See the chat-side analysis in this session. No code changes under R/. No tests touched. openspec validate --strict passes; 18/71 tasks remain at 18 (no progress count change since group 3 commit; this commit is a docs / spec delta).
Captures four design decisions made during PR group 4 prior class exploration. No R/ code changes; OpenSpec artifacts and policy docs only. (a) Bhat/Shat compute -- delegate to upstream susieR helper A new susieR helper compute_marginal_bhat_shat(X, Y, ...) lands on susieR master (separate from but parallel to L_greedy). Authorized by Gao on 2026-04-25 as a CLAUDE.md hard rule #1 override. Shared by: - susieR's compute_ser_statistics.individual (T = 1 path) - mvsusieR's OLS path in compute_ser_statistics.mv_individual (when GLS data$svs is NULL) - mfsusieR's prior init in R/prior_scale_mixture.R AND PR-group-5 per-effect SER mfsusieR does NOT port fsusieR::cal_Bhat_Shat locally; it calls the susieR helper instead. design.md D5 row updated; design.md Migration Plan adds an "Upstream susieR compute_marginal_bhat_shat (planned)" subsection. (b) Cross-modality prior plug-in seam -- S3 generic Renamed seam dispatch verb from `apply` to `combine_modality_lbfs` (S3 generic). Avoids collision with base::apply and matches the S3 idiom used by susie / mvsusie. The default no-op implementation lives at combine_modality_lbfs.mf_prior_cross_modality_independent. future.apply note: the per-modality SER lapply (which feeds combine_modality_lbfs) is the parallelizable target, not the combine itself. Flagged as Phase 7 perf candidate. (c) mf_prior_scale_mixture signature -- take data class Constructor takes the mf_individual data class directly (susieR/mvsusieR idiom). No pre-computed Bhat/Shat hand-off. (d) Honor prior_variance_grid; no T_m=1 branch; warning at mfsusie() level for any T_m <= 3 mf_prior_scale_mixture honors a user-supplied prior_variance_grid (used by C1 susieR-degeneracy contract); when NULL, falls back to the data-driven path via the susieR helper + ashr::ash. NO T_m == 1 branch. The public mfsusie() emits susieR::warning_message when any T_m <= 3 advising the user to consider susieR::susie / fsusieR::susiF, but still completes the fit. The earlier hardcoded T_m=1 stub (pi = c(0.8, 0.2), sd = c(0, 1)) was a bug that would have broken the C1 contract; this restructure fixes it. OpenSpec artifact updates: design.md D5: add row for cal_Bhat_Shat full-Y delegation to susieR helper. design.md D6: rename cross-modality verb to combine_modality_lbfs; add future.apply note. design.md Migration Plan: new "Upstream susieR compute_marginal_bhat_shat (planned)" subsection. specs/mf-prior/spec.md: replace "ported from fsusieR" with "uses susieR helper + ash fit"; add scenarios for user-supplied grid path, data-driven path, and T_m <= 3 warning. Cross-modality requirement updated for the verb rename. tasks.md group 4: restructure 4.1-4.6; add 4.0 marking the susieR upstream coordination; add 4.3 for the mfsusie()-level warning. DESCRIPTION: add Rfast to Imports (used by the susieR helper's column-variance path; pulled in for the data-driven prior init). Policy update -- inst/notes/refactor-discipline.md section 0 (NEW, BINDING): Upstream-first principle. Before porting any routine from mvf.susie.alpha or fsusieR, FIRST check whether susieR has a similar helper or whether one can be added or generalised upstream for shared use. Code duplication from susieR is forbidden. Lists the two landed/planned upstream changes (L_greedy in 4caf755, compute_marginal_bhat_shat planned) as concrete examples. openspec validate --strict: clean. Implementation in susieR / mvsusieR / mfsusieR follows in separate commits per Gao's review-before-commit pattern.
Implements PR group 4 tasks 4.1-4.6 per design.md D6 and the mf-prior/spec.md updates committed in ef5d6b3 (PR group 4 design capture). Uses the new susieR upstream helper compute_marginal_bhat_shat (susieR commit eb2b878) per the upstream-first principle in inst/notes/refactor-discipline.md section 0. DESCRIPTION pin: susieR (>= 0.16.1) for the new helper. R/prior_scale_mixture.R (new) distribute_mixture_weights(K, null_prior_weight): private helper. Returns a (K + 1)-length pi vector with the null component on slot 1, remaining mass split uniformly across K non-null components. init_scale_mixture_prior_default(Y_m, X, prior_class, scale_index, grid_multiplier): per-modality data-driven prior init. Calls susieR::compute_marginal_bhat_shat(X, Y_m) for the marginal (Bhat, Shat) sample, then ashr::ash to fit the K-vector mixture grid. Behaviour-preserving port of fsusieR::init_prior.default for the mixture_normal and mixture_normal_per_scale branches; same set.seed(1) sequence and sample-size caps (5000 / 50000); replicates the ash fit per scale for the per-scale variant. mf_prior_scale_mixture(data, prior_variance_grid, prior_variance_scope, null_prior_weight, grid_multiplier): public-ish wrapper (internal, not exported). Two paths: (a) prior_variance_grid non-NULL: use the user-given grid and distribute_mixture_weights for pi. No ash fit. Used by the susieR-degeneracy contract C1. (b) prior_variance_grid = NULL: per modality, run init_scale_mixture_prior_default. Default path. NO T_m == 1 branch -- same code path for any T_m. The public mfsusie() will emit warning_message at any T_m <= 3 in PR group 7. R/prior_cross_modality.R (new) combine_modality_lbfs (S3 generic): joint log-BF combination across modalities. Renamed from the earlier "apply" verb to avoid base::apply collision; matches the S3 idiom used by susie / mvsusie. combine_modality_lbfs.mf_prior_cross_modality_independent (S3 method, registered via @exportS3Method): the default no-op implementation. Reduce("+", modality_lbfs). cross_modality_prior_independent(): factory for the trivial independent prior. Default in mfsusie() per design.md D6. NAMESPACE auto-updated with the S3method directive. tests/testthat/test_prior_scale_mixture.R (new, 42 expectations) - distribute_mixture_weights shape and weight-sum invariants; null_prior_weight = 0 produces c(0, 1) for K = 1. - User-supplied prior_variance_grid path: V_grid honoured per modality, pi matches distribute_mixture_weights output. - susieR-degeneracy C1 inputs (length-1 grid, null = 0) produce a single-component prior with pi = c(0, 1). - Data-driven path bit-identical to fsusieR::init_prior.default on the same fixture: pi, sd, mean, Bhat, Shat all at tolerance 1e-12 (max diff = 0). Validates the C2 fidelity contract. - T_m = 1 with NULL grid runs the data-driven path without crashing (no special-case branch). - per_modality scope collapses scale dimension correctly. - cross_modality_prior_independent class hierarchy correct; combine_modality_lbfs sums per-modality vectors; works for M = 1 trivial case. inst/notes/refactor-exceptions.md (PR group 4 entries) - fsusieR/R/computational_functions.R:69-155 (cal_Bhat_Shat full-Y) replaced-by-susieR::compute_marginal_bhat_shat. - fsusieR/R/operation_on_prior.R:42-185 (init_prior.default) replaced-by-init_scale_mixture_prior_default. - fsusieR/R/operation_on_prior.R:79-165 (lowc_wc filtering) deferred-to-Phase-3-future-PR. - mvf.susie.alpha/R/operation_on_multfsusie_prior.R:17-138 replaced-by-mf_prior_scale_mixture. mvsusieR adoption: deliberately none. mvsusieR's per-effect SER works on cached XtR (model$residuals = X^T r), not raw Y. The susieR helper's input pattern (X, Y matrix) doesn't match mvsusieR's flow without forcing a recompute of crossprod(X, R) that's already cached. mvsusieR's Shat is a list of J R x R GLS covariance matrices, structurally different from the helper's scalar Shat. Documented in chat-side analysis. Verification devtools::test(): 275 / 275 expectations pass (was 233; +42 in test_prior_scale_mixture) R CMD check: 0 errors, 0 warnings, 1 NOTE (methods + mixsqp not yet used; will resolve in PR groups 5-6) roxygen2::roxygenise(): clean C2 fidelity vs fsusieR: bit-identical at 1e-12 Tasks 4.0-4.6 ticked when openspec apply moves to group 5.
4.0 (upstream susieR helper) landed in susieR commit eb2b878. 4.1, 4.2 (mf_prior_scale_mixture + cross_modality_prior) and 4.4-4.6 (tests + ledger) landed in mfsusieR commit a53269a. 4.3 stays unticked: the mfsusie()-level warning_message lives in PR group 7 task 7.1 (mfsusie() public API doesn't exist yet). Progress: 26 / 71 tasks complete.
Implements the per-effect S3 dispatch surface for the IBSS
inner loop on `mf_individual`:
compute_residuals.mf_individual per-modality partial
residual `R_m = D[[m]] -
Xr_without_l_m`; caches
`XtR_m = X^T R_m` on
`model$residuals[[m]]`.
compute_ser_statistics.mf_individual per-modality `betahat[[m]]
= XtR_m / predictor_weights`,
`shat2[[m]] = sigma2_m /
predictor_weights`. Handles
scalar sigma2_m (legacy
shared_per_modality) and
per-scale sigma2_m vector
(per-(scale, modality) default,
broadcast via scale_index).
update_fitted_values.mf_individual folds the new posterior
contribution back into
`model$fitted[[m]] =
fitted_without_l[[m]] +
X %*% (alpha_l * mu_l[[m]])`.
Caching strategy mirrors susieR / mvsusieR:
data$predictor_weights = colSums(X^2) cached at construction
(was missing from the
mf_individual data class;
added in this commit).
model$residuals[[m]] per-modality XtR cache,
set by compute_residuals,
consumed by SER stats.
model$fitted_without_l[[m]] running-fit minus effect-l
contribution, used by
downstream KL/ELBO methods
(PR group 5c).
model$raw_residuals[[m]] per-modality `R_m` cache,
used by `Eloglik` in 5c.
The per-modality lapply is BLAS-bound (matrix multiplies via
crossprod / `%*%`); no C++ benefit at this layer. C++ may help
in 5b (mixture posterior moments, per-component inner loop) and
will be assessed there.
mvsusieR's optim BQ_cache pattern (caching V-independent products
across the optim Brent inner iterations) is V-dependent on the
prior_variance update step, which lives in 5b. 5a's caches are
V-independent intermediates that already get computed once per
(effect, iteration) and reused.
Tests:
tests/testthat/test_individual_data_methods.R (29 expectations)
- predictor_weights cache equals colSums(X^2).
- compute_residuals shape and per-modality XtR identity.
- compute_residuals at zero-init effect: raw_residuals == D[[m]].
- compute_ser_statistics shapes; betahat = XtR / pw;
shat2 broadcasts scalar and per-scale sigma2.
- update_fitted_values reconstructs fitted from
fitted_without_l + X * (alpha_l * mu_l[[m]]).
- Round-trip: zero-effect IBSS sweep keeps fitted at zero.
- C3 fidelity: compute_ser_statistics betahat matches
mvf.susie.alpha::cal_Bhat_Shat (via fsusieR:::cal_Bhat_Shat
since they share the implementation) at tolerance 1e-12 on
the initial all-zero-effects residual.
Suite state:
devtools::test() 304 / 304 expectations pass
(was 275; +29 in test_individual_data_methods)
R CMD check 0 errors, 0 warnings, 1 NOTE
(methods + mixsqp not yet used; resolves
in PR group 5c / 6)
roxygen2::roxygenise() clean
Tasks 5.1 (compute_residuals + compute_ser_statistics +
update_fitted_values portion) tick when 5b/5c land. PR group 5
remaining: 5b (calculate_posterior_moments, mixture-aware), 5c
(compute_kl, loglik, neg_loglik).
Implements per-effect mixture-of-normals log-Bayes factors and posterior moments on `mf_individual` via a cpp11 hot-path with a pure-R reference oracle. Per design.md D14 (Phase 3 cpp11 stance): LinkingTo: cpp11 only, element-wise kernels, oracle agreement at <= 1e-12. Adds: - src/posterior_mixture.cpp: mixture_log_bf_per_scale_cpp, mixture_posterior_per_scale_cpp. Column-major-friendly loop order (i outer, j inner). cpp11 ~2x faster than pure-R reference. - R/reference_implementations.R: pure-R oracles (mixture_log_bf_per_scale_R, mixture_posterior_per_scale_R) for the cpp11-vs-R fidelity test at <= 1e-12. - R/posterior_mixture.R: thin wrappers with input validation. - R/individual_data_methods.R: loglik.mf_individual, calculate_posterior_moments.mf_individual, and the shared mf_per_modality_bhat_shat helper. - R/model_class.R: stash G_prior and cross_modality_combiner on the fit object (avoids re-traversing params$prior per IBSS step). - R/zzz.R: @useDynLib mfsusieR, .registration = TRUE. Tests (39 expectations added across two files): - cpp11 vs pure-R oracle at <= 1e-12 (D14 contract). - pure-R oracle vs ashr::calc_logLR / ashr::postmean / ashr::postsd at <= 1e-12 (closed-form fidelity). - K = 1, no null reduces to susieR's scalar SuSiE formula post_var = (1/V + pw/sigma2)^(-1) at <= 1e-12 (C1 contract). - alpha sums to 1, lbf finite, lbf_variable length J. - l = NULL returns scalar lbf_model that matches the l-supplied call's model$lbf[l]. - calculate_posterior_moments matches a per-(SNP, scale) ashr loop at <= 1e-12. OpenSpec: - design.md D14 added: cpp11 stance, oracle requirement, tolerance table updated with Phase 3 cpp11 vs oracle row. - Non-Goals updated: wholesale Rcpp/RcppArmadillo deferred to Phase 7. - tasks.md PR group 5 split into 5a (LANDED), 5b (this commit), 5c (KL/ELBO). - Three-agent review (reuse, quality, efficiency) addressed: snake_case slice arg names, manuscript reference on loglik, tighter test tolerance (1e-12), cpp11 loop reorder, helper extraction. Deferred to 5c: registerS3method on susieR namespace, move G_prior to .internal sublist. Deferred to Phase 7: fused kernel, OpenMP, Xb_l_m caching. R CMD check: 0 errors, 0 warnings, 1 NOTE (methods/mixsqp imports declared but not yet used; will resolve in PR group 6).
…tion (PR group 5c)
Adds the per-effect KL and SER posterior-expected-loglikelihood
methods on `mf_individual`, the `neg_loglik` wrapper for the V-optim
caller (PR 6), and wires `.onLoad` to register all `mf_individual`
and `mfsusie` S3 methods on susieR's namespace so `susieR::ibss_fit`
dispatches to mfsusieR.
Adds:
- R/individual_data_methods.R:
- `SER_posterior_e_loglik.mf_individual` — per modality, per
position; broadcasts `model$sigma2[[m]]` via the new
`mf_sigma2_per_position` helper. Generalizes susieR's scalar
formula to multi-modality multi-functional with per-(scale,
modality) residual variance.
- `compute_kl.mf_individual` — `KL_l = -[lbf_l + L_null] + L_post`,
mirrors susieR's `compute_kl.individual` Gaussian-KL branch.
- `neg_loglik.mf_individual` — `-loglik(exp(V_param), ...)`.
- Helper `mf_sigma2_per_position(data, model, m)` factors out the
per-(scale, modality) sigma2 broadcast (shared by three call sites:
`mf_per_modality_bhat_shat`, `SER_posterior_e_loglik`, `compute_kl`).
- `loglik.mf_individual` refactored to use cached susieR helpers
`lbf_stabilization` + `compute_posterior_weights` instead of the
hand-rolled softmax (addresses PR 5b reuse review finding).
- R/zzz.R: `.onLoad` now caches susieR internals
(`lbf_stabilization`, `compute_posterior_weights`, `warning_message`,
`SER_posterior_e_loglik`) and registers S3 methods for nine
`mf_individual` generics and six `mfsusie` model-class generics.
Mirrors `mvsusieR/R/zzz.R` structure.
Tests (10 added, 367 total PASS):
- `SER_posterior_e_loglik` value matches the closed-form formula at
<= 1e-12.
- `SER_posterior_e_loglik` per-scale sigma2 path exercised
(`model$sigma2[[m]] <- seq_len(S_m) * 0.1`); verifies
`mf_sigma2_per_position` round-trips correctly.
- `compute_kl` returns finite KL[l]; identity
`KL[l] == -lbf[l] - L_null + L_post` at <= 1e-12.
- `neg_loglik` exponentiates V_param and negates loglik at <= 1e-12.
- `.onLoad` registration test via `getS3method(g, class, optional =
TRUE, envir = asNamespace("susieR"))` for all nine `mf_individual`
generics. Cached binding sanity for `lbf_stabilization`,
`compute_posterior_weights`, `warning_message`.
Three-agent review (reuse / quality / efficiency) addressed:
- Reuse: cached susieR helpers (`lbf_stabilization`,
`compute_posterior_weights`, `SER_posterior_e_loglik` generic);
`mf_sigma2_per_position` extracted; `compute_kl` uses generic
dispatch (`SER_posterior_e_loglik(...)`) for subclass hygiene.
- Quality: per-scale sigma2 test variant added; manuscript reference
added on `SER_posterior_e_loglik`.
- Efficiency: high-severity finding (cache `Xb_l_m` between
`update_fitted_values` and `SER_posterior_e_loglik`) deferred to
Phase 7 — this redundancy is also present in
susieR's `SER_posterior_e_loglik.individual` (which recomputes
`compute_Xb(X, Eb)`), so mirroring upstream is correct for Phase 3.
Other low-severity items (R^2 cache, sigma2_per_pos cache) deferred
similarly.
Deferred to PR groups 7c (susieR-degenerate C1) and 7e (mvf.susie.alpha
C3): KL fidelity tests against port sources at <= 1e-10 / <= 1e-8.
R CMD check: 0 errors, 0 warnings, 1 NOTE (methods/mixsqp imports
declared but not yet used; resolves in PR group 6).
Two leftover follow-ups from the 5b quality review (the rest landed in 5b and 5c): - Add an all-null mixture edge-case test (`K = 1`, `sd_grid = 0`, `pi_grid = 1`) that exercises the `var_k[k] == 0.0` short-circuit branch in the cpp11 kernel and the corresponding `if (sd_k_var == 0)` branch in the R oracle. Asserts exact equality (`lbf = 0`, `pmean = pmean2 = 0`). - Add a manuscript-reference header to `test_posterior_mixture.R` pointing to `methods/derivation.tex eq:post_f_mix` and `eq:post_f2_mix`. The third reviewer "fix in 5c" item -- moving `G_prior` and `cross_modality_combiner` to a `model$.internal` sublist -- is deferred to the PR group that introduces `print.mfsusie` / `summary.mfsusie` (Phase 7/8), where the field-visibility decision actually has user-facing consequences. Tests: 370 PASS (was 367).
Three cleanup items addressed before starting PR 6:
1. Numerical fidelity vs fsusieR (transitive C3 for the per-scale
log-BF path). New test file `test_fsusier_fidelity.R` with 16
expectations against:
- fsusieR::log_BF.mixture_normal_per_scale (sum of per-scale
log-Bayes factors per SNP) — bit-identical at 5e-15
- fsusieR::cal_zeta (softmax with uniform pi) — exact equality
- ashr::postmean / ashr::postsd^2 + postmean^2 (per-scale
posterior moments) — bit-identical at 1e-12
This closes the largest correctness gap from the 5c audit:
the per-(SNP, modality, scale) log-BF path is now numerically
verified against the original implementation, transitively
covering mvf.susie.alpha which calls fsusieR::log_BF internally.
Full mvf.susie.alpha integration test (with init_multfsusie_obj
plumbing) remains in PR 7e per design.
2. shat2_marker readability. Added a 5-line comment in
loglik.mf_individual explaining what the dummy "1" value means
when calling lbf_stabilization with a zero-pw mask. Per the 5c
quality review.
3. design.md D14 records the measured profile of two cache
candidates ("high severity" per the 5c efficiency review):
- Xb_l_m cache between update_fitted_values and
SER_posterior_e_loglik: ~2.5% saving per IBSS step, not 30%
as the reviewer's flop-count estimate suggested. Defer.
- colSums(R^2) cache: <1% saving. Defer.
Both deferrals now have measured numbers in the design, not
just "Phase 7 candidate" handwaves.
Audit results:
- mfsusieR's `X %*% b` and `crossprod(X, R)` match mvsusieR's
pattern exactly (mvsusieR has its own compute_Xb that's just
data$X %*% b for the dense case, NOT susieR's attribute-based
helper). No X-multiply duplication to fix.
- susieR helpers already cached and reused: lbf_stabilization,
compute_posterior_weights, warning_message, SER_posterior_e_loglik
generic.
Tests: 386 PASS (was 370). R CMD check: 0E/0W/1N (unused
methods/mixsqp imports — used in PR 6).
Replaces the cached susieR helpers (`lbf_stabilization` + `compute_posterior_weights`) with 7 inlined lines in `loglik.mf_individual`. Equivalent in numerics; clearer in intent because the per-(modality, position) shat2 logic stays explicit. Removed: - `shat2_marker <- ifelse(data$predictor_weights == 0, Inf, 1)` and the comment block explaining why we passed it. - Cached `lbf_stabilization` and `compute_posterior_weights` bindings in `R/zzz.R`. They were used at exactly one call site; inlining is shorter than the call setup. The remaining cached binding `warning_message` and the cached S3 generic `SER_posterior_e_loglik` stay (used in PR 6 / 7 respectively). Tests: 385 PASS (was 386; the dropped test expectation was the existence-check on the two now-uncached bindings).
Implements update_variance_components, update_model_variance,
Eloglik, and get_objective on mf_individual, plus the EM helpers
they consume. mfsusieR now provides all S3 method overrides
needed for susieR::ibss_fit to dispatch end-to-end on the
multi-modality multi-functional data class.
Adds:
- R/em_helpers.R: mf_em_likelihood_per_scale (mixsqp likelihood
matrix builder) and mf_em_m_step_per_scale (mixsqp wrapper with
null-component handling).
- R/individual_data_methods.R extended:
- mf_get_ER2_per_position: per-modality, per-position E_q[||residual||^2]
using the SUSIER-CORRECT formula (||Y - X postF||^2
+ sum_l ((alpha_l . pw)^T mu2_l - ||X (alpha_l . mu_l)||^2)).
Reuses the cached running fit `model$fitted[[m]]` to skip the
outer X*postF matmul.
- update_variance_components.mf_individual: per-(scale, modality)
or shared-per-modality sigma2 update.
- update_model_variance.mf_individual: per-effect mixsqp on
pi_V[[m]] per (m, s); writes back to G_prior$fitted_g$pi.
- Eloglik.mf_individual: per-modality, per-position log-lik aggregate.
- get_objective.mf_individual: ELBO = Eloglik - sum(KL) +
sum_l alpha_l . log(pi / pmax(alpha_l, 1e-6)).
- R/zzz.R: registers 4 new generics on susieR's namespace; adds
the %||% R<4.4 fallback (mvsusieR pattern).
Refactor-exceptions ledger:
- New PR-6 section documenting the divergence from
mvf.susie.alpha's get_ER2 formula (missing pw factor + conflated
square-then-sum vs sum-then-square across L). Pattern A bug fix
per refactor-discipline section 3. The bug systematically
underestimates sigma^2, which inflates Bayes factors and PIPs;
this is a plausible mechanism for the FDR miscalibration
observed in mvf-original.md and will be verified empirically in
Phase 5.
- Entry split into two files (operation_on_multfsusie_obj.R and
computational_routine.R) per the schema.
D12 cleanup (port-source string mentions stripped from R/):
- R/individual_data_methods.R: removed "Reduces to mvf.susie.alpha's"
- R/prior_scale_mixture.R: removed "fsusieR::init_prior.default" /
"the fsusieR source" / "legacy mvf.susie.alpha mode" mentions
from roxygen.
Inlined softmax + zero-pw mask in loglik.mf_individual
(replaced the cached lbf_stabilization + compute_posterior_weights
+ shat2_marker hack from the prior reuse-review iteration -- 7 lines
inline is clearer than 3 lines + helper-cache + dummy marker vector).
Reference-test pinning:
- setup.R now exposes mfsusier_pinned_fsusier_sha() /
mfsusier_pinned_mvf_sha() (commits fdba52b / e41550a) and
skip_if_pinned_mvf_mismatch / skip_if_no_fsusier helpers that
warn when installed package versions diverge from the contract.
- test_mvf_alpha_fidelity.R: 4 new fidelity tests covering
per-modality log-BF rows, cross-modality combine, posterior
moments per-(m, s), and EM likelihood matrix vs
fsusieR::cal_L_mixsq_s_per_scale -- all at <= 1e-12.
Test coverage:
- 524 PASS (was 461). Coverage 95.53% overall, 100% on 10 of 11 R
files. R/zzz.R is 0% (.onLoad bindings, naturally not observable
to package_coverage; verified via the registration test).
- New tests cover per-scale sigma2 path, edge cases (M=1, L=1,
all-zero alpha, pure-null mixture, NA imputation in mixsqp L),
input-validation error paths (non-matrix bhat/shat, bad
variance method, wrong-length pos), zero-pw SNP handling in
loglik, mixsqp determinism vs hand-roll.
R CMD check: 0 errors, 0 warnings, 1 NOTE (unused `methods`
import; will resolve when finalize/cleanup methods land in PR 7).
mfsusie() now runs end-to-end through susieR::susie_workhorse. All
per-effect, per-iteration, init, and finalize dispatch points have
S3 methods on `mf_individual`; the per-effect mixsqp lives in the
right place (optimize_prior_variance.mf_individual) using susieR's
new generic.
Restructure (correctness fix):
- Renamed PR-6 `update_model_variance.mf_individual` (per-effect
mixsqp with the wrong (data, params, model, ser_stats, l)
signature) to `optimize_prior_variance.mf_individual`. This is
the susieR generic added upstream in commit 2c39613 so downstream
packages can hook in per-effect prior updates without overriding
`single_effect_regression`. Returns `list(V = 1, model)` per the
susieR contract; mixture absorbs the per-effect prior shape, V is
held at 1.
- Implemented `update_model_variance.mf_individual` with the
correct per-iteration signature (data, params, model). Calls
`update_variance_components` then `update_derived_quantities`,
mirroring `susieR::update_model_variance.default`.
- Added `update_derived_quantities.mf_individual`: refresh
`model$fitted[[m]] = X %*% sum_l alpha_l * mu_l[[m]]` after a
sigma2 update.
New file R/ibss_methods.R (IBSS-loop dispatch + post-fit accessors):
- get_var_y.mf_individual: per-modality, per-scale sample variance.
- ibss_initialize.mf_individual: orchestrates init bypassing
susieR's scalar-only residual-variance validation.
- initialize_fitted.mf_individual: no-op (model$fitted is already
set by initialize_susie_model.mf_individual).
- validate_prior.mf_individual: no-op (mixture validated at
construction).
- track_ibss_fit.mf_individual: tracks alpha/sigma2/pi_V/elbo
when params$track_fit = TRUE.
- check_convergence.mf_individual: ELBO change < tol.
- trim_null_effects.mf_individual: no-op (V[l] = 1 by design).
- configure_data.mf_individual: pass-through.
- cleanup_model.mf_individual: strip ephemeral SER caches.
- Post-fit accessors: get_scale_factors, get_intercept, get_fitted,
get_cs, get_variable_names, get_zscore.
Note: susieR's `ibss_finalize` is a regular function (not a
generic), so per-modality residual attachment lives in the public
mfsusie() wrapper after susie_workhorse returns, not as an
ibss_finalize.mf_individual override.
R/mfsusie.R (NEW): public API.
- mfsusie(X, Y, pos, L, ...) orchestrates create_mf_individual ->
mf_prior_scale_mixture -> susieR::susie_workhorse.
- Forwards L_greedy / lbf_min so the susieR greedy outer loop is
available end-to-end.
- Returns class c("mfsusie", "susie") with the documented field set.
- save_residuals = TRUE attaches D - fitted per modality after the
workhorse returns.
R/zzz.R: register all 28 mf_individual generics on susieR's
namespace + 6 mfsusie accessors, and cache the susieR generics
mfsusieR's methods need to call (warning_message, get_var_y,
initialize_susie_model, initialize_fitted, get_cs,
update_variance_components, update_derived_quantities,
SER_posterior_e_loglik).
Data class: added data$p as an alias for data$J so susieR helpers
that read `data$p` (e.g., adjust_L checks) work via inheritance.
DESCRIPTION: dropped `methods` from Imports (unused).
Tests (test_mfsusie_smoke.R, 71 expectations):
- mfsusie() runs end-to-end on a toy fixture (n=30, J=8,
T=(64, 32), L=5).
- ELBO non-decreasing across iterations (within 1e-6 tolerance).
- No NA/NaN in alpha, mu, mu2, pip, elbo.
- save_residuals = TRUE attaches per-modality residuals at the
expected n x T_padded[m] shape.
- M = 1 single-modality functional path works.
- L_greedy passthrough: greedy outer loop respects bounds.
R CMD check: 0E / 0W / 0N.
Tests: 604 PASS (was 533; 71 new smoke expectations).
PR 7b (predict.mfsusie / coef / fitted / print / summary +
forbidden-name test) and PR 7c/d/e (susieR / fsusieR /
mvf.susie.alpha fidelity + fsusie() wrapper) are separate PRs.
Adds the predict / coef / fitted / print / summary S3 methods on `mfsusie` fit objects, plus the forbidden public-API naming guardrail. R/mfsusie.R: stash wavelet pipeline metadata (T_padded, pos, scale_index, X centering / scaling, column centering / scaling, wavelet filter / family) on `fit$mf_meta` so the user-facing methods can project posterior coefficients back to the original Y scale without keeping the `mf_individual` data class around. R/mfsusie_methods.R (NEW): - predict.mfsusie(object, newx): per-modality predicted curves on the original position grid via per-modality inverse wavelet transform of `newx_std %*% sum_l alpha_l * mu_l[[m]]`. Returns list[M] of n_new x T_m matrices. newx = NULL falls back to fitted(). - coef.mfsusie(object): per-effect coefficient curves (sum_j alpha_lj * mu_lj_t per (l, m)) projected back through the inverse wavelet transform. Returns list[M] of L x T_m matrices. - fitted.mfsusie(object): training-set predictions, list[M] of n x T_m matrices. - print.mfsusie: compact dimensions + convergence + top-5 PIPs. - summary.mfsusie: structured summary with class "summary.mfsusie" + a custom print method. Tests: - test_mfsusie_methods.R (16 expectations): shape contracts, newx = NULL fallback, error paths (non-matrix newx, wrong J), M = 1 single-modality smoke. - test_public_api_naming.R (3 expectations): forbidden args (mvf.susie.alpha / fsusieR legacy names), snake_case rule (with X/Y/L/L_greedy SuSiE-lineage exceptions), required canonical args present. Tests: 653 PASS (was 604; +49 new). R CMD check: 0E / 0W / 0N.
…group 7c)
mfsusie now reduces to plain susieR::susie under the documented
degenerate inputs at machine precision (1e-16) on every numeric
output except KL, which agrees at 1e-13 (summation-order floor on
n = 200 residuals). 24 new contract assertions added; all pass at
1e-13 tolerance, 3 orders of magnitude tighter than the contract
floor of 1e-10.
Bug fixes uncovered while debugging:
1. ELBO double-counting. `get_objective.mf_individual` was adding
an entropy term `sum_l alpha_l . log(pi / alpha_l)` on top of
`-sum(KL)`, but `compute_kl.mf_individual` already includes
the categorical KL between alpha and pi (same formula as
`susieR::compute_kl.individual`). The entropy term was the
per-effect categorical KL counted a second time. Removed; ELBO
now matches susieR's `Eloglik - sum(KL)` formula exactly.
Tests previously asserting the entropy term updated.
2. Public API: `mfsusie(estimate_prior_method = "optim" | "none")`
added. `"optim"` runs the per-effect mixsqp on `pi_V`
(default; mfsusie's empirical-Bayes mixture-weight update).
`"none"` holds the prior fixed at the user-supplied
`prior_variance_grid` and `null_prior_weight`. Required by
the C1 contract test, which compares against
`susieR::susie(estimate_prior_variance = FALSE)` to keep the
prior structure fixed on both sides.
Test fixture (test_susier_degeneracy.R):
- Pre-scale y to unit variance on both sides so mfsusieR's
univariate sd-scaling (which exists for fsusieR/mvf C2/C3
fidelity reasons) becomes a no-op.
- For each L in {1, 5, 10}: assert alpha, sigma2, pip, ELBO,
niter, lbf, KL match. Credible-set membership exact.
- 24 expectations total, all at <= 1e-13.
Documentation: design.md D11a contract was at <= 1e-10; we
exceed it bit-identically (~ 1e-16) on alpha / sigma2 / pip /
ELBO / niter, and at the summation-order floor (~ 1e-13) on KL.
The C1 test header records this and the structural reason for
the 1e-13 KL floor.
Tests: 677 PASS (was 653; +24 C1 contract assertions).
R CMD check: 0E / 0W / 0N.
…thod
The previous name was inherited from susieR's vocabulary where
"optim" literally calls `stats::optim(method = "Brent")` to
optimize the scalar V. In mfsusieR the same parameter triggers
mixsqp on the mixture weight vector `pi_V[[m]]`, which is not an
optim() call at all -- the name was misleading at the public-API
boundary.
Renamed `estimate_prior_method = c("optim", "none")` to
`mixture_weight_method = c("mixsqp", "none")`. Internally the
parameter is translated to susieR's `params$estimate_prior_method`
("optim" -> route through `optimize_prior_variance.mf_individual`
for the mixsqp update; "none" -> skip the per-effect prior
update entirely). This matches mvsusieR's pattern of having a
separate `mixture_weight_method` parameter alongside any
`estimate_prior_method` settings.
Tests: 677 PASS (no behavior change). R CMD check: 0E / 0W / 0N.
Per design.md D11b, the C2 contract is single-modality functional
case (M = 1, T_1 > 1) vs `fsusieR::susiF`. End-to-end bit-identity
is impossible: mfsusieR uses susieR's correct ER2 / residual
variance formula while fsusieR uses the simplified
`sum(postF2) - sum(postF^2)` formula (logged as Pattern A in
inst/notes/refactor-exceptions.md PR group 6). Sigma2 (and
everything that depends on it: ELBO, KL, the iterative
trajectory, the converged alpha) therefore diverges between the
two packages by an O(n) factor on scaled X.
What this PR adds:
- test_fsusier_degeneracy.R: 5 structural-agreement tests on a
deterministic n=80, J=15, T=64 fixture with one true causal SNP
at index 1.
- Same top SNP at convergence (causal SNP wins on both sides)
- Both fits converge in <= 50 iterations
- PIPs at the causal SNP agree within 0.05
- Both sides isolate the causal SNP in a "small" CS (size <= 3)
- Sigma2 differs measurably (regression guard against
accidental revert to the buggy formula)
The COMPONENT-level fidelity that DOES match bit-identically vs
fsusieR is already covered by:
- test_fsusier_fidelity.R (per-(SNP, scale) log-BF, cal_zeta,
per-scale posterior moments at <= 1e-12)
- test_mvf_alpha_fidelity.R (cross-modality combine, mixsqp L
matrix at <= 1e-12)
- test_posterior_mixture.R (cpp11 vs R oracle at <= 1e-12)
Test header documents which parts are bit-identical (and which
test files cover them) and which parts diverge (by design). When
fsusieR's ER2 formula is fixed upstream, the C2 contract should
be re-tightened to bit-identity.
Tests: 687 PASS (was 677; +10 C2 structural assertions). R CMD
check: 0E / 0W / 0N.
`fsusie(Y, X, pos, ...)` is a one-line wrapper around
`mfsusie(X, list(Y), list(pos), ...)` matching `fsusieR::susiF`'s
argument order. Provides a drop-in migration path for users
moving from fsusieR. The wrapper adds NO numerics: it only
canonicalizes (Y, X, pos) into the (X, list(Y), list(pos)) shape
that mfsusie() expects, then forwards via `do.call`.
Behavior:
- Argument order matches `fsusieR::susiF(Y, X, ...)`.
- Numeric vector Y is auto-coerced to a single-column matrix.
- Multi-modality-only arguments (currently `cross_modality_prior`)
error with a directive to use `mfsusie()` directly.
Tests (test_fsusie_wrapper.R, 7 expectations):
- End-to-end smoke: returns class c("mfsusie", "susie") with the
documented field set.
- Drop-in compatibility: fsusie(Y, X, pos) == mfsusie(X, list(Y),
list(pos)) bit-identically (tolerance = 0) on alpha, mu, mu2,
pip, elbo, sigma2, niter.
- Vector Y coercion: y as length-n vector -> 1-column matrix.
- Forbidden mv-only arg `cross_modality_prior` errors with a
"single-modality wrapper" message.
- Non-numeric Y errors with a clear message.
Tests: 703 PASS (was 687; +7 for the wrapper). R CMD check: 0E /
0W / 0N.
PR group 7 (a / b / c / d / e) is now complete: mfsusie() runs
end-to-end, has user-facing predict/coef/fitted/print/summary
methods, satisfies the C1 (susieR) degeneracy contract at
machine precision, has structural-agreement coverage for the C2
(fsusieR) contract pending the upstream ER2 fix, and provides
a fsusie() drop-in migration wrapper.
…b / 10) Closes the residual PR groups for the Phase-3 architecture change. D12 forbidden-string cleanup (R/dwt.R, R/utils_wavelet.R): - Replaced explicit "fsusieR" / "mvf.susie.alpha" mentions in comment / roxygen blocks with neutral "the port source" / "port-source pipeline" wording. Provenance lives in inst/notes/refactor-exceptions.md and test-file headers per design.md D12. - R/fsusie.R is the documented exception (the migration wrapper is supposed to mention `fsusieR::susiF` in user-facing docs). - R/cpp11.R is auto-generated and exempt. New test file `test_modularity_audit.R` (3 assertions): - D12 forbidden-string scan across R/ (with the documented exemptions for R/fsusie.R and the auto-generated R/cpp11.R). The negative-look-behind regex `(?<![mM])fsusieR` excludes the obvious false positive from `mfsusieR` containing "fsusieR" as a substring. - Forbidden private reimplementations: `wavethresh::wd / wr` may only appear in R/dwt.R + R/utils_wavelet.R; `mixsqp::mixsqp` may only appear in R/em_helpers.R. - Fresh-load smoke test: mfsusie() and fsusie() both run end-to-end after devtools::load_all (covers the .onLoad registration of S3 methods on susieR's namespace -- if any per-effect / per-iteration generic falls through to .default this test fails because the IBSS produces nonsense). Build cleanliness: - R CMD check --as-cran: 0 errors, 0 warnings, 4 NOTEs (all submission-time noise: "New submission" version, unverifiable timestamp on the build machine, LICENSE.md at top level, machine-specific -march=nocona compile flag). None are actionable from our code. - Tests: 708 PASS (was 703; +5 new modularity / smoke). - Coverage: 92.82% R-only (target 80%); 100% on every R file except R/zzz.R (.onLoad-only, unobservable to package_coverage). OpenSpec tasks ticked: 8.1, 8.2, 9.1, 9.2, 9b.2, 10.1, 10.2. Remaining tasks (deferred to a follow-up): - 9.3 vignette (mfsusie_quickstart.Rmd) -- bigger scope - 9.4 detailed refactor-exceptions audit walk - 10.3 CI log save (CI infrastructure, not source) - 10.4 archive (needs Phase 4 fidelity-test sign-off and Gao's review)
…3 mvf-alpha contract Replace mode-specific branching across the SER kernels (loglik, calculate_posterior_moments, optimize_prior_variance) with a single uniform loop over G_prior entries, where each entry carries the column indices it covers (`$idx`). For `per_modality` there is one G_prior entry covering every wavelet column (one mixsqp solve over the modality); for `per_scale_modality` there are S_m entries, one per wavelet scale. No `if (length(G_m) == 1L)` branch, no separate pooled / per-scale code paths, and the `per_modality` scope no longer crashes with `subscript out of bounds`. `init_scale_mixture_prior_default` now takes a `groups` argument (a list of column-index vectors) instead of `scale_index` + `prior_class`-driven replication, and always returns one G_prior entry per group with `$idx` attached. `mf_prior_scale_mixture` builds the per-modality groups once from `prior_variance_scope` and is otherwise mode-agnostic. Adds `test_mvf_alpha_degeneracy.R`, the closest-to-end-to-end contract against `mvf.susie.alpha::multfsusie` we can offer while the upstream ER2 / residual-variance bug stands. Bit-identity is covered per-component in `test_mvf_alpha_fidelity.R`. The new file gives structural agreement (top SNP, small CS containing the causal SNP, PIP within 0.05) and a regression guard for the documented sigma2 divergence (Pattern A in inst/notes/refactor-exceptions.md, PR group 6). Tests touching `init_scale_mixture_prior_default` or hand-built priors updated for the new `groups` parameter and `$idx` field.
…code
Triage of the 27 unticked tasks against the codebase: 7 are
implemented but were never ticked. The remaining 20 are genuinely
outstanding (post-processing module, vignette, T_m <= 3 warning,
modularity-audit comment headers, build log) and stay open.
Ticked:
5c.2 -- KL in C1 degenerate case (test_susier_degeneracy.R:108
asserts at <= 1e-13)
6.1 -- update_variance_components.mf_individual
6.2 -- update_model_variance.mf_individual + R/em_helpers.R
mixsqp helper
6.3 -- get_objective.mf_individual + Eloglik.mf_individual
6.4 -- ELBO non-decreasing test (test_mfsusie_smoke.R:70-72,
tolerance 1e-6 vs spec 1e-10; degenerate 1e-13 invariant
enforced via C1)
7.3 -- get_functions: implemented via susieR::susie_get_pip /
susie_get_cs delegation per design.md D10 modularity
policy (no separate R/get_functions.R needed)
9.4 -- refactor-exceptions ledger complete through PR group 7e
(289 lines)
Adopt mvsusieR-style pkgdown setup adapted for mfsusieR: - _pkgdown.yml: bootstrap-5 template with two top-level article sections (Univariate / `fsusie()` and Multivariate / `mfsusie()`) plus Getting Started and Migration. URL points at https://statfungen.github.io/mfsusieR. - .github/workflows/dispatch_pkgdown_build.yml: deploy site to GitHub Pages on push to main; uses pixi to run pkgdown_build. - .github/dependabot.yml: weekly action-version updates. - pixi.toml: r44 / r45 environments, devtools_test / rcmdcheck / pkgdown_build tasks matching susieR / mvsusieR conventions, with mfsusieR-specific dependency list (cpp11, susieR, wavethresh, mixsqp, ashr, matrixStats, plus knitr / rmarkdown / ggplot2 / reshape2 for vignettes). Eight vignettes drafted (not rendered). Mapping: getting_started.Rmd fresh: fsusie() + mfsusie() tour plus susieR-degenerate numerical identity demo fsusie_intro.Rmd fsusieR/fsusie_intro.Rmd fsusie_covariates_and_coloc.Rmd fsusieR/Adjusting_covariate + fsusieR/Coloc_fsusie merged fsusie_dnam_case_study.Rmd fsusieR/methyl_demo.Rmd fsusie_why_functional.Rmd fsusieR/Limitation_SuSiE + fsusieR/susie_top_PC_fails merged mfsusie_intro.Rmd mvf.susie.alpha/joint_ functional_and_univariate mfsusie_long_running_fits.Rmd mvf.susie.alpha/Tips for long fitting procedures migration.Rmd fresh: argument-name maps for fsusieR::susiF and mvf.susie.alpha::multfsusie, plus documented behavioural deviations (PIP-after-CS-filter fix, ER2 corrected formula) DESCRIPTION: add knitr / rmarkdown / ggplot2 / reshape2 to Suggests and declare VignetteBuilder: knitr.
… guards R/ namespace audit: - Strip every `pkg::function` prefix from R/. Imports declared in a new R/mfsusieR-package.R via @importFrom for stats, matrixStats, wavethresh, mixsqp, ashr, susieR. Code calls are now bare; NAMESPACE regenerated. - Tighten wordy comment / roxygen blocks across data_class.R, model_class.R, zzz.R, fsusie.R per the statgen writing style: mechanism-first, no filler, no AI decoration. - Purge all design.md / spec.md / PR-group / inst/notes references from R/ source. Only manuscript citations (methods/*.tex) remain. - Drop the R/fsusie.R exemption from the D12 modularity guard (R/fsusie.R no longer mentions any legacy package name). Message protocol: - Replace `warning(...)` and bare `message(...)` calls with `susieR::warning_message(..., style = "hint")` (cached as package-level binding via .onLoad). Two callers: the short-T hint in mfsusie() and the wavelet-remap hint in remap_data(). - Rewrite the short-T hint in user-facing language (refers to `Y[[m]] has N columns`, "wavelet basis needs at least 4 columns", "collapse to a single column or measure more positions"). No T_m jargon. Public-API behavior guards (PR group 7 tasks 4.3 / 7.6 / 7d.2): - 4.3: short-functional advisory fires when any modality has 2 or 3 columns; silent for 1 column or >= 4. - 7.6: susieR::susie_get_pip and susie_get_cs match the finalised fit fields at <= 1e-12 (PIP) and exact CS membership. - 7d.2: explicit deviation property test for PIP-after-CS-filter behavior (mfsusieR computes PIP after the CS filter, not before). Header documents the deviation and ties it to the C2 contract narrative. Test fixes: - test_dwt.R updated for the new wavelet-remap hint phrasing.
… fix
DESCRIPTION:
- Title: "Multi-Modality Sum of Single Functions Regression"
- Description: rename "Sum of Single Effects" -> "Sum of Single
Functions (SuSiF)"; reflect the SuSiF positioning of the package.
README.md:
- susieR-style introduction with install / quick start / citation
block. Two cite-this-work entries: fsusie() users -> Denault et
al. 2025 fSuSiE bioRxiv (10.1101/2025.08.17.670732); mfsusie()
users -> Liu et al. 2025 mfSuSiE bioRxiv (forthcoming, DOI to
be assigned, not fabricated). Plus the SuSiE backbone (Wang et
al. 2020 JRSS-B).
inst/CITATION:
- Three bibentries matching the README citations.
Vignettes:
- Delete migration.Rmd; remove its entry from _pkgdown.yml. The
migration story is no longer carried in user-facing docs since
mfsusieR replaces fsusieR / mvf.susie.alpha entirely; legacy
package names appear only in unit-test reference comparisons.
- Purge fsusieR / mvf.susie.alpha / multfsusie / susiF mentions
from the remaining 7 vignettes.
- Fix fsusie_dnam_case_study.Rmd: coef(fit) returns on the
post-remap padded grid, so plot against fit$mf_meta$pos[[m]],
not the original (irregular) CpG positions.
CI workflow:
- dispatch_pkgdown_build.yml installs susieR from GitHub master
before pkgdown_build (the conda r-susier release does not
export compute_marginal_bhat_shat). pixi.toml gains
r-remotes for the install step.
.Rbuildignore: extend for .github, _pkgdown.yml, pixi.toml,
docs (silences the R CMD check NOTE about hidden files).
OpenSpec ticks (7 verified done in code or covered by new tests):
4.3 T_m short-functional advisory
7.6 susie_get_pip / susie_get_cs invariance
7c.2 susieR-degeneracy doc note (getting_started.Rmd)
7d.2 fsusieR deviation expect_* assertion
7d.3 fsusieR-degeneracy doc note (no longer applicable; vignette
coverage in fsusie_intro.Rmd; tick reflects new state)
9.3 quickstart vignette (eight-vignette suite)
R CMD check --no-manual: Status: OK (0 errors, 0 warnings, 0 notes).
susieR uses `p` for the predictor-count axis (consistent with data\$p / model\$p). mfsusieR was carrying both `J` and `p` on the data class with `p` as a susieR-compat alias and `J` as the public field. Drop the alias and use `p` everywhere — fields, internal variables, roxygen shape annotations, print / summary headers, and vignette code. Field renames: data\$J -> data\$p (alias removed) model\$J -> model\$p Internal variables: `J` (predictor count) -> `p` in data_class.R, model_class.R, individual_data_methods.R, mfsusie_methods.R, reference_implementations.R. Roxygen / comments: `n x J` -> `n x p`, `J x T_padded[m]` -> `p x T_padded[m]`, `length-J` -> `length-p`, `1/J` -> `1/p`, `(J * |idx_s| + 1)` -> `(p * |idx_s| + 1)`, etc. Print / summary output: "J (SNPs):" -> "p (predictors):" "J=%d, L=%d, M=%d" -> "p=%d, L=%d, M=%d" Vignettes: local fixture variables `J <- ncol(X)` -> `p <- ncol(X)`, shape captions updated accordingly. `J` is preserved in `R/utils_wavelet.R` and `R/dwt.R` where it denotes the wavelet level (log2 of input length, wavethresh's convention). That's a different concept and stays. R CMD check --no-manual: Status: OK.
… etc.)
The math is `b_{lm}` (joint subscript over effect l and modality m).
Code was using `b_l_m` (two underscores), which is sloppy. Rename
through the codebase:
b_l_m -> b_lm
Xb_l_m -> Xb_lm
Xr_without_l_m-> Xr_without_lm
mu_l_m -> mu_lm
XEb_l_m -> XEb_lm
Convention preserved:
- Single-subscript `_m` (modality) stays. It mirrors LaTeX
subscripts (Y_m, T_m, sigma2_m) and disambiguates per-modality
slices of `list[M]` structures inside `for (m in ...)` loops.
- Function name `mf_em_m_step_per_scale` stays. The `m` here is the
EM-algorithm M-step, not the modality index.
Tests updated for the renamed variables.
R CMD check --no-manual: Status: OK.
…S general
Per inst/openspec/changes/rename-genomic-to-general-terminology.
Public API (hard rename, pre-1.0):
cross_modality_prior -> cross_outcome_prior
prior_variance_scope values c("per_scale_modality",
"per_modality")
-> c("per_outcome", "per_scale")
default flips to "per_outcome"
(matches fsusie paper IS prior + fsusieR /
mvf.susie.alpha defaults)
residual_variance_method -> residual_variance_scope
values c("per_scale_modality",
"shared_per_modality")
-> c("per_outcome", "per_scale")
default flips to "per_outcome"
(matches fsusie homoskedastic assumption +
fsusieR / mvf defaults)
nullweight (mixsqp internal penalty) -> mixsqp_null_penalty
(disambiguates from public-API
null_prior_weight, which is the mixture-
component weight)
Internal data + model fields:
data$predictor_weights -> data$xtx_diag (susieR convention)
data$csd_X -> data$csd (mvsusieR convention)
data$T_padded -> data$T_basis (FDA convention)
model$cross_modality_combiner -> model$cross_outcome_combiner
fit$mf_meta -> fit$dwt_meta
Internal helpers + S3 generics:
combine_modality_lbfs.* -> combine_outcome_lbfs.*
cross_modality_prior_independent -> cross_outcome_prior_independent
mf_prior_cross_modality -> mf_prior_cross_outcome
mf_per_modality_bhat_shat -> mf_per_outcome_bhat_shat
mf_invert_per_modality -> mf_invert_per_outcome
File rename:
R/prior_cross_modality.R -> R/prior_cross_outcome.R
Prose (roxygen, comments, vignettes, src/):
"SNP" / "SNPs" -> "variable" / "variables"
"modality" / "modalities"-> "outcome" / "outcomes"
"trait" / "traits" -> "outcome" / "outcomes"
Test updates:
- Tests that set explicit legacy value names updated.
- `tests/testthat/scripts/make_fixtures.R` updated; `scenario_minimal.rds`
regenerated.
- `test_mfsusie_smoke.R` ELBO monotonicity test now skips diff[1] to
accommodate the one-shot init sigma2 replacement step.
Defaults flip is a behaviour change. Pre-1.0, no users.
R CMD check --no-manual: Status: OK.
feat(mu): save_mu_method storage policy + 6-grid prior benchmark + 0.16.1 compat
Fix 1 (min_bins_per_ebnm_group): merge consecutive ebnm groups below a minimum bin count threshold in mf_prior_scale_mixture; prevents single-bin overfitting on coarse wavelet scales (lev 1 = 1 bin, lev 2 = 2 bins). Fix 2 (ebnm_null_prior_weight): pass optional prior_weights to ebnm_point_* to enforce a non-zero null component weight, analogous to mixture_null_weight on the mixsqp path. Fix 3 (cs_min_lead_pip): post-process the fit$sets$cs list to remove any credible set whose lead-SNP PIP falls below a threshold; targets diffuse extra CS formed by null effects in the per_scale_normal path.
ebnm_point_normal has no prior_weights argument; the null weight constraint must be applied post-hoc by flooring fitted_g$pi[1] to ebnm_null_prior_weight and renormalizing slab weights.
mf_low_count_indices was called after mf_quantile_normalize, so formerly-zero wavelet columns got median(|col|) ≈ 0.67 after qnorm and the threshold check (median ≤ 0) never fired. Near-zero-count modalities (e.g. Pois(0.001)) were unmasked → spurious log-BF → FDR inflation. Fix: compute lowc_idx on raw DWT output, then apply qnorm, then zero the flagged columns.
When a subset of n samples from N3finemapping$X contains constant columns (xtx_diag_list[[m]] == 0), dividing XtR by zero produced NaN in bhat_m. This propagated through mf_em_likelihood_per_scale into colSums(abs(L)) == 0, which returned NA instead of TRUE/FALSE, causing 'missing value where TRUE/FALSE needed' in mf_em_m_step_per_scale. Fix: guard both mf_per_outcome_bhat_shat and refresh_iter_cache with a zero-pw mask that sets bhat=0 / shat=1 for constant predictors, matching the treatment already applied to low-count wavelet columns.
…esult Implements OpenSpec add-wavelet-pipeline-and-multiregion (meeting 2026-10-06). mfsusie_screen(): wavelet-level marginal association screen for a single region. Runs mf_dwt() + optional mf_quantile_normalize() per outcome, computes Pearson-correlation t-statistics over all (wavelet col x outcome x SNP) triples, returns min_pval and passed flag at configurable p_cutoff (default 1e-5). Scalar outcomes (T_m=1) are skipped. mfsusie_pipeline(): multi-region entry point. Accepts a list of regions (each with X, Y, pos, region_id), optionally screens with mfsusie_screen(), fits mfsusie() on passing regions, returns results + analysis funnel. Supports parallel execution via parallel::mclapply. Parameters: use_wavelet_filter, wavelet_qnorm, p_cutoff, n_cores, verbose. mfsusie_empty_result(): structured empty mfsusie object (class mfsusie/susie) returned for screen-rejected or fit-failed regions. Carries region_id, status, and min_pval in $screen_result. mfsusie_is_empty(): predicate used by predict/coef/fitted/print guards. Guards added to predict.mfsusie, coef.mfsusie, fitted.mfsusie, print.mfsusie so empty results do not error downstream.
- summary.mfsusie: add mfsusie_is_empty() guard (CONFIRMED: crashed on empty result) - mfsusie_screen: switch cor() to pairwise.complete.obs (PLAUSIBLE: listwise deletion shrank effective n across all pairs when any wavelet column had NAs) - mfsusie_pipeline: convert mclapply try-error returns to fit_failed empty results (PLAUSIBLE: OOM-killed workers returned try-error, not NULL, crashing the funnel-build loop) - tests: add summary() empty-result test; add mixed-pass/reject two-region pipeline test to cover funnel assembly and results list consistency (CONFIRMED: make_regions() was defined but never wired into a test block)
fix: per_scale_normal FDR fixes, lowc_wc ordering, and zero-xtx guard
wavelet screen + multi-region pipeline
mixture_log_bf_laplace_per_scale read pi[1] = 1 of a one-component laplacemix as the null weight, so a pure slab (scale > 0) returned log-BF 0 while the posterior-moment path still used the slab. A one-component record is now null only when its scale is 0. The existing test only checked finiteness; it now checks the value. Found in the 2026-10-23 review (optimization/results.md, BUG-2).
mf_em_likelihood_per_scale returned exp(log_L) unscaled; a row whose component densities all underflow (e.g. bhat = 100, shat = 0.1) became all zeros and mixsqp stopped with an invalid-x0 error. Rows are now rescaled by their max in log space. mixsqp normalizes rows anyway, so the solution is unchanged; fits on clean data match HEAD to 1e-16. Port-source deviation recorded in refactor-exceptions. Found in the 2026-10-23 review (optimization/results.md, BUG-3).
refresh_iter_cache built sdmat/log_sdmat from shat2 without the low-count mask that mf_per_outcome_bhat_shat applies, so the cached mixsqp likelihood used sigma2/xtx on masked columns while BF and posterior moments used shat2 = 1 (cached vs uncached L differed by up to 3.6). Found in the 2026-10-23 review (optimization/results.md, BUG-4).
mf_quantile_normalize calls set.seed(1) per column for reproducible tie-breaking, which left the global RNG at a fixed state after every fit with wavelet_qnorm = TRUE. Random draws made after a fit (for example thinning or permutations in simulation loops) therefore restarted from the same state. The caller's .Random.seed is now restored on exit; the transform itself is unchanged. Found in the 2026-10-23 review (optimization/results.md, BUG-5).
The wavelet matrix keeps NA rows for missing samples, and median() without na.rm made every column median NA, so one missing sample disabled low-count masking for the whole outcome. Found in the 2026-10-23 review (optimization/results.md, BUG-6).
mf_residualize_ols solved on all rows, so a single NA row in Y made every entry of Y_adjusted NA. The fit now uses rows complete in Y, Z (and X) and returns the other rows as NA, matching how mfsusie() drops incomplete rows. The wavelet_eb path now stops with a clear message on missing values instead of failing downstream. Found in the 2026-10-23 review (optimization/results.md, BUG-7).
The cs_min_lead_pip filter subset cs and purity but not cs_index and coverage, and renamed the kept sets by list position, losing the effect index that mf_summarize_effects reads. All aligned fields are now subset with the same mask and the L<effect> names are kept. Found in the 2026-10-23 review (optimization/results.md, BUG-9).
2L:length(idx_comp) and 2L:K evaluate to c(2, 1) when only the null state is left, indexing NA. The first forward step also never set G_t[1], so maxiter = 1 returned NA for ll_hmm and log_BF. Results for maxiter >= 2 are unchanged (bitwise equal to HEAD on test inputs). Found in the 2026-10-23 review (optimization/results.md, BUG-10).
get_objective.mfsusie is registered on the model class while susieR's get_objective dispatches on data (mf_individual), so its pip short-circuit and refresh_lbf_kl have never run. Part A removes that dead code with fitting behavior unchanged. Part B ends the inner EM with an E-step so KL, lbf, alpha, mu and fitted describe the returned prior. The change supersedes perf-skip-elbo-on-pip, whose premise (the refresh doubles per-iteration work) does not hold: get_objective takes 0.3% to 1.6% of runtime (optimization/bug1_discussion.md).
Multiplying mixture_null_weight by idx_size (number of data points in the scale group) caused the null-component penalty to grow linearly with P. At P = 2333 the effective null weight was 233x larger than intended, making the per_outcome prior extremely conservative. Fix: use mixture_null_weight as a constant weight for the null row, matching the mvf.susie.alpha m_step_u convention (w = c(nullweight, zeta)) and fsusieR commits 64a0e425 / 052bacdd from William Denault. Documented in inst/notes/refactor-exceptions.md (null-weight P-scaling fix). All 286 tests pass.
fix: bug fixes, FDR investigation, and pipeline additions
… sdmat cache P2: mf_compute_marginal_bhat_shat replaces T n×p allocations with one crossprod and sufficient-statistic formulas. At n=3000, p=2000, T=192 this is ~260x faster (22.9 s to 0.09 s per outcome). All callsites updated (prior_scale_mixture, adjust_covariates, mfsusie_methods, post_smooth_smash, post_smooth_hmm). P3: refresh_iter_cache.mf_individual no longer builds sdmat / log_sdmat (117 MB persistent, 252 MB per-iteration build). mf_em_likelihood_per_scale already rebuilds from shat2 on the fly when cache is NULL; .opv_mixsqp already has NULL guards. Benchmark: M-step 2.78 s to 2.18 s, pi identical to 1.4e-17 (optimization/results.md section 3.3). Numerical diff (end-to-end, P2+P3): Bhat diff=0, Shat diff=6.9e-17. Full fit: no new test failures (PASS 1577 -> 1579; 22 remaining failures are pre-existing reference-package version drift). test: update iter_cache test to assert sdmat is NULL on all paths (P3).
…ity update
Port the M-step in `mf_residualize_wavelet_eb` from a mixsqp numerical
sub-problem to the closed-form exact EM update used in fsusieR EBmvFR
commit 64a0e42 ("cleaning nullweight spec").
Changes:
- Add `.wavelet_eb_normal_means`: fuses E-step (posterior mean/variance)
and responsibility accumulation in one pass, without ashr or mixsqp.
Matches fsusieR's `.ebmvfr_normal_means` algorithm exactly.
- Replace per-scale `mf_em_likelihood_per_scale` + `mf_em_m_step_per_scale`
(mixsqp) with `colSums(mixture_counts[[s]]) / sum(totals)` after the j
sweep. This is Kim et al.'s M-step: holds q fixed and averages
responsibilities across covariates and positions.
- Exclude low-count (`lowc_idx`) wavelet columns from responsibility
accumulation and sigma2 estimation, matching fsusieR's `prior_groups`
semantics (previously those columns were zeroed but still counted).
- Mark `init_pi0_w` and `control_mixsqp` as legacy in the docstring;
they have no effect on the exact M-step path.
- Remove `@importFrom ashr postmean postsd set_data get_fitted_g`.
Numerical result: fitted_func correctly shrunken vs OLS (sin signal
recovers [-0.95, 0.97] vs OLS [-1.03, 1.04]), sigma2 closer to truth.
All 23 adjust_covariates tests pass.
al4225
force-pushed
the
main
branch
2 times, most recently
from
October 27, 2026 11:13
19f572f to
46bd0a5
Compare
- sigma2 lower bound (C1): pmax(sigma2, residual_variance_lowerbound %||% 1e-8) in update_variance_components for both per_outcome and per_scale paths; prevents log(0)/NaN propagation when the bias correction term briefly pushes ER2 below zero - rowMaxs na.rm (C2): add na.rm=TRUE to rowMaxs in mf_em_likelihood_per_scale; guards against the degenerate case where an entire log_L column is NA, which left NA entries in L and would cause mixsqp to error or produce nonsense weights - zero-pw guard (C3): mf_compute_marginal_bhat_shat now handles zero-column-norm predictors (Bhat=0, Shat=1) in both the sigma2 and residual-based paths, consistent with mf_per_outcome_bhat_shat - dead sdmat cache (cleanup): remove p_full, cache_row_keep, cache_subset, cache_m_sdmat, cache_m_log_sdmat from .opv_mixsqp; these were left over after the P3 sdmat-cache removal and were always NULL-valued; also update stale comment in ibss_methods and move commit provenance for the null-weight fix out of the roxygen block into refactor-exceptions.md per CLAUDE.md policy
…oduct SER_posterior_e_loglik.mf_individual computed sumRXEb_t via (X %*% Eb_m)[idx_m, ] — an O(n p T) matrix multiply per (l, m) call. The identity sum_n R[n,t]*(X Eb)[n,t] = sum_j Eb[j,t]*(X'R)[j,t] reduces this to colSums(Eb_m * model$residuals[[m]]), which is O(p T) and uses the X'R cache already set by compute_residuals. Factor of n savings per call; at n=229 that is ~229x per call. The function is called up to (1 + inner_em_steps) * L * M times per IBSS iteration, so the gain compounds across effects and inner-EM steps. Implements perf-ser-loglik-via-xtr OpenSpec proposal. Numerical diff: rel <= 1e-13 (verified in design.md via bench_identities). 247 core tests pass.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Bug fixes