Two independent defects in Horseshoe::GlobalUpdate()
have been fixed. Both affected every fit using
prior_type = "horseshoe" (including
HorseTrees() and CausalHorseForest()). Results
from earlier versions will not reproduce under this release.
global_hp had no effect. The scale
entered the auxiliary draw as
1.0 / alpha_global * alpha_global, which by operator
precedence is the constant 1.0. The user’s
global_hp was silently discarded, so every fit behaved as
though global_hp = 1 regardless of what was supplied. This
was reported by a referee of the accompanying manuscript.
The error scale and the forest scale were
swapped. EtaPrior::GlobalUpdate() declares
(..., sigma, omega, ...) but was being called with
(omega, sigma). Both are const double&, so
the swap compiled silently. The tau draw therefore divided
by sigma where it should have divided by
omega. The effect is small where omega = 1
(HorseTrees(), ShrinkageTrees()) and
substantial where omega = 1/2 (the causal forests), where
it grew with the number of covariates.
The old default of k = 0.1 was chosen while defect (1)
was present, so it was tuned for a model in which the global scale was
pinned at 1. With the defect corrected, that value shrinks far too
aggressively: in our simulations it drove pointwise coverage of the
conditional treatment effect from roughly 0.95 down to about 0.4.
Re-running existing scripts against this release without
updating k is therefore worse than not upgrading at
all.
The defaults are now calibrated for the corrected sampler:
| function | old | new |
|---|---|---|
HorseTrees() |
k = 0.1 |
k = 1.0 |
CausalHorseForest() |
k = 0.1 |
k = 1.5 |
ShrinkageTrees() |
local_hp, global_hp required |
default to 1.0 / sqrt(number_of_trees) |
CausalShrinkageForest() |
four scales required | default to 1.5 / sqrt(number_of_trees_*) |
Values of k between roughly 0.5 and 1.5 (single-forest
models) or 1 and 2 (causal models) worked well across a range of
simulated settings, with smaller values shrinking more aggressively and
larger values being more conservative. The defaults sit in the middle of
each range. The causal range is higher because the treatment forest
enters as b * tau(x) with b = +/- 1/2, halving
its contribution to the response; measured on the scale of that
contribution the two recommendations nearly coincide. k
remains exposed on all four functions and is a natural target for
cross-validation.
The local_hp / global_hp arguments are also
no longer mandatory for prior_type = "horseshoe"; either
may be supplied alone, with the other taking its default. Only their
product identifies the prior, so setting them equal loses nothing.
prior_type = "horseshoe_fw" is unchanged and still requires
both to be given explicitly.
ovarian
rebuilt, ovarian_truth added (breaking)left_time and
right_time columns and a duplicated
year_of_diagnosis.1; the frame is now 357 x 1004. Old fits
leaked the outcome, and positional column selection now picks different
covariates.ovarian_truth: the prognostic surface, treatment
effect, propensity score and uncensored event time behind the
semi-synthetic cohort. Its times are in months;
ovarian$OS_time is in days.year_of_diagnosis was documented as an instrumental
variable; it is a confounder. Wording only.predict()
passed the training design in the wrong memory orderpredict() on a ShrinkageTrees object handed
the stored training design matrix straight to the C++ sampler. The C++
layer reads designs as a flat row-major buffer, but R
coerces a matrix to a vector column-major, so the
training data was transposed inside the sampler on every call whenever
n != p. The test design on the same call was flattened
correctly, so the two disagreed. It ran and returned plausible
numbers.
predict() on a CausalShrinkageForest object
was unaffected: it already flattened both designs correctly.
posterior_projection()Lower-dimensional posterior summarisation following Woody, Carvalho and Murray (2021). Every posterior draw of the fitted function, or of the treatment-effect and prognostic surfaces for causal fits, is projected onto a simpler summary model: a linear model (unpenalised, ridge, lasso, or elastic net), a spline-additive model, or a shallow CART tree. Returns a posterior over summaries with credible intervals, plus the posterior of the summary R-squared.
The penalised linear families are the exception: they project the posterior mean of the surface once and report point estimates, with no intervals. The spread of a penalised solution across draws has no coverage guarantee for the projection of the true surface, and the implied inclusion frequencies read as selection probabilities while being driven by the penalty.
glmnet and rpart are new suggested
packages.
right_time onto left_time because the C++
layer needs a finite boundary, but that collapse happened
before the mean and scale were estimated. Those come from
survival::Surv(left_time, right_time, type = "interval2"),
under which left == right denotes an exact event, so every
right-censored subject was treated as observed. The centring mean and
the sigma prior were biased downward, by more the heavier the censoring.
All four interval-censored paths were affected:
HorseTrees(), ShrinkageTrees(),
CausalHorseForest() and
CausalShrinkageForest(). SurvivalBCF() and
SurvivalShrinkageBCF() were not. Interval-censored fits
containing right-censored observations will change.predict() silently dropped the Dirichlet splitting rule
and re-ran under the standard prior, affecting every
SurvivalDART() and SurvivalShrinkageBCF() fit.
Predictions from those models will change.HorseTrees()
base (0.95, not 1.05) and q (0.90, not 1.00);
CausalShrinkageForest() listed three of seven
prior_type values; plot(type = "path") shows
the cross-validation curve, not the summary R-squared.stats, utils and
splines in Imports.There is no setting that reproduces pre-2.1.0 results exactly,
because the sampler itself has changed. The closest correspondence is
that the old code was equivalent to the corrected code with
global_hp = 1, since the global scale was inert.
summary(), plot(type = "ate"), and
predict() for CausalShrinkageForest (and
CausalHorseForest) now default to a Bayesian-bootstrap
posterior for the average treatment effect: at each MCMC iteration the
per-observation CATEs are reweighted with Dirichlet(1, …, 1) weights
before being summed, giving a draw from the posterior of the
population ATE (PATE). Credible intervals are correspondingly
wider than before because they now propagate uncertainty in the
covariate distribution, not only in tau(x).
bayesian_bootstrap = TRUE on
summary.CausalShrinkageForest(),
plot.CausalShrinkageForest() (for
type = "ate"), and
predict.CausalShrinkageForest(). Set to FALSE
to recover the previous equal-weight mixed ATE (MATE).predict() now also returns an ate summary
for the new data and retains the posterior CATE sample matrix as
cate_samples.bayesian_bootstrap_ate() returns
PATE and MATE summaries (means, CIs, and full posterior draws) from
either a fitted CausalShrinkageForest or a
CausalShrinkageForestPrediction.This is a breaking change for printed/plotted numerics:
existing scripts will report wider CIs than before. Use
bayesian_bootstrap = FALSE to reproduce previous
output.
ovarian dataset
restructuredThe ovarian dataset is now a single data frame
(previously a list with $clinical and $X
elements). Clinical columns (OS_time,
OS_event, treatment, age,
figo_stage, tumor_grade) and the 2000 gene
expression columns are combined into one data frame with 2006 columns.
This simplifies data access and aligns the format with the
pdac dataset.
Code that previously used ovarian$clinical or
ovarian$X must be updated — see ?ovarian for
the new structure.
plot(fit, type = "ate") and
plot(fit, type = "cate") for
CausalShrinkageForest models incorrectly subtracting the
control forest predictions from the treatment forest predictions.ovarian)Added the ovarian dataset: a processed TCGA-OV cohort (n
= 357) for high-dimensional survival prediction and causal
inference.
See ?ovarian and examples/test-ovarian.R
for a full worked example covering survival prediction (SurvivalBART,
SurvivalDART, HorseTrees) and causal inference (SurvivalBCF,
SurvivalShrinkageBCF, CausalHorseForest).
treatment_coding)All causal model functions — CausalHorseForest(),
CausalShrinkageForest(), SurvivalBCF(), and
SurvivalShrinkageBCF() — now accept a
treatment_coding argument controlling how the treatment
indicator enters the BCF decomposition y = f(x) + b * tau(x) + epsilon.
Four options are available:
"centered" (default): b_i in {-1/2, 1/2}. This is the
original behaviour."binary": b_i in {0, 1}. Standard binary coding."adaptive": b_i = z_i - e_hat(x_i), where e_hat(x_i) is
the estimated propensity score. This follows Hahn, Murray & Carvalho
(2020) and is implemented in the bcf package. Requires a
propensity vector."invariant": Parameter-expanded (invariant) treatment
coding. The coding parameters b_0 and b_1 are assigned N(0, 1/2) priors
and estimated within the Gibbs sampler via conjugate normal updates,
yielding a parameterisation that is invariant to the coding of the
treatment indicator (Hahn et al., 2020, Section 5.2). The treatment
effect is tau(x) = (b_1 - b_0) * tau_tilde(x). Posterior draws of b_0
and b_1 are returned in the fitted object.The predict() method for
CausalShrinkageForest objects automatically carries forward
the treatment coding used at training time. A
propensity_test argument is available for supplying
test-set propensity scores (defaults to 0.5).
All survival-capable functions now support
interval-censored data in addition to right-censored
data. Supply left_time and right_time vectors
(with outcome_type = "interval-censored") instead of
y and status. Three censoring types are
distinguished:
left_time == right_time.left_time < right_time.right_time = Inf.This convention follows
survival::Surv(type = "interval2"). Censored event times
are augmented within the AFT Gibbs sampler. The following functions are
affected:
HorseTrees(), ShrinkageTrees()
(single-forest models)CausalHorseForest(),
CausalShrinkageForest() (causal models)SurvivalBART(), SurvivalDART(),
SurvivalBCF(), SurvivalShrinkageBCF()
(survival wrappers)n_chains)All four primary model-fitting functions —
ShrinkageTrees, HorseTrees,
CausalHorseForest, and CausalShrinkageForest —
now accept an n_chains argument (default 1).
When n_chains > 1:
parallel::mclapply (falls back to sequential execution on
Windows).min(n_chains, parallel::detectCores()).N_post * n_chains total draws.ShrinkageTrees or
CausalShrinkageForest instance, so all existing
print, summary, and predict
methods work without modification.SurvivalBART,
SurvivalDART, SurvivalBCF, and
SurvivalShrinkageBCF inherit n_chains support
through ....print and summary output adapts
automatically: single-chain models show Posterior draws,
multi-chain models show Chains and Draws per chain,
with per-chain acceptance ratios listed separately.ShrinkageTrees and
CausalShrinkageForest with constructors in
constructors.R.print methods for both classes, displaying model
specification, MCMC settings, acceptance ratio, and posterior mean
sigma.summary methods for both classes, returning an
inspectable object with posterior sigma (mean, SD, 95% CI), prediction
summaries, variable importance (posterior inclusion probabilities), and
— for causal models — ATE with credible interval (when
store_posterior_sample = TRUE) and CATE heterogeneity.predict method for ShrinkageTrees,
enabling posterior predictive inference on new data by re-running the
sampler with stored training data and hyperparameters. Returns a
ShrinkageTreesPrediction object with posterior mean and
credible interval bounds. For survival models, the prediction object
additionally stores predictions_sample (full posterior
draws on the original scale) and sigma (posterior draws on
the log-time scale), enabling posterior predictive survival curve
plotting.predict method for
CausalShrinkageForest, returning a
CausalShrinkageForestPrediction object with three
components: prognostic (\(\mu(X)\)), cate (\(\tau(X)\)), and total (\(\mu(X) + \tau(X)\)), each with posterior
mean and credible interval bounds. For survival models with
timescale = "time", predictions are back-transformed to the
original time scale and the CATE becomes a multiplicative time
ratio.print and summary methods for
ShrinkageTreesPrediction and
CausalShrinkageForestPrediction.coda)as.mcmc.list() S3 method for
ShrinkageTrees objects, converting the sigma posterior
(split by chain) into a coda::mcmc.list. This enables all
standard coda diagnostics: Gelman–Rubin R-hat, effective
sample size, Geweke test, Heidelberger–Welch test, autocorrelation
plots, and more.summary() now automatically reports effective
sample size (ESS) and — for multi-chain fits — the
Gelman–Rubin R-hat when the suggested package
coda is installed.coda to Suggests in
DESCRIPTION.plot)S3 plot() methods added for ShrinkageTrees,
CausalShrinkageForest, and
ShrinkageTreesPrediction. Requires the suggested package
ggplot2.
plot(fit, type = "trace") — sigma traceplot; one line
per chain, useful for assessing mixing.plot(fit, type = "density") — overlaid posterior
density of sigma, one curve per chain.plot(fit, type = "vi") — posterior credible intervals
for variable inclusion probabilities (top n_vi
predictors).plot(fit, type = "ate") — posterior density of the ATE
with 95 % credible region (causal models only; requires
store_posterior_sample = TRUE).plot(fit, type = "cate") — point estimates and 95 %
credible intervals for the CATE of each training observation, sorted by
posterior mean (causal models only; requires
store_posterior_sample = TRUE).plot(fit, type = "vi", forest = "both") — side-by-side
VI for the control and treatment forests (causal models
only).plot(fit, type = "survival") — posterior survival
curves \(S(t | x_i) = 1 - \Phi((\log t -
\mu_i) / \sigma)\) derived from the AFT log-normal model
(survival outcomes only).
obs = NULL): computes \(\bar{S}(t) = n^{-1} \sum_i S(t | x_i)\) at
each MCMC iteration with pointwise credible bands.obs = c(1, 5, ...)): one curve per selected training
observation with its own credible band.level controls the credible band width (default
0.95).t_grid allows a custom time grid; auto-generated if
NULL.km = TRUE overlays the Kaplan–Meier estimate as a
dashed black step function (population-averaged plot only; requires
survival package). Ignored with a message when
obs is not NULL.plot(pred, type = "survival") — posterior
predictive survival curves for new (out-of-sample) data
from predict(). Same obs, t_grid,
and level arguments as above. The KM overlay is not
available for prediction objects.HorseTrees,
ShrinkageTrees, SurvivalBART,
SurvivalDART, SurvivalBCF,
SurvivalShrinkageBCF, CausalHorseForest,
CausalShrinkageForest), all S3 methods (print,
summary, predict, plot),
multi-chain MCMC, and a full TCGA PAAD case study.SurvivalBART(), SurvivalDART(),
SurvivalBCF(), and SurvivalShrinkageBCF() now
accept store_posterior_sample as an explicit parameter
(default TRUE), avoiding a “matched by multiple actual
arguments” error when passing it via ....CausalHorseForest and
CausalShrinkageForest failing with
"argument 'y_train' is missing" when called directly
(broken constructor call introduced in 2.0.0 S3 refactor).plot(..., type = "vi") crashing with
"argument must be coercible to non-negative integer" in all
four model functions: covariate matrices were being stored as flat
numeric vectors instead of matrices, making ncol() return
NULL.HorseTrees and
ShrinkageTrees where sigma_hat,
y_mean, and lambda were not initialised in the
binary (probit) branch.SurvivalBCF wrapper for AFT-based Bayesian Causal
Forests.SurvivalShrinkageBCF with Dirichlet structural
sparsity.SurvivalDART for sparse high-dimensional survival
modeling.SurvivalBART for survival modeling using standard
BART.....std::map structure with a more
efficient vector-based lookup, improving overall computational speed by
approximately 30%.standard prior type option for
ShrinkageTrees, corresponding to the conventional BART
implementation without reversible jumps.CONTRIBUTING.md file with guidelines for
contributors🎉 First CRAN release of ShrinkageTrees!
This package provides Bayesian regression tree models with shrinkage priors, supporting:
It includes four core functions:
HorseTrees(): fits a single regression tree with a
standard Horseshoe prior.ShrinkageTrees(): fits a single tree with customizable
shrinkage priors.CausalHorseForest(): fits a causal forest using the
standard Horseshoe prior.CausalShrinkageForest(): fits a flexible causal forest
with user-defined shrinkage priors and tuning options.The ...Trees functions use a single learner to estimate
the outcome model directly. In contrast, the
Causal...Forest variants fit separate models for the
treated and control regression function. This enables estimation of
conditional average treatment effects (CATEs).