--- title: "dbarts as a Component" author: "Vincent Dorie" date: "2026-08-15" output: pdf_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) library(dbarts, quietly = TRUE) ``` A `dbartsSampler` is a mutable object that runs one sweep at a time, so BART can be a conditional model inside a larger Gibbs or Metropolis-Hastings scheme rather than the whole of it. Seven recipes follow. The first six are complete, runnable outer samplers; each answers the same three questions, which are the ones a composition gets wrong: * **which channel** the outer block writes -- response, offset, weights, sigma, or a predictor column; * **which fit** the outer block reads back, since the six reporting methods do not report the same quantity (`?dbartsSampler`, "Reading the fit"); * **whose prior** is in force, since a swapped response does not re-anchor the leaf prior the sampler was built with. `?dbarts-embedding` is the reference page for the channels and the reads; this vignette is the worked set. Recipes 1 through 6 are each also a cell in the package's test suite, because a composition that runs to completion can still be sampling the wrong posterior. The seventh recipe is different in kind: it drives `bart()` itself through a compiled per-draw callback rather than a `dbartsSampler` loop, and is not evaluated here -- it needs a C++ compiler `dbarts` does not depend on -- so its own check is a plain-C copy of the same callback under `inst/tinytest`, not a chunk in this file. Throughout, `run(0L, 1L)` advances the sampler exactly one sweep, and every recipe builds its sampler the same way: ```{r common} n <- 200L nBurn <- 100L nDraws <- 100L nTrees <- 25L # a single seeded chain that serializes no state on a per-sweep mutation recipeControl <- function(seed) { dbartsControl(n.chains = 1L, n.threads = 1L, updateState = FALSE, seed = seed) } ``` ## 1. Data augmentation written in R A binary response is fit by drawing a per-observation latent and regressing on it. `dbarts` will do this internally with `family = "probit"`, but a host that needs the latent for its own blocks -- a selection equation, a shared factor -- can run the Albert-Chib step itself with `dbartsDrawLatents` and `dbartsWorkingResponse` over an ordinary gaussian sampler. Two details make it correct rather than merely runnable. The sampler's residual standard deviation is pinned at 1 with `family = gaussian(sigma = fixed(1))`, because the augmented model has no free residual scale. And the leaf prior is stated with `setLeafPrior`, because a gaussian sampler built on a cold-start vector would otherwise inherit its prior scale from that vector's range -- an accident of the starting values, not a modeling choice. ```{r recipe1} set.seed(101L) x <- matrix(runif(n * 3L), n, 3L) fTrue <- 2 * (x[, 1L] - 0.5) + sin(2 * pi * x[, 2L]) o <- rep(0.25, n) y <- rbinom(n, 1L, pnorm(fTrue + o)) control <- recipeControl(31L) # the host: a gaussian sampler on the working response, sigma pinned at 1 host <- dbarts(x, 2 * y - 1 - o, control = control, n.trees = nTrees, family = gaussian(sigma = fixed(1))) host$setLeafPrior(normal(sd = 1)) composed <- matrix(0, n, nDraws) for (i in seq_len(nBurn + nDraws)) { # the fit WITHOUT the offset is what the augmentation conditions on; both # helpers form the linear predictor themselves as fit + offset z <- dbartsDrawLatents("probit", host$getFitsWithoutOffset(), y, offset = o) host$setResponse(dbartsWorkingResponse("probit", z, y, offset = o)) host$run(0L, 1L) if (i > nBurn) composed[, i - nBurn] <- host$getFitsWithoutOffset() } # a native probit sampler, given the SAME leaf prior, targets the same posterior native <- dbarts(x, y, offset = o, family = "probit", control = control, n.trees = nTrees, leaf.prior = normal(sd = 1)) nativeFit <- rowMeans(native$run(nBurn, nDraws)$train) - o c(agreement = cor(rowMeans(composed), nativeFit), signal = cor(rowMeans(composed), fTrue)) ``` Agreement, not identity: these are two chains of one posterior, so they differ by Monte Carlo error. The package's tests state the tolerance. ## 2. An outer block on the offset channel The offset is the channel for anything additive that another block owns. Here the model is partially linear, $y_i = f(x_i) + \alpha w_i + \epsilon_i$: BART holds $f$, a conjugate normal block holds $\alpha$, and each writes its own contribution into the other's conditioning through `setOffset`. The **seeding contract** is the whole of the correctness argument. The offset in force during a sweep must be the value the outer block currently holds, so the sampler is created with the initial offset installed and `setOffset` is called before `run`, never after. A sampler run once against a stale offset has drawn $f$ from the wrong conditional, and nothing later reports it. ```{r recipe2} alphaTrue <- 1.25 tau <- 2 set.seed(11L) x <- matrix(runif(n * 3L), n, 3L) w <- rnorm(n) fTrue <- 3 * (x[, 1L] - 0.5)^2 + sin(2 * pi * x[, 2L]) y <- rnorm(n, fTrue + alphaTrue * w, 0.5) alpha <- 0 control <- recipeControl(8L) # created WITH the offset the loop's first sweep will run against sampler <- dbarts(x, y, offset = alpha * w, control = control, n.trees = nTrees) keepAlpha <- rep(NA_real_, nDraws) for (i in seq_len(nBurn + nDraws)) { sampler$setOffset(alpha * w) # install, then sweep sampler$run(0L, 1L) f <- as.numeric(sampler$getFitsWithoutOffset()) sigma <- as.numeric(sampler$getSigmas()) v <- 1 / (sum(w * w) / sigma^2 + 1 / tau^2) alpha <- rnorm(1L, v * sum(w * (y - f)) / sigma^2, sqrt(v)) if (i > nBurn) keepAlpha[i - nBurn] <- alpha } c(mean = mean(keepAlpha), sd = sd(keepAlpha), truth = alphaTrue) ``` ### Checking the composition Recovering the truth on one data set is not evidence that a sampler targets its posterior. `dbartsValidateComposition` runs simulation-based calibration over the same one-sweep step: it draws $\theta_0$ from the prior, simulates data at it, starts the chain there, and ranks each functional's value at $\theta_0$ among the retained draws. Uniform ranks are what a correct kernel produces. Stating the prior draw here means drawing $f$ from BART's own prior, which `sampleTreesFromPrior` and `sampleLeafParametersFromPrior` do on the sampler itself. The problem is deliberately small, since every replication runs a full chain. ```{r recipe2_sbc} nSbc <- 60L sigma0 <- 1 tauSbc <- 1 set.seed(4L) xSbc <- matrix(runif(nSbc * 2L), nSbc, 2L) wSbc <- rnorm(nSbc) sbcControl <- dbartsControl(n.chains = 1L, n.threads = 1L, updateState = FALSE) # the response the sampler is BUILT on fixes the response transform, and the # leaf prior is then stated outright, so every replication shares one prior newSampler <- function() { s <- dbarts(xSbc, seq(-3, 3, length.out = nSbc), control = sbcControl, n.trees = 10L, family = gaussian(sigma = fixed(sigma0^2))) s$setLeafPrior(normal(sd = 1)) s } drawPrior <- function() { s <- newSampler() s$sampleTreesFromPrior() s$sampleLeafParametersFromPrior() list(sampler = s, f = as.numeric(s$getFitsWithoutOffset()), alpha = rnorm(1L, 0, tauSbc)) } simulateData <- function(theta0) rnorm(nSbc, theta0$f + theta0$alpha * wSbc, sigma0) initState <- function(theta0, data) { theta0$sampler$setResponse(data) list(sampler = theta0$sampler, y = data, f = theta0$f, alpha = theta0$alpha) } step <- function(state) { state$sampler$setOffset(state$alpha * wSbc) state$sampler$run(0L, 1L) state$f <- as.numeric(state$sampler$getFitsWithoutOffset()) v <- 1 / (sum(wSbc * wSbc) / sigma0^2 + 1 / tauSbc^2) state$alpha <- rnorm(1L, v * sum(wSbc * (state$y - state$f)) / sigma0^2, sqrt(v)) state } functionals <- function(state) c(alpha = state$alpha, f.mean = mean(state$f)) dbartsValidateComposition(drawPrior, simulateData, initState, step, functionals, n.replications = 40L, n.draws = 50L, n.thin = 5L, n.burn = 20L, seed = 12L) ``` Reporting `v * sum(wSbc * (state$y - state$f)) / sigma0^2` in place of the draw -- the defect this diagnostic exists to catch -- flags `alpha` here, and is a test cell of its own. ## 3. Several forests, one prior budget An additive decomposition $f(x) = \sum_{k=1}^K g_k(x)$ over $K$ samplers is ordinary backfitting: each sampler regresses on the partial residual through `setResponse`. What does not come for free is the prior. Each sampler was built on the full response, so each carries a leaf prior sized for the *whole* of $f$; their sum then has $\sqrt{K}$ times the prior standard deviation a single sampler would have, and nothing warns about it. `setLeafPrior(normal(sd = base / sqrt(K)))` restates each forest's share of one budget. The check below draws from the composed prior directly -- `sampleTreesFromPrior` acts on the engine, so it sees what `setLeafPrior` wrote -- and compares it with a single sampler's. ```{r recipe3} K <- 2L set.seed(21L) x <- matrix(runif(n * 4L), n, 4L) colnames(x) <- paste0("x", 1:4) fTrue <- 2 * sin(pi * x[, 1L]) + x[, 2L] + 3 * (x[, 3L] - 0.5)^2 - 2 * x[, 4L] y <- rnorm(n, fTrue, 0.5) control <- recipeControl(5L) blocks <- list(x[, 1:2, drop = FALSE], x[, 3:4, drop = FALSE]) samplers <- lapply(blocks, function(xb) dbarts(xb, y, control = control, n.trees = nTrees)) base <- samplers[[1L]]$getLeafPrior()$k.scale / samplers[[1L]]$getK()[1L] for (s in samplers) s$setLeafPrior(normal(sd = base / sqrt(K))) # the composed prior on the TOTAL, against a single sampler's priorSd <- function(samplers) { draws <- replicate(200L, { total <- 0 for (s in samplers) { s$sampleTreesFromPrior() s$sampleLeafParametersFromPrior() total <- total + as.numeric(s$getFitsWithoutOffset()) } total }) mean(apply(draws, 1L, sd)) } single <- dbarts(x, y, control = control, n.trees = nTrees) c(composed = priorSd(samplers), single = priorSd(list(single))) ``` ```{r recipe3_fit} fits <- rep(list(rep(0, n)), K) total <- matrix(0, n, nDraws) for (i in seq_len(nBurn + nDraws)) { for (k in seq_len(K)) { samplers[[k]]$setResponse(y - Reduce(`+`, fits[-k])) samplers[[k]]$run(0L, 1L) fits[[k]] <- as.numeric(samplers[[k]]$getFitsWithoutOffset()) } if (i > nBurn) total[, i - nBurn] <- Reduce(`+`, fits) } c(composed = cor(rowMeans(total), fTrue), single = cor(rowMeans(single$run(nBurn, nDraws)$train), fTrue)) ``` `setResponse` leaves the response transform where creation put it (`updateScale = FALSE`, the default), which is what makes a per-sweep swap coherent: every sweep's partial residual is read on one scale. This is the case where every forest enters the response with a multiplier of one. `vignette("gibbs_sampler_mixture_model", package = "dbarts")` carries the general form, where forest $f$ enters through a known per-observation $m_{f,i}$, along with its measured price and what it costs to reach the same posterior a `forests = ` sampler targets. ## 4. A latent predictor shared by two samplers When the quantity being sampled is a *predictor* rather than a response -- a latent ability, a propensity, an unobserved dose -- it moves through `updatePredictorPerObservationJointly`, which installs each observation in every sampler or in none. An observation is declined when its new value would empty a leaf in any tree of any sampler. Score each proposal *before* installing it, with `predict` while the trees are held: row $i$'s fit depends on row $i$'s predictors alone, so the per-row ratio factorizes. Decide accept or reject per row, then install only the accepted rows. The returned mask is **part of the Metropolis accept decision, not a diagnostic**: a row is declined when it is the only row in a leaf of some tree, so moving it would empty that leaf. An empty leaf has zero prior, so the decline is the correct outcome and counts as a rejection. The random-walk proposal is symmetric, so the ratio needs no proposal correction. Keeping the host's copy equal to what the samplers hold means taking the mask, not the proposal. Nothing is reverted. ```{r recipe4} set.seed(404L) nJ <- 80L thetaTrue <- rnorm(nJ) d1 <- data.frame(theta = thetaTrue, u = runif(nJ)) d1$y <- 1.5 * thetaTrue + rnorm(nJ, 0, 0.3) d2 <- data.frame(v = runif(nJ), theta = thetaTrue) d2$y <- rbinom(nJ, 1L, pnorm(thetaTrue)) control <- recipeControl(44L) sA <- dbarts(y ~ theta + u, d1, control = control, n.trees = nTrees, power = 0.5) sB <- dbarts(y ~ v + theta, d2, family = "probit", control = control, n.trees = nTrees, power = 0.5) invisible(sA$run(0L, 20L)) invisible(sB$run(0L, 20L)) # per-observation log likelihood of a latent theta, scored with the trees held # fixed; nothing is installed logLik <- function(theta) { fA <- as.numeric(sA$predict(transform(d1[c("theta", "u")], theta = theta))) fB <- as.numeric(sB$predict(transform(d2[c("v", "theta")], theta = theta))) dnorm(d1$y, fA, as.numeric(sA$getSigmas()), log = TRUE) + ifelse(d2$y == 1, pnorm(fB, log.p = TRUE), pnorm(-fB, log.p = TRUE)) + dnorm(theta, 0, 1, log = TRUE) } current <- sA$data@x[, "theta"] proposal <- current + rnorm(nJ, 0, 0.5) wanted <- log(runif(nJ)) < logLik(proposal) - logLik(current) # install only the wanted rows; a declined install is a rejection installed <- wanted installed[wanted] <- updatePredictorPerObservationJointly( list(sA, sB), ifelse(wanted, proposal, current), column = "theta" )[wanted] theta <- ifelse(installed, proposal, current) c(wanted = sum(wanted), accepted = sum(installed), n = nJ) ``` ## 5. An outer sampler that owns sigma `family = gaussian(sigma = fixed(value))` holds the residual variance at `value` and suppresses the sampler's own draw, which is exactly what a host wants when sigma is drawn elsewhere -- from a hierarchical prior, jointly with another outcome, or from a scale shared across several samplers. `setSigma` then writes it, and it stays written. ```{r recipe5} set.seed(505L) x <- matrix(runif(n * 3L), n, 3L) fTrue <- 2 * (x[, 1L] - 0.5) + sin(2 * pi * x[, 2L]) y <- rnorm(n, fTrue, 0.75) sigma <- 1 control <- recipeControl(55L) sampler <- dbarts(x, y, control = control, n.trees = nTrees, family = gaussian(sigma = fixed(sigma^2))) keepSigma <- rep(NA_real_, nDraws) for (i in seq_len(nBurn + nDraws)) { sampler$setSigma(sigma) # the sampler will not redraw it sampler$run(0L, 1L) f <- as.numeric(sampler$getFitsWithoutOffset()) sigma <- sqrt(sum((y - f)^2) / rchisq(1L, n)) if (i > nBurn) keepSigma[i - nBurn] <- sigma } c(mean = mean(keepSigma), truth = 0.75) ``` The channel exists only where the model has a residual scale to own. A family that fixes it by definition refuses the write rather than accepting it and quietly rescaling every leaf posterior: ```{r recipe5_guard, error=TRUE} probit <- dbarts(x, rbinom(n, 1L, 0.5), family = "probit", control = control, n.trees = nTrees) probit$setSigma(1) ``` ## 6. Reading the fit an incremental scheme needs The last recipe is about a read, not a write. A host driving the sampler one sweep at a time conditions on $f(x_i)$, and `$getFitsWithoutOffset()` is the method that reports it. `run()$train` reports the same location with the installed offset folded in, and `$getLatents()` reports the augmentation variable, which carries no offset at all. Differencing the last two is the trap: `getLatents() - run()$train` is short by exactly the offset, so an outer block that regresses it on the offset's own covariate recovers a coefficient shrunk toward zero. ```{r recipe6} gammaTrue <- 1.5 set.seed(314L) x <- matrix(runif(n * 3L), n, 3L) v <- rnorm(n) fTrue <- 2 * (x[, 1L] - 0.5) + sin(2 * pi * x[, 2L]) y <- rbinom(n, 1L, pnorm(fTrue + gammaTrue * v)) gamma <- 0 control <- recipeControl(20L) sampler <- dbarts(x, y, offset = gamma * v, family = "probit", control = control, n.trees = nTrees) keepGamma <- rep(NA_real_, nDraws) for (i in seq_len(nBurn + nDraws)) { offset <- gamma * v sampler$setOffset(offset) train <- as.numeric(sampler$run(0L, 1L)$train) z <- as.numeric(sampler$getLatents()) f <- as.numeric(sampler$getFitsWithoutOffset()) # z - f, NOT z - train: the training channel carries the offset already vv <- 1 / (sum(v * v) + 1 / 4) gamma <- rnorm(1L, vv * sum(v * (z - f)), sqrt(vv)) if (i > nBurn) keepGamma[i - nBurn] <- gamma } c(mean = mean(keepGamma), truth = gammaTrue, shortfall = max(abs((z - train) - ((z - f) - offset)))) ``` The last column is zero by construction, which is the point: `z - train` and `z - f` differ by the installed offset exactly. Substituting `train` for `f` above pulls the posterior mean of $\gamma$ from `r round(mean(keepGamma), 2)` to about 0.29, roughly seventeen oracle standard errors below the truth, on a run that reports no error at all. ## 7. A per-draw callback that never returns the array Every recipe above drives a `dbartsSampler` from R, one sweep at a time. This one drives `bart()` itself and reads no returned fit back at all. `callback` fires a *compiled* function once per saved draw per chain, over `const` pointers into the engine's own buffers, on the thread that owns that chain -- so a reduction is taken as each draw is produced instead of over the `n` x `S` x `C` array `bart()` would otherwise build only to have the caller reduce it afterward. The shape is the one a compiled consumer wants: raw pointers and an opaque `void*` context, no `SEXP` crossing into the callback and nothing blocking on R. `Rcpp` is not a dependency of `dbarts`, so the code below does not run as part of building this vignette; `?bart`'s `callback` and `keepFits` entries document the contract in full, and `inst/tinytest/capi/consumer.c` carries a plain-C copy of the same callback, checked against a fit's own `yhat.train.mean`. The task: a running posterior mean of the training fit, over `x`/`y` as in the recipes above, without ever holding the training draws at all. ```{r recipe7_r, eval=FALSE} Rcpp::sourceCpp("drawMean.cpp") # makeMeanContext, meanCallbackPtr nChains <- 4L acc <- numeric(nrow(x) * nChains) # the ONLY per-observation allocation ctx <- makeMeanContext(acc, nrow(x), nChains) fit <- bart(x, y, n.chains = nChains, callback = list(fn = meanCallbackPtr(), context = ctx)) means <- rowMeans(matrix(acc, nrow(x), nChains)) stopifnot(contextStatus(ctx) == 0L) ``` ```cpp #include #include // dbarts_draw; PKG_CPPFLAGS -I the include dir struct MeanContext { double* out; // REAL() of the R vector, taken BEFORE the run std::size_t n, numChains; std::size_t* counts; // per chain: no chain reads another's counter int status; // set here, read from R after the run }; extern "C" int meanCallback(void* context, const dbarts_draw* draw) { MeanContext* ctx = static_cast(context); if (draw->train == nullptr || draw->numObservations != ctx->n || draw->numFittedValuesPerObservation != 1 || // train is n x F; F = 1 here draw->chainIndex >= ctx->numChains) { ctx->status = 1; // record and continue; nonzero would ABORT the run return 0; } double* out = ctx->out + draw->chainIndex * ctx->n; // disjoint slice const double m = static_cast(++ctx->counts[draw->chainIndex]); for (std::size_t i = 0; i < ctx->n; ++i) out[i] += (draw->train[i] - out[i]) / m; // running mean return 0; } // freeMeanContext (not shown) deletes counts and the context // [[Rcpp::export]] SEXP meanCallbackPtr() { // ...PtrFn: casting a function pointer to void* is // what -Wpedantic flags under R CMD check return R_MakeExternalPtrFn(reinterpret_cast(&meanCallback), R_NilValue, R_NilValue); } // [[Rcpp::export]] SEXP makeMeanContext(Rcpp::NumericVector acc, std::size_t n, std::size_t C) { MeanContext* ctx = new MeanContext{REAL(acc), n, C, new std::size_t[C](), 0}; // acc goes in the PROTECTED slot (third argument): that is what keeps the // array alive as long as the context is SEXP ptr = PROTECT(R_MakeExternalPtr(ctx, R_NilValue, acc)); R_RegisterCFinalizerEx(ptr, &freeMeanContext, TRUE); UNPROTECT(1); return ptr; } ``` Four rules make this recipe correct rather than merely running, and every one is the caller's to keep, not the engine's. `REAL(acc)` is read exactly ONCE, before the fit starts: a later reassignment of `acc` in R builds a DIFFERENT vector, and the callback goes on writing the old one, so nothing between `makeMeanContext` and `bart()` may touch `acc`. Chains write DISJOINT slices -- `draw->chainIndex * ctx->n` onward -- which is what lets every chain's calls run concurrently with no lock at all; two chains sharing a slice would need an atomic or a lock, giving up the bitwise reproducibility `dbarts` otherwise guarantees, for nothing this recipe needs. A context is built PER RUN, because `drawIndex` restarts at 0 on every `bart()` call (and every `$run`, for a `dbartsSampler` driven directly) -- reusing one context across two fits overwrites the first fit's slots rather than extending them. And the return value is an ABORT switch, not an error channel: nonzero stops the whole run at that draw and discards every chain's results immediately, so a callback that merely disagrees with a draw's shape -- the `if` above -- records a `status` and returns 0, leaving the caller to check `contextStatus(ctx)` (not shown; it reads `ctx->status` back) once the fit returns. Getting either of the last two wrong does not raise an R condition: it crashes the session outright, with nothing to `tryCatch`. Under all four sits the rule the header states in these terms: NO R API inside the callback, ever -- not `Rf_allocVector`, not `PROTECT`, not `Rf_error`, and in C++ not the construction or destruction of an `Rcpp` proxy type, since `Rcpp::NumericVector` and its siblings touch R's protection stack on both, allocation visible or not. R's evaluator, allocator and protection stack belong to the main thread, and `Rf_error` would longjmp out of a context that thread established. That is the whole reason `REAL(acc)` is taken in `makeMeanContext` and only a raw `double*` crosses into `meanCallback`. Supplying `callback` defaults `keepFits` to `FALSE` automatically, so this fit returns no `yhat.train` at all -- `means` above, not a component of `fit`, is the only fitted value there is; passing `keepFits = TRUE` alongside `callback` keeps the array too, at the memory cost the callback exists to avoid. A strided variant that keeps every draw rather than reducing it on the fly is the same shape with a `numDraws` field added to `MeanContext` and the loop body replaced by one `memcpy` into `out + (draw->chainIndex * numDraws + draw->drawIndex) * ctx->n`. ## Where to look next `?dbarts-embedding` collects the channels, the reads and the rules the recipes above apply. `?dbartsSampler` documents each method, including which quantity each reporting method returns; `?dbartsAugmentation` documents the two augmentation helpers; `?dbartsValidateComposition` documents the calibration check; and `?bart`'s `callback` and `keepFits` entries document recipe 7's contract, including the crash and lost-fitted-value warnings the recipe states in prose.