Bayesian Image Regression with Soft-thresholded Conditional Autoregressive Prior
Authors: Yuliang Xu
Jian Kang
Venue: ICLR 2025 Poster (2025)
Keywords: Brain fMRI Image
Variational Inference
Scalar-on-Image Regression
Image-on-Scalar Regression
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We propose a novel bayesian prior for efficient and accurate brain imaging analysis with high-Dimensional fMRI data.

Abstract: In the analysis of brain functional MRI (fMRI) data using regression models, Bayesian methods are highly valued for their flexibility and ability to quantify uncertainty. However, these methods face computational challenges in high-dimensional settings typical of brain imaging, and the often pre-specified correlation structures may not accurately capture the true spatial relationships within the brain. To address these issues, we develop a general prior specifically designed for regression models with large-scale imaging data. We introduce the Soft-Thresholded Conditional AutoRegressive (ST-CAR) prior, which reduces instability to pre-fixed correlation structures and provides inclusion probabilities to account for the uncertainty in choosing active voxels in the brain. We apply the ST-CAR prior to scalar-on-image (SonI) and image-on-scalar (IonS) regression models—both critical in brain imaging studies—and develop efficient computational algorithms using variational inference (VI) and stochastic subsampling techniques. Simulation studies demonstrate that the ST-CAR prior outperforms existing methods in identifying active brain regions with complex correlation patterns, while our VI algorithms offer superior computational performance. We further validate our approach by applying the ST-CAR to working memory fMRI data from the Adolescent Brain Cognitive Development (ABCD) study, highlighting its effectiveness in practical brain imaging applications.

Published as a conference paper at ICLR 2025
IMAGE REGRESSION WITH SOFT- THRESHOLDED CONDITIONAL AUTOREGRESSIVE PRIOR
BAYESIAN
Yuliang Xu
Jian Kang
Department of Statistical Science Duke University yul lang. xu@duke. edu
Department of Biostatisties University of Michigan j jiankang@umi umi ch edu
ABSTRACT
In the analysis of brain functional MRI (fMRI) data using regression models, Bayesian methods are highly valued for their flexibility and ability to quantify uncertainty. However, these methods face computational challenges in high- dimensional settings typical of brain imaging, and the often pre-specified cor- relation structures may not accurately capture the true spatial relationships within the brain. To address these issues, develop general prior specifically de- signed for regression models with large-scale imaging data. We introduce the Soft- Tettshhreshlled Conditional AutoRegressive (ST-CAR) priof" which reduces instability to pre-fixed correlation structures and provides inclusion probabilities account for the uncertainty in choosing active voxels in the brain. We apply the ST-CAR prior to scalar-on-image (SonI) and image-on -nage-r (IonS) regres- sion models-bott critical in brain imaging studies-aae
develop efficient com- putational algorithms using variational inference (VI) and stochastic subsampling techniques. Simulation studies demonstrate that the ST-CAR prior outperforms existing methods in identifying active brain regions with complex correlation pat- terns, while our VI algorithms offer superior computational performance. We further validate our approach by applying the ST-CAR to working memory fMRI data from the Adolescent Brain Cognitive Development (ABCD) study, highlight- ing its effectiveness in practical brain imaging applications.
INTRODUCTION
Regression problems with high-dimensional components have wide-ranging applications. For indi- vidual i, let a random function Ms(s), s E Rd denote the high-dimensional component, and it can be time-series (d = 1), 2D images (d 2), in particular interest, the functional Magnetic Resonance Imaging (IMRI) data (d = 3) of the human brain. In the context of the Adolescent Brain Cognitive Development (ABCD) study (Casey et al. 2018). fMRI data reflects children's brain development. Two scientific problems of interest are (i) the impact of brain development on children's general cogmitive ability, formulated scalar-on- on-image (SonI) regression; and (ii) the impact of the parental education level on different areas of the children's 's brain development, for- mulated as an image-on-scalar (IonS) regression. However, the complex anatomical structure of the human brain, the sparse signal in fMRI data, and the computational burden associated with handling high-resolution 3D images all pose challenges to image regressions.
Recent Development in Image Regressions To address these issues, state-of-the-art methods mostly impose sparsity and spatial dependence on the functional parameters. edged in the neuroimaging community that the fMRI data signals assumed to be sparse and piecewise smooth functions (Kang et al.] 2018 Taylor & Worsley (2077; [Smthh Fahrmeir] 22077 Mao al. 2017). This includes frequentist penalization via total variation distance (Wang et al. 2017). Bayesian Ising priors to control sparsity, and Gaussian Markov Random Field (MRF) or other
well acknowl-
"To whom correspondence should be addressed: Jian Kang (jinnkang@umich.c.m
STCAR is available R package Github httppp:://gittuuu....../

---

Published as a conference paper at ICLR 2025
latent Gaussian Processes to control spatial-correlation (Goldsmith et al. 2014 Huang al. 2013
Li et al. 2015; Kang et al. 2018: Zeng et al., 2022). A common issue with these approaches is that the spatial smoothness is pre-determined, either through tuning penalty parameters imposing prior structures on MRF. The data-adaptive way of tuning these smoothness parameters resorts to cross-validation, but with high-resolution 3D fMRI data, cross-validation can be computationally expensive. The T-LoHo method (Lee et al., 2021) proposed a tree-based graph partition prior that can learn spatial correlation treating voxels in 2D space as nodes in a graph, thus avoiding pre- specified smoothness assumption, but this approach is still computationally expensive for 3D fMRI data. In IonS, the popular approach is to use low-rank mapping (Ramsay & Silverman, 2005 Reiss et al. 2010, Yu et al.) 202 Li et al. 2021: Zhu et al., 2014). Again, the choice of the low-rank mapping using basis functions requires fine-tuning to reflect the true correlation structure of the high-dimensional data. /Zhang et al. (2023a) proposed to learn the high-dimensional effect through neural networks. Although tuning-free, this only provides point estimation with no uncertainty quantification the active areas.
Varitional Inference (VI) Posterior sampling for high-dimensional large-scale data has been chal- lenging for traditional MCMC methods, especially for imaging applications. As discussed in [Beii ett al. (0177). VI approximates the posterior sampling problem by an optimization problem, to min- imize the Kullback-Leibler (KI.) divergence between the posterior density and the candidate den- sity function over a family of densities. This allows us to borrow optimization techniques such as stochastic optimization with subsampling, and develop scalable algorithms for massive imaging data. Under the mean- mean-field assumption, our proposed ST-CAR prior can be easily implemented us- ing Coordinate Ascent VI (CAVI). In addition, inspired by the mini-batch Langevin Dynamic (Wul et al.) 2022 algorithm, we also provide simple modification on CAVI using stochastic subsampling to compute the posterior gradient, which we refer to as Stochastic Subsampling VI (SSVI). One crit- icism of VI algorithms is that, although VI can provide accurate point estimation, its ability to provide inference and uncertainty quantification on the point estimates is less desirable compared to traditional MCMC algorithms such as the Gibbs sampler. However, by designing sparsity inducing prior such as the ST-CAR prior, we can at least get the Posterior Inclusion Probability (PIP) using Vls, which provides statistical confidence whether the effect on certain voxels is nonzero. Down- stream FDR control methods (Meyer et al. 2015) can be applied using PIP.
sparsity-iduucing
Motivated by the above challenges and computational algorithms, we propose the Soft-Thresholded Conditional AutoRegressive (ST-CAR) prior, which uses a latent variable with CAR structure that is not sensitive the pre-specified prior correlation and provides shrinkage. Through extensive simu- lation and real data examples using the ABCD data on both Sonl and lonS regressions, we find that ST-CAR with VI algorithms can achieve the best computational performance among all comparable existing methods, and provide top-tier estimation and selection accuracy with state-of-the-an meth- ods. In summary, the proposed ST-CAR prior and its VI algorithms are (i) computationally efficient; (ii) prone to misspecified prior correlation and adaptive to data-driven smoothness; (iii) suitable for active voxels selection and downstream FDR control using PIP; (iv) generalizable not only to Sonl and IonS regressions, but also to logistic regression with imaging predictor, or other generalized regression models.
2 ST-CAR PRIOR
2.1 GENERAL NOTATIONS
Let N(,, 0)) represent a normal distribution with mean i and variance o², For the index set (I....,P). let [-j] denote the set (1,...,p} \ {(j). For a square matrix A, let Amin(A) and Xmax ( (A) be the smallest and the largest eigenvalues of A respectively. Let C+((A) denote the half Cauchy distribution with density function f(x) = (2/TA)/(x > 0)/(1 + =²/A²). Let I2 denote the q- dimensional Euclidean space. I, is the q by dimensional identity matrix. IG b) stands for the inverse- inverse-gamma gamma distribution. We use (() to denote neighborhood index set for j-th component, and (VG) to denote the cardinality of (f).. For a vector a. diag {a} is used to denote the diagonal matrix with diagonal vector a. We use sg)( = T(x > 0) Txx < 0) to denote the sign of x..
2

---

Published as a conference paper at ICLR 2025
2.2 ST-CAR PRIOR
The soft-thresoolding operator, defined as (2)) :- {i { [[ > v) for any v 2. 0, is continuous function that shrinks small values below v towards zero. The soft-thresholded Gaussian Process (Kang et al. 20188] STGP) prior applies T, to a latent GP. Thus, STGP allows the estimation of functional parameters within the Reproducing Kernel Hilbert Space (RKHS) of the late GP kernel. As the GP kernel and the corresponding RKHS change, the smoothness of the functional parameter space also changes. Without knowing the true signal pattern, srecippii the GP kernel can harm STGP's ability to learn the functional parameters (shown in Figure 2 Motivated by this idea, we relax the prior assumption on the smoothness the functional parameters, and propose novel noise-relaxed version of the STGP, referred to the Soft-thresholded conditional auto-regressive (ST-CAR) prior.
Definition 1. A sparse, spatially correlated parameter B(s) on a fixed grid ,1y..., Sp follows the ST-CAR prior if
B(s;) ind ((Tolu)). 03),, j= i=1....,p
P3 (NGG) ((NGG) T2s), BNG) = Pj kEN()
Σ b3.
We use B ~ ST-CAR(V, B) to denote 33 = (3(81)..... ((5)))))..(((()))) follow the ST-CAR prior with thresholding parameter and the neighborhood matrix B where (B)j,k bj.k
Here, 3()) is the target spatially varying parameter. The prior mean of B(s;) is the soft-thresholded μj. where Hy latent parameter for location 8y. The prior mean of (4 is determined by the correlation coefficient Pj. the neighborhood set N((). and the neighborhood weights bj.a-
The variance parameter ag is not identifiable when the ST-CAR prior is applied to Sonl or IonS regression. To impose sparsity on 3, we use the annealing idea on ag. and let as decay to O as the iteration increases. This will force the spatially independent B3 to converge to the sparse and spatially correlated μ (Ju...., (kp) The decay rate of ag has an impact on the variable selection accuracy especially for low signal-to-noiss ratio data. A general rule of thumb is to set as relatively large the beginning to allow for more flexibility and decays to a small value at the end.
The correlation parameter Pj can either be pre-fixed all locations, or updated by
6Gp. 85 Ber(ps),
where ūj is pre-fixed, and the binary indicator dy help to adaptively determine whether the value of My is strongly correlated with its neighborhood mean. We find these two approaches of updating or fixing Pj have similar performances in variable selection accuracy, and the ST-CAR prior is not sensitive to the choice of p or the bandwidth [W(j).. However, taking the second approach to adap- tively update p can give us extra information on the correlation structure of the high-dimensional coefficient.
The spatially-correlattd correlated structure of the latent μ is imposed by the Conditional Auto- Regreeegeee (CAR) covariance structure (Gelfand & Vounatsou 2003). Define matrix B and diagonal matrix Do, as
(B) bj,k
Wj.k
(Dow)ss Tu.j
Wi+
wi+
where Wiek are the kj. index of symmetric matrix W, and Wy+ Ek= Wy,k is the row(column) summation. The matrix W represents the correlation structure, and in practice we set Wj,k exp (-d(s)..s)) s.)). exponentially negatively associated with the distance between s and 8y. In addition to the CAR structure, we set a bandwidth [W(i) on the number of components in- cluded in the neighborhood. N(j,, such that for each j. bj.k is nonzero only if Wj,k is within the first the joint density of μu takes the
[W())) largest values among (w form
Denote = (pi....,p)
1 - DD-1 (T diag (p) B)
f(1la)) x exp
3

---

Published as a conference paper at ICLR 2025
If we denote Ep 1 {I diag (p) B), it can be shown that n is symmetric, positive definite matrix when each pj E (Amin(13)
(B) Note that B is not a symmetric matrix in general. In the construction in equation [ Amax
and Amin (B) < Hence we choose Pj E [0, 1) for any j. and the joint density of guaranteed to be non-degenerative. One caveat of doing so that the value of μú and the neighborhood mean [N(()) only allowed to be either positively correlated or independent (pj 2 0), but the negative correlation is not taken into consideration. This constraint makes sense in brain imaging applications, because the true signal is assumed to be sparse and piecewise smooth, which excludes the where the signal across neighboring voxels has a sharp drop from positive to negative values. In general, the ST-CAR prior is suitable for the case where the positive and negative areas do not share boundaries.
The proposed ST-CAR prior enjoys good computational properties as it has a conditional conjugate posterior when applied to spatial parameter in a regression problem. The main challenge in updat- ing a thresholded parameter is that the thresholding function such as T is non-linear. But we can show that μy conditional on all other M[-j] and has a mixture of truncated normal distribution as its posterior.
Proposition 1. Within the ST-CAR prior, the posterior of can be expressed as a mixture of three trancated normal distributions.
<(u (, uu 0z)
Pst N Not+ (1g Vz) Pf. ((-N)((NN)) Vo) + Pj -(--0,,-(15.V))
The expression for P P j. 1r M, Vg. Vo can be found in the proof of Proposition I in the I
MH
Appendix. Proposition provides the elosed-form posterior density for the sparse-mean latent pa- rameter (1z) as mixture of truncated normal distributions in the ST-CAR prior. Due to the and the prior
normal mixture conjugacy, the Gibbs sampler can be directly applied. Proposition [ property (Bylm; = 0) = will allow us to approximate the posterior inclusion probability using 19) given data.
The proposed ST-CAR prior is general prior that can be applied to many nigh-dimmessinal regres- sion settings, where the coefficient is assumed be smooth and sparse across their spatial domain. SonI and lonS regressions in Section ].3are are two examples to illustrate the power of ST-CAR prior. Other potential applications include logistic regression with -igh-dimensional exposure and other types of generalized linear models.
ST ST-CAR prior
Dardon PrCcHEIRE
Recression
OF
M
Obsorved Image Data
Spatial conficiants
0.0
0.5
Active Signal
location
Figure 1: Illustration to use ST-CAR prior for regression models with imaging data.
2.3 APPLICATION TO IMAGE REGRESSION MODELS
The ST-CAR prior can be applied to various models with sparse and spatially varying functional parameters as shown in Figure [ This section demonstrates its advantage using Sonl and IonS models. We focus on the SonI model and defer IonS model details to the Appendix.
Scalar-on-Image (SonI) Regression Let Mi(8;) denote the fMRI signal intensity at location 8y for individual X; E 1S be vector-valued confounder variables. Let Y, denote the scalar-valued outcome for subject i I,, I....pp.
Y = Σ
((s,) + TX,, xxiter te 6 = N(,,0))
A=1
β ST-CAR(U, B)
(3)
4

---

Published as a conference paper at ICLR 2025
where B IR is the nigh-dimensional spatially-varying coefficient of interest, and 7 R is the vector- valued coefficient for the confounders X,
The remaining parameters C+ (1). For all of the Half- Cauchy parameters, C+ (1) is equivalent to o² (1/22, 1/ay), ay ~ (G 1). Because B ~ ST-CAR B) essentially assigns spatially independent prior / with latent spatial-correlated mean function, we can use singular value decomposition (SVD) on the design matrix M E R=xp to boost further the computation speed
assigned the following priors 7 o0, oy (1),
their equivalent conjugate form to update:
Image-on-Scalar (IonS) Regression In IonS regression, let the spatially varying outcome be the fMRI data (sy). The exposure of interest is denoted X,, and the confounders are denoted as C, E R The lonS model is as follows
M,(s;)= = a(s;) E Ex(s))Cax + m/(ss) + a.s. as iid N(0.OM)
Rm
a ~ ST-CAR(V, B),
(4)
with other prior specification follows: & iid GP(O, 1,..., Y)i GP(O, or оM (1),, og ~ C+))), on (1)... Here, we only assign ST-CAR to for selecting active regions for the exposure. We defer details of IonS to Appendix A.l
idd
VARIATIONAL INFERENCE ALGORITHMS
This section uses the Sonl equation as an example, and introduces variational inference algorithms: Coordinate Ascent Variational Inference (CAVI) and Stochastic subsampling version of variational inference (SS More computational details can be found in A.2 and A.3
The variational inference methods (CAVI, SSVI) are based on the mean-field assumption (Blei al. 2017). If we denote f = (B, Y n, Ty ..) as the collection of all parameters. The mean- field variational inference minimizes the evidence lower bound
min E (KL(q(0) p( M,X) s.l. q(0) =(
Coordinate Ascent Variational Inference (CAVI) The classic CAVI algorithm iteratively refines the approximated density g by updating each parameter in successive iterations by the following density approximation,
log 90 (3) x (gt-N()saav.. flogp(3 X, Y,
Because each prior in SonI has closed-form posterior density, directly apply this iterative approach.
One issue with the conventional CAVI is that although it give good point estimation an optimization algorithm, but cannot directly make inferences such as the credible interval, com- pared with MCMC sampling methods. The novelty in our proposed ST-CAR prior is that can use the mixing probability in Proposition [ as the uncertainty quantification measure for selecting significant regions, circumventing the requirement for credible interval based on MCMC samples, while leveraging the computational efficiency provided by CAVI. Proposition T gives the posterior probability of μú belonging to the positive group [2, oo), zero group [-v, v),, and negative group (-20, -v). When using CAVI, we can directly compute the Posterior Inclusion Probability (PIP) under density as (P) + P,)) in equation 2 as a measure of coefficient significance.
Stochastic Subsampling Variational Inference (SSVI) The idea of SSVI is similar to the Stochas- tic Gradient Langevin Dynamics (Welling & Teh 2011 SGLD), in the spirit that instead of using the entire sample, we take a subsample of data, indexed by 7 C {1,...,th n.), gradient updating algorithm. In the context of VI update of Bj. let 8t be the step size (-th iteration,
and apply stochastic
Details this derivation be found in the Appendix
5

---

Published as a conference paper at ICLR 2025
let be the prior density of Bj.
Ege [3,] 4 1) {B,} + St
V 11- log Σ P(YY M, X 0) g log π (3,)
We defer details of SSVI to Section A.3 SSVI particularry useful for large-scale problems such as IonS where the likelihood computation is expensive due to the high diheensinnii outcome, whereas the performances of SSVI and CAVI for Sonl are similar.
4 NUMERICAL EXAMPLES
In this section, we focus on simulation for Sonl equation 3 and defer simulation results on IonS qquation@eegressiinn A.4 Our primary goal is to compare the proposed prior ST-CAR with other existing methods using CAVI. Because CAVI offers a balance between computational efficiency and selection accuracy in Sonl setting. Empirical comparisons of the posterior computation algorithms including Gibbs sampler, SSVI, and CAVI are included in A.10
Simulation Scalar-on-image regression with CAVI
For SonI model equation 3 we compare ST-CAR with 3 other methods: (1) Soft-thresholding Gaus-
sian Process prior (STGP) Kang et al. 2018), (2) T-LoHo Lee al. (0221), (3) Elastic Net Zou &
Hastie 22005).
For the elastic net implemented in the glmnet package (Friedman al. 2010a). the mixing parameter is set to 0.5, and the penalty parameter À is chosen using tross-validation.
The STGP prior is based on soft-thresholding on the latent Gaussian Process. When (s)) ~ STGP(v, k). there exists a corresponding latent Gaussian Process (s)) ~ GP(O, such that ((s) = (3).. This method requires pre-specified kernel function K, and the posterior sam- pling algorithm is the Metropolis-adjusted Langevin algorithm (MALA). In this simulation, we use the exponential square kernel
<(s, s'; b) = (Os), = exp(-a(s2 + 5 - b(s - 8')3}
(5)
where a = 0.01, b = 10. The implementation is based on BIMA package F (Xu & Kang 2023)
Note that this implementation of STGP allows the users to specify different regions in the image and specify a region-wise independent kernel in order to speed up the computation in high dimensions and boost selection accuracy in each region. Hence for the simulation pattern shown in Figure 2 we evenly split the entire 2D region into 4 sub-regions, and use the modified exponential square kernel on each sub-region. The basis function is generated using Kang (022) with 10 degrees of Hermite polynomials for each sub-region. We use the elastic net result as the initial value for B, and run total of 10 iterations with the last 20% the converged MCMC sample. The thresholding parameter v is set to be 0.2. For the variable selection accuracy, we use the Posterior Inclusion
Probability (PIP) based on the MCMC sample of B, defined as PIP; = Ein T(3 + 0)/T for the location j with MCMC sample.
The T-LoHo method is designed for clustering nodes in graph models into finite discrete values, and it shows great performance for this purpose especially under low SNR. However, it is not well suited to capture continuous functional parameters. In addition, T-LoHo does not impose sparsity. We use the 95% credible interval to select active voxels. T-LoHo package does not provide the confounder coefficients estimation, hence we set true O. We use total of 50000 MCMC iterations and take the last 10000 as the converged sample.
For ST-CAR prior updated using the CAVI algorithm, value for B, and set the initial value for μu to be all 0. The thresholding parameter v is set to be the largest marginal value in estimated from ridge regression. This is because setting be large value can reduce false discoveries, and is still able to recover the true signal pattern even when starting from all 0 initial values. This algorithm much less sensitive to the thresholding parameter
the ridge regression result as the initial
http://iinieeea//tooooooo//o//aa *TLOHO httpt: htt.
T1.HO

---

Published as a conference paper at ICLR 2025
compared to STGP. The decay rate of o3 is set to be 0.5(1 + -)7 where represents the number of iterations. Since we use annealing on of instead of fully conjugate update, we can longer iteration, compute the difference of ( and
ELBO stopping rule. Instead, at the g(+++), defined as ( ) /p, to determine whether the optimization has converged. The tolerance is set to be 10-10, For the neighboring matrix B, we set the number of neighbors as and the correlation parameter ý is set 0.9. The variance parameter 0y fixed at 1 for the CAVI update.
Sonl:
Sonl:
Sonl:
Sonl:
STGP
ca
0O
-10
(a) Case 1. 71 600, 1600, o
1.
Sonl: beta
Soni: ST-CAR
Sonl: ElastiNNel
Sonl. T-LoHo
Sonl: STGP
Valle
10
V
2
0.
T
48-
40
10
20
0.0
(b) Case 5. 6000, 10000, of 5.
Figure 2: Sonl result illustration for all competing methods. The first figure in each row is the true signal.
Figure and 2b provide visual comparison under two simulation settings. The true 3 image is
designed to include several challenging patterns where the active area can decay smoothly to almost 0, to have complex correlation structure such as the M-shape on the top-left corner, and to include both positive and negative patterns. Case 1 (Figure 2a is the low-resolution and low SNR setting, and Case 5 (Figure 2b is the high-resolution and high SNR setting.
From this visual comparison, STGP is good at estimating smooth function patterns such as the bottom-left circle, but without further tuning the Gaussian Process kernel, estimating more complex patterns such the M-shape would be difficult. Here, STGP already considers the region partition. we were to generate a Gaussian Kernel over the entire support, the result would be even smoother (appear more blurry effects, but as the true signal decays smoothly towards O, T-LoHo can ignore some small non-zero effects, resulting in lower statistical power, as shown in the bottom-left circle in Figure 2b Elastic Net can identify the spatial pattern to a certain extent, but the power is limited as it does not leverage spatial information. Consequently, it may yield a noisy estimation in case 1. Even in case where the point estimation is favorable, Elastic Net can still introduce background noises. ST-CAR can estimate each pattern relatively well without specifying any region partition or tuning the correlation matrix and it automatically adapts to different signal patterns. Although some small effects such as the bottom tip of the T-shape can still be missed, ST-CAR provides the best overall performance compared to other priors across different settings without any tuning procedure.
sharp edges). T-LoHo
clustering method is good at grouping larger
Table I provides a detailed numerical comparison. The evaluation criteria for estimation accuracy include (i) Selection accuracy: false discovery rate (FDR), true positive rate (TPR), and overall ac- curacy (ACC); (ii) Point estimation: root mean squared error (RMSE); (iii) Goodness-of dictive mean squared eтror the outcome Y, using training and testing data (train and test pMSE). We also include the computational time comparison averaged over 100 replications. Because CAVI an optimization algorithm, we can set stopping rule, whereas for MCMC sampling algorithms
the pre-

---

Published as a conference paper at ICLR 2025
Table I: Numeric result for Sonl simulation, under 100 replications.
(a) Sonl: Comparison of estimation accuracy. The evaluation criteria for estimation includes false discovery rate (FDR), positive rate (TPR), overall accuracy (ACC), and root mean squared error (RMSE), all multiplied by 100. The evaluation criteria predictive performance includes training and testing predictive MSE, denoted as Train Test pMSE respectively
Case 1. n = 600, 1600, of
Case 2. nI 600. = 900, o.
ST-CAR
ElasNet
STGP
T-LoHo
ST-CAR
ElasNet
STGP
T-LoHo
FDR
9.52
9.54
10.42
4.41
FDR
9,40
1.82
9.77
3.92
TPR
97.47
44.85
95.37
97.32
TPR
99.81
98.45
94.06
99.34
ACC
98.05
90.87
97.60
98.62
ACC
98.45
99.50
97.63
98.65
RMSE
11.00
30.45
13.44
9.46
RMSE
5.51
7.12
13.00
6,56
Train PMSE
2.29
1.03
5.97
3.61
Train pMSE
1.10
0.33
3.40
1.55
Test pMSE
7.18
60.31
12.54
6.21
Test PMSE
2.02
3.04
6,36
2.24
Case 3.
1000.
1600, o.
Case 4.
600,
1600, o. =
ST-CAR
ElasNet
STGP
T-LoHo
ST-CAR
ElasNet
STGP
T-LoHo
FDR
9.39
0.39
9.89
1.10
FDR
9.64
9.42
9.83
7.09
TPR
100.00
99.08
97.93
99.79
TPR
90.06
37.11
92.00
94.19
ACC
98.42
99.80
98.05
99.76
ACC
97.02
89.82
97.25
97.09
RMSE
4.27
6.53
12.04
5.83
RMSE
14.65
33.31
14.72
12.51
Train p.MEE
1.08
0.48
6.47
1.92
Train PMSE
4.85
2.55
9.17
8.56
Test pMSE
2.03
3.79
9.92
2.78
Test pMSE
17.76
76.84
19.31
14.37
Case
6000,
10000, o.
Case 6.
6000, p 10000, o.
10
ST-CAR
ElasNet
STGP
T-LoHo
ST-CAR
ElasNet
STGP
T-LoHo
FDR
0.46
0.27
17.16
6.84
FDR
1.18
2.47
17.15
4.78
TPR
99.98
99.39
98.52
99.67
TPR
99.91
98.21
98.53
99.52
ACC
99.92
99.86
96.57
95.96
ACC
99.80
99.33
96.58
97.71
RMSE
3.73
6.73
13.78
5.25
RMSE
3.14
3.94
6.44
2.48
Train PMSE
5.48
3.12
80.90
12.78
Train PMSE
10.62
3.81
85.91
18.60
Test pMSE
9.37
22.94
86.88
15.54
Test PMSE
18.43
41.06
92.25
22.44
(b) Computation time for Sonl simulation, averaged 100 replications.
Computation time
Total time (seconds)
Number of iteratios per second
Case
ST-CAR
STGP
T-LoHo
ST-CAR
STGP
T-LoHo
Case 1, 600,p p= 1600. 05 =
103.0
503.0
306.5
11.4
208.0
262.6
Case 1000.p 1600 of =1 600,
900, oy
24.2
250.8
205.0
42.3
420.5
111.2
866.2
426.2
10.2
122.9
189.5
Case 4, 600.p 1600, a3 5
108.5
486.0
312.5
11.1
212.4
259.8
Cuse 5. 6000,p 10000,o0 o
8034.9
40658. 8
11141.1
0.2
2.6
7,3
Cuse 6. n 6000, 10000, of
10
7811.3
40839.1
11297.8
0.2
2.6
7.2
for STGP and T-LoHo, a lot more iterations are required. Hence we report both the total time and the number of iterations per second.
For the variable selection results for ST-CAR, Elastic Net, and STGP, we use tuning procedure find cutoff such that the FDR can be controlled below 10% within fixed tuning window. For STGP and ST-CAR, the PIP is used to control FDR. For elastic net, effect size of 3 is used to control FDR. For T-LoHo, the 95% is used.
Based on Table I we can see that ST-CAR has the lowest testing pMSE in relatively high SNR cases (Case 2,3.5). For Case 1 and 4 with relatively low SNR, ST-CAR has the second-best perfor- mance next to ToLoHo. For the computation time in Table Tb, ST-CAR has the best computational efficiency across all cases.
APPLICATION TO ABCD STUDY
In this section, we apply ST-CAR to analyze the Adolescent Brain Cognitive Development (ABCD) study release 1 data (Casey et al. 2018). The ABCD study is long-term study on the brain de- velopment of children in the United States. In this real data analysis, we use the 2-back 3mm task fMRI contrast data (Sripada al. 2020a). The preprocessing and the general cognitive score) follow from the previous study Sripada et al. (2020b). The scientific questions of interest are: (i) whether the brain signals in different regions have different impacts the children's IQ score
(g-
8

---

Published as a conference paper at ICLR 2025
(Sonl); (ii) whether parents with higher education degrees has impact on the children's 's cognitive ability development (lonS). For the task fMRI data, after preprocessing (Sripada et al.) 2020b) and removing subjects with missing covariates, we have p 47636 voxels and n 1861 subjects in total.
To answer (i) with Sonl model equation we use the children's IQ score the scalar outcome Y;, and use the task fMRI data as the high- high-dimensional predictor M[ys;), where S, stands for voxel locations in the brain. The confounders include parental education level (binary, if the parent has bachelor's ' degree or higher), age, gender, race and ethnicity (Asian, Black, Hispanic, Other, White), and household income (less than 50k, between 50k and 100k, greater than 100k). The coefficient of interest is 3 in equation 3 We expect / to be very sparse and have small effects, since the interpretation for 3(s,) = b is that one unit increase in the brain signal in location 83 is associated with b amount of change in the children's IQ score, and the range of the standardized IQ score is (-2.84, 3.26), a small range compared to large number of predictors p = 47636.
To answer (ii) with the IonS model equation 4 we use the task fMRI data as the outcome, and the parental education level as the exposure. The confounders include age, gender, race and ethnicity, and household income. For the IonS model, for &k, 7/ that are assigned GP priors. The interpretation for a(s,) = a in equation 4 is that parents with bachelor's 's degrees or higher are associated with a amount of change in the brain signal at location sj. Hence expect the effect size of to be relatively larger than that of (3.
In ST-CAR prior, the two most important tuning parameters the thresholding parameter v. and the initial value for oă which controls how close / is to the latent sparse μ.. In theory (Kang et al., 2018), the choice v does not have a huge impact as long as the initial values are close enough to the truth, the MCMC sampling algorithm can run long enough to fully explore the parameter space. Because we are using VI algorithms, it is important to start with a good initial value. Hence we perform sensitivity analysis select the best v and initial o3 in terms of the smallest testing pMSE. The entire data set is split into 70% training data and 30% testing data. Based on the sensitivity analysis results in Table A2 and Appendix Table A4 and A3 we choose v = 0.007, the initial value for of to be 10 10-5, bandwidth and decay rate y = 0.35 in the decay rate function of oî for SonI, and choose v = 0.005, the initial value for o to be 0.1, bandwidth 26 and decay rate ς = 0.45 for lonS. Although varying v. bandwidth and decay rate have little influence the results. Table A2 also reflects that our method has better testing pMSE compared to the competing methods Elastic Net, STGP, and MUA. Due the computational limitation of other methods, we chose not to run all competing methods on the real data.
We use CAVI SonI, which takes 1.7 hours to run, and SSVI lonS, which takes 7.3 hours to run. Due to the vast sparsity and low SNR in 3, the computational time of SonI is similar to STGP (1.6 hours). However, the IonS model with the SSVI algorithm shows huge computational improvement compared to STGP (85.9 hours).
We present the final data analysis result in both visual illustrations in Figure 3 and numeric values in Table Figure is a visualization of the positive significant voxels in SonI and IonS. The color range for the plots is between (a, b). where only voxels with values greater than a are shown, and voxels with values greater than b are shown in the brightest color. From Figure 3a due to the low SNR in SonI, both the effect size and PIP are small, and only a small amount of voxels with large effect size aligns with the mapping of PIP greater than 0.25, In comparison, a in IonS has larger effect size, and as shown in Figure 3b the large effect areas align well with the mapping of PIP greater than 0.98.
In Table 2 we show the region-level numeric result. Note that, although both Sonl and IonS have small amount of negative effects, they are very close to 0 compared to the positive effeett scale, hence we only report the positive effect here. From Table 2 for Sonl, Precuneus I. is the region with the largest positive effect, which means brain development in this region can have the most positive effect on the children's 's IQ score. This aligns with the previous study in Xu & Kang (2023) and scientific findings (Wallentin et al.) 2006) that Precuneus is related to memory tasks. For IonS, Frontal Mid region in both the left and right hemispheres have the largest positive effect, and have been shown to play a key role in the development of literacy (left Prontal_Mid) and numeracy (right Frontal Mid) in previous findings (El- & Schury 2020).

---

Published as a conference paper at ICLR 2025
-79
-66
(a) SonI: values of / (Red) w. color range [0.0005,0.001], and values of PIP (overlaying blue) w. color range [0.25, 0.5].
(b) IonS: values of a (Red) w. color range [0.05,0.1], I], and values of PIP (overlaying blue) w. color range [0.95,1].
Figure 3: Visual illustration of 3 in Sonl and a in lonS.
Table 2: Numeric result for the top 10 regions sorted by number of significant positive voxels in Sonl and lonS. For Sonl, sig count the number significant voxels (PIP, 2 0.25) in each region, pos_sig count is the number of significant voxels with ((s,) 2 0.0005, and pos sum is Ejes. (Bss))(((s,) > 0), the sum of positive effect for all voxels in region r. The IonS result has the same interpretation, except the cutoff for significant voxels is PIP, 2 0.95, and the cutoff for positive effect sis_si count is 0.05.
SonI
lonS region code
region
region code
sig
pos_sig count
pos
region
sig
pos_sig
pos
name
count
sum
name
count
count
sum
Precuneus
12
0.25
Parictal Inf
61
382
357
38.94
Temporal Sup R
12
9
0.16
Precuneus
67
377
312
37.26
Temporal Inf R
0.18
Precentral
305
293
33.55
Precuneus.R
0.14
Precuneus
68
316
285
36,24
Temporal.Int Inf
0.15
Frontal Mid_R
322
270
43.12
Occipital _Md__MidLL. I.
0.12
Frontal Mid_L
272
244
43.06
Parietal Inf T.
8
0.15
Supp Motor Area
215
205
23.82
Frontal Sup Orb R
6
6
0.08
Parietal Sup L
59
224
167
18.63
Frontal_Mid L
7
11
0.15
Temporal_Mid MMddR
86
175
154
27.18
Frontal_Mid Orbb..
10
5
0.08
Frontal.Sup Sup_LL
155
147
21.81
6 DISCUSSION AND CONCLUSION
In this work, we have proposed the ST-CAR prior, which is a general and flexible prior that could be applied to any regression problems with imaging components. Variational inference algorithms are proposed for the ST-CAR prior. Especially, we implemented the coordinate ascent variational inference (CAVI) as a baseline VI algorithm that can provide good estimation accuracy in low SNR settings, and we proposed novel stochastic subsampling variational inference (SSVI) algorithm that is more computationally efficient. We demonstrated the use of the ST-CAR prior in both scalar- on-image and image-on-sealle regression models. Through comparisons in numeric studies, we find proposed method has better performance in terms of estimation and computation, compared with existing methods such as T-LoHo, STGP, and SBIOS. The proposed method applied to the significant region to
ABCD study with task fMRI image data and identifies the left Precuneus the children's IQ development, and the development of the middle frontal gyrus as the significant region that can be most positively impacted by parental education level.
10

---

Published as a conference paper at ICLR 2025
REFERENCES
Tracy Packiam Alloway and Ross Alloway. Working memory: Is it the new iq? Nature Precedings, pp. 1-1, 2008.
Sunil Arya, David Mount, Samuel E. Kemp, and Gregory Jefferis, RANN: Fast Nearest Neighbour Search ( (Wraps ANN Library) Using L2 Metric, 2019. URL https::/cRAA, R org/package-RANy RANN
package version 2.6.1.
Yoav Benjamini and Yosef Hochberg. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal statistical society: series B (Methodological), 57(1):289-300 300, 1995.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: review for statisti- cians. J. Am. Stat. Assoc., 112(518):859-8888 April 2017.
BJ Casey, Tariq Cannonier, May Conley, Alexandra O Cohen, Deanna M Barch, Mary M Heitzeg, Mary Soules, Theresa Teslovich, Danielle V Dellarco, Hugh Garavan, al. The adolescent brain cognitive development (abed) study: imaging acquisition across 21 sites. Developmental cognitive neuroscience, 32:43 54, 2018.
Pavla Cermakova, Adam Chlapeëka, Zsófia Csajbók, Lenka Andrýsková, Milan Brázdil, and Klára Marečková. Parental education, cognition and functional connectivity of the salience network. Scientific Reports, 13(1):2761, 2023.
Dirk Eddelbuettel and Conrad Sanderson. Repparmadillo: Accelerating r with high-performance C++ linear algebra. Computational Statistics and Data Analysis, 71:1054-1063, March 2014. URL http:://ss> doi wwtpp/..gg/...111////1.11113
2013. 005
Rami El-Baba and Mark Schury. Neuroanatomy, frontal cortex. 2020.
Jerome Friedman, Trevor Hastie, and Rob Tibshirani. Regularization paths for generalized linear models via coordinate descent. Journal of statistical software, 33(1):1, 2010a.
Jerome Friedman, Robert Tibshirani, and Trevor Hastie. Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 33((): 2010b. doi: 10.18637/
vs3..003.
Alan E Gelfand and Penelope Vounatsou. Proper multivariate conditional autoregressive models for spatial data analysis. Biostatistics, 4(1):11-25 January 2003.
Jeff Goldsmith, Lei Huang, and Ciprian M Crainiceanu. Smooth Scalar-on-Image regression via spatial bayesian variable selection. J. Comput. Graph. Stat., 23(1):46-64 January 2014.
Olivia M Halabicky, Jennifer A Pinto-Martin, Peggy Compton, and Jianghong Liu. Low level lead exposure in early childhood and parental education on adolescent iq and working memory: a cohort study. Journal of exposure science & environmental epidemiology. 33(2):168 1166)66,66 2023.
Lei Huang, Jeff Goldsmith, Philip T Reiss, Daniel S Reich, and Ciprian M Crainiceanu. Bayesian scalar-on-image regression with application to association between intracranial DTI and cognitive outcomes. Neuroimage, 83:210-223 December 2013.
Jian Kang. BayesGPfit: Fast Bayesian Gaussian Process Regression Fitting, 2022. URL https: / R-BANN.-projecee wwwpee.trt....ccaaaaaaeeeeeeeeee
R package version 1.1.0.
Jian Kang, Brian Reich, and Ana-Maria Staicu. Scalar-on-image regression via the soft- thresholded gaussian process. Biometrika, (05())115--184, 2018.
Changwoo Lee, Zhao Tang Luo, and Huiyan Sang. T-loho: A bayesian regularization model for structured sparsity and smoothness on graphs. Advances in Neural Information Processing Sys- fems, 34:598 (.05996699.... 2021.
Fan Li, Tingting Zhang, Quanli Wang, Marlen Z Gonzalez, Erin L Maresh, and James A Coan. Spatial bayesian variable selection and grouping for high- nigh-dimensional scalar-on one--oo-iiin -miae regression. 2015.
11

---

Published as a conference paper at ICLR 2025
Xinyi Li, Li Wang. Huixia Judy Wang, and Alzheimer's Disease Neuroimaging Initiative. Sparse learning and structure identification for nal of the American Statistical Association, 1166536):1994-2088 2021.
image-on-seed regression. Jour-
Dun Mao, Steffen Kandler, Bruce L McNaughton, and Vincent Bonin. Sparse orthogonal population representation of spatial context in the retrosplenial cortex. Nature communications, 8(1):243, 2017.
Mark J Meyer, Brent A Coull, Francesco Versace, Paul Cinciripini, and Jeffrey S Morris. Bayesian function-on-fuuccion regression for multilevel functional data. Biometrics, 71(3):563-574, 2015.
Emily L Morris, Kevin He, and Jian Kang. Scalar on network regression via boosting. The annals of applied statistics, 16(4):2755, 2022.
James O Ramsay and Bemarrd W Silverman. Firting differential equations to functional data: Prin- cipal differential analysis. Springer, 2005.
Carl Edward Rasmussen and Christopher K Williams. Gaussian processes for machine learning. Adaptive Computation and Machine Learning series. MIT Press, London, England, November 2005.
Philip T Reiss, Lei Huang, and Maarten Mennes. Fast function-on-scatt regression with penalized basis expansions. The international journal of biostatistics, 6(1), 2010.
Michael Smith and Ludwig Fahrmeir. Spatial bayesian variable selection with application to func- tional magnetic resonance imaging. Journal of the American Statistical Association, 102(478): 417-431, 2007.
Chandra Sripada, Mike Angstadt, Saige Rutherford, Aman Taxali, and Kerby Shedden. Toward a "treadmill test" for cogmition: Improved prediction of general cognitive ability from the task activated brain. Human brain mapping, 1166:377633 2020a.
Chandra Sripada, Saige Rutherford, Mike Angstadt, Wesley K Thompson, Monica Luciana, Alexan- der Weigard, Luke H Hyde, and Mary Heitzeg. Prediction neurocognition in youth from resting state IMRI. Mol. Psychiatry, 25(12):34133422) December 2020b.
Jonathan E Taylor and Keith Worsley. Detecting sparse signals in random fields, with application to brain mapping. Journal of the American Statistical Association, 102(479):913-928, 2007.
Mikkel Wallentin, Andreas Roepstorff, Rebecca Glover, and Neil Burgess. Parallel memory systems for talking about location and age in precuneus, caudate and broca's region. Neuroinage, 32(4): 1850-1864, 2006.
Xiao Wang, Hongtu Zhu, and Alzheimer's ' Disease Neuroimaging Initiative. Generalized scalar-on- image regression models via total variation. Journal of the American Statistical Association, 112 2017.
Max Welling and Yee W Teh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th international conference on machine learning (ICML-11), pp. 681-688, 2011.
Tung-Yo Yu Wu, YX Rachel Wang, and Wing H Wong. Mini-batch metsopolis-haatingg with reversible sgid proposal. Journal of the American Statistical Association, 117(537):36-334,, 2022.
Yuliang Xu and Jian Kang. Bayesian image mediation analysis. arXiv preprint arXiv:2310.16284 2023.
Yuliang Xu, Timothy D. Johnson, Thomas E. Nichols, and Jian Kang. Scalable bayesian image-on- scalar regression for population-scale neuroimaging data analysis, 2024.
Shan Yu, Guannan Wang, Li Wang, and Lijian Yang. Multivariate spline estimation and inference for image-on-scalar regression. Statistica Sinica, 31(3):1463-1487, 2021.
Zijian Zeng, Meng Li, and Marina Vannucci. Bayesian image-on-scalle regression with a spatial
global-local spike-and-slab prior. Bayesian Analysis, 101):1226, 2022.
12

---

Published as a conference paper at ICLR 2025
Daiwei Zhang, Lexin Li, Chandra Sripada, and Jian Kang. Image response regression via deep neural networks. Journal of the Royal Statistical Society Series B: Statistical Methodology, pp. qkad073, 07 2023a. ISSN 1369-7412. doi: 10. 10.1093/jssssb(kka007.. URL https://doi,
Daiwei Zhang, Lexin Li, Chandra Sripada, and Jian Kang. Image response regression via deep neural networks. Journal of the Royal Statistical Society Series B: Statistical Methodology, 85 (5):1599-1614, 1614, 2023b.
Hongtu Zhu, Jianqing Fan, and Linglong Kong. Spatially varying coefficient model for neuroimag- ing data with jump discontinuities. Journal of the American Statistical Association, 109(507): 1084-1098, 2014.
Hongtu Zhu, Tengfei Li, and Bingxin Zhao. Statistical learning methods for neuroimaging data analysis with applications. Annual Review of Biomedical Data Science, 6(1):73-104, 2023.
Hui Zou and Trevor Hastie. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society Series B: Statistical Methodology, 67(2):301-320, 2005.
A APPENDIX
A.1 DETAILS ON IONS REGRESSION
For confounder coefficients & and the individual effects 7h, we assign the Gaussian Process prior with the same kernel function x for computational convenience. The individual effect parameter f separates the spatially correlated noise from the noise term 6, and avoids setting a dense correlation matrix for the noise term which speeds up the computation. This is similar to the correlated noise model in //hu et al. al.] (2014). The identifiability of model equation has been shown in Zhang (2023a) under the following sufficient conditions: (1) the design matrix X = (X, C) Bax FT) is full rank matrix, (2) for any and any 8y, denote m)s)) = (1h(s)) (..) Th (sy)) R",
(s,) = 0. The first condition is easily satisfied when the design matrix X is not linearly dependent.
For the Gaussian Process prior update of ε and Th. we use the basis decomposition approach. Leveraging Mercer's 's theorem, which asserts that for any function g(s) following Gaussian Process with mean zero and covariance function x²x(s,)) )), we can utilize the following basis decomposition.
((x) = Σ (g)))()),, 0g ind (0,, 1=1
where d is the 1-th eigen-value, and di is the I-th eigen-function (see Section 4.2 in Rasmussen & Williams (2005)). In practice, we choose a finite L as the cutoff on the number of bases, and approximate g(8) by E/= (s). The number of basis I is chosen such that the summation Et d over 90% of EP di. The choice of the kernel function includes exponential square kernel, Matérn kernel, and other kernel functions. For the simulation section, we use the modified exponential square kernel, x's,s':a,) analysis with ABCD data, the kernel is a pre-tuned Matérn kernel with region region-specific smoothness parameters that can best align with the empirical correlation of the observed image data, same as in Xu & Kang (2023).
exp[-a(s² s²2) b(s z¹)2). For the real data
A.2 COMPUTATIONAL DETAILS
For the neighborhood matrix B ST-CAR B), in order to speed up the computation, we use sparse matrix structure in ReppArmadillo (Eddelbuettel & Sanderson 2014), and set a fixed band- width [[] for all j. For a given fixed gnd [s.,..., $p in Ra, we use RANN package (Arya et al. 2019 to efficiently search for the nearest neighbors in the high-dimensional setting.
13

---

Published as a conference paper at ICLR 2025
A.3 STOCHASTIC SUBSAMPLING VARIATIONAL INFERENCE (SSVI)
To make the variational inference method scalable for large data sets, we propose a stochastic sub- sampling version of CAVI, referred to SSVI. The main computational bottleneck of CAVI is to update B. which is a high-dimensional parameter, and the latent variable I further requires complex computation of mixed truncated normal densities. Hence given μt, when updating B3, we randomly select a subsample of data, indexed by I C {I..... n). and apply a stochastic gradient update sim- ilar to the Stochastic Gradient Langevine Dynamics (SGLD) (Welling Teh, (201T) Let 8t be the step size at t-th iteration, be the total number of observations, ns be the subsample size, and be the prior density of B3 at the jth voxel.
1gg0 (3/) < B) {B,} + 9t
-
(.yutn log Ep(Y,, 0) + VI log π (8))
'EI
This is because under the mean-field assumption, the optimum density q (B(3) has a closed-form solution: normal density with mean and variance
By (8,) Var, (3))x
N Σ
Bg' (o 2)
Mij
Y; Ey
TX,
Σ KE[-5l
Eg 3g + Eg foz (o32To(us)}}
1=1
Varg' (8,)= Bag (oy²)
Σ
Mag+ + E.g (07²)
And Egg (By) is also the maximizer
N log PY, M, X, 0) + Eg- log π (82).
Egg
We require the step size 8t to decrease to 0 as t + oc. In practice, we use the decay function s, a(b + ? suggested in Welling & Teh 2011).
In practice, we find that in low signal-to-noie toiise ratio (SNR) settings, the CAVI algorithm gives better accuracy. Hence we recommend using CAVI for the Sonl model, where the SNR be very low, especially in brain imaging data, and using SSVI for the lonS model since the lonS model has a much higher SNR for the coefficient at each voxel. The implementation is available as R package STCAR in the supplementary materia
A.4 SIMULATION II: IMAGE-ON-SCALAN REGRESSION WITH SSVI
For IonS model equation we compare ST-CAR with other methods: (1) STGP prior, (2) Scal- able Bayesian Maage-on-Scalaat regression (SBIOS) (u et al. 2024). (3) Mass Univariate Analysis (MUA). For the lonS regression equation 4 estimation of a has a larger SNR compared to esti- mating in SonI equation 3 hence we impose GP prior for the confounder parameters &a and individual effect nj. the GP kernels used in this simulation use region-wis wise independent kernels for the GP priors in equation The GP kernel is the same as equation 5 with a 0.01 and b 10.
SSVI for this application for ST-CAR prior. Because
all the same for STGP, ST-CAR, and SBIOS for fair comparison. We also
The mass univariate analysis (MUA) is one of the most commonly used methods for IonS regression. MUA analyzes lonS spatially independent problem and treats the IonS regression as p indepen- dent linear regression problems with exposure X; and confounders C, To select active voxels, we use the Benjamini-Hosee Hochberr adjusted p- p-values (Benjamini & Hochberg 1995) to control the false discovery rate. The active voxels selected by MUA have an adjusted p-value below 0.05.
The STGP method is similar to what has been discussed in the Sonl regression. For the IonS regres- sion, we use total of 2 x 10 iterations and take the last 10% as the converged MCMC sample. The
SSTCAR is also available Github https: thub. http://ggthhh...como/o//aa/g////
14

---

Published as a conference paper at ICLR 2025
thresholding parameter y is set to be 0.2. We use the point estimates of a and & from MUA as the initial value for the MALA algorithm.
The Scalable Bayesian Image-on-Scalat regression (SBIOS) (Xu et al. 2024) is another Bayesian approach where the parameter of interest can be expressed as a(s) = a (s))(s). The latent spatially smooth function á is assigned GP prior, and the binary selection variable ō(s,) is assigned an independent prior Ber(p(s;)) for each location 5j- SBIOS is designed to analyze a large-scale data set by using batch updates with stochastic gradient Langevin dynamics algorithm (SGLD). Hence it is more appropriate to be compared with the SSVI implementation of ST-CAR, since both methods are based on stochastic gradient updates of small random sample drawn from the entire observed data. Different from the idea of SSVI where we simply use stochastic gradient update for an optimization problem, SGLD gives a smooth transition from optimization to MCMC sampling as the step size decays to O WWeling & Teh 2011). Similar to STGP, we can use the MCMC sample of 8(s;) to determine the PIP at location j. PIP, Σi_ o(s,) # 0/T for T MCMC sample. In the simulation, we use 5000 SGLD iterations, with the decay function of the step size set as 8, 0.0001 (10 + +1) We 200 subsample in each iteration. The prior for 0(8;) is set to be (er().55) for all locations. The last 20% of iterations is used to compute the point estimation of a and PIP.
The ST-CAR method implemented using the SSVI algorithm requires a stochastic gradient update of a. We use step of 10 and subsample of 100 for the SGD optimization, The decay rate function for oa is (1 + -).4 We use C+ (1) as the prior for oμ in equation
Because of the randomness in the SGD update, we cannot use the difference between (1 and ( or ELBO as a stopping rule. For the simulation, we simply run 10a iterations. In practice, the convergence of SSVI can be roughly determined by the convergence of of For the point estimation and inference of a, we use the averaged values over the last 20% iterations as the posterior mean of a and PIP to avoid the randomness from SGD.
Note that updating the individual effects 1/s for i = 1,..., n is computationally challenging for all Bayesian methods. We choose to update 77 every 100 iterations for ST-CAR, SBIOS, and every 1000 iterations for STGP.
STRCARI
CA-CARCA
12-
10
O
0 02-
3 DU
Z
7
7
42
41
24
69
16
40
09
91
10
100
Valun
1.
MUA
9B/O8
STGP
11-
0
10-
10-
-10
05
3
80
00
7
Z
85
10
Figure A1: Point estimation result of IonS regression for all competing methods, n 600, p
6400, 0h = 5. The top-left figure is the true a signal.
Figure Al provides visualization of the point estimation for each method. MUA has the most noisy point estimation since it does not consider the spatial correlation, and there is no sparsity constraint directly imposed other than using the adjusted p-value to determine the level of significance for each voxel location. STGP suffers from the same issue as in Figure 2 where the pre-specified kernel is too smooth for the Z-shape and recycle shape (top-left). SBIOS uses the same kernel, but the binary selection parameter ō(s,) has spatially independent prior, and can get clearer edge compared to STGP and better selection, but the latent GP kernel still too smooth that the edge of the recycle shape and Z-shape tends to be O. In the ST-CAR plot, although the functional parameters εx, 1)
15

---

Published as a conference paper at ICLR 2025
are all assigned GP prior with the same kernel as STGP and SBIOS, still that ST-CAR can give a very clear edge for all 4 shapes. This demonstrates that ST-CAR prior is very flexible to different correlation patterns without much tuning the neighborhood matrix B or correlation coefficient p. especially for the high SNR cases.
Table Al: Numeric result for IonS simulation, under 100 replications.
(a) IonS: Comparison of estimation accuracy. The evaluation criteria include false discovery rate (FDR), true positive rate (TPR). overall accuracy (ACC), and root mean squared error (RMSE), all multiplied by 100.
Case 11 = 600, == 1600, 03
5
Case 2. f = 600, p = 900, 03. =
Criteria
ST-CAR
MUA
STGP
SBIOS
Criteria
ST-CAR
MUA
STGP
SBIOS
FDR
5,8
7.98
16.88
6.3
FDR
4.95
8.23
12.75
4.69
TPR
95.41
93.2
94.18
94.94
TPR
84.34
81.42
85,01
84.61
ACC
97.86
96.95
94.89
97.66
ACC
96.27
95.19
94.9
96.36
RMSE
7.86
9.35
10.53
10.79
RMSE
7.88
9.4
6,85
7.36
Case = 1000, p = 1600, d3 5
Case n = 60D, P = 1600, T7 = 10
Criteria
ST-CAR
MUA
STGP
SBIOS
Criteria
ST-CAR
MUA
STGP
SBIOS
FDR
7.1
8.1
19.24
7.97
FDR
5.07
8.06
14.87
4.12
TPR
98.38
97.32
95.57
97.53
TPR
85.31
82.61
90.84
86.19
ACC
98.13
97.69
94.43
97.76
ACC
96.06
94.96
94.88
96.42
RMSE
6.44
7.21
10.16
10.35
RMSE
10.32
13.14
11.21
11.74
Case n1
600, = 6400, oh
5
Case 6, 1000, p
6400, d3
20
Criteria
ST-CAR
MUA
STGP
SBIOS
Criteria
ST-CAR
MUA
STGP
SBIOS
FDR
5.97
7.84
20.78
6.04
FDR
5.93
7.93
19.94
2.95
TPR
93.64
91.78
97.19
93.62
TPR
81.83
80.09
96.64
84.46
ACC
97.35
96.55
93.9
97.33
ACC
95.06
94.32
94.18
96.16
RMSE
8.52
9.19
9.79
10.25
RMSE
13.84
14.19
9.76
10.72
(b) Computation time for IonS simulation, averaged over replications.
Computation time
Total time (saconds)
Number iterations per second
Case
ST-CAR
STGP
SBIOS
ST-CAR
STGP
SBIOS
Case I. 600. p 1600, of.
73.7
588.4
717.4
137.5
3.4
7.3
Case 2. 600, p 900, ou
55
381.6
874.5
186
5.3
6.3
Case 3. 1000.p 1600 di 5
117.3
1062.
1968.4
88.9
1.9
3.1
Case 4. 600. p 1600, z 10
82.9
621.8
1214.1
122.7
3.2
4.9
Case 5. 600, 64400,
dM 5
409.2
2190.3
1049.8
24.7
0.9
5.3
Case 6. 600, 6400, OM
20
596.7
5090.2
1733.2
17
0.4
3.2
provides paaooidss the numerical result on the IonS model based on 100 replications in six different settings. Because the predictive MSE on the outcome averaged over all voxel locations very close for all methods, we do not report it here. Instead, we focus the estimation of the coefficient a. The proposed ST-CAR prior with the SSVI algorithm gives the lowest RMSE except for Case and 6, for which the STGP has the lowest RMSE, although STGP has much larger FDR in both cases. To control the FDR below 10%, we use the Henjamini-Hochbeee (BH) adjusted p-values on MUA and set threshold such that (s;)) with the adjusted p-value below 0.1 are selected active voxels. For ST-CAR, STGP, and SBIOS, we compute the proportion of the active voxels selected by MUA, and apply the same proportion to get the cutoff on PIP. In this way, we select roughly the same proportion of voxels as active. Based on the result in Table Ala this selection method can control the FDR for ST-CAR and SBIOS to be below 10%, whereas for STGP, the FDR still over 10%. The MUA has the worst power (TPR) in all scenarios after controlling for FDR. The total running time shown in Table(AID]also shows a great improvement in the computational speed for the SSVI algorithm when compared with other MCMC sampling algorithms. On average, STGP takes 7.6 times longer compared to SSVI, and SBIOS takes 10.4 times longer compared to SSVI. In the Supplemental Material, also provide additional results of using CAVI under ST-CAR prior and compare the performance with SSVI. SSVI still slightly outperforms CAVI the IonS regression in terms of both estimation and computation speed.
Table
For a more comprehensive comparison between different estimation algorithms (Gibbs, CAVI, SSVI) under ST-CAR prior, we include small -ineenneniooaal comparison for the Sonl regres- sion in the Supplemental Material. The result suggests that CAVI tends to have better estimation accuracy in Sonl regression where the signal-to-noise ratio is low.
16

---

Published as a conference paper at ICLR 2025
A.5 POSTERIOR DERIVATION
For the fully conjugate posterior derivations, the hierarchical model of applying the sparse-mean prior on the scalar-on-image regression can be written as
iid Yi = Σ B(s)) + + X, t ~ j=1
N(0,0,},
~
B(5)) ind (TT,())),3))
oB 38~0t(1)
M P[-) ~ N (PNO)) +f), o o2 ( o0 of (1)
Based on this hierarchical model, we can derive the posterior distributions of each parameter. Note that the posterior for most of the parameters is straightforward, except for Mjj which involves the soft-threshodding operator Ty. The posterior μaj can be expressed in terms of a mixture of trun- cated normal distributions with 3 components.
A.5.1
PROOF OF PROPOSITION
Proof. The posterior of μy can be expressed as
Witt (M3 PNGG)) log B. n. oB) x 20 (3(,) - Tx(1)))²
(6)
20
<(u ....) ) PT (u,t V2) + PP N ((- Vo) P N4- Vz)
where Na (H, d²) is notation for the truncated normal distribution supported on [a, with mean j and variance o². The middle component is just the truncated normal on [-v, v/ with the prior mean ANG) and variance Vo-
Lit. For the other two components,
Vn =
Wi+
1r = Vi
(B(s)) +v + ANG)) j
Wi+
7 3
,
03
(3(s)) v) + Wist PNG)
1 = Vi
of
The density of this 3 component mixture can be expressed
*(u ....) 2 Z
where
normalizing 19.5 represent the density functions of 3 truncated Vo), No (N.-No...u) V2) respectively.
24,29.27-11 Z3
represent different
constant, and normal distributions Hence the mixing
N No,+ (uf, N probabilities can be represented
z+ Z,
7 PP
7 Z
p
PJ
Z
Now denote f,,* {-, 0, +} as the RHS in equation -v), v), (v, +00) respectively.
supported
log(2 = log f - log f
x E (v, +00)
= 203 (8(x)) + v 11)²
Wi+
(n3 - -
20
log VV,
1
(4g 11)2 2 log
1-
1/ - 11)
2V,
VV,
17

---

Published as a conference paper at ICLR 2025
(og( = log /2] log [g:
x [-v, +u]
203 (B(s)))²
Wjt (1g iNNG))
20 /
{ log 1
2Vo (1² AN())² - loog 1 V (N( (
VV
VVo
log( (,) ) log fj-
log fj. I E
20 (B(s)) M4 v)²
wit
(es (NNG))) )
20
log
1
(nj - -²) - log 1- 1
Hi
2
VV,
Hence the entire density function is complete.
A.5.2 VARIATIONAL INFERENCE: Q-DENSITIES FOR SCALAR-ON-IMAGE REGRESSION
In the following derivation, we denote the vector Y E R", matrix M E R X E R"x to denote the outcome and design matrices.
Q-density for 33r using SVD
First, we use Singular Value Decomposition (SVD) on M and re-express the scalar-on-image re- gression model as follows.
Let the compact SVD of M E RRX be M UDVT where U E R=xn VT R"*p, and U = UU In, VTV In B- Tx(w), = MTx(H) X = M3+e.
Now apply the rotation matrix U on both sides, Y UT DVT 3 c.c The q-density for 3 is now normal density with mean and variance
1
VVT
Var, (31~~ ~) ~ Egg
Lp
Ey
033
1 VDTT.
(8l ~) = Var, (3) ~2)
En
Tr
Note that Ey( ~) can be further simplified,
Eg VVD VD
1 c
In + Dz
y.
Eg
where 22=
Then Eg (3) ~) = [3] ~) + Tp(p).
The g-density for og is as follows. Note that we use the hierarchical expression to sample half- Cauchy prior 0g ~ C+ (1): oz ~ 1 d ag ~ 1G 1).
p log(o - (5 (1+1 log ag oz 2
1
( ~)
(B(s)) Tx(1)))
20
Hence the q-density for o3 follows IG
p+1 1 Σラ= E. (B(s)) - ((5))) + 2)) 2 :2
Note that - T, (1) and the marginal variance Var, (B(,)))))))) =
Egg 11B - T,(11 = (3)) + [3,(3)
E. ( Er M2 + Eg
The q-density for 0, is IG 1, 1 1
18

---

Published as a conference paper at ICLR 2025
Q-density for y The q-density of follows the multivariate normal distribution with mean and variance
Var,(1| ()) = By ox.
X X + By
04 (1
E(7l () = (1l ~)
B
Σ (Y; Mf X;
To speed up the computation, we use eigen-decomposition x = QAxQ", and the variance update can written as
Egg ( Ax + Egg
1
Varg(1) ())(0))~~--)0))))) ~) Qdiag
Q².
0.
Similarly, the q-density for of follows IG
B 17 2 + Eg (2m)
2
The q-density for ay is IG 1, E
oo })
Q-density for oy
111- 1
1
o₂ 2 IG
Belly
MB X7|12
E,g
ax
ay IG
1,1+ Bg
oy
ELBO derivation
ELBO -R. flog 7 (Y M, 3, 23.
-E, E. flog (B) + log ((7) log (og)(o))) log ((o)) log {(}))}
-F, (log Y| r)]
+ E. flog T[3] ~) - log ((3) + E. (log n (n) ~) - log (1x))))
+ E, flog nT ( ~) - -log -logg(o3)) + E. flog T Tapl ~) - logg
F. flog ( ~) log + F.y flog n (a ~) {o } + B, flog (ay) log {(a)))
E. (log n ( ~)) - logg(o})} + E., flog π (ayl ~) log ((r)))
In the implementation, we separately compute each part of the ELBO and add them together.
ELBO logi -E,y (log <(Y
M, B, o7 o o
1
= -11=
Eg
Ey
Eg Y -M3 ~xy12
2
2
o
Here, denote E,SSE = EggyY - MB X
SSSSE = || - MEy3 XE + Tr { _ M Var, (9)) + Tr {x" X (7)
With the eigen decomposition on XTX,
Ax + ) +)
Tr Trfx"xVarz (7)) = Tr Axdiag By (
Eg
A.6 SENSITIVITY ANALYSIS CHOOSE HYPERPARAMETERS
In addition, Table A3 provides additional sensitivity analysis results on Sonl when the bandwidth is 26. Table Table provides additional sensitivity analysis results on IonS when the bandwidth is 9 and the decay rate parameter 0.35.
19

---

Published as a conference paper at ICLR 2025
Table A2: Sensitivity Analysis on SonI and IonS regressions.
(a) Sonl for varying and initial value for oz. The Elastic Net and STGP results shown as comparison. Bandwidth=9 in the ST-CAR model. Additional sensitivity analysis where bandwidth=26 and varying decay rate y for oß is available in the Appendix.
10-5
10-5 5
10-5
10-5
5 10
10-5
5 10-6
ElasNet
STGP
0.003
0.005
0.007
0.01
0.005
0.005
0.005
test PMSE
0.58
0.5
0.48
0.48
0.61
0.5
0.55
0.53
0.5
Test R
0.16
0.28
0.30
0.30
0.12
0.28
0.20
0.23
0.28
train PMSE
0.17
0.24
0.27
0.29
0.11
0.24
0.21
0.45
0.49
Train R
0.75
0.65
0.61
0.58
0.84
0.65
0.70
0.35
0.29
(b) IonS for varying v, initial value for o². and decay rate y for a. The total test pMSE is the summation of all voxel-level PMSE. Bandwidth is 26. Additional sensitivity results where bandwidth=9 are available in the appendix.
initial af
0.1
0.01
0.1
0.1
0.1
0,1
0.1
0.1
0.005
0.005
0.005
0.05
0.005
0.005
0.005
MUA
Decay rate
0.35
0.35
0.35
0.35
0.35
0.35
0,25
0.45
0.55
total test PMSE
47357.74
47351.78
47354.06
47354.13
47351
47354.1
47354.12
47350.7
47352.68
47487.05
Table A3: Additional sensitivity analysis for Sonl when the bandwidth is 26, on three parameters: (i) the initial value of o.. (ii) the thresholding parameter in ST-CAR prior, (iii) the decay rate y in the decay rate function for o² where (az) a(b + 17
c3
6
x 10-5
10-5
1 10-4
x 10
1 10-5
10-5
1 10-5
x 10-5
0.007
0.007
0.007
0.007
0.005
0.01
0.012
0.007
0.007
0.35
0.35
0.35
0.35
0.35
0.35
0.35
0.25
0.45
test PMSE
0.54
0.51
0.56
0.64
0.52
0.52
0.51
0.52
0.51
train pMSF
0.29
0.37
0.31
0.22
0.33
0.4
0.41
0.36
0.37
A.7 ABCD DATA ANALYSIS USING ELASNET AND MUA
In this section, we perform the baseline methods ElasNet for Sonl regression and MUA for IonS regression as sanity check the ABCD data analysis.
For SonI using Elastic Net with a = 0.5 and À chosen with cross-validation, we visualize the point estimation 3 in Figure A2a The Elastic Net has straight forward uncertainty quantification method such as p-values, hence we can only compare the point estimates. As expected, the scale of using Elastic Net is also very small, and there are only few small areas with near or greater than 0.001, as highlighted by the red areas.
For lonS using MUA, we note that after using BH-adjusted p-values in the MUA, none of the p- values can pass below the 0.05 threshold, which means MUA cannot identify any signal with a reasonable adjusted p-value. Nonetheless, we still visualize the point estimation of on the whole brain. Based on Figure A2b we can see that the MUA estimates of a with large positive effects pick up lots of noises, and compared with A the ST-CAR's estimate of active a subarea of the MUA's estimate and is validated by high PIP (> 0.95).
In addition, we present the estimation by ST-CAR after thresholding by PIP 0.25 in Figure A3 In fact, both Figure/A2aland A2a]and Figure A3 demonstrate the true signals are very sparse in the Sonl problem. ST-CAR still tends to identify more active voxels than ElasNet, and some active areas still overlap the sagittal plane (first plot from the right) with ElasNet result.
To address the concern that the identified voxels in SonI might be due to random noises, we have conducted stability check: we rerun the SonI on slightly reduced subset of the ABCD data, by randomly removing 10% (or 30%) of the data to re- re-train the model, and check if the identified activa- tion voxels in the full sample can still be recovered in each random subsample. We have conducted a total of 50 splits, and report the results in Table A55 where we combine the identified voxels by ST- CAR and ElasticNet in the full sample analysis to as the "ground truth".. The numbers in each
20

---

Published as a conference paper at ICLR 2025
Table A4: Additional sensitivity analysis for IonS when the bandwidth is > 0.35 in the decay rate function, on two parameters: (i) the initial value of d (ii) the thresholding parameter v in ST-CAR prior
1
0.1
0.01
0.1
0.1
À
0.01
0.01
0.01
0.005
0.05
total test PMSE
47362.8
47351.85
47356.9
47354.97
47353.58
(a) Sonl-ElasNct: values of B (Red) color range [0.0005,0.001]]
(b) lonS-MUA: point estimates of (Red) w. color range [0.05,0.
Figure A2: Visual illustration of 33 in Sonl and using MUA in lonS.
row represent the proportion of voxels being repeatedly selected in 50 splits. For example, in the top half of Table A5 0.89 for row "ST-CAR: PIP>-0.25" means that 89% of the voxels are selected by at least 10% of the 50 subsampled studies, with the activation criteria "PIP:>0.25", similarly, 0.21
means 21% % of voxels are repeatedly selected by at least 70% of the 50 subsampled studies.
Table A5: Result of stability check on Sonl
When randomly removing 10% of the data
threshold ST-CAR: PIP50.25
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.89
0.83
0.72
0.59
0.42
0.32
0.21
0.13
ST-CAR: [8(s)] > 10-4
0.89
0.83
0.72
0.59
0.32
0.21
0.13
ElasticNet: [8(s)
> x 10
0.23
0.21
0.18
0.15
0.14
0.12
0.1
0.08
When randomly removing 30% of the data:
threshold ST-CAR: PIP>:0255
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.99
0.83
0.63
0.46
0.31
0.12
0.08
0.06
ST-CAR: [B(s)] > 5 10-4
0.77
0.58
0.39
0.27
0.13
0.07
0.05
0.05
ElasticNet: [3(x) > 5 x
0.38
0.32
0.24
0.18
0.14
0.09
0.04
0.03
-79
Figure A3: Sonl: values of 3 T(PIP> 0.25) (Red) color range [0.0005,0.001].
21

---

Published as a conference paper at ICLR 2025
The results show that, in both cases, ST-CAR consistently identifies more voxels than ElasticNet. This indicates that ST-CAR effectively borrows information across spatial locations, reducing sus- ceptibility from noises compared to ElasticNet. The stability results indicate more robust perfor- mance of ST-CAR, particularly in scenarios with low SNR and unknown ground truth. This stability check approach has been widely applied in other neuroimaging studies where the underlying ground truth is not available and performance evaluation is challenging (Zhang et al. 2023b Morris et al. 2022
Zhu al. 2023).
A.8 SCIENTIFIC EVIDENCE FOR THE MOTIVATION FOR ABCD STUDY
For the two scientific problems mentioned in Section 5, we provide more scientific evidence on the importance of this Sonl and IonS problem in this section.
Recent studies (Cermakova et al. 2023 Halabicky et al. 2023) have demonstrated that parental education levels are significantly associated with children's 's cognitive abilities, including specific cognitive functions such as working memory. In particular, children's 's cognitive functions are often associated with their brain development, which fMRI captures during tasks like the working memory task in the ABCD study.
We chose to use general cognitive ability the outcome because it provides comprehensive as- sessment of child's overall cognitive function, which encompasses not only working memory but also other critical skills such as reading. spelling, and math abilities (Alloway & Alloway| 2008). often correlated and influenced by common neural processes, making general
These abilities cognitive ability a robust and representative measure for examining broader cognitive development. By using this summary outcome, we can account for the cumulative effect of parental education on multiple facets cognitive performance, providing a more holistic view of the relationship between socioeconomic factors and brain development.
A.9 TABLES OF STANDARD DEVIATION FOR TABLE AND A
In this section, we present the tables of standard deviations over 100 replicates for the mean values presented in Table and Ala respectively. shown in Table A6
A. ADDITIONAL SIMULATION RESULTS
In the first simulation AI, we provide low dimensional comparison on the Sonl regression between three different implementation of the ST-CAR model, the Gibbs sampler, the CAVI algorithm, and the SSVI algorithm.
In the second simulation A2, provide the additional to results to the simulation II IonS with further comparison between CAVI and SSVI in high-dimensional settings.
A.10.1 SIMULATION A1: Low DIMENSIONAL COMPARISON (SONI)
We compare the proposed Gibbs, CAVI, SSVI with ST-CAR prior the classical penalized gression method glmeet, and Bayesian method where 3 is assigned a Soft-thresholded Gaussian Process prior (Kang et al. 2018) implemented in the IMAA package.
The Frequentist penalized regression is implemented using R package glmnet (Friedman et al. 2010b) with lasso penalty (a 1), using 10-fold cross-validation.
The BIMA method requires a pre-specified kernel function, and the posterior sampling algorithm is Metropolis-adjusted Langevin algorithm (MALA). In this simulation we sue the exponential square kernel
k(s, s'; a, cor ((())))(())) exp(-a(é + 2²2) b(s -s))²)
where a = 0.01, b - 10, and used L = 66 basis functions.
For the four Bayesian methods (Gibbs, CAVI, SSVI. BIMA), set the thresholding parameter = 0.1. To evaluate the variable selection accuracy, for the variational inference ST-CAR meth- the mixing probabilities shown in define the posterior inclusion
ods (CAVT, SSVI),
22

---

Published as a conference paper at ICLR 2025
Table A6: Standard deviation result for Sonl simulation, under 100 replications.
(a) Sonl: standard deviations over 100 replicates. The evaluation criteria for estimation includes fälse discovery rate (FDR), true positive rate (TPR), overall accuracy (ACC), and root squared emror (RMSE), all multi- plied by 100. The evaluation criteria predictive performance includes training and testing predictive MSE, denoted Train Test pMSE respectively
Case 1. n. 600,
1600, of
Case 2. nI 600, = 900,
ST-CAR
ElasNet
STGP
T-LoHo
ST-CAR
ElasNet
STGP
T-LoHo
FDR
0.41
0.47
1.37
9.34
FDR
0.40
1.38
1.32
12.86
TPR
1.37
7.69
1.53
1.23
TPR
0.47
0.95
1.42
0.70
ACC
0.20
1.04
0,36
3.08
ACC
0.10
0.26
0.28
5.10
RMSE
11.00
30.45
13.44
9.46
RMSE
5.51
7.12
13.00
6.56
Test PMSE
0.37
1.19
0.49
0.68
Test PMSE
0.10
0.05
0.24
0.29
Train p.MSE
1.28
9.32
1.09
1.29
Train PMSE
0.20
0.33
0.56
0.45
Case 3.
1000.
1600, o.
Case 4. 600, p 1600, o. =
ST-CAR
ElasNet
STGP
T-LoHo
ST-CAR
ElasNet
STGP
T-LoHo
FDR
0.95
0.46
0.89
4.40
FDR
0.29
0.50
1.13
12.86
TPR
0.00
0.73
0.60
0.31
TPR
3,02
7.01
2.61
2.05
ACC
0.17
0.13
0.19
0.91
ACC
0.40
0.95
0.43
7.18
RMSE
4.27
6.53
12.04
5.83
RMSE
14.65
33.31
14.72
12.51
Test PMSE
0.07
0.06
0.32
0.23
Test PMSE
0.62
2.60
0.89
0.97
Train p.MEE
0.17
0.50
0.47
0.36
Train PMSE
2.47
10.40
1.75
1.74
Case
6000,
10000, o.
Case 6.
6000, P
10000, o..
10
ST-CAR
ElasNet
STGP
T-LoHo
ST-CAR
ElasNet
STGP
T-LoHo
FDR
0.24
0.15
0.57
12.85
FDR
0.80
0.62
0.61
14.36
TPR
0.07
0.25
0.09
0.17
TPR
0.10
0.36
0.09
0.19
ACC
0.04
0.05
0.13
9.76
ACC
0.13
0.12
0.14
10.27
RMSE
2.26
3.00
6.43
2.30
RMSE
3.14
3.94
6.44
2.48
Test PMSE
0.18
0.23
1.40
1.54
Test PMSE
0.31
0.40
1.76
1.42
Train PMSE
0.32
1.33
1.58
1.88
Train PMSE
0.50
2.17
1.67
1.92
(b) IonS: standard deviations 100 replicates. The evaluation criteria for estimation includes false discovery rate (FDR), true positive rate (TPR), overall accuracy (ACC), and root mean squared error (RMSE), all multi- plied by 100. The evaluation criteria for predictive performance includes training and testing predictive MSE, denoted as Train and Test PMSE respectively
Case 1, n 600, p 1600, o
5
Case n 600, p = 900, TM
Criteria
ST-CAR
MUA
STGP
SBIOS
Criteria
ST-CAR
MUA
STGP
SBIOS
FDR
1.49
1.55
1.42
1.49
FDR
2.16
2.44
2.35
2.12
TPR
1.31
1.45
0.88
1.25
TPR
3.26
3.11
2.01
3.04
ACC
0.36
0.44
0.42
0.35
ACC
0.66
0.68
0.63
0.57
RMSE
0.24
0.32
0.12
0.15
RMSE
0.28
0.35
0.25
0.22
Case 3. n 1000, p
1600, of
Case 4. n 600, p = 1600, 0M
Criteria
ST-CAR
MUA
STGP
SBIOS
Criteria
ST-CAR
MUA
STGP
SBIOS
FDR
1.42
1.44
1.27
1.47
FDR
1.48
1.76
1.55
1.35
TPR
0.55
0.65
0.56
0.61
TPR
2.16
1.95
1.29
2.15
ACC
0.35
0.38
0.40
0.38
ACC
0.45
0.47
0.43
0.41
RMSE
0.16
0.20
0.07
0.08
RMSE
0.33
0.46
0.17
0.21
Case 5, n 600, p = 6400, of
Case
1000, p = 6400, 03
20
Criteria
ST-CAR
MUA
STGP
SBIOS
Criteria
ST-CAR
MUA
STGP
SBIOS
FDR
0.73
0.74
0.79
0.77
FDR
0.88
0.68
0.90
0.65
TPR
0.80
0.89
0.28
0.75
TPR
1.22
1.23
0.41
1.19
ACC
0.21
0.25
0.26
0.21
ACC
0.29
0.26
0.28
0.22
RMSE
0.26
0.29
0.04
0.08
RMSE
0.28
0.32
0.06
0.12
probability(PIP)
PIP(3(s,)) =1 PM
where both CAVT and SSVT can trace the mixing probability Pp. We use the converged value at the last iteration of P in CAVI to compute PIP. Since SSVI is stochastic method, use the averaged po over the last 2000 iterations compute its PIP. For the MCMC methods (Gi BTM,
23

---

Published as a conference paper at ICLR 2025
directly the posterior sample of (13s) (for Gibbs) or ((s,) (for BIMA) being the last 20% of iterations as the posterior inclusion probability. For the final selection reported in Table A7 we the true generating image 3, and set threshold (5,)) = o, otherwise (85) equals the posterior sample mean or the variational mean. By tuning t, control the FDR to be below 10%.
PIP: if PIP(((s;)) <
true:beta
BiMa
gibbs:beta
cavitbeta
sgd:beta
035
040
(a) Case
true:bete
BIMArbeta
gibbs:beta
cavitetta
agd:beta
035
040
(b) Case 2
Figure A4: Illustration of estimated / under 2 different cases.
Simulation I provides a relatively low-dimensional small-scale example, where n 200, p 400. We simulate two testing image cases for 33 in equation 3 as shown in Fig A4b In case 1, the true image intensity has a smooth transition from to and the voxels around the edge of the signal tend to have low signal-to-noise ratio, but the signal region is smooth round shape that can be easily estimated by smooth Gaussian process. In case 2, the true image of B is a sharp triangular shape, but the edge voxels of the signal has sharp contrast to O, with higher signal noise ratio compared to case 1.
Table A7 provides the simulation results of the posterior mean estimates of 8. with mean and standard deviation computed across 100 replications. SSVI has the best time efficiency across Bayesian methods.
A.10.2 SIMULATION A2: (IONS)
COMPARISON BETWEEN CAVI AND SSVI
In the second simulation A2, we provide the additional to results to the simulation II lonS with a further comparison between CAVI and SSVI in high-dimensional settings. Table A8 provides the additional result on the performance of CAVI compared to SSVI.
24

---

Published as a conference paper at ICLR 2025
Table A7: Simulation results based on 100 replications, with standard deviation in the bracket. All values are timed by 100 except for time (in seconds). FDR (false discovery rate) is the proportion of times that zero coefficients are identified as nonzero among all identified nonzero coefficients. Power the proportion of times that nonzero coefficients are identified as nonzero among all nonzero coefficients. Accuracy is the proportion of times the prediction is correct. RMSE is the root mean square error over all voxels.
Casel FDR
Gibbs
CAVI
SSVI
BIMA
glmnet
5.4(()
(2)
(2)
11)5
3.0 (3)
Power
80.0 (7)
(3)
84)9
100.0
24.1 (7)
Accuracy
92.6 (2)
91)
91)
95)
72)
RMSE
(2)
(1)
(1.0
0.5
19.7
time
97.7 (5)
(322(10)
12.7 (0)
21.4
1.2 (0)
(a) Case
Case2
Gibbs (3)
CAVI
SSVI
BIMA
glmnet
FDR
3.7 (0)
(2)
16.6
0.0
Power
100.0 (0)
100.0 (0)
97.
100.0
94.7 (4)
Accuracy
98.8 (0)
99.5
99.3 (0)
(1)4
99.3
RMSE
(1)
109
71)
2.2
118
time
101.4 (11)
15.9
(1)
21)
1.2 (0)
(b) Case
Table A8: Additional Simulation results to Simulation II. Comparison between CAVI and SSVI for ST-CAR prior, based on 100 replications.
FDR
TPR
ACC
Case
SSVI
CAVI
SSVI
CAVI
CAVI
Case n = 600, p 1600, o =
5.8
6.89
95.41
95.15
97.86
96.94
Case 2. 600. p 900, 0 5
4.95
6.23
84.34
84.05
96.27
95.32
Case 1000. p 1600, 07, 5
7.1
7.91
98.38
98.31
98.13
97.32
Case n 600. p 1600. 2 = 10
5.64
85.31
84.79
96.06
95.85
Case n 600, p 6400, 0 =5 5
5.97
11.72
93.64
94.18
97.35
91.3
Case 6. 600, p 6400,
0 M 20
5.93
8.57
81.83
83.58
95.06
91.73
Total time (seconds)
Number of iteratios
RMSE
per second
Case
SSVI
CAVI
SSVI
CAVI
SSVI
CAVI
Case 600, p 1600, o3 5
7.86
9.13
73.7
239.3
137.5
41.3
Case n 600. p 900, 0
9.17
55
155.8
186
63.8
Case n = 1000, p = 1600, 0 = 5
6.44
7.11
117.3
429.8
88.9
23.1
Case 4. n 600, p 1600, o = 10
10.32
12.7
82.9
255.4
122.7
39,2
Case 5. 600, p 6400, CM 5
8.52
8.99
409.2
1282.8
24.7
6.7
Case 600, P 6400,
20
13.84
13.72
596.7
2641.5
17
3.4
25


Reviews:
[{"review_id": "9yqntOxfmJ", "reviewer": "Reviewer_RCFm", "summary": "The authors propose a new Bayesian method for sparse regression. To solve the sparse regression settings, they introduce a method based on stochastic subsampling and variational inference. They verify their model using simulation data, and finally apply their model on fMRI task activation maps to predict the effect of brain regions on IQ development and how the task activation map is affected by parental education level. Throughout the paper the authors consider two regression settings; Image on Scalar (IonS) regression and Scalar on Image regression (SonI).", "questions": "1. L102: I don’t quite understand why the authors say that the function on L101 shrinks values below $v$ to zero, since $I(\\cdot)$ is the indicator function isn’t any value below $v$ zero?\n2. Based on Figure 1 I was wondering whether it’s not just possible to use a different kernel to get better results both on the simulated data, i.e. reducing the variance in the exponential kernel or using a kernel like the Jump Gaussian process [1]? The authors select hyperparameters for the STGP model on L288, but do not provide any reason for these hyperparameters. I would strongly encourage the authors to do a thorough hyperparameter search on a validation set for all of the methods they compare against.\n3. On Line 213/214 the authors mention the design matrix M, but isn’t M the task activation map? The design matrix normally contains regressors of interest and nuisance regressors.\n4. Have the authors looked at how initialization affects each model? I understand there are differences that may make this hard for some methods, but ST-CAR is initialized differently than the other methods.\n5. Why did the authors choose to look at IQ and parental education level effects in the brain? And why did the authors use the 2-back task activation map?\n\n[1] Park, C. (2022). Jump Gaussian process model for estimating piecewise continuous regression functions. Journal of Machine Learning Research, 23(278), 1-37.", "limitations": "", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "I think the authors set up an important and interesting problem both in their simulated setting and for fMRI data. Specifically, I believe that encoding spatial relationships in both regression settings is important, and I think the significance of the authors’ method can therefore be high. The authors compare their method to good baselines, and I think the comparisons are relatively fair. The idea is original and builds on the soft-thresholded Gaussian process, which the authors cite and compare against. I think the authors also do a relatively good job of explaining their method, and their results on simulated data are good. However, the clarity of the method can be improved.", "weaknesses": "Major weakness(es): \\\nI think the two main major weaknesses of this work are 1) that the paper lacks a thorough hyperparameter search and the method is very sensitive to hyperparameters and 2) that the fMRI analysis lacks depth.\n1. First, authors describe a lot of hyperparameters on Pages 6 and 9, but none of these hyperparameters are explained or searched over. Specifically on Page 6 for the STGP model: a=0.01, b=10, 10 degrees of Hermite polynomials per sub-region, 10^4 iterations, thresholding parameter is 0.2. For ST-CAR: thresholding parameter, decay rate, number of neighbors, correlation parameter, and fixed variance parameter. The same thing applies for Page 9. It is of course not necessary to search over all of these hyperparameters, but it’s good to explain why you set them and search over the most important hyperparameters for competing methods. Especially since hyperparameters may make a big difference. Second, the method is very sensitive to its hyperparameters, as can be seen Table A2 &A3, and the training and test R^2 values do not correlate with each other. Specifically, the best training R^2 in Table A.2 (0.84 R^2) achieves the worst test R^2 (0.12). If the authors had thus used the best training model, they would have performed worse than ElasticNet and STGP on the held-out test set (0.23 and 0.28 R^2 on the test set respectively). The same behavior is found in Table A3, where the best test R^2 (0.64) is achieved at the worst training R^2 (0.22) and vice-versa (0.51 vs 0.41). In fact this to me puts into question the simulated results, and I would strongly encourage the authors to perform a hyperparameter search, pick the best hyperparameter based on validation set results (i.e. separate the data into a training, validation, and test set), and then evaluate each model (including their own) with those hyperparameters on the test set. I specifically mention this point as well because the authors select the best v and $\\sigma_\\beta^2$ based on the best testing performance (L443), which is bad practice, since it can lead to overfitting. Especially since there is almost an inverse relationship between the training and test performance (as described previously).\n2. The authors mention throughout the paper that many of the assumptions they make in their model hold for fMRI data (L162 for example). Hence, I would assume this paper to have more of a focus on actually showing why this model works well for fMRI data. However, I think the fMRI analysis lacks depth, and the results in Figure 3a do not look convincing because of how noisy the significant regions look. The noisiness is underlined by the results in Table A2: the IQ predictions do not seem to generalize well from the training to the test set. Additionally, the authors do not show whether their method actually identifies different areas than any of the baseline methods, which would help verify that their model is a helpful tool for the neuroimaging field compared to the baseline methods. Although I understand that the authors choose not to go too in-depth about the results since this is a machine learning conference, I think given the focus on fMRI in this paper the authors should expand a lot on their current discussion of what these results indicate (i.e. why does it make sense that these regions are related to IQ and parental education level?), atleast in the Appendix. Lastly, and this is not as important to address, but I think there are additional potentially more interesting questions one can ask using the ABCD dataset, i.e. regressions to mental health assessments [1].\n\nMinor weaknesses:\n1. The authors mention that there is a sparse signal in fMRI data, the authors should expand on this claim. Similarly, on Line 162 the authors claim that the signal is assumed to be sparse and piecewise smooth, but the authors should expand on this claim and potentially provide references. I agree generally with the statement, but I do think it’s important for the authors to explain why.\n2. The notation the authors use can be confusing, the authors use $I_d$ for an identity matrix of dxd dimensions, but also use $I(\\cdot)$ as the indicator function. I would encourage the authors to use $\\mathbb{1}_{\\\\{\\cdot\\\\}}$ as notation for the indicator function instead. Similarly, instead of using $\\mathcal{N(j)}$ for a neighborhood and $\\mathbf{N}(\\cdot)}$ for a Normal distribution, I would use $\\mathcal{N}$ for a Normal distribution and other notation for the neighborhood. Lastly, $p_j$ is not defined in the text on Line 131.\n3. The authors do not compare their model with T-LoHo on the data that is used in the T-LoHo paper. I assume because ST-CAR cannot handle non-smooth spatial flips from positive to negative values well (i.e. Line 161), but the authors do not mention this in the limitations section or show any results on these types of edge cases.\n4. The authors should also and more importantly report FLOPs instead of computation time because computation time is machine-dependent.\n5. The authors mention fMRI pre-processing  on Line 422, but do not explain in the paper what types of pre-processing they perform. Moreover, they do not explain how they select the subjects they use and how they split the data into a training and test set.\n6. The authors only show the positive effects in their fMRI results (Figure 3) and Table 2, which they should do, since it isn’t clear why only positive effects are expected for the hypotheses they test.\n7. The quality of Figure 3 is not up to par, it is hard to read the color bar numbers.\n\nSpelling/grammar:\n- L36-37: “…, and it can be time-series …, or in our particular interest, the functional Magnetic Resonance Imaging data … of human brain signal.” -> .This high-dimensional component can be a time-series, …, or in our case, functional magnetic resonance imaging (fMRI) data of the human brain.\n- L39/40: “…, formulated as a scalar-on-image (SonI) regression…” -> formulated as scalar-on-image (SonI) regression\n- L54/55: “… is still computationally expansive...” -> is still computationally expensive\n- L109/109: “…, referred to as the…” -> called the\n- L247: “… algorithm, but cannot directly give inference results …” -> algorithm, but cannot directly infer\n- L271/272: “… selection accuracy in SonI setting.” -> selection accuracy in the SonI setting.\n- L305: “…, we use ridge regression result …” -> we use the ridge regression result\n- L371/372: “… blow 10% …” -> below 10%\n\n[1] Chen, J., Tam, A., Kebets, V., Orban, C., Ooi, L. Q. R., Asplund, C. L., ... & Yeo, B. T. (2022). Shared and unique brain network features predict cognitive, personality, and mental health scores in the ABCD study. Nature communications, 13(1), 1-17.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "BIrHNGNLew", "reviewer": "Reviewer_85HH", "summary": "This paper introduces a novel prior, the Soft-Thresholded Conditional AutoRegressive (ST-CAR) prior, for handling high-dimensional imaging data in regression models, specifically for brain imaging studies. The ST-CAR prior adapts to data-driven spatial smoothness, controlling sparsity, and incorporating inclusion probabilities for selecting active voxels.\n\nThe ST-CAR prior is applied to two common regression models in neuroimaging: scalar-on-image (SonI) and image-on-scalar (IonS). In SonI, brain imaging data are used to predict scalar outcomes, such as cognitive scores, while in IonS, scalar predictors (e.g., parental education level) are used to predict imaging outcomes. The paper introduces variational inference (VI) algorithms, including Coordinate Ascent Variational Inference (CAVI) and Stochastic Subsampling VI (SSVI), to make the approach computationally feasible for large-scale datasets. Results from simulation studies and a real-data analysis using the Adolescent Brain Cognitive Development (ABCD) study demonstrate that the ST-CAR prior outperforms traditional methods in identifying active brain regions and offers improved computational efficiency.", "questions": "### Questions \n\n- Can you discuss how this method could be applied to other imaging modalities, such as structural MRI, DTI, and PET imaging? Would any assumptions need to be adjusted when handling these modalities instead of fMRI data?\n\n- Is it necessary to apply IonS and SonI regressions in such high-dimensional spaces? Are there atlas-based methods that could produce ROI-level features that are more interpretable? By using such features, one might make inferences at a more aggregated level. So, my question is: is it essential to perform voxel-level inferences to determine the relationships between brain ROIs and clinical outcomes in such high-dimensional data? I would like the authors to elaborate on the need to work in this high-dimensional space.\n\n- Clinical study data are often curated and carefully processed. Do you believe that ST-CAR could handle real imaging and clinical data that are frequently noisy and of lower quality?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors have made a comprehensive presentation of their method and the presented results seem convincing. They have implemented both synthetic experiments and also a real ones, which strengthens their method. Overall, I believe that this paper can have an impact on the neuroimaging domain and making inferences on high-dimensional imaging data in curated clinical studies.", "weaknesses": "### Weaknesses \n1. Dependence on Thresholding Parameter v. \n- While the paper claims that the ST-CAR prior is relatively insensitive to the choice of the thresholding parameter v, the results could still be impacted by this setting. A suboptimal threshold could lead to either excessive sparsity (losing important signal) or insufficient sparsity (resulting in noise). Without clear guidelines for tuning v, users may struggle with finding the right balance, especially across different datasets. I would suggest to implement a sensitivity analysis showing how different values of v affect the results, or propose a data-driven method for choosing v.\n\n2. Limited Applicability for High Signal-to-Noise Ratios (SNR)\n - The paper notes that the ST-CAR prior performs well in low-SNR settings, particularly with SonI models. However, in settings where SNR is high, the method might not be optimal compared to other methods that handle high-SNR data more effectively. This could limit the generalizability of the model to other types of imaging data or modalities with higher SNR.  Can you discuss potential challenges or modifications needed to apply ST-CAR to noisier, lower quality data. Have you tested or plan to test the method on less processed datasets? Can you propose specific preprocessing steps or model adjustments that might be necessary for such data?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "b7Uh8H5m3k", "reviewer": "Reviewer_1CfH", "summary": "This paper presents a Soft-Thresholded Conditional AutoRegressive (ST-CAR) prior for Bayesian regression in high-dimensional brain imaging data, demonstrating its advantages in capturing complex spatial dependencies and quantifying uncertainty. The authors validate the ST-CAR model's effectiveness through simulations and real data, while employing variational inference and stochastic subsampling to enhance computational efficiency.", "questions": "1. Given the model’s focus on identifying significant brain regions, could it be extended to analyze inter-regional connections or brain circuits in cognitive functions?\n\n﻿2. How sensitive is the ST-CAR model to different choices of the number of neighbors, the thresholding parameter and the decay rate? Do those hyperparameters stay robust in other datasets (e.g., HCP datasets)\n\n3. Would it be possible to consider structural connectivity among regions of interest (extracted from DTI) as the exposures or confounders？", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The ST-CAR prior introduces a novel method to handle spatial dependencies in high-dimensional imaging data, where the coefficient is assumed to be smooth and sparse across their spatial domain.\n\n2. The proposed framework for solving the problem achieves substantial improvements in computational speed, making it more feasible for high-dimensional applications.\n\n3. The framework is validated on both simulated data and real fMRI data from the ABCD study, demonstrating its effectiveness in identifying significant brain regions associated with cognitive development.", "weaknesses": "1. The ST-CAR model involves several tuning parameters (e.g., the number of neighbors, thresholding parameter and decay rate), it seems that they require careful designed to achieve optimal performance. \n\n2. The ST-CAR model effectively identifies significant regions of WM tasks in the ABCD dataset, but it seems that it cannot provide insights into broader brain circuits or networks, which are increasingly recognized as critical for understanding complex brain functions and disorders.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Pe1zTZN9NY", "reviewer": "Reviewer_NeTL", "summary": "In this work, the authors present a novel prior for use in bayesian image regression from functional magnetic resonance imaging (fMRI). Motivated by limitations of previous priors which require the prior specification of a gaussian process kernel, which can lead to bias and other issues in the regression unless the true signal is known beforehand, the authors provide a novel prior which relaxes smoothness assumptions on functional parameters. The authors demonstrate that this novel prior autoperforms previous priors by a significant margin.", "questions": "All of the questions I have regard the performance over multiple repetitions and would be addressed simply by providing those standard deviations. Every other question I had was substantially addressed in the appendix.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "I find this work to be altogether very strong work with some omissions in the experimentation that should be easily addressed. At its core this work is well-motivated by intuition regarding the shortcomings of previous priors - the relaxation of the smoothness assumption is easy to grasp, and provides a firm foundation for the reader to follow the following, more dense theoretical work. I found no issues with the theoretical work, and appreciated the rigor taken by the authors.\n\nThe inclusion of substantial numerical and real-data experiments is a huge boon to this work. The authors do a substantial amount of experimental validation here, and demonstrate that this novel prior outperforms existing priors. \n\nThe authors also include a substantial appendix with additional work which I think further solidifies this as a strong paper. \n\nWell done!", "weaknesses": "The only minor quibble I have with this work is the lack of standard deviations in most of the results. For example, in table 1 you average over 100 replications, but don't provide any standard deviation measures. These should be provided to fully solidify the performance of your approach over other methods.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)