The reviewed record of science sign in
Pith

arxiv: 2606.29114 · v1 · pith:2XXEUGTB · submitted 2026-06-27 · stat.ME · math.ST· stat.ML· stat.TH

Multivariate Varying-Coefficient BART with Graphical Horseshoe Priors

Soham Ghosh , Sameer K. Deshpande This is my paper

Reviewed by Pith T0 review T1 audit T2 compute T3 formal T4 kernel 2026-06-30 08:13 UTCgrok-4.3pith:2XXEUGTBrecord.jsonopen to challenge →

classification stat.ME math.STstat.MLstat.TH
keywords multivariate regressionBayesian additive regression treesvarying coefficient modelsgraphical horseshoe priorposterior contraction ratesresidual dependenceBayesian nonparametric methodshigh-dimensional data
0
0 comments X

The pith

A new multivariate BART model with independent ensembles per coefficient and a graphical horseshoe prior on residuals achieves the first posterior contraction rates for jointly estimated dependence.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper develops multiVCBART to address multivariate regression settings where multiple outcomes have nonlinear, outcome-specific effects that vary with modifiers and where residual dependencies among outcomes carry scientific information. It represents each entry of the coefficient matrix with its own BART ensemble so that tree structures need not be shared across outcomes, while a graphical horseshoe prior on the precision matrix encourages sparse conditional dependence. An auxiliary-variable sampler converts the joint multivariate Gaussian likelihood into a sequence of independent scalar pseudo-response regressions, allowing separate updates for the trees and the precision matrix. The central theoretical result is the first set of posterior contraction rates for any multivariate BART procedure that estimates residual dependence jointly, with the rates adapting near-minimax to the unknown smoothness of the coefficient surfaces and the sparsity pattern of the precision matrix.

Core claim

multiVCBART represents each entry of the coefficient matrix B(x) by an independent BART ensemble and places a Graphical Horseshoe prior on the residual precision matrix Ω; the accompanying sampler reduces the multivariate likelihood to a sequence of scalar pseudo-response updates that decouple tree backfitting from the graphical horseshoe step. This construction yields the first posterior contraction rates for a multivariate BART model with jointly estimated residual dependence, establishing near-minimax adaptation to underlying smoothness and structural sparsity.

What carries the argument

Independent BART ensembles for each entry of the coefficient matrix B(x), paired with a Graphical Horseshoe prior on the residual precision matrix Ω, made tractable by a pseudo-response sampler that converts the joint multivariate Gaussian likelihood into separate scalar regressions.

If this is right

  • The model outperforms existing multivariate tree-based methods and Bayesian seemingly unrelated regression competitors on sparse high-dimensional data.
  • Application to the Genomics of Drug Sensitivity in Cancer data recovers distinct biomarker signals per outcome together with a coherent residual pharmacologic network.
  • Outcome-specific coefficient surfaces can be estimated without forcing shared tree architecture across responses.
  • Posterior inference jointly adapts to smoothness of the mean functions and sparsity of the residual precision matrix.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The pseudo-response trick may generalize to other non-Gaussian multivariate likelihoods if the auxiliary-variable representation can be maintained.
  • The separation of mean modeling from dependence modeling suggests the framework could be combined with other sparse precision priors without altering the tree-update step.
  • Empirical gains on genomics data indicate potential utility in any domain with multiple correlated responses whose predictors differ across outcomes.

Load-bearing premise

The reduction of the multivariate Gaussian likelihood to a sequence of scalar pseudo-response updates preserves the correct joint posterior without introducing bias or approximations that change the contraction rates.

What would settle it

A simulation on a small multivariate problem in which the full joint posterior (sampled directly) differs systematically from the posterior produced by the pseudo-response updates, or a data set where the empirical posterior fails to contract at the claimed near-minimax rate.

Figures

Figures reproduced from arXiv: 2606.29114 by Sameer K. Deshpande, Soham Ghosh.

Figure 1
Figure 1. Figure 1: Schematic comparison of tree partitions. [PITH_FULL_IMAGE:figures/full_fig_p006_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Schematic comparison of the classical GGV sieve argument and our two-sieve shelling [PITH_FULL_IMAGE:figures/full_fig_p015_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: High-dimensional Friedman–SUR benchmark ( [PITH_FULL_IMAGE:figures/full_fig_p017_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Estimated residual conditional dependency network among the seven GDSC drug re [PITH_FULL_IMAGE:figures/full_fig_p020_4.png] view at source ↗
read the original abstract

Modern multivariate regression problems involve several related outcomes whose regression effects are not only nonlinear, heterogeneous, and outcome-specific, but also where the residual dependence among outcomes is scientifically meaningful. Existing multivariate Bayesian tree-based methods typically address only part of this problem: some impose substantial sharing of tree architecture across outcomes, which is overly restrictive when responses depend on distinct predictors or effect modifiers, while others accommodate residual dependence but retain simpler mean structures. This paper develops multiVCBART, a multivariate varying-coefficient Bayesian additive regression tree framework that jointly models flexible outcome-specific coefficient surfaces and a sparse residual precision matrix. Each entry of the coefficient matrix $B(x)$ is represented by an independent BART ensemble, allowing predictor effects to vary nonlinearly with modifiers $x$ across outcomes, while a Graphical Horseshoe prior on the precision matrix $\Omega$ captures parsimonious residual conditional dependence. To permit efficient computation, we introduce a sampler that reduces the multivariate Gaussian likelihood to a sequence of scalar pseudo-response updates, decoupling the tree backfitting from the Graphical Horseshoe step. Theoretically, we establish the first posterior contraction rates for a multivariate BART model with jointly estimated residual dependence, proving near-minimax adaptation to underlying smoothness and structural sparsity. Empirically, multiVCBART outperforms existing multivariate tree models and Bayesian SUR competitors on sparse, high-dimensional datasets. Finally, in a re-analysis of the Genomics of Drug Sensitivity in Cancer dataset, our method identifies distinct biomarker signals and recovers a coherent residual pharmacologic network.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 1 minor

Summary. The manuscript proposes multiVCBART, a multivariate varying-coefficient BART model in which each entry of the coefficient matrix B(x) is modeled by an independent BART ensemble to allow outcome-specific nonlinear surfaces, paired with a Graphical Horseshoe prior on the residual precision matrix Ω to capture sparse conditional dependence. A new sampler is introduced that reduces the multivariate Gaussian likelihood to scalar pseudo-response updates, decoupling tree backfitting from the precision-matrix step. The central theoretical contribution is the first posterior contraction rates for a multivariate BART model with jointly estimated residual dependence, claimed to achieve near-minimax adaptation to smoothness and structural sparsity; empirical results on sparse high-dimensional data and a re-analysis of the GDSC dataset are also reported.

Significance. If the sampler targets the exact joint posterior, the contraction-rate results would constitute a meaningful advance by supplying the first such guarantees for multivariate BART with joint residual-dependence estimation and by establishing near-minimax adaptation to both smoothness and sparsity. The empirical outperformance and real-data application would then provide supporting practical evidence.

major comments (2)
  1. [Abstract and methodological sampler section] Abstract and sampler description: the reduction of the multivariate Gaussian likelihood to a sequence of scalar pseudo-response updates is presented as enabling efficient computation while preserving the target posterior, yet no explicit verification is supplied that the pseudo-responses are exact conditional draws (rather than conditional expectations or fixed-point approximations). Because the contraction rates are asserted for the joint posterior p(B, Ω | data), any bias or approximation in this step would render the rates inapplicable to the claimed target.
  2. [Theoretical results section] Theoretical results section: the near-minimax adaptation claim for both smoothness and structural sparsity is load-bearing for the paper's primary contribution, but the argument must explicitly incorporate the Graphical Horseshoe prior on Ω and confirm that the independent BART ensembles for the entries of B(x) do not alter the rate derivations; without this linkage the adaptation statement remains unanchored.
minor comments (1)
  1. [Abstract] The abstract states that the method 'outperforms existing multivariate tree models and Bayesian SUR competitors' but does not name the specific competitors or report the quantitative metrics and error-bar summaries used in those comparisons.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for their detailed and constructive feedback on our manuscript. Below we provide point-by-point responses to the major comments and outline the revisions we plan to make.

read point-by-point responses

  1. Referee: [Abstract and methodological sampler section] Abstract and sampler description: the reduction of the multivariate Gaussian likelihood to a sequence of scalar pseudo-response updates is presented as enabling efficient computation while preserving the target posterior, yet no explicit verification is supplied that the pseudo-responses are exact conditional draws (rather than conditional expectations or fixed-point approximations). Because the contraction rates are asserted for the joint posterior p(B, Ω | data), any bias or approximation in this step would render the rates inapplicable to the claimed target.

    Authors: We appreciate the referee pointing out the need for explicit verification. Upon review, the current manuscript describes the sampler as preserving the target posterior but does not include a formal proof of exactness for the pseudo-response step. In the revised version, we will add a dedicated subsection or lemma proving that the scalar pseudo-response updates correspond to exact conditional draws from the joint posterior, ensuring the contraction rates apply to the correct target distribution. This will involve deriving the conditional posterior for each component under the multivariate normal model. revision: yes

  2. Referee: [Theoretical results section] Theoretical results section: the near-minimax adaptation claim for both smoothness and structural sparsity is load-bearing for the paper's primary contribution, but the argument must explicitly incorporate the Graphical Horseshoe prior on Ω and confirm that the independent BART ensembles for the entries of B(x) do not alter the rate derivations; without this linkage the adaptation statement remains unanchored.

    Authors: We agree that the theoretical contribution would be strengthened by making these linkages explicit. The current proof sketch relies on the properties of the Graphical Horseshoe for sparsity and univariate BART rates for the mean functions, but does not spell out the combination in detail. In the revision, we will expand the theoretical results section to explicitly show how the Graphical Horseshoe prior enables the structural sparsity adaptation in the joint model and to confirm that modeling the entries of B(x) with independent BART ensembles preserves the contraction rates derived from the univariate case, as the dependence is captured solely through Ω. revision: yes

Circularity Check

0 steps flagged

No significant circularity; contraction rates derived from standard BART theory applied to new model

full rationale

The paper introduces multiVCBART with independent BART ensembles per coefficient entry and a Graphical Horseshoe prior on Ω, then presents a sampler reducing the multivariate likelihood to scalar pseudo-response updates. The claimed posterior contraction rates are established as the first such results for this class of model, adapting near-minimax rates to smoothness and sparsity. No step reduces a prediction or rate to a fitted parameter by construction, nor does any load-bearing premise rest solely on self-citation. The sampler is asserted to target the exact joint posterior p(B, Ω | data), allowing the theoretical analysis to apply directly without auxiliary approximations that would invalidate the rates. This is a standard non-circular theoretical development.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

The central claims rest on standard BART assumptions (tree priors, Gaussian errors) plus the new Graphical Horseshoe prior and the validity of the pseudo-response reduction; no invented entities are introduced.

axioms (2)
  • domain assumption The multivariate responses follow a Gaussian likelihood conditional on the coefficient surfaces and precision matrix.
    Invoked to justify the likelihood and the pseudo-response sampler.
  • domain assumption The Graphical Horseshoe prior induces the desired sparsity on the precision matrix.
    Standard property of the prior used to capture parsimonious residual dependence.

pith-pipeline@v0.9.1-grok · 5805 in / 1477 out tokens · 26378 ms · 2026-06-30T08:13:23.803914+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

18 extracted references · 5 canonical work pages · 1 internal anchor

  1. [1]

    Agapiou, S., Castillo, I., and Egels, P. (2026). Heavy-tailed and Horseshoe priors for regression and sparse Besov rates.arXiv:2505.15543. Ando, T. and Zellner, A. (2010). Hierarchical Bayesian analysis of the seemingly unrelated re- gression and simultaneous equations models using a combination of direct Monte Carlo and importance sampling techniques.Bay...

    work page arXiv 2026

  2. [2]

    R., and Chen, Y

    Bai, R., Boland, M. R., and Chen, Y. (2023). Scalable high-dimensional Bayesian varying coeffi- cient models with unknown within-subject covariance.Journal of Machine Learning Research, 24(259):1–49. Baldwin, S. A., Bauer, D. J., Stice, E., and Rohde, P. (2014). Analyzing multiple outcomes in clinical research using multivariate multilevel models.Journal ...

    2023

  3. [3]

    J., Vannucci, M., and Fearn, T

    Brown, P. J., Vannucci, M., and Fearn, T. (1998). Multivariate bayesian variable selection and prediction.Journal of the Royal Statistical Society: Series B, 60(3):627–641. Cai, T. T., Ren, Z., and Zhou, H. H. (2016). Estimating structured high-dimensional covariance and precision matrices: Optimal rates and adaptive estimation.Electronic Journal of Stati...

    1998

  4. [4]

    A., George, E

    21 Chipman, H. A., George, E. I., and McCulloch, R. E. (2010). BART: Bayesian Additive Regression Trees.Annals of Applied Statistics, 4(1):266–298. Clark, T. E., Huber, F., Koop, G., Marcellino, M., and Pfarrhofer, M. (2023). Tail forecasting with multivariate Bayesian additive regression trees.International Economic Review, 64(3):979–1022. Deshpande, S. ...

    2010

  5. [5]

    Garnett, M., Edelman, E., Heidorn, S., et al. (2012). Systematic identification of genomic markers of drug sensitivity in cancer cells.Nature, 483(7391):570–575. Ghosal, S., Ghosh, J. K., and van der Vaart, A. W. (2000). Convergence rates of posterior distri- butions.Annals of Statistics, 28(2):500–531. Ghosh, S., Bhogale, S., and Deshpande, S. K. (2025)....

    work page arXiv 2012

  6. [6]

    Semiparametric Mixed-Scale Models Using Shared Bayesian Forests

    Hastie, T. and Tibshirani, R. (1993). Varying-coefficient models.Journal of the Royal Statistical Society: Series B, 55(4):757–796. Hoover, D. R., Rice, J. A., Wu, C. O., and Yang, L.-P. (1998). Nonparametric smoothing estimates of time-varying coefficient models with longitudinal data.Biometrika, 85(4):809–822. Huang, A. and Wand, M. P. (2013). Simple ma...

    work page internal anchor Pith review Pith/arXiv arXiv 1993

  7. [7]

    Yeung, C. C. S., Egan, D., and Radich, J. P. (2016). Molecular monitoring of chronic myeloid leukemia: present and future.Expert Review of Molecular Diagnostics, 16(10):1083–1091. Zellner, A. (1962). An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias.Journal of the American Statistical Association, 57(298):348...

    2016

  8. [8]

    For Gaussian measures, the squared Hellinger distance admits the factorization (Pardo, 2005, Chap- ter 1): h2 Nq(η1,i,Σ 1),N q(η2,i,Σ

    in eFn, and writeδ i :=η 1,i −η2,i and ¯Σ := 1 2(Σ1+Σ2). For Gaussian measures, the squared Hellinger distance admits the factorization (Pardo, 2005, Chap- ter 1): h2 Nq(η1,i,Σ 1),N q(η2,i,Σ

    2005

  9. [9]

    The mapF(Σ) = log det(Σ) has Hessian∇ 2F(Σ)[H, H] =−tr(Σ −1HΣ −1H), hence ∇2F(Σ)[H, H] ≤ ∥Σ −1∥2 op ∥H∥2 F

    Averaging overiyields 1 n nX i=1 Tµ,i ≤ 1 8 RΩ,n ∥η1 −η 2∥2 2,n.(S1.10) (ii) Determinant term.Write Adet := |Σ1|1/4|Σ2|1/4 |¯Σ|1/2 = exp n 1 4 log|Σ 1|+ 1 4 log|Σ 2| − 1 2 log| ¯Σ| o . The mapF(Σ) = log det(Σ) has Hessian∇ 2F(Σ)[H, H] =−tr(Σ −1HΣ −1H), hence ∇2F(Σ)[H, H] ≤ ∥Σ −1∥2 op ∥H∥2 F . On eFn,∥Σ −1∥op =∥Ω∥ op ≤R Ω,n, so for all Σ on the line segmen...

    1977

  10. [10]

    By A1,B 0,jr(x) =B 0,jr(xJ0,jr ) andB 0,jr ∈ H αjr ([0,1] d0,jr ;K)

    Let us fixr∈[q] andj∈S B0,r. By A1,B 0,jr(x) =B 0,jr(xJ0,jr ) andB 0,jr ∈ H αjr ([0,1] d0,jr ;K). Partition [0,1] d0,jr intom d0,jr jr congruent cubes and letB jr,δ be the cell-average step function. Standard H¨ older approximation following the arguments in Roˇ ckov´ a and van der Pas (2019, Lemma 3.2) gives ∥Bjr,δ −B 0,jr ∥∞ ≤C m −αjr jr , for a constan...

    2019

  11. [11]

    Consider any noise leafb jrtℓ withj /∈S B0,r. Marginalizing overλ jr (and noting truncatingλ jr ≤n AB only decreases tails), we have P |bjrtℓ |> u noise |τ B =τ ≤C HS τ unoise √ M log 1 + unoise √ M τ ≤C HS t0 unoise √ M log 1 + unoise √ M t0 . With the choice oft 0 in our lemma, we haveu noise √ M /t0 =C 0QBM L0 log(eQB), hence log(1 + unoise √ M /t0)≤2 ...

    2019

  12. [12]

    Therefore: ∥ηnoise(xi)∥2 2 ≤pD 2M X (j,r)/∈bSB MX t=1 LjrtX ℓ=1 b2 jrtℓ =pD 2M EB,noise(b)

    Furthermore,B jr(xi)2 = (P t gjrt(xi))2 ≤ MP t g2 jrt(xi)≤M P t,ℓ b2 jrtℓ. Therefore: ∥ηnoise(xi)∥2 2 ≤pD 2M X (j,r)/∈bSB MX t=1 LjrtX ℓ=1 b2 jrtℓ =pD 2M EB,noise(b). 53 By the sieve constraint,E B,noise(b)≤∆ B,n = r2 n 16D2p2qM . Hence: ∥ηnoise(xi)∥2 2 ≤pD 2M r2 n 16D2p2qM ≤ r2 n 16 . Averaging overnyields∥η noise∥2,n ≤r n/4. Thus, by covering onlyη sig ...

    2000

  13. [13]

    Upper Bounding the Truncated NumeratorN F n .To evaluate the numerator strictly over the sieve, we utilize the existence of a global test functionϕ n from Ghosal et al

    Thus, on the eventE n, the truncated denominator satisfies DF n ≥exp(−C Dn(ε† n)2) forC D =C prior +C. Upper Bounding the Truncated NumeratorN F n .To evaluate the numerator strictly over the sieve, we utilize the existence of a global test functionϕ n from Ghosal et al. (2000, Section 7), which separates Θ 0 from{Θ∈ F n :H(Θ,Θ 0)> M ε † n}. We bound the ...

    2000

  14. [14]

    This establishes the Hellinger contraction rateε † n for the posterior distribution under the truncated prior on the effective sieveF n

    Consequently, the ratio ΠF n (Un |Y) =N F n /DF n →0 inP Θ0-probability asn→ ∞. This establishes the Hellinger contraction rateε † n for the posterior distribution under the truncated prior on the effective sieveF n. To establish testing power, we must prove that as the noise indices (k, m) grow, the Hellinger dis- tance from the true parameter Θ0 strictl...

    2000

  15. [15]

    Throughout, we suppress MCMC iteration superscripts for readability. S2.1 Outcome-wise pseudo-response representation The key computational simplification is that, conditional on the remaining outcomes and on Ω, the multivariate Gaussian likelihood reduces to a scalar Gaussian working likelihood. Fixr∈ {1, . . . , q}. WritingE i = (Ei1, . . . , Eiq)⊤, we ...

    2016

  16. [16]

    Similar to the first setting, again we takez i =1.The mean vectorη 0(xi) = (η 0,1(xi), . . . , η0,10(xi))⊤ is 65 defined by ten nonlinear component functions: η0,1(xi) = 4 sin πxi1xi2 , η 0,2(xi) = 3 cos πxi3 , η0,3(xi) = 4(x2 i4 −0.33), η 0,4(xi) = 3xi5, η0,5(xi) = 4 exp −2x 2 i6 , η 0,6(xi) = 3xi7xi8, η0,7(xi) = 4 |xi9| −0.5 , η 0,8(xi) = 3 sin πxi10 , ...

    work page arXiv 2023

  17. [17]

    The first five predictors are continuous withX 1, . . . , X5 iid∼ U(0,1),the next three are binary,X 6, X7, X8 iid∼Bern(0.5),the next two are ordinal categorical, X9, X10 ∈ {0,1,2,3,4},and the remaining 40 predictors are independent noise covariates generated fromU(0,1). Treatment assignment followsT|X∼Bern(X 4),so thatX 4 acts as a confounder. The two po...

    2023

  18. [18]

    σB RMSEtest CRPStest PI Coveragetest Time (s) 0.22.353 1.2620.937 126.124 1 2.493 1.341 0.94077.974 5 2.490 1.3480.942125.093 70

    work page arXiv