Smoothing effect curves and credible bands

Why post-process

mfsusie() and fsusie() fit the model in the wavelet domain. The per-effect posterior is summarized by coef(fit), which inverts the wavelet representation back to position space. The raw inverse curve is correct but visually wiggly: small fine-scale wavelet coefficients contribute pointwise noise that obscures the underlying shape, and the fit returns no pointwise credible band by default.

The post-processor mf_post_smooth(fit, method = ...) runs the chosen smoother on the per-effect wavelet posterior mean and stores the result on the fit at fit$smoothed[[method]]. The slot carries:

fit$smoothed[[method]]$effect_curves[[m]][[l]] — smoothed position-space curve for outcome m, effect l.
fit$smoothed[[method]]$credible_bands[[m]][[l]] — T_basis x 2 matrix [lower, upper].
fit$smoothed$HMM$lfsr_curves[[m]][[l]] — per-position local false sign rate (HMM only).

The fit accumulates per-method results: every call to mf_post_smooth(fit, method = ...) adds (or, with overwrite_previous = TRUE, replaces) one entry in fit$smoothed. The raw inverse-DWT curves remain accessible via coef(fit); pass coef(fit, smooth_method = "TI") to get the TI-smoothed curves on the same list[M] of L x T_basis matrix shape.

mfsusie_plot() reads from fit$smoothed and picks a smoother in the priority order TI > smash > HMM > scalewise when the user does not name one explicitly. When more than one smoother has been applied, the plot emits a hint listing the alternatives so the user can compare. After method = "HMM" smoothing the per-position lfsr curves are also overlaid on a secondary axis with a threshold line at lfsr_threshold = 0.05. Affected regions (where each CS’s credible band excludes zero) are marked with thick bars at the bottom of each effect panel.

A small example

library(mfsusieR)
library(susieR)
set.seed(1)

Simulate one functional outcome with a smooth bell-shaped effect at each of two causal variables:

data(N3finemapping)
X <- N3finemapping$X[, seq_len(150)]
n <- nrow(X); p <- ncol(X); T_m <- 64

positions <- seq_len(T_m)
shape <- exp(-((positions - T_m / 2)^2) / (2 * (T_m / 8)^2))
beta <- matrix(0, p, T_m)
beta[42, ] <-  1.5 * shape
beta[97, ] <- -1.0 * shape
Y <- X %*% beta + matrix(rnorm(n * T_m, sd = 0.5), n)

fit <- fsusie(Y, X, pos = positions, L = 15, L_greedy = 5, verbose = TRUE)
#> iter          ELBO       delta   sigma2      mem      V  extras
#>    1   -51400.7084           -   [0.998, 0.998, 0.999]   0.17 GB  [1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00]  pi_null=[1.00, 1.00]
#> iter   2: max|dPIP|=1.75e-12, V=[1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00] -- converged [mem: 0.17 GB]
#> [L_greedy] 1 round, greedy_lbf_cutoff=0.100, final L=5
#> round  L     min(lbf)   action
#> 1      5     0.000      saturated

Before smoothing

mfsusie_plot() on the unsmoothed fit shows curves recovered from coef() (raw inverse DWT). The curves track the true bell shape, but fine wavelet scales contribute high-frequency noise, and the fit returns no credible bands by default.

plot_dims_b <- mfsusie_plot_dimensions(fit)

mfsusie_plot(fit)

After smoothing

fit_s <- mf_post_smooth(fit)            # default method = "TI"
names(fit_s$smoothed)                   # which smoothings are on the fit
#> [1] "TI"

The same plot call now reads fit_s$smoothed$TI and adds 95% pointwise credible bands as transparent ribbons:

plot_dims_a <- mfsusie_plot_dimensions(fit_s)

mfsusie_plot(fit_s)

The shape is preserved, high-frequency noise is suppressed, and the ribbon marks the positions where the band excludes zero.

Stacking multiple smoothers on one fit

Every call to mf_post_smooth(fit, method = ...) adds an entry to fit$smoothed rather than overwriting the previous result. Comparing methods is a few lines:

fit_s <- mf_post_smooth(fit_s, method = "HMM")
fit_s <- mf_post_smooth(fit_s, method = "scalewise")
names(fit_s$smoothed)
#> [1] "TI"        "HMM"       "scalewise"

coef() exposes the same multi-method state via the smooth_method argument:

str(coef(fit_s),                          max.level = 1L)  # raw inverse DWT
#> List of 1
#>  $ : num [1:5, 1:64] 0.0429 -0.1049 0 0 0 ...
str(coef(fit_s, smooth_method = "TI"),    max.level = 1L)
#> List of 1
#>  $ : num [1:5, 1:64] 0.00601 0.01743 0.03055 0.03055 0.03055 ...
#>  - attr(*, "smooth_method")= chr "TI"
str(coef(fit_s, smooth_method = "HMM"),   max.level = 1L)
#> List of 1
#>  $ : num [1:5, 1:64] -5.47e-08 1.53e-06 -9.83e-04 -9.83e-04 -9.83e-04 ...
#>  - attr(*, "smooth_method")= chr "HMM"

Because coef(fit) defaults to the raw inverse DWT, smoothing is always an explicit opt-in. Pass smooth_method = "<name>" when you want the post-processed curves.

Choosing a smoother

mf_post_smooth(method = ...) supports four methods:

"TI" (default) — translation-invariant wavelet denoising via cycle spinning. Best when the per-position band is the visual focus.
"scalewise" — per-scale soft-thresholding of the wavelet posterior mean at threshold_factor * sd_s * sqrt(2 log T), followed by inverse DWT. Cheapest, no extra dependencies. Lowers fine-scale jitter while preserving the overall shape.
"HMM" — hidden Markov model on per-position regression estimates with ashr mixture refinement. Adds a per-position local false sign rate (lfsr_curves); mfsusie_plot() overlays it on a secondary axis with the threshold line.
"smash" — smashr::smash.gaus empirical-Bayes wavelet shrinkage on the per-position regression estimate. Requires the smashr package (Suggests).

All four smoothers can coexist on the same fit. mfsusie_plot() picks one by priority (TI > smash > HMM > scalewise) and emits a hint listing the others; pass smooth_method = "<name>" to plot a specific one.

fit_scalewise <- mf_post_smooth(fit, method = "scalewise",
                                threshold_factor = 0.5)
plot_dims_sc <- mfsusie_plot_dimensions(fit_scalewise,
                                        smooth_method = "scalewise")

mfsusie_plot(fit_scalewise, smooth_method = "scalewise")

fit_hmm <- mf_post_smooth(fit, method = "HMM")
plot_dims_hm <- mfsusie_plot_dimensions(fit_hmm)

mfsusie_plot(fit_hmm)   # picks HMM (the only smoother on the fit)

fit_smash <- mf_post_smooth(fit, method = "smash")
plot_dims_sm <- mfsusie_plot_dimensions(fit_smash)

mfsusie_plot(fit_smash)

When several smoothers are stacked, plotting picks the highest-priority one and lists the others as a hint:

fit_multi <- mf_post_smooth(fit, method = "TI")
fit_multi <- mf_post_smooth(fit_multi, method = "HMM")
fit_multi <- mf_post_smooth(fit_multi, method = "scalewise")
plot_dims_ml <- mfsusie_plot_dimensions(fit_multi)
#> HINT: Plotting smoothing 'TI'. Other smoothings on this fit: 'HMM', 'scalewise'. Pass `smooth_method = '<name>'` to plot a different one.

mfsusie_plot(fit_multi)                              # picks TI; hint prints
#> HINT: Plotting smoothing 'TI'. Other smoothings on this fit: 'HMM', 'scalewise'. Pass `smooth_method = '<name>'` to plot a different one.

mfsusie_plot(fit_multi, smooth_method = "HMM")       # plot HMM directly

Multi-outcome usage

mf_post_smooth() works on mfsusie() fits the same way; the chosen entry of fit$smoothed[[method]] carries effect_curves[[m]] and credible_bands[[m]] for every functional outcome. Scalar outcomes (T_m = 1) are passed through (no smoothing needed) and get a band derived from the wavelet posterior SD.

# Two functional outcomes at different sampling resolutions
# (T = 64 and T = 32). Same two causal variables, different
# per-outcome bell shapes.
T_per <- c(64L, 32L)
shape64 <- exp(-((seq_len(64) - 32)^2) / (2 * 8^2))
shape32 <- exp(-((seq_len(32) - 16)^2) / (2 * 4^2))
beta1 <- matrix(0, p, 64L); beta1[42, ] <- 1.5 * shape64
                            beta1[97, ] <- -1.0 * shape64
beta2 <- matrix(0, p, 32L); beta2[42, ] <- 1.0 * shape32
                            beta2[97, ] <- -0.6 * shape32
Y_list <- list(
  X %*% beta1 + matrix(rnorm(n * 64, sd = 0.5), n),
  X %*% beta2 + matrix(rnorm(n * 32, sd = 0.4), n)
)

fit_m   <- mfsusie(X, Y_list, L = 15, L_greedy = 5, verbose = TRUE)
#> iter          ELBO       delta   sigma2      mem      V  extras
#>    1   -77079.1854           -   [0.998, 0.998, 0.999]   0.18 GB  [1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00]  pi_null=[1.00, 1.00]
#> iter   2: max|dPIP|=1.60e-18, V=[1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00, 1.00e+00] -- converged [mem: 0.18 GB]
#> [L_greedy] 1 round, greedy_lbf_cutoff=0.100, final L=5
#> round  L     min(lbf)   action
#> 1      5     0.000      saturated
fit_m_s <- mf_post_smooth(fit_m, method = "TI")

plot_dims_ms <- mfsusie_plot_dimensions(fit_m_s)

mfsusie_plot(fit_m_s)

Method

Each method takes the per-effect wavelet-domain posterior (mean_w, var_w) = (alpha[l, ] %*% mu[[l]][[m]], alpha[l, ] %*% mu2[[l]][[m]] - mean_w^2) and produces a smoothed position-space curve plus a pointwise credible band:

"TI" — cycle spinning: at every shift, run a stationary DWT with squared filter coefficients to compute pointwise variance, apply ash to the wavelet coefficients per scale, and average the inverse reconstructions across shifts. Credible bands come from the average-of-basis variance kernel.
"scalewise" — soft-threshold mean_w[scale_index[[s]]] per wavelet scale s at threshold_factor * sd_s * sqrt(2 log T), then invert via mf_invert_dwt(). The pointwise variance is the position-space variance from the wavelet posterior: Var(pos[t]) = sum_k W^T_{tk}^2 * var_w[k] where W^T is the inverse-DWT matrix. The band is `mean_pos +/- z(level)
- sqrt(Var(pos))`.
"HMM" — discrete-state hidden Markov model on per-position regression estimates with hub-and-spoke transitions and ash refinement at each EM iteration. The smoothed curve is the posterior mean. The pointwise SD is computed via the law of total variance over the HMM state mixture: `Var(pos[t]) = sum_k prob[t, k] * (sd_kt^2 + mean_kt^2)
- (sum_k prob[t, k] * mean_kt)^2, whereprob[t, k]is the forward-backward state probability and(mean_kt, sd_kt)are the per-state ash posterior moments. Band asmean_pos +/- z(level) * sqrt(Var(pos))`.
"smash" — smashr::smash.gaus empirical-Bayes wavelet shrinkage on the per-position regression estimate of the effect against the lead variable. The credible band uses smash’s post.var = TRUE posterior variance.

Session info

This is the version of R and the packages that were used to generate these results.

sessionInfo()
#> R version 4.4.3 (2025-02-28)
#> Platform: x86_64-conda-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS/LAPACK: /home/runner/work/mfsusieR/mfsusieR/.pixi/envs/r44/lib/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] susieR_0.16.1  mfsusieR_0.0.2
#> 
#> loaded via a namespace (and not attached):
#>  [1] generics_0.1.4        sass_0.4.10           smashr_1.2-7         
#>  [4] bitops_1.0-9          ashr_2.2-63           lattice_0.22-9       
#>  [7] magrittr_2.0.5        digest_0.6.39         caTools_1.18.3       
#> [10] evaluate_1.0.5        grid_4.4.3            RColorBrewer_1.1-3   
#> [13] fastmap_1.2.0         plyr_1.8.9            jsonlite_2.0.0       
#> [16] Matrix_1.7-5          reshape_0.8.10        mixsqp_0.3-54        
#> [19] scales_1.4.0          truncnorm_1.0-9       invgamma_1.2         
#> [22] textshaping_1.0.5     jquerylib_0.1.4       cli_3.6.6            
#> [25] rlang_1.2.0           zigg_0.0.2            crayon_1.5.3         
#> [28] LaplacesDemon_16.1.8  cachem_1.1.0          yaml_2.3.12          
#> [31] otel_0.2.0            tools_4.4.3           SQUAREM_2026.1       
#> [34] parallel_4.4.3        dplyr_1.2.1           ggplot2_4.0.3        
#> [37] wavethresh_4.7.3      Rfast_2.1.5.2         vctrs_0.7.3          
#> [40] R6_2.6.1              matrixStats_1.5.0     lifecycle_1.0.5      
#> [43] fs_2.1.0              htmlwidgets_1.6.4     MASS_7.3-65          
#> [46] ragg_1.5.2            irlba_2.3.7           pkgconfig_2.0.3      
#> [49] desc_1.4.3            pillar_1.11.1         pkgdown_2.2.0        
#> [52] RcppParallel_5.1.11-2 bslib_0.10.0          gtable_0.3.6         
#> [55] data.table_1.17.8     glue_1.8.1            Rcpp_1.1.1-1.1       
#> [58] systemfonts_1.3.2     tidyselect_1.2.1      tibble_3.3.1         
#> [61] xfun_0.57             dichromat_2.0-0.1     knitr_1.51           
#> [64] farver_2.1.2          htmltools_0.5.9       rmarkdown_2.31       
#> [67] compiler_4.4.3        S7_0.2.2

Gao Wang, Anjing Liu and William Denault