---
title: "Checking Convergence and Fit"
author: "Noah Greifer"
date: "`r Sys.Date()`"
output:
  html_vignette:
    toc: yes
vignette: >
  %\VignetteIndexEntry{Checking Convergence and Fit}
  %\VignetteEngine{knitr::rmarkdown_notangle}
  %\VignetteEncoding{UTF-8}
bibliography: references.bib
link-citations: true
---

```{=html}
<style>
pre {overflow-x: auto;}
pre code {word-wrap: normal; white-space: pre;}
</style>
```

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6.5,
                      fig.height = 3.0, dpi = 110, fig.align = "center",
                      message = FALSE, warning = FALSE)
options(digits = 3)

has_post <- rlang::is_installed("posterior")
# marginaleffects summarizes posterior draws with collapse, so both are needed.
has_me <- rlang::is_installed(c("marginaleffects", "collapse"))
has_bp <- rlang::is_installed("bayesplot")
has_loo <- rlang::is_installed("loo") && rlang::is_installed("rstantools")
has_rst <- rlang::is_installed("rstantools")
run <- TRUE

# The fits below take minutes, so their output is saved ahead of time and
# replayed where the code is not run, as on CRAN; see saved-results.R.
source("saved-results.R", local = knitr::knit_global())
knitr::opts_chunk$set(saved = TRUE)
```

## Introduction

Checking a fit involves questions that often get asked together and are worth
separating. Has the sampler converged, meaning has it explored the posterior
properly? And does the model fit, meaning does it describe the data?

The first is about the algorithm and the second is about the model; a fit can
converge beautifully on a badly chosen family, and a well chosen family can be
fitted by a chain that has not run long enough.

In this guide, we will take those questions in that order. First we'll fit a
model with four chains, hand it to `diagnose()`, and then unpack what that
summary reports and why, including the one row (the additive predictor `eta`)
whose conventional threshold does not apply to a forest and the one quantity
that is deliberately left out of the table. Next we'll ask whether the model
describes the data, using posterior predictive checks, a calibration plot for
the binary case in which those checks are weak, residuals, and a Bayesian $R^2$.
A short checklist at the end collects what we recommend running before reporting
anything from a fit.

```{r libs, eval = run}
library(bartisan)

# For parallelization; optional
if (rlang::is_installed("future")) {
  future::plan(future::multisession)
}

data("rhc")

set.seed(2026)

fit <- bartisan(death ~ rhc + age + sex + race + edu + aps + meanbp + resp +
                  hema + pafi + paco2 + crea + surv2m + card,
                data = rhc, family = binomial(), chains = 4)
```

## Convergence

### Every Statistic in One Call (`diagnose()`)

`diagnose()` computes every convergence and mixing statistic in this section,
says which of them fall short, and says what to do about each. Everything it
reports comes from the stored draws, so it needs no other package.

```{r diagnose, eval = run}
diagnose(fit)
```

The rest of this section explains what it reports and why, which is worth
reading once; when the summary itself is enough, `vignette("bartisan")` is the shorter
tour of the whole workflow.

### Running More Than One Chain (`chains`)

Setting `chains = 4` above is doing the work. The default is one chain, which produces
estimates but no way to check most of what matters, because R-hat compares chains
to each other; with one chain it is computed by splitting that chain, which
catches drift but cannot catch two chains settling in different places. Running
four costs four times as much sampling, which for most fits is a few seconds, and
less than that under a *future* plan.

More chains is not, however, a way to improve `rhat`, and it is worth knowing
which way it cuts before reaching for it. Effective sample size depends on the
total number of draws and not on how they are divided, so twice as many chains
and twice as long a chain buy the same amount of it. `rhat` is not symmetric that
way: because it compares chains against each other, its value under a sampler
with nothing wrong is about `1 + chains / ess`, so at 400 effective draws four
chains sit at 1.010 and sixteen chains at 1.040. Adding chains raises the bar it
has to clear. When a fit fails on `rhat`, longer chains are the fix; when it
fails only on effective sample size, either will do.

The table is in `diagnose(fit)$table` when it is wanted as a data frame, and the
checks and the advice are in `$checks` and `$advice`.

`rhat` compares the variance between chains to the variance within them
[@vehtari2021]; if the chains have found the same posterior the two agree and
the ratio is near 1. `ess_bulk` and `ess_tail` are effective sample sizes,
meaning how many independent draws the correlated ones are worth, in the middle
of the distribution (i.e., the bulk) and in the tails respectively.

The conventional thresholds are `rhat` below 1.01 and effective sample sizes
above about 400 for a quantity we intend to report precisely. Those come from
the general MCMC literature and are a reasonable default here, with one
exception described below. `diagnose()` applies them, and takes both as
arguments so they can be moved.

Some of the columns are worth naming separately. `rhat_late` is the same
statistic computed on the second half of the retained draws alone, which tells a
warmup that ended too early (i.e., too small a `num_burn`) from chains that have
each settled somewhere different: throwing away the early draws has the same
effect as a larger `num_burn`, so if that fixes `rhat`, warmup was the problem.
The `splits.*` row is the total number of splitting rules in the forest at
each draw, the one quantity here that is about the trees rather than about the
fitted values, and which no general-purpose MCMC diagnostic would think to look
at.

### The Rows of the Table

`loglik` is the log likelihood of the whole dataset at each draw. It is a useful
scalar summary of the fit, and usually the slowest-mixing row in the table,
since it moves with every observation's fitted value at once.

`aux.*` are the nuisance parameters of the family (i.e., the parameters that are
not part of the additive predictor), and are usually the best behaved rows in
the table. The binomial likelihood has none of its own, so none appear here,
though a `vc()` coding contributes its coefficients to `aux`, which is where the
`aux.b.*` rows in the `bcf()` table below come from; a
Gaussian fit would show `aux.sigma`, and a survival fit its scale or baseline
hazard.

`eta.*` is the additive predictor, and it gets two rows. One summarizes the
worst 5% of observations rather than the single worst, because the worst of a
thousand values is extreme even when every chain has converged; the other is the
average over observations. Reading them against each other is the point of having
both, because they routinely differ by orders of magnitude: a forest settles the
level of the fitted function within a sweep or two and takes far longer to settle
which observation gets which share of it. A prediction for one observation is
governed by the worst-5% row.

Neither row settles an effect, and the average row in particular should not be
read as though it did. An effect is a *contrast*, and a contrast can mix badly
where the function it is a contrast of mixes well, so what
`estimate_effect()` or `avg_comparisons()` reports needs diagnosing on its own
draws rather than inferring from this table. The section below does that.

One quantity is deliberately not in the table. The leaf prior scale, in
`fit$sigma_mu`, mixes badly in every BART implementation: on one dataset its
counterparts in *dbarts* and in *stochtree* both come out above `rhat` 1.1 with
effective sample sizes in the tens out of 3000 draws, and *stochtree*
draws it exactly from its full conditional, so no sampler can do better. It is a hyperparameter
whose disagreement between chains does not reach the fitted function, which on
those same fits has `rhat` 1.00 and thousands of effective draws, so reporting it
beside `eta` only produces an alarming row with nothing to act on. It is out of
this table and only out of this table: the draws are in `fit$sigma_mu` and reach
`as_draws()`, so it can still be diagnosed by anyone who wants to.

### The Exception for `eta`

The `eta` row is a high percentile over every observation, so it is conservative
by construction, and forests mix slowly on their fitted values. Two chains can visit
quite different collections of trees, and a forest settles the level of the
fitted function in a sweep or two while taking far longer to settle which
observation gets which share of it, which is why the per-observation R-hats run
high without meaning the two chains disagree
about anything we would report, so whether two chains found the same trees is
not a question the `eta` rows can answer; `splits.*` is the row that asks it.

On harder problems this shows up clearly. Below we fit the Friedman function,
first with a chain far too short and then at the defaults.

```{r friedmanbad}
set.seed(3)
fr <- as.data.frame(matrix(runif(400 * 10), 400, 10))
names(fr) <- paste0("x", 1:10)
fr$y <- 10 * sin(pi * fr$x1 * fr$x2) + 20 * (fr$x3 - 0.5)^2 +
  10 * fr$x4 + 5 * fr$x5 + rnorm(400)

too_short <- bartisan(y ~ ., fr, family = gaussian(), chains = 4,
                      control = bartisan_control(num_trees = 20, num_burn = 50,
                                                 num_draws = 50))

diagnose(too_short)$table
```

Almost everything is bad. Every `rhat` but the average over observations is
far above the 1.01 threshold, and the
effective sample sizes are a fraction of the 200 draws kept (4 chains x 50 kept
draws per chain) on every row but that one. `rhat_late` is no
better than `rhat` on any row but the residual standard deviation, which is the
informative part: discarding the
early draws does not fix it, so this is not a warmup that ended too soon but
chains that have each settled somewhere different. Nothing from this fit should
be used.

```{r friedmangood}
long_enough <- bartisan(y ~ ., fr, family = gaussian(), chains = 4)

diagnose(long_enough)$table
```

Everything has improved by a large factor, and the residual standard deviation
is nearly at the threshold. But the log likelihood and `eta` are both still above
1.01, and `eta`'s effective sample size is still in the tens. That residual is
the phenomenon described above rather than a chain that merely needs to be
longer, and adding draws does not reliably remove it.

When the `eta` row is elevated, the distribution behind it is more informative
than the maximum:

```{r perobs, eval = run && has_post}
per_observation_rhat <- function(fit) {
  draws <- fit$eta$eta
  per <- nrow(draws) / fit$chains
  index <- matrix(seq_len(per * fit$chains), per, fit$chains)
  apply(draws, 2, function(column) {
    posterior::rhat(matrix(column[index], per, fit$chains))
  })
}

round(quantile(per_observation_rhat(fit), c(0.5, 0.9, 0.99, 1)), 3)
```

The median is a little above 1 and so is most of the distribution, so the
worst-5% row in the table is not one badly behaved observation: the fitted
function as
a whole mixes slowly here. That is the phenomenon described above rather than a
broken chain, and the two are told apart by what happens to the quantities we
report. If the log likelihood has converged, along with any nuisance parameters
the family has, and an `avg_comparisons()` estimate is stable across chains and
across a longer run, the fit is usable.

If instead the log likelihood were badly elevated too, or the estimate moved when
the chain was lengthened, that would be a chain that has not converged, and the
answer is more draws or a simpler model.

### Looking at the Chains (`as_draws()`)

`as_draws()` hands the fit to *posterior* and *bayesplot*. It carries the
scalar parameters and a representative set of `eta` columns.

```{r summarise, eval = run && has_post}
library(posterior)

as_draws(fit) |>
  summarise_draws()
```

```{r trace, eval = run && has_bp, fig.height = 2.6}
library(bayesplot)

as_draws(fit, eta = 1) |>
  mcmc_trace(pars = c("loglik", "eta[1]"))
```

What we want to see is four chains overlapping, wandering around the same level,
with no drift and no long excursions. Setting `eta = 1` selects the first
observation, `eta = TRUE` a spread of ten, and `eta = FALSE` none.

### Diagnosing the Quantity Being Reported

Everything above is about the parameters the sampler draws. What gets reported is
usually not one of them: it is an effect, a contrast, a predicted probability at
some covariate value, a difference between two subgroups. Those need their own
diagnostics, because a contrast can mix badly where the function it is a contrast
of mixes well, and nothing in the table above will say so.

When estimating a treatment effect, the quantity of interest is produced by `estimate_effect()`, the output of which is accepted by `diagnose()` directly. Below we fit a Bayesian causal forest model to estimate the ATE:

```{r diagnose_effect, eval = run}
bcf_fit <- bcf(death ~ age + sex + race + edu + aps + meanbp + resp +
                 hema + pafi + paco2 + crea + surv2m + card,
               treat = ~rhc, data = rhc,
               family = binomial("probit"), chains = 4)

estimate_effect(bcf_fit, estimand = "ATE") |>
  diagnose()
```

The output has a row for the contrast and one for each of the two average
potential outcomes it is a contrast of. Here the contrast clears the thresholds
and one of the averages, `Y[1]`, falls just short on bulk effective sample
size. It does not always come out that way round: a difference can mix worse
than either of its parts, so all three rows are worth reading rather than the
contrast alone.

Now compare what the fit's own table says about the same sampler:

```{r diagnose_bcf_fit, eval = run}
diagnose(bcf_fit)$table
```

The two forests are far worse than the estimand. Their R-hats run from 1.2 to
1.6 with effective sample sizes in the tens, where the effect clears 400. That
is not a contradiction: a varying
coefficient model is $f_0(x) + z\,f_1(x)$, and on the treated units the two
forests are only jointly identified, so the sampler can move level between them
while their sum, and much of what is built from it, stays where it was. The
effect inherits part of that wandering and not all of it.

So the estimand's diagnosis cannot be deduced from the fit's in either
direction. Read only the `eta` rows here and the fit looks unusable; read only
the effect's own table and the sampler looks close to adequate; neither reading
gives the other. Both are worth
running, and it is the second that decides whether the number being reported can
be reported.

### Diagnosing Anything Else (`summarise_draws()`)

For a quantity this package has no method for, *posterior* is the general tool:
anything that can be arranged as draws by chains can be handed to
`summarise_draws()`, which computes the same statistics `diagnose()` reports.

The arranging is the only step with a decision in it. The sampler stacks its
chains, so the draws of any quantity are one long vector in chain order, and
folding it into an iterations-by-chains matrix recovers the structure R-hat
needs. Below, the quantity is an average comparison computed by
*marginaleffects*, whose draws come out of `posterior_draws()`:

```{r me_draws, eval = run && has_me && has_post}
library(marginaleffects)

drawn <- avg_comparisons(fit, variables = "rhc") |>
  posterior_draws()

chains <- fit$chains
per_chain <- nrow(drawn) / chains

folded <- array(drawn$draw, dim = c(per_chain, chains, 1L),
                dimnames = list(NULL, NULL, "ATE"))

folded |>
  as_draws_array() |>
  summarise_draws("mean", "rhat", "ess_bulk", "ess_tail")
```

The same statistics `diagnose()` reports, computed by a different route on a
quantity it has no method for,
which is the point: any posterior quantity can be diagnosed this way, including
ones neither package knows about. Draw it, fold it, summarize it.

The fold has to match the order the draws are stored in, which here is all of
chain one, then all of chain two, and so on; specifying
`dim = c(per_chain, chains, 1)` fills the array in that order. And the
quantity has to be the one being reported: diagnosing a prediction is not
diagnosing the contrast of two predictions, for the reason the section above
gives.

### Remedies for Poor Mixing

There are a few remedies worth trying, given here in the order in which they
usually help. **More draws**, from increasing `num_burn` and `num_draws`, fix
most cases. **A smaller forest** helps next, since reducing `num_trees` leaves
fewer ways to represent the same function, so the sampler mixes faster. **A
different family** can be needed last, because a likelihood that fits the data
badly can produce a posterior that is hard to explore, so the family is worth
checking when more draws and a smaller forest have not helped;
`vignette("families")` covers the alternatives.

Note these latter two options change the model itself, so only use them when required, not as a routine fix.

The first step is to rule out a false alarm. `rhat` above 1.01 on a quantity
carrying only a handful of effective draws is not yet evidence that the chains disagree, since
that is roughly where `rhat` sits when nothing is wrong, and `diagnose()` says so
rather than reporting a disagreement it cannot support. In that case the
effective sample size needs fixing, and more draws is the whole of the
remedy.

### A Worked Example: More Draws, Same Model

The first remedy is the only one that leaves the model alone, so it is worth
seeing how far it goes on its own. Below, a fit that fails and then the same
model run longer; the formula, the family, the tree count and every prior setting
are identical between the two, and only `num_burn` and `num_draws` move.

```{r mixmodel, eval = run}
model <- death ~ rhc + age + aps + meanbp + surv2m

set.seed(2026)
short <- bartisan(model, data = rhc, family = binomial(), chains = 4)

diagnose(short)
```

Several things are flagged. `rhat` is above the threshold, and the check under it
says that the R-hat in question rests on too few effective draws for the
threshold to mean anything, which is the case the previous section says to rule
out first: the chains are not demonstrably disagreeing, they have not run long
enough for the statistic to be readable. Both effective sample sizes are short of
400 as well. Everything here points the same way, and none of it points at the
model.

So we run it longer, increasing the burn-in to 2000 and the number of draws per chain to 8000. Nothing else changes, so we can call `update()` on the previous model with our new arguments supplied:

```{r mixfixed, eval = run}
set.seed(2026)
longer <- update(short, num_burn = 2000, num_draws = 8000)

diagnose(longer)
```

Most of it is cleared. The split-count row and both `eta` rows are well past
the thresholds, and every tail effective sample size is above 900. What is left
sits on the log likelihood alone: its bulk effective sample size rose from 85 to
208 rather than tenfold, and its R-hat of 1.02 is again the kind the check calls
unreadable, since four chains average 1.019 at that count even when they agree.
The log likelihood is usually the slowest-mixing row in a forest, because it
moves with every fitted value at once, and it is not a quantity we report; the
fitted function and its contrasts are, and those are fine here. This fit can be
reported from, with the log-likelihood row noted rather than chased.

The arithmetic of that fix carries over to other fits, with one qualification.
Effective sample size grows roughly in proportion to the draws, so the shortfall
tells you the factor you need: a fit an order of magnitude short of 400 needs
about an order of magnitude more draws, not a little more. The `eta` rows and
the split count grew five- to ninefold here on ten times the draws.

The qualification is that the estimate of effective sample size is itself
noisy, and can fall as the chain lengthens: a short chain cannot see
autocorrelation at long lags, so it reports an efficiency the chain does not
have. The log likelihood above did exactly that: its efficiency fell from 2.6
to 0.6 effective draws per hundred kept once the longer chain could see its
slow component. Reading a single ESS figure as exact invites chasing
it; the factor should be read instead.

Increasing draws does not always suffice. It did most of the work here because
the flagged R-hat was the unreadable kind, resting on too few effective draws.
Where R-hat instead stays above the threshold however many
draws are added, as it does for
the Friedman fit above, the chains genuinely disagree and more of them will not
help. That is when the remedies that change the model come up, and the order in
the list above is the order to try them in.

## Fit

Convergence says the sampler did its job; it says nothing at all about whether
the model is right. The checks below assess the fit of the model to the data.

### Reading the Prior Back (`prior_summary()`)

Everything else in this section compares the model to the data it was fitted to.
Some checks are worth running before any of that information has been spent, and
the first is simply to read the prior back.

```{r priorsummary, eval = run && has_rst}
rstantools::prior_summary(fit)
```

This is the prior the engine was given rather than the one we wrote, and the two
are not the same object: we named no priors at all, and the defaults are a
function of the response and of the number of trees. Note in particular that the
leaf scale is not a number anyone typed. It is `3 * s / (k * sqrt(num_trees))`,
so halving the number of trees widens each leaf's prior by a factor of
$\sqrt 2$ and leaves the prior on their sum alone, which is why `num_trees` is
not the tuning knob it looks like.

### Before the Data: The Prior Predictive (`prior_only`)

Reading the prior back is not the same as knowing what it implies. The settings
in `bartisan_control()` are statements about trees and leaves, and nobody has
intuition for what `k = 2` implies about a patient's chance of dying. A prior
predictive check puts the prior on the outcome's own scale, where it can be
judged.

Setting `prior_only = TRUE` fits the same model to no data. Every observation is given a
weight of zero, and since the weight multiplies that observation's contribution
to the likelihood, the likelihood goes flat and the sampler draws from the prior.

```{r prioronly, eval = run}
set.seed(2026)

prior <- bartisan(death ~ rhc + age + sex + race + edu + aps + meanbp + resp +
                    hema + pafi + paco2 + crea + surv2m + card,
                  data = rhc, family = binomial(), prior_only = TRUE)
```

What comes back is an ordinary fit whose draws are prior draws, so everything
that reads a fit reads this one. The question to put to it is whether the risks
the prior considers plausible are ones a critically ill patient could actually
have.

```{r priorquant, eval = run}
p <- c(.01, .1, .5, .9, .99)

quantile(predict(prior, type = "response", draws = TRUE), p)

quantile(predict(fit, type = "response", draws = TRUE), p)
```

The prior runs from about .01 to about .99, as we would want of a prior on a
probability: it rules almost nothing out and asserts almost nothing. The fitted
model runs from about .18 to about .95 over the same quantiles, so the data have
brought it a long way in from where it started. A prior that had come back
concentrated near zero and one, or confined to a narrow band in the middle, would
be worth changing before going any further.

There is one thing not to read into this. The median sits close to the observed
death rate because the additive predictor is anchored at an intercept-only fit,
which is how a BART prior is specified and not a leak. The replicates take their
location and scale from the response and everything else from the prior, so read
them for shape and spread rather than for level.

Every family supports this except `ordinal()` and `ordbeta()`, whose cutpoints
are drawn from the likelihood alone and so have nothing to fall back on when the
likelihood goes flat; `?bartisan` has the details.

### Posterior Predictive Checks (`pp_check()`)

A posterior predictive check simulates outcomes from the fitted model and
compares their distribution to the observed one.

```{r ppc, eval = run && has_bp, fig.height = 3.0}
pp_check(fit)
```

For a continuous outcome this is the workhorse check, and we look for systematic
differences: replicates that are too narrow, that miss a second mode, or
that put mass where the outcome cannot go. Simulating negative values for an
outcome that cannot be negative says the family is wrong, and
`vignette("families")` covers the alternatives.

For a binary outcome it is a weak check, which is worth knowing before
reading too much into it. There are only two values the replicates can take, so
they will match the observed proportion unless the model has gone badly wrong;
passing this check tells us almost nothing. The calibration check below should
be read instead.

#### Checks Worth Reaching For

The `type` argument names any of *bayesplot*'s checks without its `ppc_` prefix, and four of
them earn their place for a forest.

The default compares whole distributions, which is a coarse question. A
test statistic is a sharper one: setting `type = "stat"` compares any summary of the
replicates with the same summary of the outcome, so the spread, a tail
quantile, or the share of zeros can each be asked about on its own. The share of
zeros finds an unmodeled point mass, the case for which
`vignette("families")` recommends `zi_poisson()` and `tweedie()`.

```{r ppcstat, eval = run && has_bp, fig.height = 2.6}
pp_check(fit, type = "stat", stat = "sd")
```

**Calibration of the whole predictive distribution** is checked by setting
`type = "loo_pit_ecdf"`. Each observation is transformed through its own
leave-one-out predictive distribution, which should leave a uniform if the model
is calibrated, and the plot compares the result with one. A curve that leaves
the band says the predictive spread is wrong even where the mean is right, which
the default check cannot see.

```{r ppcloo, eval = run && has_bp && has_loo, fig.height = 2.8}
pp_check(fit, type = "loo_pit_ecdf")
```

**The location of a residual**, rather than its size, is shown by setting
`type = "error_scatter_avg_vs_x"`: the average residual against a
predictor, which is where a missing interaction appears as structure. And
setting `type = "intervals"` draws a predictive interval per observation, so systematic
under-coverage is visible as a run of points outside their intervals.

The rest apply where the response allows, and say so when it does not: a
rootogram wants counts, a bar plot wants discrete values, a calibration plot
wants a binary outcome. The one written for censored data is
`type = "km_overlay"`, which overlays the replicate survival curves on the
observed Kaplan-Meier curve; it takes the censoring indicator from the fit's own response, and needs
the *ggfortify* package installed alongside *bayesplot*.

### Calibration

The useful question for a binary outcome is whether the predicted probabilities
mean what they say: among the patients the model gave a 30% chance of dying, did
about 30% die?

```{r calibration, eval = run && has_bp && has_loo, fig.height = 3.4}
pp_check(fit, type = "loo_calibration")
```

A line on the diagonal means the probabilities are calibrated, and this one sits
close to it across the whole range. The dots along the bottom show where the
predictions fall, so a departure over a stretch holding few patients can be
discounted.

A line above the diagonal at the left and below it at the right means the
predictions are too extreme, which is the usual failure of an overfitted model;
the opposite pattern means they are too timid, which heavy shrinkage
produces.

Each patient is judged here against a probability estimated without them, which
separates this from the in-sample version, `type = "calibration"`. That
one reads optimistically, by however much the model has fitted the individual
patients rather than the pattern. `fitted()` returns the probabilities
themselves, for binning them some other way.

### Residuals

For a continuous outcome, residuals plotted against fitted values read the way
they would for a linear model, with one difference. A forest shrinks its
predictions toward the overall mean, so a mild positive trend is expected even
when the model is correct: the highest fitted values are pulled down and the
lowest pulled up. A strong slope is not expected, and a fan shape means the
spread of the outcome depends on the predictors, which `gaussian_ls()` models
directly.

For a binary outcome, the residuals take two values for any given fitted
probability and the plot is not informative; the calibration check above should
be used instead.

### The Bayesian R-Squared (`r2()`)

```{r r2, eval = run && rlang::is_installed("performance")}
performance::r2(fit)
```

This is the Bayesian $R^2$ of @gelman2019, computed from the posterior rather
than from a single fit, so it comes with an interval and cannot exceed 1. For a
binary outcome it is bounded well below 1 by the outcome's own randomness, so it
is best read as a relative measure rather than against any absolute standard.

## A Checklist

Before fitting in earnest we read the prior back with `prior_summary()` and run
the model once with `prior_only = TRUE`, checking that the replicates fall where
the outcome can actually go, since a prior worth changing is cheapest to change
before any of the data's information has been spent. Then, before reporting anything from a fit, we run it with
`chains = 4` and pass it to
`diagnose()`, reading what that says to change: it applies the thresholds, names
the rows that fall short, says when an R-hat rests on too few effective draws to
be read at all, and separates a warmup that ended too early
from chains that have each genuinely settled somewhere different. We then run
`pp_check()` and look for systematic differences, checking in particular that
the replicates respect any bound the outcome has (e.g., a count that cannot go
below zero), and for a binary outcome we read a calibration plot instead, since
the predictive check has little to say there.

## Further Reading

`vignette("comparison")` covers choosing between models once each of them fits,
and `vignette("families")` covers what to do when the posterior predictive check
says the family is wrong.

## References