Code
?bace
BACE: Bayesian Augmentation by Chained Equations for Comparative Biology
1 Introduction
We provide a brief introduction to the BACE package for Bayesian phylogenetic and correlated trait imputation. The BACE package allows users to impute missing data in comparative datasets using Bayesian methods that account for phylogenetic relationships among species.
2 When to use BACE?
General-purpose multiple-imputation tools such as mice, Amelia and missForest assume that rows are independent: the missing value for row i is imputed using the marginal distribution of the column and the observed covariates of row i, with no structural relationship between rows. For comparative biology datasets this assumption is wrong. Closely related species tend to share trait values because they share evolutionary history.
BACE keeps that information by fitting every per-variable imputation model as a Bayesian phylogenetic mixed model via MCMCglmm, with a phylogenetic random effect that borrows information between related taxa. Imputation is done via chained equations — one conditional model per variable with missing data — and final inference uses multiple independent imputations whose posterior draws are combined (see Section 6.6).
BACE is the right tool when:
- Your rows are species (or populations nested within species) and you have a phylogenetic tree that covers them.
- Your traits are likely to show phylogenetic signal — i.e. closely related species have correlated values. If you are unsure, fit a single-variable phylogenetic model to a complete trait first (e.g.
phytools::phylosig()or anMCMCglmmfit with aphylorandom effect) and check whether the phylogenetic variance component is non-trivial. We also provide aphylo_signal_summary()function for this very check. - You want posterior inference on downstream models rather than point imputations followed by ad-hoc variance adjustments.
- Your missing-data pattern spans mixed variable types (continuous, binary, categorical, ordinal, count), but it of course can be used when variables are all of one type too (e.g., continuous).
BACE is probably the wrong tool when:
- Rows are independent observational units with no phylogenetic structure — generic MI tools are faster and simpler.
- The trait of interest has no detectable phylogenetic signal. The phylogenetic random effect is then noise that slows down MCMC without helping.
- Missingness is rare and plausibly MCAR — a complete-case analysis may be defensible and much easier.
- You need fast imputations. Each chained-equations iteration runs an MCMC fit, so runtimes are measured in minutes, not milliseconds.
3 Installation
You can install the BACE package from GitHub using the following command:
If you encounter any issues during installation, please ensure that you have the required dependencies installed.
4 Documentation
The main user-facing function is bace(), which runs the whole imputation workflow end-to-end. Under the hood it chains four steps — bace_imp(), assess_convergence(), bace_final_imp() and pool_posteriors() — and each of those can also be called directly when you need finer control (see Section 8); for almost all analyses, though, one call to bace() is all you need. You can access the documentation using the following commands:
As you will see, the function allows you to specify fixed effects, random effects, the phylogenetic tree, and other parameters for the imputation process. At the bare minimum, you need to provide a formula for the fixed effects, a formula for the random effects (including phylogeny), a phylogenetic tree, and the dataset with missing values.
At present, BACE only works with intercept-only phylogenetic and species-level random effects, but we are working on adding support for other types of random effects in the future (e.g., random slopes). For most comparative datasets, this will be sufficient for modeling the phylogenetic and species-level variation in the data, but not always. BACE is particularly powerful for when you have non-gaussian missingness as it can handle variables of all kinds of different distribution families, unlike other packages.
5 The BACE workflow at a glance
Fig. 1 maps the package functions onto the workflow. bace() is the top-level wrapper function. Internally it glues together four steps that you can also run individually when you need finer control:
- Initial chained imputation —
bace_imp(). For each variable with missing data, fit aMCMCglmmphylogenetic mixed model, impute the missing cells, move on to the next variable, and repeat forrunsiterations. At each iteration the working dataset is updated on the raw scale so that later draws condition on the most recent imputations for every other variable. The first iteration warm-starts predictor NAs with marginal means (for continuous variables) or empirical draws (for categorical variables); subsequent iterations use the chained draws. Crucially, imputations in this phase are deterministic point predictions (posterior means / modal classes) — that is what lets the chained loop settle to a stable state. - Convergence assessment of the chained loop —
assess_convergence(). Checks whether the sequence of imputed datasets has stabilised. This is not within-chain MCMC convergence — see Section 9.2. It’s useful to know whether the imputed values for each variable have settled, because the final imputation runs start from this converged dataset. - Final imputation runs —
bace_final_imp(). Starting from the converged dataset, launchn_finalindependent MCMCglmm fits per variable (defaultn_final = 50; parallelisable vian_cores). Each run produces one complete imputed dataset and one per-variable posterior distribution. In contrast to phase 1, missing cells are re-drawn from the posterior predictive distribution (not a point estimate), so the returned imputed datasets differ from run to run. This between-imputation variability is what makes BACE proper multiple imputation — without it, downstream intervals would ignore missing-data uncertainty. When more than one variable has missing data, each run additionally cycles through the variables several times (thesweepsargument, default 3), so that every variable — including the first one in your formula — is imputed conditional on fresh random draws of the other variables rather than on fixed point estimates; only the last pass is kept. With a single incomplete variable one pass is exact, and BACE skips the extra passes automatically. - Posterior pooling —
pool_posteriors(). Concatenate the per-imputation MCMC chains into a single set of posterior draws per parameter. These stacked draws approximate the marginal posterior of the analysis model integrating over the imputation distribution.
bace() runs all four steps and additionally retries step 1 with more chained iterations if step 2 fails (up to max_attempts times). As Fig. 1 shows, there is also a second, model-agnostic route from step 3’s imputed datasets to pooled inference: fit any downstream model to each imputed dataset with with_imputations() and combine the estimates with Rubin’s rules via pool_mi() — see Section 7. The rest of this vignette walks through the workflow first with bace() and then shows how to run each step individually.
bace_imp()), chained-loop convergence checking (assess_convergence()), posterior-predictive final imputations (bace_final_imp()), and posterior pooling (pool_posteriors()), with accessors (get_pooled_model(), get_imputed_data(), get_imputation_intervals()) reading the finished object. The light band marks the steps orchestrated by the bace() wrapper, including its automatic retry of the chained phase on non-convergence; every downstream box takes the finished bace() result directly, and bace_imp()/bace_final_imp() are called individually only in the two-phase workflow for very large datasets. The lower-right branch is the alternative Rubin’s-rules pathway: with_imputations() fits any downstream model to each imputed dataset and pool_mi() combines the estimates. A finished run can be validated with bace_validate(), which hides known cells, re-imputes them with the run’s own data, tree and settings, and scores the recovery. phylo_signal_summary() screens phylogenetic signal per variable before imputation, and sim_bace() simulates test data with known truth for benchmarking. Dashed boxes and arrows show optional advanced controls: per-model MCMC control and diagnostics, and the convergence-plot family.
6 Example Usage
6.1 Simulating Data
We’ll first simulate a dataset with missing values and a phylogenetic tree. We can do this using the sim_bace() function.
Code
# Set seed for reproducibility
set.seed(123)
# Simulate data. Note that threshold4 indicates an ordinal variable with 4 levels, and multinomial3 indicates a categorical variable with 3 levels. You can replace the numbers
sim_data <- sim_bace(response_type = "poisson",
predictor_types = c("gaussian", "binary", "threshold4", "multinomial3"),
phylo_signal = c(0.55, 0.75, 0.8, 0.7, 0.7),
missingness = c(0.2, 0.3, 0.4, 0.2, 0.5),
n_cases = 100,
n_species = 100)Generated default beta_matrix with pre-set sparsityGenerated default beta_resp: -0.043, 0.457, -0.047, 0.178We can then view the simulated data and phylogenetic tree because the sim_bace() function returns a list containing both objects. It also has a lot of other useful information about the simulation users can explore.
Code
# Object contains both the simulated data and phylogeny
data <- sim_data$data
tree <- sim_data$tree
# View first few rows of the data
head(data) Species y x1 x2 x3 x4
1 pJUhd NA -2.72232649 <NA> 3 <NA>
2 h62z2 NA 0.07569633 <NA> <NA> <NA>
3 69Haj 0 -3.04527643 <NA> <NA> B
4 gv77g 3 NA <NA> 3 C
5 Arhl9 2 -1.78548552 <NA> 3 A
6 wLS5r 0 NA 0 2 <NA>Here, we can see that the dataset contains lots of missing values across different variable types. Of course, normally, we would need to exclude these cases to use a complete dataset. However, with BACE we simply feed in the entire dataset, specify the models we want to use for imputation (or the real model of interest), and impute the missing data.
6.2 Pre-flight: check phylogenetic signal
Before committing to a full BACE run (which can take minutes to hours for larger datasets), it’s worth checking whether the variables carry enough phylogenetic information for the phylogenetic random effect to help imputation. BACE provides phylo_signal_summary(), which fits one univariate phylogenetic mixed model per variable and reports:
- H² — phylogenetic heritability (fraction of variance attributable to the phylogeny) on the latent scale that BACE’s per-variable models actually operate on. Posterior mean + 95% HPD. Works for all variable types (gaussian, poisson, binary, ordered, multinomial).
- λ, K (continuous variables only) — Pagel’s (Pagel 1999) and Blomberg’s (Blomberg, Garland, and Ives 2003) classical statistics, via
phytools::phylosig(). - D (binary variables only) — Fritz-Purvis D (Fritz and Purvis 2010), via
caper::phylo.d(). - Convergence flags — every fit is checked with
coda::effectiveSize()andcoda::geweke.diag(); rows with low effective sample size or failed Geweke tests are flaggedunreliableand marked with an asterisk in the printout. Do not interpret flagged rows as data-driven estimates.
The same H² formulas that BACE’s own per-variable models imply are used here, so the number is directly informative about the imputation model’s structure. The function automatically updates these statistics depending on BACE’s two-RE modes: species = FALSE fits a single phylogenetic random effect; species = TRUE fits the dual structure (phylogenetic + non-phylogenetic species) which BACE uses when there are within-species replicates.
Code
# Quick first-pass signal check using quick = TRUE (halves iterations):
signal <- phylo_signal_summary(
data = data,
tree = tree,
species_col = "Species",
species = FALSE,
quick = TRUE
)
print(signal)=============================================================
Phylogenetic signal summary
=============================================================
Random-effect structure: single (phylo only)
Variables : 5
--- Warnings ---
[quick_mode] quick = TRUE: nitt and burnin halved. Expect more low_ess flags.
--- Table ---
variable type n_obs H2_mean H2_lo H2_hi lambda K D
y poisson 80 0.125* 5.89e-08 0.509 7.35e-05 0.0669 NA
x1 gaussian 70 0.713* 4.89e-01 0.904 6.93e-01 0.0665 NA
x2 threshold 60 0.728* 4.34e-02 0.999 NA NA 0.746
x3 threshold 80 0.998* 9.95e-01 1.000 NA NA NA
x4 categorical 50 0.764* 3.07e-01 0.982 NA NA NA
min_ess interpretation flag
178.60 low low_ess;geweke_fail;unreliable
750.00 high low_ess;unreliable
3.51 high low_ess;unreliable
3.73 high low_ess;unreliable
46.43 high low_ess;geweke_fail;unreliable
Interpretation bands: H2 < 0.2 = low, 0.2-0.5 = moderate, >0.5 = high.
Flag meanings:
low_ess -- effective sample size below threshold (default 1000)
geweke_fail -- Geweke diagnostic |z| > 2 on any variance component
unreliable -- H2 cell should NOT be trusted (low_ess or geweke_fail)
low_n -- n_obs < 20; posterior is data-sparse
fit_error -- MCMCglmm fit failed (check warnings)
Note: phylogenetic signal is necessary but NOT sufficient for good
imputation. A high-H2 trait can still impute poorly under MNAR missingness;
a low-H2 trait can still impute well if coupled to an observed predictor.Or, more concisely, call bace() with phylo_signal = TRUE to get the same diagnostic without running the imputation pipeline:
Code
# Signal-only preview: returns the same phylo_signal table as above and
# does NOT run the imputation pipeline (so we don't evaluate it twice here).
preview <- bace(
fixformula = y ~ x1 + x2 + x3 + x4,
ran_phylo_form = ~ 1 | Species,
phylo = tree,
data = data,
phylo_signal = TRUE
)How to interpret the table:
- H² < 0.2 = low signal — the phylogenetic random effect contributes little for this variable; imputation will rely mostly on observed predictors.
- 0.2 ≤ H² < 0.5 = moderate signal.
- H² ≥ 0.5 = high signal — the phylogenetic random effect carries substantial information; BACE’s phylogenetic structure is expected to help.
Important
Phylogenetic signal is necessary but not sufficient for good imputation. A high-H² trait can still impute poorly if the data are missing-not-at-random (MNAR); a low-H² trait can still impute well if it is tightly coupled to an observed predictor. The signal table sets expectations; it does not predict imputation accuracy.
If you have weak phylogenetic signal in some variables all is not lost so long as you have auxillary variables you can include in the imputation model.
6.3 Using BACE for Imputation
Now that we have the data and the phylogenetic tree we can run the imputation. For almost all users this means a single call to bace(); the sub-functions such as bace_imp() exist for finer stepwise control and are covered in Section 8.
Caution
It’s very important to ensure that all variables are correctly classified in the dataset (e.g., factors for categorical variables, numeric for continuous variables, etc.) before running the imputation. BACE will fit different models for different variable types, and so, if they have not been correctly classified you will be modeling the variable incorrectly.
BACE decides the imputation family for each variable deterministically from how the column is stored:
| Stored as | Detected type | Imputation model |
|---|---|---|
| numeric (double) | gaussian | Gaussian phylogenetic mixed model |
| integer, all values ≥ 0 | count | Poisson (log link), count-specific priors |
| integer with negative values | gaussian | treated as integer-coded continuous |
| factor with 2 levels | binary | binary threshold (probit) model |
| ordered factor with ≥ 3 levels | ordinal | ordinal threshold model |
| unordered factor with ≥ 3 levels | categorical | multinomial probit model |
Two consequences of these rules deserve emphasis:
- Integer storage means counts. A column stored as
integerwhose values are all non-negative is imputed as a Poisson count. If you have an integer-coded variable that is really continuous (e.g. body mass recorded to the nearest gram), store it as a double (as.numeric()) before running BACE. Conversely, true counts must be stored as integers to get the Poisson model. - Unordered factors with ≥ 3 levels get a true multinomial probit model by default (
ovr_categorical = FALSE). Settingovr_categorical = TRUEinstead fits one-vs-rest binary threshold models per category — this is retained for backwards compatibility, but the multinomial parameterisation is both better calibrated and substantially faster in our simulations, so we recommend the default.
If you want to see what BACE has classified each variable as, you can inspect the types element of any fitted bace_imp() or bace() object. For example, after running the simple workflow below you can do:
Code
# After running bace() (see next section), inspect the inferred variable types
bace_final1$initial_results$typesImportantly, users also need to think carefully about the distributions of the variables themselves and transform them if necessary prior to running bace_imp (e.g., log-transforming skewed continuous variables).
There are two ways you can use bace/bace_imp. The first is to specify a list of all the formulas for each variable you would like to fit. It is important here that, if you provide a list, that any variable with missing data is found as a response variable. This must be the case because BACE will not automatically build models for other variables if you provide a list of formulas (but see below).
For now, we’ll focus on the simplified bace workflow because it does the full process for you, but later we will show how different components of the workflow can be run separately for more control over the imputation process.
6.4 Simple BACE Workflow
The bace() function is the main function to conduct a BACE analysis running each of the core steps: 1) initial imputations to reach convergence in imputed data; 2) convergence assessment to ensure that imputed data has reached a stable state; 3) imputing missing data n_final times (50 by default) by drawing from the posterior predictive distribution of the final models, re-fitting and estimating parameters and then pooling posteriors across the final runs.
If convergence is not achieved, it will automatically restart the initial steps to allow chains to converge before running the final set of models. You can also specify the number of runs, MCMC parameters, and whether to plot the results. The final output will include pooled models that account for imputation uncertainty.
Code
# Run the full bace workflow. Note: we set n_final = 10 here only to keep this
# vignette quick to build; for real analyses the default (n_final = 50) is a
# better choice (see the n_final bullet below).
bace_final1 <- bace(fixformula = "y ~ x1 + x2",
ran_phylo_form = "~ 1 | Species",
phylo = tree,
data = data,
species = FALSE,
runs = 15, nitt = 100000, burnin = 10000, thin = 10, n_final = 10, plot = TRUE, skip_conv = FALSE, sample_size = 1000)
=======================================================
BACE: Bayesian Augmentation by Chained Equations
=======================================================
Step 1: Initial imputation for convergence assessment
-------------------------------------------------------
Run 1, Imputed variable: y
Run 1, Imputed variable: x1
Run 1, Imputed variable: x2
Run 2, Imputed variable: y
Run 2, Imputed variable: x1
Run 2, Imputed variable: x2
Run 3, Imputed variable: y
Run 3, Imputed variable: x1
Run 3, Imputed variable: x2
Run 4, Imputed variable: y
Run 4, Imputed variable: x1
Run 4, Imputed variable: x2
Run 5, Imputed variable: y
Run 5, Imputed variable: x1
Run 5, Imputed variable: x2
Run 6, Imputed variable: y
Run 6, Imputed variable: x1
Run 6, Imputed variable: x2
Run 7, Imputed variable: y
Run 7, Imputed variable: x1
Run 7, Imputed variable: x2
Run 8, Imputed variable: y
Run 8, Imputed variable: x1
Run 8, Imputed variable: x2
Run 9, Imputed variable: y
Run 9, Imputed variable: x1
Run 9, Imputed variable: x2
Run 10, Imputed variable: y
Run 10, Imputed variable: x1
Run 10, Imputed variable: x2
Run 11, Imputed variable: y
Run 11, Imputed variable: x1
Run 11, Imputed variable: x2
Run 12, Imputed variable: y
Run 12, Imputed variable: x1
Run 12, Imputed variable: x2
Run 13, Imputed variable: y
Run 13, Imputed variable: x1
Run 13, Imputed variable: x2
Run 14, Imputed variable: y
Run 14, Imputed variable: x1
Run 14, Imputed variable: x2
Run 15, Imputed variable: y
Run 15, Imputed variable: x1
Run 15, Imputed variable: x2
========================================
BACE MCMC Diagnostics Summary
========================================
Number of models: 3
Total parameters: 15
Fixed effects: 9
Variance components: 6
Fixed units (threshold/categorical): 1
Mean Effective Sample Size:
Fixed effects: 4277.4
Variance components: 2948.3
Convergence Assessment:
Acceptable : 2 parameters (13.3%)
Check : 6 parameters (40.0%)
Fixed : 1 parameters (6.7%)
Good : 6 parameters (40.0%)
Overall: CHECK - Many parameters show convergence issues
Consider increasing nitt, adjusting priors, or checking model specification
Convergence by Response Variable:
----------------------------------------
[~] y : 3/ 5 parameters good
[CH] x1 : 3/ 5 parameters good
Check convergence: Species, units
[CH] x2 : 0/ 4 parameters good (1 fixed)
Check convergence: (Intercept), x1, y, Species
========================================
Step 2: Convergence assessment
-------------------------------------------------------
Step 3: Final imputation runs for posterior pooling
-------------------------------------------------------
=== Running 10 final imputation iterations ===
Run 1/10 (sweep 1/3) - Imputed variable: y
Run 1/10 (sweep 1/3) - Imputed variable: x1
Run 1/10 (sweep 1/3) - Imputed variable: x2
Run 1/10 (sweep 2/3) - Imputed variable: y
Run 1/10 (sweep 2/3) - Imputed variable: x1
Run 1/10 (sweep 2/3) - Imputed variable: x2
Run 1/10 (sweep 3/3) - Imputed variable: y
Run 1/10 (sweep 3/3) - Imputed variable: x1
Run 1/10 (sweep 3/3) - Imputed variable: x2
Run 2/10 (sweep 1/3) - Imputed variable: y
Run 2/10 (sweep 1/3) - Imputed variable: x1
Run 2/10 (sweep 1/3) - Imputed variable: x2
Run 2/10 (sweep 2/3) - Imputed variable: y
Run 2/10 (sweep 2/3) - Imputed variable: x1
Run 2/10 (sweep 2/3) - Imputed variable: x2
Run 2/10 (sweep 3/3) - Imputed variable: y
Run 2/10 (sweep 3/3) - Imputed variable: x1
Run 2/10 (sweep 3/3) - Imputed variable: x2
Run 3/10 (sweep 1/3) - Imputed variable: y
Run 3/10 (sweep 1/3) - Imputed variable: x1
Run 3/10 (sweep 1/3) - Imputed variable: x2
Run 3/10 (sweep 2/3) - Imputed variable: y
Run 3/10 (sweep 2/3) - Imputed variable: x1
Run 3/10 (sweep 2/3) - Imputed variable: x2
Run 3/10 (sweep 3/3) - Imputed variable: y
Run 3/10 (sweep 3/3) - Imputed variable: x1
Run 3/10 (sweep 3/3) - Imputed variable: x2
Run 4/10 (sweep 1/3) - Imputed variable: y
Run 4/10 (sweep 1/3) - Imputed variable: x1
Run 4/10 (sweep 1/3) - Imputed variable: x2
Run 4/10 (sweep 2/3) - Imputed variable: y
Run 4/10 (sweep 2/3) - Imputed variable: x1
Run 4/10 (sweep 2/3) - Imputed variable: x2
Run 4/10 (sweep 3/3) - Imputed variable: y
Run 4/10 (sweep 3/3) - Imputed variable: x1
Run 4/10 (sweep 3/3) - Imputed variable: x2
Run 5/10 (sweep 1/3) - Imputed variable: y
Run 5/10 (sweep 1/3) - Imputed variable: x1
Run 5/10 (sweep 1/3) - Imputed variable: x2
Run 5/10 (sweep 2/3) - Imputed variable: y
Run 5/10 (sweep 2/3) - Imputed variable: x1
Run 5/10 (sweep 2/3) - Imputed variable: x2
Run 5/10 (sweep 3/3) - Imputed variable: y
Run 5/10 (sweep 3/3) - Imputed variable: x1
Run 5/10 (sweep 3/3) - Imputed variable: x2
Run 6/10 (sweep 1/3) - Imputed variable: y
Run 6/10 (sweep 1/3) - Imputed variable: x1
Run 6/10 (sweep 1/3) - Imputed variable: x2
Run 6/10 (sweep 2/3) - Imputed variable: y
Run 6/10 (sweep 2/3) - Imputed variable: x1
Run 6/10 (sweep 2/3) - Imputed variable: x2
Run 6/10 (sweep 3/3) - Imputed variable: y
Run 6/10 (sweep 3/3) - Imputed variable: x1
Run 6/10 (sweep 3/3) - Imputed variable: x2
Run 7/10 (sweep 1/3) - Imputed variable: y
Run 7/10 (sweep 1/3) - Imputed variable: x1
Run 7/10 (sweep 1/3) - Imputed variable: x2
Run 7/10 (sweep 2/3) - Imputed variable: y
Run 7/10 (sweep 2/3) - Imputed variable: x1
Run 7/10 (sweep 2/3) - Imputed variable: x2
Run 7/10 (sweep 3/3) - Imputed variable: y
Run 7/10 (sweep 3/3) - Imputed variable: x1
Run 7/10 (sweep 3/3) - Imputed variable: x2
Run 8/10 (sweep 1/3) - Imputed variable: y
Run 8/10 (sweep 1/3) - Imputed variable: x1
Run 8/10 (sweep 1/3) - Imputed variable: x2
Run 8/10 (sweep 2/3) - Imputed variable: y
Run 8/10 (sweep 2/3) - Imputed variable: x1
Run 8/10 (sweep 2/3) - Imputed variable: x2
Run 8/10 (sweep 3/3) - Imputed variable: y
Run 8/10 (sweep 3/3) - Imputed variable: x1
Run 8/10 (sweep 3/3) - Imputed variable: x2
Run 9/10 (sweep 1/3) - Imputed variable: y
Run 9/10 (sweep 1/3) - Imputed variable: x1
Run 9/10 (sweep 1/3) - Imputed variable: x2
Run 9/10 (sweep 2/3) - Imputed variable: y
Run 9/10 (sweep 2/3) - Imputed variable: x1
Run 9/10 (sweep 2/3) - Imputed variable: x2
Run 9/10 (sweep 3/3) - Imputed variable: y
Run 9/10 (sweep 3/3) - Imputed variable: x1
Run 9/10 (sweep 3/3) - Imputed variable: x2
Run 10/10 (sweep 1/3) - Imputed variable: y
Run 10/10 (sweep 1/3) - Imputed variable: x1
Run 10/10 (sweep 1/3) - Imputed variable: x2
Run 10/10 (sweep 2/3) - Imputed variable: y
Run 10/10 (sweep 2/3) - Imputed variable: x1
Run 10/10 (sweep 2/3) - Imputed variable: x2
Run 10/10 (sweep 3/3) - Imputed variable: y
Run 10/10 (sweep 3/3) - Imputed variable: x1
Run 10/10 (sweep 3/3) - Imputed variable: x2
=== Final imputation complete ===
Saved 10 imputed datasets with corresponding models
--- Final-phase MCMC diagnostics (10 runs; worst case across runs) ---
Flagged parameters (min ESS < 200 or > 25% of runs flagged):
x2 Species [Variance] min ESS 2 (run 10), 100% runs flagged
x2 x1 [Fixed] min ESS 16 (run 10), 100% runs flagged
x2 y [Fixed] min ESS 18 (run 5), 90% runs flagged
x2 (Intercept) [Fixed] min ESS 251 (run 9), 40% runs flagged
Consider longer chains (nitt/nitt_cat_mult) or a per-variable prior
(see ?bace section 'Priors') for the variables above.
Full table: $diagnostics$by_param
Step 4: Pooling posteriors across imputations
-------------------------------------------------------
Using posterior sampling: 1000 samples per imputation
Posterior pooling complete!
Pooled 3 variable model(s) across 10 imputations
=======================================================
BACE Analysis Complete!
=======================================================
Convergence achieved: TRUE
Number of attempts: 1
Final imputations: 10
Access results via:
- get_pooled_model(result, "y"): pooled posterior model for a variable ($pooled_models)
- get_imputed_data(result): list of 10 imputed datasets ($imputed_datasets)
- $convergence: Convergence diagnostics
- $final_results: Full final imputation results (per-run models and datasets)We can now check what BACE classified each variable as to just make sure they have been modelled correctly:
Code
# The variable types inferred by BACE for our simulated data
bace_final1$initial_results$types$y
[1] "poisson"
$x1
[1] "gaussian"
$x2
[1] "threshold"As we can see bace provides a comprehensive overview of the imputation process including each of the core steps. This includes the initial set of runs to allow data to converge, followed by convergence assessment checks both across runs and within parameters of models.
Note
If convergence is not achieved, the function will automatically restart the initial steps to allow chains to converge before running the final set of models. However, users should be aware that the check can sometimes fail even when convergence has likely been achieved. If you think this is the case then you can do the initial runs and skip the convergence assessment (see below). It is also not always the case that parameter checks signal problems with the models. Users will need to explore final models carefully before adjusting MCMC parameters or the number of runs. More on this below.
It’s worth highlighting a few key arguments within bace.
- First, the
skip_convargument allows users to skip the convergence checks and just run the imputation for the number of runs specified. This is not recommended, but it can be useful if you have a large dataset and want to run the imputation for a long time without checking convergence. - Second, the
sample_sizeargument allows users to specify the number of posterior samples to pool across final runs. This can be useful if you have a large dataset and want to reduce the size of the finalbaceobject. - Third, the
speciesargument allows users to specify whether to include a non-phylogenetic random effect for species identity in addition to the phylogenetic random effect. This is important if multiple replicates of species are included in the dataset (e.g., multiple populations of the same species). - Fourth, the
n_finalargument sets how many independent posterior-predictive imputations are drawn in the final phase (default 50). More imputations give a better Monte Carlo approximation to the marginal posterior and tighter calibration of per-cell prediction intervals, at a linear cost in compute. The final runs are completely independent of one another, so they can be run simultaneously: on multi-core machines setn_coresto spread them across cores. Supplyingseedgives each run its own random-number stream, making the results reproducible and identical.
Once the imputation is complete, you can access the final pooled models and evaluate them as you would any MCMCglmm model for each of the models fit. This is important for checking model fit and convergence along with downstream inference. For example, you can check the summary of the model for the response variable y as follows:
Code
# Extract the pooled model for the response variable using get_pooled_model()
y_model <- get_pooled_model(bace_final1, variable = "y")
summary(y_model)
+----------------------------------------------------------------+
| BACE Pooled MCMCglmm summary (stacked per-imputation draws) |
+----------------------------------------------------------------+
Stacked from 10 imputations
Total posterior draws: 10000
(= 1000 sampled per imputation x 10 imputations)
Original samples per imputation: 9000
Posterior summaries below are computed from stacked draws that
approximate p(theta | Y_obs) (Monte Carlo integration over imputed
values). Validity requires each per-imputation chain to have
converged on its own; check a subset via
plot(final$all_models[[1]]$<var>) and coda::effectiveSize().
----------------------------------------------------------------
Iterations = 1:10000
Thinning interval = 1
Sample size = 10000
DIC: 257.6983
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 0.1096 1.106e-08 0.4492 9548
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.5067 0.07734 0.9672 1688
Location effects: y ~ x1 + x2
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 0.06241 -0.38314 0.49030 1645.96 0.720
x1 -0.02012 -0.32049 0.28775 287.36 0.886
x2.L 0.52301 0.06323 0.97196 90.42 0.019 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
plot(y_model) 
Code
# We can now extract posteriors in the same way we would use a normal MCMCglmm model, but now we have the pooled results that account for imputation uncertainty. For example, we can extract the posterior samples for the fixed effects as follows:
hist(y_model$Sol[, "x1"], main = "Posterior distribution for x1", xlab = "x1 coefficient")
Code
# The same accessor without the 'variable' argument returns the full named
# list of pooled models (one per imputed variable) - so all_models$y is
# identical to y_model above.
all_models <- get_pooled_model(bace_final1)
names(all_models)[1] "y" "x1" "x2"The banner above reports whether the chained loop passed the convergence check (and, if not, how many retry attempts bace() made). Whatever the banner says, it is worth looking at the trace and density plots of the pooled MCMCglmm models: variance components in particular can mix slowly, and if they look poorly converged it’s best to adjust the bace MCMC parameters (nitt, burnin, thin) to allow for longer runtimes and better mixing before trusting downstream inference.
Another important feature is that users can also access all the final datasets with imputed values for each run. This can be useful for checking the imputed values (i.e., cross-validation) and exploring the imputation process in more detail. You can access the imputed datasets using the get_imputed_data() function, which returns them as either a list or a single stacked data frame:
Code
# Access the imputed datasets as a list (one data frame per imputation run)
imputed_datasets <- get_imputed_data(bace_final1, format = "list")
length(imputed_datasets)[1] 10Code
head(imputed_datasets[[1]]) y x1 x2 Species
1 1 -2.72232649 0 pJUhd
2 4 0.07569633 0 h62z2
3 0 -3.04527643 1 69Haj
4 3 -4.14199517 1 gv77g
5 2 -1.78548552 1 Arhl9
6 0 0.55217341 0 wLS5rCode
# Or stack all imputed datasets into one data frame with an .imputation column
imputed_df <- get_imputed_data(bace_final1, format = "data.frame")
table(imputed_df$.imputation)
1 2 3 4 5 6 7 8 9 10
100 100 100 100 100 100 100 100 100 100 For reporting individual imputed values with honest uncertainty, use get_imputation_intervals(). For every imputed cell it gives the best guess (the mean across the n_final imputations) and a calibrated 95% interval built with the standard multiple-imputation formula [Rubin’s t-interval; Rubin (1987)].
Code
head(get_imputation_intervals(bace_final1)) row species variable mean sd lower upper n_imputations
1 1 pJUhd y 1.3 1.059350 -1.213382 3.813382 10
2 2 h62z2 y 3.0 2.708013 -3.424951 9.424951 10
3 14 4MBXV y 0.6 0.843274 -1.400727 2.600727 10
4 17 jdT90 y 0.6 0.843274 -1.400727 2.600727 10
5 20 cwSPe y 1.9 1.595131 -1.884562 5.684562 10
6 38 JE8PP y 1.7 2.057507 -3.181579 6.581579 10These are exactly the intervals to put in a table of imputed trait values, or to draw as error bars when plotting imputed alongside observed values.
6.5 Propagating phylogenetic (tree) uncertainty
Nobody knows the true phylogeny — there is always a multitude of plausible trees. Imputing with a single tree makes your results look more certain than they are. BACE handles this with one argument change: pass the whole tree set as phylo, and each of the n_final final imputation runs is fit under its own tree. The spread between imputations — the quantity every pooling route uses — then automatically includes tree uncertainty. This is standard multiple-imputation logic with the tree treated as missing information, following Nakagawa and de Villemereuil (2019), who also show a modest tree sample (~50) is typically enough.
Code
# trees: a multiPhylo, e.g. read with ape::read.tree()/read.nexus()
res <- bace(
fixformula = "y ~ x1 + x2",
ran_phylo_form = "~ 1 | Species",
phylo = trees, # the only change: a tree SET, not one tree
data = mydata
)
# which tree each imputation run used:
res$final_results$tree_indexImportantly, you do not need as many trees as imputations. When there are fewer trees than runs they are recycled in a balanced rotation. With more trees than runs an evenly spaced subset is used. The convergence phase runs on the first tree only (its job is to let point imputations settle, which needs a fixed target). When trees differ across runs, pool_posteriors() aligns species-level effects by name and drops the internal-node effects (node identities differ between trees) — fixed effects, variance components and cutpoints pool exactly as usual.
6.6 Interpreting the pooled output
It is worth being precise about what the pooled posterior returned by bace() actually is, because the terminology in the multiple-imputation literature can be confusing.
pool_posteriors() concatenates the MCMC draws from the n_final per-imputation fits into a single set of draws per parameter. Each per-imputation chain targets the conditional posterior \(p(\theta \mid Y_{\text{obs}}, Y_{\text{mis}}^{(m)})\), where \(Y_{\text{mis}}^{(m)}\) is the m-th complete imputed dataset. Stacking these chains is an approximation to the marginal posterior \[
p(\theta \mid Y_{\text{obs}}) = \int p(\theta \mid Y_{\text{obs}}, Y_{\text{mis}})\, p(Y_{\text{mis}} \mid Y_{\text{obs}})\, dY_{\text{mis}},
\] which integrates over uncertainty in the imputed values as well as over parameter uncertainty. Posterior summaries — means, quantiles, HPD intervals — computed from the stacked draws therefore reflect both sources of uncertainty simultaneously. This is the standard Bayesian-MI combiner; see Zhou and Reiter (2010) and Chapter 18 of Gelman et al. (2013).
Two practical notes follow from this:
- This is not Rubin’s rules. Rubin’s rules combine scalar point estimates and variances from each imputation using a specific within/between-imputation variance decomposition. BACE does provide a Rubin’s-rules pathway if desired —
with_imputations()+pool_mi()— see Section 7. - Validity depends on per-imputation MCMC convergence. Because each per-imputation chain is supposed to target its own conditional posterior, the stacked draws inherit any bias from chains that have not converged. Before trusting pooled summaries, spot-check a handful of the per-imputation fits with standard diagnostics. For example:
Code
# Inspect MCMC diagnostics on the first per-imputation fit for 'y'.
# Sol also stores a sample for every species-level effect (pr = TRUE),
# so we restrict the output to the fixed effects (the first nfl columns);
# drop the subsetting to inspect every parameter.
first_fit <- bace_final1$final_results$all_models[[1]]$y
fixed_sols <- coda::as.mcmc(first_fit$Sol[, 1:first_fit$Fixed$nfl, drop = FALSE])
plot(fixed_sols)
Code
coda::effectiveSize(fixed_sols)(Intercept) x1 x2.L
3407.018 5199.212 3874.261 Code
coda::autocorr.diag(fixed_sols) (Intercept) x1 x2.L
Lag 0 1.00000000 1.000000000 1.000000e+00
Lag 10 0.27232137 0.223192668 3.102074e-01
Lag 50 0.04301813 -0.002202535 1.445306e-02
Lag 100 0.01242483 0.014935373 3.519735e-02
Lag 500 0.01843817 0.014438562 7.525545e-05How to read this output. The trace plots should look like fuzzy caterpillars with no drift or sticky patches, and the density plots should be smooth and unimodal. effectiveSize() reports how many independent draws the autocorrelated chain is worth: values close to the number of stored samples indicate good mixing, while values below a few hundred mean the posterior summaries for that parameter are unstable and the chain should be run longer. autocorr.diag() should drop towards zero within a few lags; values that stay high at large lags are the same warning in a different form. If these look poor, increase nitt and/or burnin and re-run bace() rather than simply trusting the pooled banner.
7 Downstream analysis with Rubin’s rules
The pooled MCMCglmm models returned by bace() are the imputation models themselves. Often, though, the model you actually care about is a different downstream analysis — a phylogenetic GLS on the completed data, a GLM with a different predictor set, a phylolm fit — and you want that analysis to propagate imputation uncertainty. This is the classical multiple-imputation setting, and BACE supports it directly with two functions (the lower-right branch of Fig. 1):
with_imputations(object, .f)applies your model-fitting function.fto each of then_finalimputed datasets and returns the list of fits. It accepts the output ofbace()orbace_final_imp()(or a plain list of completed data frames), shows progress, and captures per-fit errors rather than aborting the whole loop.pool_mi(fits)combines the per-imputation coefficient estimates with Rubin’s rules (Rubin 1987): the pooled variance is the within-imputation variance plus (1 + 1/M) times the between-imputation variance, degrees of freedom use the Barnard and Rubin (1999) small-sample correction when the complete-data residual df is available, and the output reports the fraction of missing information (fmi) and relative increase in variance (riv) per coefficient. The implementation matchesmice::poolto machine precision on identical inputs.
Continuing with the bace_final1 object fitted above (recall y is a count):
Code
# 1. Fit a downstream GLM to each imputed dataset
fits <- with_imputations(bace_final1,
function(d) glm(y ~ x1 + x2, family = poisson, data = d))
# 2. Pool with Rubin's rules
pool_mi(fits)Pooled estimates from 10 imputations (Rubin's rules)
Confidence level: 95%
term estimate std.error df statistic p.value conf.low conf.high
(Intercept) 0.3110 0.1447 52.69 2.148 0.03631 0.02059 0.6013
x1 -0.01887 0.04952 25.52 -0.3811 0.7063 -0.1208 0.08301
x2.L 0.4787 0.1387 273.1 3.451 0.0006462 0.2056 0.7517
fmi riv
0.4344 0.7044
0.6223 1.462
0.1874 0.2218
Note: max fmi = 0.62. Consider more imputations for stable SEs.Each row gives the pooled estimate, standard error, Barnard–Rubin degrees of freedom, confidence interval, and the fmi/riv diagnostics. High fmi for a coefficient means much of its uncertainty comes from the missing data itself — a useful, standard MI diagnostic that the stacked-posterior pathway does not report.
Because with_imputations() is model-agnostic, phylogenetic downstream models work too. If your function declares a tree argument, the phylogeny supplied to with_imputations() is passed through to every fit:
Code
# Phylogenetic GLS on each imputed dataset, pooled with Rubin's rules
fits_gls <- with_imputations(
bace_final1,
function(d, tree) {
nlme::gls(y ~ x1 + x2, data = d,
correlation = ape::corBrownian(phy = tree, form = ~Species))
},
tree = tree
)
pool_mi(fits_gls)When to use which combiner:
pool_posteriors()(stacking) — when the per-imputation analysis is Bayesian (MCMCglmm) and you want full marginal posteriors. No normality assumption; the right choice for variance components, phylogenetic signal, and heritability, whose posteriors are skewed and bounded.pool_mi()(Rubin’s rules) — when the per-imputation analysis is frequentist, or when you want the standard MI reporting apparatus (pooled SEs, df,fmi). It assumes the per-imputation estimates are approximately normal, which is reasonable for fixed-effect coefficients but not for variance components — never pool variance components, phylogenetic signal, or heritability withpool_mi(); usepool_posteriors()for those.
pool_mi() will also accept MCMCglmm fits (it takes the posterior mean and covariance of the fixed effects as the per-imputation estimate and variance), which is convenient for computing fmi on a Bayesian analysis — but for coefficient inference from MCMCglmm fits, stacking remains the principled combiner. For the general theory see Rubin (1987) and Chapter 2 of van Buuren (2018).
8 Advanced: running the workflow step by step
For most users bace() is the right entry point — it runs the four-step workflow end-to-end and retries if the chained loop has not converged. But there are cases where you want finer control: you may want to skip convergence assessment, inspect intermediate objects, run more or fewer final imputations than the default, or substitute your own pooling step. Here we show how to call bace_imp(), bace_final_imp() and pool_posteriors() directly.
8.1 Controlling the imputation models with a list of formulas
The first thing you might want to control is which model is fit for each variable with missing data. bace() and bace_imp() both accept a list of formulas, one per variable with missing data. This is useful when you want different predictors for different variables (for example, when a variable should only be imputed from a subset of the others, or when you want to include an interaction in one model but not another).
Note that if you pass a list, every variable with missing data must appear as a response somewhere in the list — BACE will not auto-generate models for variables you have left out.
Code
# First check that all variables are correctly classified. If using sim_bace() this is done automatically, but always worth checking.
str(data)
# Perform BACE imputation. Here we specify models for each variable in the dataset.
bace_impute <- bace_imp(fixformula = list("y ~ x1 + x2",
"x1 ~ x2",
"x2 ~ x4",
"x3 ~ x1 + x2",
"x4 ~ x1 + x2"),
ran_phylo_form = "~ 1 | Species",
phylo = tree,
data = data,
runs = 5, verbose = FALSE)What BACE will do here is fit a MCMCglmm model for each of the response variables specified in the formula list. At the first iteration, it will fill in missing values for any predictors using either their mean (for continuous variables) or randomly sampled values (for categorical variables). Models will then be fit using this complete predictor set and predictions will be made for missing values of the response variable.
In the next run, model predictions will be used to ‘fill’ in missing values in any predictors instead. This cycle will continue until the predicted values converge. This process is repeated for the number of iterations specified by the runs argument. All the models will include the phylogenetic random effect specified in the ran_phylo_form argument.
The second way to use BACE is to just provide a model with your main variables. This may be because you have an a priori model in mind. As such, you can also simply feed in a single formula as below. Here, BACE will automatically build models for the other variables based on their types and the variables present in the dataset. It will impute missing values in the same way as described above, but using all other variables as predictors. You can also specify interactions if you wish. We can run BACE in this way as follows (this time using verbose = TRUE to see more output during the imputation process):
Code
# Fit BACE with a single formula
bace_impute2 <- bace_imp(fixformula = "y ~ x1 + x2",
ran_phylo_form = "~ 1 | Species",
phylo = tree,
data = data,
runs = 5, nitt = 20000, burnin = 5000, thin = 10) Run 1, Imputed variable: y
Run 1, Imputed variable: x1
Run 1, Imputed variable: x2
Run 2, Imputed variable: y
Run 2, Imputed variable: x1
Run 2, Imputed variable: x2
Run 3, Imputed variable: y
Run 3, Imputed variable: x1
Run 3, Imputed variable: x2
Run 4, Imputed variable: y
Run 4, Imputed variable: x1
Run 4, Imputed variable: x2
Run 5, Imputed variable: y
Run 5, Imputed variable: x1
Run 5, Imputed variable: x2
========================================
BACE MCMC Diagnostics Summary
========================================
Number of models: 3
Total parameters: 15
Fixed effects: 9
Variance components: 6
Fixed units (threshold/categorical): 1
Mean Effective Sample Size:
Fixed effects: 780.5
Variance components: 477.8
Convergence Assessment:
Check : 6 parameters (40.0%)
Fixed : 1 parameters (6.7%)
Good : 8 parameters (53.3%)
Overall: CHECK - Many parameters show convergence issues
Consider increasing nitt, adjusting priors, or checking model specification
Convergence by Response Variable:
----------------------------------------
[CH] y : 3/ 5 parameters good
Check convergence: Species, units
[OK] x1 : 5/ 5 parameters good
[CH] x2 : 0/ 4 parameters good (1 fixed)
Check convergence: (Intercept), x1, y, Species
========================================As you can see from the printout above, BACE provides a summary of the imputation process, including the number of models fitted, parameters estimated, and effective sample sizes which evaluate how well the MCMC chains for parameters are mixing. This is important to check to ensure that the imputation is performing well. Note, it is not always the case that these checks signal problems with the models. Users will need to explore final models carefully before adjusting MCMC parameters or the number of runs. More on this below.
8.2 Evaluating Final Models
BACE returns a list of models for each variable that was imputed. You can access these models and evaluate them as you would any MCMCglmm model. This is important for checking model fit and convergence along with downstream inference. For example, you can check the summary of the model for the response variable y as follows:
Code
# Access the model for the response variable 'y' fit on the LAST chained
# iteration. bace_imp() stores these as $models_last_run (one MCMCglmm
# object per response variable with a formula).
model_y <- bace_impute2$models_last_run[["y"]]
# View summary of the model
summary(model_y)
Iterations = 5001:19991
Thinning interval = 10
Sample size = 1500
DIC: 256.6774
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 0.1178 4.842e-11 0.4324 150.8
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.4366 0.03206 0.8267 128.5
Location effects: y ~ x1 + x2
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 0.05625 -0.40965 0.46382 375.4 0.73867
x1 -0.01813 -0.31315 0.24886 827.0 0.90133
x2.L 0.60916 0.24815 0.98636 576.2 0.00267 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
# Plot trace and density plots for model parameters
plot(model_y)
At short chain lengths like these you will often see the phylogenetic variance component (the Species random effect) mixing relatively poorly. If so, we should go back and re-run the imputation with more iterations or possibly adjust the prior to ensure better mixing for this parameter — the priors argument of bace() accepts a per-variable MCMCglmm prior (see the Priors section of ?bace).
9 Phylogenetic and Non-Phylogenetic Random Effects
BACE allows users to include both phylogenetic and non-phylogenetic random effects in the imputation models. This is important if multiple replicates of species are included in the dataset (e.g., multiple populations of the same species). In this case, users can specify a non-phylogenetic random effect for species identity in addition to the phylogenetic random effect. This is done by adjusting the species = FALSE argument in bace_imp, setting it to TRUE.
Code
# Simulate data with multiple replicates per species
set.seed(123)
sim_data2 <- sim_bace(response_type = "poisson",
predictor_types = c("gaussian", "binary", "threshold4", "multinomial3"),
phylo_signal = c(0.55, 0.75, 0.8, 0.7, 0.7),
missingness = c(0.2, 0.3, 0.4, 0.2, 0.5),
rr = TRUE,
n_cases = 100,
n_species = 20,
rr_form = list(species = c("x1")))
# Object contains both the simulated data and phylogeny
data2 <- sim_data2$data
tree2 <- sim_data2$tree
# View first few rows of the data
head(data2) Species y x1 x2 x3 x4
1 P9vI3 3 -0.7731097 1 3 C
2 eJBm5 NA NA <NA> 3 <NA>
3 P9vI3 2 -0.6260633 <NA> 3 <NA>
4 JE8PP 1 -2.9886121 1 3 <NA>
5 eJBm5 NA -4.8716976 <NA> 3 <NA>
6 sA8eL 2 NA <NA> <NA> <NA>Code
# We can also make a table of the number of cases per species to confirm that there are multiple replicates per species
dim(data2)[1] 100 6Code
length(unique(data2$Species))[1] 20We now have a dataset with 100 cases but only 20 species, meaning that there are multiple replicates per species. We can run BACE on this dataset while including a non-phylogenetic random effect for species identity as follows:
Code
bace_impute3 <- bace_imp(fixformula = "y ~ x1 + x2",
ran_phylo_form = "~ 1 | Species",
phylo = tree2,
data = data2,
runs = 5, nitt = 20000, burnin = 5000, thin = 10, species = TRUE) Species effect decomposition enabled: 20 out of 20 species (100%) have replicated observations.Run 1, Imputed variable: y
Run 1, Imputed variable: x1
Run 1, Imputed variable: x2
Run 2, Imputed variable: y
Run 2, Imputed variable: x1
Run 2, Imputed variable: x2
Run 3, Imputed variable: y
Run 3, Imputed variable: x1
Run 3, Imputed variable: x2
Run 4, Imputed variable: y
Run 4, Imputed variable: x1
Run 4, Imputed variable: x2
Run 5, Imputed variable: y
Run 5, Imputed variable: x1
Run 5, Imputed variable: x2
========================================
BACE MCMC Diagnostics Summary
========================================
Number of models: 3
Total parameters: 18
Fixed effects: 9
Variance components: 9
Fixed units (threshold/categorical): 1
Mean Effective Sample Size:
Fixed effects: 978.8
Variance components: 384.9
Convergence Assessment:
Acceptable : 2 parameters (11.1%)
Check : 5 parameters (27.8%)
Fixed : 1 parameters (5.6%)
Good : 10 parameters (55.6%)
Overall: CHECK - Many parameters show convergence issues
Consider increasing nitt, adjusting priors, or checking model specification
Convergence by Response Variable:
----------------------------------------
[~] y : 3/ 6 parameters good
Check convergence: units
[OK] x1 : 5/ 6 parameters good
Check convergence: Species
[CH] x2 : 2/ 5 parameters good (1 fixed)
Check convergence: x1, Species, Species2
========================================Code
# We can view the final model for y to check that both random effects are included
model_y2 <- bace_impute3$models_last_run[["y"]]
summary(model_y2)
Iterations = 5001:19991
Thinning interval = 10
Sample size = 1500
DIC: 268.9279
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 0.164 8.144e-08 0.7422 354.6
~Species2
post.mean l-95% CI u-95% CI eff.samp
Species2 0.2315 4.361e-08 0.6828 210.9
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.3354 0.02608 0.6998 197.5
Location effects: y ~ x1 + x2
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 0.12982 -0.45789 0.74299 995.8 0.6040
x1 -0.06196 -0.44556 0.30930 754.0 0.7120
x2.L 0.41793 -0.03392 0.87136 623.3 0.0693 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note
Why isn’t species = TRUE automatic? BACE could detect replicated species and switch this on by itself, but the dual random-effect structure changes the imputation model, and separating the phylogenetic from the non-phylogenetic species variance is only reliable when a reasonable fraction of species are replicated — with sparse replication the two components are weakly identified. BACE therefore makes this an explicit modelling choice: it errors if you request the decomposition with fewer than two replicated species, and warns when under 30% of species carry replicates. Silent auto-switching could change your model between datasets without you noticing.
Note that in this case we get a message indicating that both phylogenetic and non-phylogenetic random effects are being included in the models because we can decompose the two variance components. If you try to fit a dataset without multiple replicates per species while setting species = TRUE you will get an error, because the two variance components are not separately estimable with one observation per species:
Code
# The first simulated dataset has one row per species (100 cases, 100 species),
# so the phylogenetic and non-phylogenetic species variances cannot be separated:
bace_impute4 <- bace_imp(fixformula = "y ~ x1 + x2",
ran_phylo_form = "~ 1 | Species",
phylo = tree,
data = data,
runs = 5, nitt = 20000, burnin = 5000, thin = 10, species = TRUE)Error in `bace_imp()`:
! Cannot decompose phylogenetic and non-phylogenetic species effects (species=TRUE) with fewer than 2 species having replicated observations. Current data has only 0 replicated species.9.1 Adjusting MCMC Parameters For Convergence
As we can see from the model summary above, some parameters may have low effective sample sizes (ESS) which indicates poor mixing of the MCMC chains. This can be a sign that the imputation is not performing well and that more iterations are needed to achieve convergence. Lets have a look:
Code
model_x1 <- bace_impute3$models_last_run[["x1"]]
summary(model_x1)
Iterations = 5001:19991
Thinning interval = 10
Sample size = 1500
DIC: 78.61307
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 1.496 1.867e-06 3.929 210.5
~Species2
post.mean l-95% CI u-95% CI eff.samp
Species2 0.6099 5.132e-06 1.469 270
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.1352 0.08621 0.1921 1500
Location effects: x1 ~ y + x2
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) -0.24515 -1.85851 0.98097 1183 0.741
y 0.02491 -0.04518 0.08989 1500 0.484
x2.L -0.27774 -0.45311 -0.10215 1500 0.004 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
model_x2 <- bace_impute3$models_last_run[["x2"]]
summary(model_x2)
Iterations = 5001:19991
Thinning interval = 10
Sample size = 1500
DIC: 74.09249
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 2.038 1.427e-08 8.305 164.1
~Species2
post.mean l-95% CI u-95% CI eff.samp
Species2 1.194 9.051e-07 5.044 171.5
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 1 1 1 0
Location effects: x2 ~ x1 + y
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 0.222283 -1.518075 2.452027 1229.3 0.727
x1 -0.690601 -1.884091 0.253856 143.5 0.124
y 0.299006 0.008017 0.617722 880.7 0.040 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1As we can see x1 and x2 models have poor mixing and low ESS for some parameters. Obviously, depending on the model complexity and data size etc, some models will need longer runtimes. This is likely because the imputation process has not yet converged. To address this, we can increase the number of iterations (i.e., nitt), increase the burn-in period (i.e., burnin), and/or adjust the thinning interval (i.e., thin) to improve mixing and convergence of the MCMC chains.
Users can adjust the MCMC parameters (i.e., nitt, thin, and burnin) for each model by providing a list of values for each parameter that corresponds to the order of the formulas provided in the fixformula argument. For example, if you have 5 formulas in your fixformula list, you can provide a list of 5 values for each MCMC parameter as follows:
Code
# Define MCMC parameters for each model
nitt_list <- list(60000, 60000, 105000)
thin_list <- list(10, 10, 10)
# Run BACE with model-specific MCMC parameters
bace_impute5 <- bace_imp(fixformula = "y ~ x1 + x2",
ran_phylo_form = "~ 1 | Species",
phylo = tree2,
data = data2,
runs = 5, nitt = nitt_list, thin = thin_list, burnin = 5000, species = TRUE)Run 1, Imputed variable: y
Run 1, Imputed variable: x1
Run 1, Imputed variable: x2
Run 2, Imputed variable: y
Run 2, Imputed variable: x1
Run 2, Imputed variable: x2
Run 3, Imputed variable: y
Run 3, Imputed variable: x1
Run 3, Imputed variable: x2
Run 4, Imputed variable: y
Run 4, Imputed variable: x1
Run 4, Imputed variable: x2
Run 5, Imputed variable: y
Run 5, Imputed variable: x1
Run 5, Imputed variable: x2
========================================
BACE MCMC Diagnostics Summary
========================================
Number of models: 3
Total parameters: 18
Fixed effects: 9
Variance components: 9
Fixed units (threshold/categorical): 1
Mean Effective Sample Size:
Fixed effects: 3583.3
Variance components: 1383.4
Convergence Assessment:
Acceptable : 8 parameters (44.4%)
Check : 1 parameters (5.6%)
Fixed : 1 parameters (5.6%)
Good : 8 parameters (44.4%)
Overall: ACCEPTABLE - Most parameters converged adequately
Convergence by Response Variable:
----------------------------------------
[~] y : 2/ 6 parameters good
Check convergence: x1
[~] x1 : 4/ 6 parameters good
[~] x2 : 2/ 5 parameters good (1 fixed)
========================================As you can imagine, longer runtimes will mean that bace_imp will need to take longer. The diagnostics still suggest possible issues with mixing for some parameters, so we may need to further increase the number of iterations or adjust the priors (via the priors argument — see ?bace) to improve mixing, however, it’s still good to look at the final models to check the model visually.
9.2 Two kinds of convergence
When working with BACE it is important to keep two distinct convergence questions separate, because they are checked in different ways and the tools that flag problems in one are silent on the other.
Within-chain MCMC convergence of each per-variable
MCMCglmmfit. This is the classical question: has the MCMC sampler explored its target posterior and mixed well? The right tools are the ones you already know — trace plots, effective sample size, autocorrelation, Geweke — applied to an individualMCMCglmmobject.bace_imp()prints a short summary of effective sample sizes and autocorrelation for the last-iteration models whenverbose = TRUE, and you can always callplot()orcoda::effectiveSize()onbace_impute$models_last_run[["y"]]or on any of the per-imputation fits stored inbace_final1$final_results$all_models[[m]][["y"]]. The final phase also checks itself:bace_final_imp()stores a worst-case summary across all runs — covering fixed effects, variance components, and ordinal cutpoints — in its$diagnosticselement and prints any flagged parameters. Failures here are fixed by increasingnitt,burnin, orthin, tightening priors (thepriorsargument), or (less often) simplifying the model.Convergence of the chained-equations loop. This is the multiple-imputation question: has the sequence of imputed datasets stabilised across chained iterations, or is it still drifting? This is not about MCMC mixing within a single fit — it is about whether the outer loop over variables has reached a stationary distribution over imputed values. The right tool for this question is
assess_convergence(), which tracks summary statistics of the imputed cells across therunschained iterations. Failures here are fixed by increasingrunsor, occasionally, by revisiting the specification of the imputation models.
bace() assesses only the second kind of convergence automatically (and retries with more runs if it fails). The first kind is your responsibility to check before trusting pooled inference — see Section 6.6.
9.3 Assessing chained-equations convergence
We need to check whether the number of runs we’ve done has led to convergence of the imputation process. We can do this by looking at the stability of the imputed values across the chained iterations. For example:
Code
# Assess whether number of runs was enough by looking at the stability of imputed points across runs. Adjust runs if needed
converge <- assess_convergence(bace_impute2, method = "summary")
plot(converge)
Code
# Now that we know that we can run BACE longer
# Fit BACE with a single formula
bace_impute2_2 <- bace_imp(fixformula = "y ~ x1 + x2",
ran_phylo_form = "~ 1 | Species",
phylo = tree,
data = data,
runs = 15, nitt = 100000, burnin = 10000, thin = 10) Run 1, Imputed variable: y
Run 1, Imputed variable: x1
Run 1, Imputed variable: x2
Run 2, Imputed variable: y
Run 2, Imputed variable: x1
Run 2, Imputed variable: x2
Run 3, Imputed variable: y
Run 3, Imputed variable: x1
Run 3, Imputed variable: x2
Run 4, Imputed variable: y
Run 4, Imputed variable: x1
Run 4, Imputed variable: x2
Run 5, Imputed variable: y
Run 5, Imputed variable: x1
Run 5, Imputed variable: x2
Run 6, Imputed variable: y
Run 6, Imputed variable: x1
Run 6, Imputed variable: x2
Run 7, Imputed variable: y
Run 7, Imputed variable: x1
Run 7, Imputed variable: x2
Run 8, Imputed variable: y
Run 8, Imputed variable: x1
Run 8, Imputed variable: x2
Run 9, Imputed variable: y
Run 9, Imputed variable: x1
Run 9, Imputed variable: x2
Run 10, Imputed variable: y
Run 10, Imputed variable: x1
Run 10, Imputed variable: x2
Run 11, Imputed variable: y
Run 11, Imputed variable: x1
Run 11, Imputed variable: x2
Run 12, Imputed variable: y
Run 12, Imputed variable: x1
Run 12, Imputed variable: x2
Run 13, Imputed variable: y
Run 13, Imputed variable: x1
Run 13, Imputed variable: x2
Run 14, Imputed variable: y
Run 14, Imputed variable: x1
Run 14, Imputed variable: x2
Run 15, Imputed variable: y
Run 15, Imputed variable: x1
Run 15, Imputed variable: x2
========================================
BACE MCMC Diagnostics Summary
========================================
Number of models: 3
Total parameters: 15
Fixed effects: 9
Variance components: 6
Fixed units (threshold/categorical): 1
Mean Effective Sample Size:
Fixed effects: 4046.7
Variance components: 2774.6
Convergence Assessment:
Acceptable : 2 parameters (13.3%)
Check : 4 parameters (26.7%)
Fixed : 1 parameters (6.7%)
Good : 8 parameters (53.3%)
Overall: CHECK - Many parameters show convergence issues
Consider increasing nitt, adjusting priors, or checking model specification
Convergence by Response Variable:
----------------------------------------
[~] y : 3/ 5 parameters good
[OK] x1 : 5/ 5 parameters good
[CH] x2 : 0/ 4 parameters good (1 fixed)
Check convergence: (Intercept), x1, y, Species
========================================Code
# Check again for convergence
converge <- assess_convergence(bace_impute2_2, method = "summary")
plot(converge)
9.4 Worked example: bace() with species random effects
Now that you are aware of the different components of the BACE workflow, it is useful to see them all combined in a second worked example — this time using the replicated-species dataset from the previous section so that both phylogenetic and non-phylogenetic species effects are estimated, and letting bace() handle the orchestration end-to-end.
Code
# Run the full bace workflow (n_final = 10 again only for vignette speed):
bace_final <- bace(fixformula = "y ~ x1 + x2",
ran_phylo_form = "~ 1 | Species",
phylo = tree2,
data = data2,
species = TRUE,
runs = 15, nitt = 100000, burnin = 10000, thin = 10, n_final = 10, plot = TRUE, skip_conv = TRUE, sample_size = 1000)
=======================================================
BACE: Bayesian Augmentation by Chained Equations
=======================================================
Step 1: Initial imputation for convergence assessment
-------------------------------------------------------Run 1, Imputed variable: y
Run 1, Imputed variable: x1
Run 1, Imputed variable: x2
Run 2, Imputed variable: y
Run 2, Imputed variable: x1
Run 2, Imputed variable: x2
Run 3, Imputed variable: y
Run 3, Imputed variable: x1
Run 3, Imputed variable: x2
Run 4, Imputed variable: y
Run 4, Imputed variable: x1
Run 4, Imputed variable: x2
Run 5, Imputed variable: y
Run 5, Imputed variable: x1
Run 5, Imputed variable: x2
Run 6, Imputed variable: y
Run 6, Imputed variable: x1
Run 6, Imputed variable: x2
Run 7, Imputed variable: y
Run 7, Imputed variable: x1
Run 7, Imputed variable: x2
Run 8, Imputed variable: y
Run 8, Imputed variable: x1
Run 8, Imputed variable: x2
Run 9, Imputed variable: y
Run 9, Imputed variable: x1
Run 9, Imputed variable: x2
Run 10, Imputed variable: y
Run 10, Imputed variable: x1
Run 10, Imputed variable: x2
Run 11, Imputed variable: y
Run 11, Imputed variable: x1
Run 11, Imputed variable: x2
Run 12, Imputed variable: y
Run 12, Imputed variable: x1
Run 12, Imputed variable: x2
Run 13, Imputed variable: y
Run 13, Imputed variable: x1
Run 13, Imputed variable: x2
Run 14, Imputed variable: y
Run 14, Imputed variable: x1
Run 14, Imputed variable: x2
Run 15, Imputed variable: y
Run 15, Imputed variable: x1
Run 15, Imputed variable: x2
========================================
BACE MCMC Diagnostics Summary
========================================
Number of models: 3
Total parameters: 18
Fixed effects: 9
Variance components: 9
Fixed units (threshold/categorical): 1
Mean Effective Sample Size:
Fixed effects: 5320.4
Variance components: 1988.9
Convergence Assessment:
Acceptable : 6 parameters (33.3%)
Check : 3 parameters (16.7%)
Fixed : 1 parameters (5.6%)
Good : 8 parameters (44.4%)
Overall: ACCEPTABLE - Most parameters converged adequately
Convergence by Response Variable:
----------------------------------------
[CH] y : 2/ 6 parameters good
Check convergence: x2.L, Species2
[~] x1 : 4/ 6 parameters good
Check convergence: Species2
[~] x2 : 2/ 5 parameters good (1 fixed)
========================================
Step 2: Convergence assessment
-------------------------------------------------------
Skipping convergence retry (skip_conv = TRUE)
Proceeding with final imputation...
-------------------------------------------------------
=== Running 10 final imputation iterations ===
Run 1/10 (sweep 1/3) - Imputed variable: y
Run 1/10 (sweep 1/3) - Imputed variable: x1
Run 1/10 (sweep 1/3) - Imputed variable: x2
Run 1/10 (sweep 2/3) - Imputed variable: y
Run 1/10 (sweep 2/3) - Imputed variable: x1
Run 1/10 (sweep 2/3) - Imputed variable: x2
Run 1/10 (sweep 3/3) - Imputed variable: y
Run 1/10 (sweep 3/3) - Imputed variable: x1
Run 1/10 (sweep 3/3) - Imputed variable: x2
Run 2/10 (sweep 1/3) - Imputed variable: y
Run 2/10 (sweep 1/3) - Imputed variable: x1
Run 2/10 (sweep 1/3) - Imputed variable: x2
Run 2/10 (sweep 2/3) - Imputed variable: y
Run 2/10 (sweep 2/3) - Imputed variable: x1
Run 2/10 (sweep 2/3) - Imputed variable: x2
Run 2/10 (sweep 3/3) - Imputed variable: y
Run 2/10 (sweep 3/3) - Imputed variable: x1
Run 2/10 (sweep 3/3) - Imputed variable: x2
Run 3/10 (sweep 1/3) - Imputed variable: y
Run 3/10 (sweep 1/3) - Imputed variable: x1
Run 3/10 (sweep 1/3) - Imputed variable: x2
Run 3/10 (sweep 2/3) - Imputed variable: y
Run 3/10 (sweep 2/3) - Imputed variable: x1
Run 3/10 (sweep 2/3) - Imputed variable: x2
Run 3/10 (sweep 3/3) - Imputed variable: y
Run 3/10 (sweep 3/3) - Imputed variable: x1
Run 3/10 (sweep 3/3) - Imputed variable: x2
Run 4/10 (sweep 1/3) - Imputed variable: y
Run 4/10 (sweep 1/3) - Imputed variable: x1
Run 4/10 (sweep 1/3) - Imputed variable: x2
Run 4/10 (sweep 2/3) - Imputed variable: y
Run 4/10 (sweep 2/3) - Imputed variable: x1
Run 4/10 (sweep 2/3) - Imputed variable: x2
Run 4/10 (sweep 3/3) - Imputed variable: y
Run 4/10 (sweep 3/3) - Imputed variable: x1
Run 4/10 (sweep 3/3) - Imputed variable: x2
Run 5/10 (sweep 1/3) - Imputed variable: y
Run 5/10 (sweep 1/3) - Imputed variable: x1
Run 5/10 (sweep 1/3) - Imputed variable: x2
Run 5/10 (sweep 2/3) - Imputed variable: y
Run 5/10 (sweep 2/3) - Imputed variable: x1
Run 5/10 (sweep 2/3) - Imputed variable: x2
Run 5/10 (sweep 3/3) - Imputed variable: y
Run 5/10 (sweep 3/3) - Imputed variable: x1
Run 5/10 (sweep 3/3) - Imputed variable: x2
Run 6/10 (sweep 1/3) - Imputed variable: y
Run 6/10 (sweep 1/3) - Imputed variable: x1
Run 6/10 (sweep 1/3) - Imputed variable: x2
Run 6/10 (sweep 2/3) - Imputed variable: y
Run 6/10 (sweep 2/3) - Imputed variable: x1
Run 6/10 (sweep 2/3) - Imputed variable: x2
Run 6/10 (sweep 3/3) - Imputed variable: y
Run 6/10 (sweep 3/3) - Imputed variable: x1
Run 6/10 (sweep 3/3) - Imputed variable: x2
Run 7/10 (sweep 1/3) - Imputed variable: y
Run 7/10 (sweep 1/3) - Imputed variable: x1
Run 7/10 (sweep 1/3) - Imputed variable: x2
Run 7/10 (sweep 2/3) - Imputed variable: y
Run 7/10 (sweep 2/3) - Imputed variable: x1
Run 7/10 (sweep 2/3) - Imputed variable: x2
Run 7/10 (sweep 3/3) - Imputed variable: y
Run 7/10 (sweep 3/3) - Imputed variable: x1
Run 7/10 (sweep 3/3) - Imputed variable: x2
Run 8/10 (sweep 1/3) - Imputed variable: y
Run 8/10 (sweep 1/3) - Imputed variable: x1
Run 8/10 (sweep 1/3) - Imputed variable: x2
Run 8/10 (sweep 2/3) - Imputed variable: y
Run 8/10 (sweep 2/3) - Imputed variable: x1
Run 8/10 (sweep 2/3) - Imputed variable: x2
Run 8/10 (sweep 3/3) - Imputed variable: y
Run 8/10 (sweep 3/3) - Imputed variable: x1
Run 8/10 (sweep 3/3) - Imputed variable: x2
Run 9/10 (sweep 1/3) - Imputed variable: y
Run 9/10 (sweep 1/3) - Imputed variable: x1
Run 9/10 (sweep 1/3) - Imputed variable: x2
Run 9/10 (sweep 2/3) - Imputed variable: y
Run 9/10 (sweep 2/3) - Imputed variable: x1
Run 9/10 (sweep 2/3) - Imputed variable: x2
Run 9/10 (sweep 3/3) - Imputed variable: y
Run 9/10 (sweep 3/3) - Imputed variable: x1
Run 9/10 (sweep 3/3) - Imputed variable: x2
Run 10/10 (sweep 1/3) - Imputed variable: y
Run 10/10 (sweep 1/3) - Imputed variable: x1
Run 10/10 (sweep 1/3) - Imputed variable: x2
Run 10/10 (sweep 2/3) - Imputed variable: y
Run 10/10 (sweep 2/3) - Imputed variable: x1
Run 10/10 (sweep 2/3) - Imputed variable: x2
Run 10/10 (sweep 3/3) - Imputed variable: y
Run 10/10 (sweep 3/3) - Imputed variable: x1
Run 10/10 (sweep 3/3) - Imputed variable: x2
=== Final imputation complete ===
Saved 10 imputed datasets with corresponding models
--- Final-phase MCMC diagnostics (10 runs; worst case across runs) ---
Flagged parameters (min ESS < 200 or > 25% of runs flagged):
y units [Variance] min ESS 706 (run 7), 30% runs flagged
Consider longer chains (nitt/nitt_cat_mult) or a per-variable prior
(see ?bace section 'Priors') for the variables above.
Full table: $diagnostics$by_param
Pooling posteriors with sampling: 1000 samples per imputation
=======================================================
BACE Analysis Complete!
=======================================================
Convergence achieved: FALSE
Number of attempts: 1
Final imputations: 10
Access results via:
- get_pooled_model(result, "y"): pooled posterior model for a variable ($pooled_models)
- get_imputed_data(result): list of 10 imputed datasets ($imputed_datasets)
- $convergence: Convergence diagnostics
- $final_results: Full final imputation results (per-run models and datasets)Code
# Extract the pooled model for y using the accessor function
y_model <- get_pooled_model(bace_final, variable = "y")
summary(y_model)
+----------------------------------------------------------------+
| BACE Pooled MCMCglmm summary (stacked per-imputation draws) |
+----------------------------------------------------------------+
Stacked from 10 imputations
Total posterior draws: 10000
(= 1000 sampled per imputation x 10 imputations)
Original samples per imputation: 9000
Posterior summaries below are computed from stacked draws that
approximate p(theta | Y_obs) (Monte Carlo integration over imputed
values). Validity requires each per-imputation chain to have
converged on its own; check a subset via
plot(final$all_models[[1]]$<var>) and coda::effectiveSize().
----------------------------------------------------------------
Iterations = 1:10000
Thinning interval = 1
Sample size = 10000
DIC: 270.209
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 0.1309 6.086e-11 0.5863 9993
~Species2
post.mean l-95% CI u-95% CI eff.samp
Species2 0.2118 3.001e-08 0.641 10000
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.3409 0.02074 0.7376 1854
Location effects: y ~ x1 + x2
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 0.1425 -0.4170 0.7070 892.62 0.544
x1 -0.1272 -0.4618 0.2291 697.13 0.444
x2.L 0.3582 -0.1938 0.8969 34.41 0.205Code
plot(y_model) 
Code
# We can now extract posteriors in the same way we would use a normal MCMCglmm model,
# but now we have draws that integrate over the imputation distribution. For example,
# we can look at the marginal posterior of the x1 coefficient:
hist(y_model$Sol[, "x1"], main = "Posterior distribution for x1", xlab = "x1 coefficient")
10 Validating the imputations
10.1 On your own data: bace_validate()
Before trusting imputed values, check how well BACE can recover values you already know. bace_validate() does the whole check in one call: it hides a share of the observed cells, runs the full pipeline on the masked data, and scores the imputations against the hidden truth — with a no-skill baseline printed next to every score, so you can see whether BACE is genuinely using the phylogenetic and trait structure.
Code
# The natural pattern: run bace(), then validate the finished run.
# bace_validate(fit) reuses the run's data, tree and settings — nothing
# is re-typed. (It re-runs the pipeline on masked data under the hood:
# a finished run cannot be scored directly, because its imputations were
# made with every observed cell visible.)
fit <- bace("y ~ x1 + x2", "~ 1 | Species", phylo = tree, data = mydata)
v <- bace_validate(
fit,
prop = 0.1, # hide 10% of each variable's observed cells
reps = 2, # repeat with a fresh mask for honest uncertainty
seed = 1
)
vThe printed report is one short block per variable. For continuous and count variables you get the zRMSE (0 = hidden values recovered exactly; 1 = no better than filling in the variable’s mean — the mean-filling baseline is printed beside it) and the 95% interval coverage (how often the interval from get_imputation_intervals() contained the hidden value; it should sit near 0.95, and for low-count variables the report prints the slightly lower level that whole-numbered data make attainable). For factors you get accuracy against the guess-the-commonest-class baseline, balanced accuracy (rare classes count equally), and the Brier score for the class probabilities (0 = perfect).
Two honest limitations, stated in ?bace_validate as well: the check costs one full BACE run per repetition, and it measures prediction where data were observed — if your data are missing not-at-random (say, small or rare species were never measured), no masking check can see that.
10.2 On simulated data with known truth
For methodological work, simulation-based benchmarking is still the gold standard: sim_bace() returns the complete dataset alongside the masked one, so the truth is known by construction rather than held out. bace_validate() automates the mask–impute–score loop there too:
Code
library(BACE)
set.seed(2026)
# A complete dataset with strong phylogenetic signal; no NAs — the masking
# is bace_validate()'s job
sim <- sim_bace(
response_type = "gaussian",
predictor_types = c("gaussian", "gaussian"),
var_names = c("y", "x1", "x2"),
phylo_signal = c(0.9, 0.9, 0.9),
n_cases = 120,
n_species = 60,
missingness = c(0, 0, 0)
)
v <- bace_validate(
"y ~ x1 + x2",
ran_phylo_form = "~ 1 | Species",
phylo = sim$tree,
data = sim$complete_data,
prop = 0.2,
runs = 4, n_final = 50,
nitt = 2500, thin = 10, burnin = 500,
skip_conv = TRUE, verbose = FALSE, seed = 2026
)
vOn a high-signal trait like this, expect a zRMSE around 0.4–0.5 (less than half the mean-filling baseline of 1) with interval coverage near 0.95. The same call handles categorical responses (swap response_type = "multinomial3"), reporting accuracy and Brier score instead.
Across the full 20-replicate sweep at default n_final = 50 settings (see dev/sim_accuracy_multirep.R), BACE is significantly better than the no-skill baseline on every response type (paired Wilcoxon p < 0.05). (That sweep scored intervals with raw empirical quantiles, which undercover by design — about 85–92% there; use get_imputation_intervals() for the calibrated t-based intervals.)
| response type | BACE metric | baseline | n_reps | paired p |
|---|---|---|---|---|
| gaussian | NRMSE 0.45 | 1.00 | 20 | <1e-4 |
| poisson | NRMSE 0.81 | 0.93 | 20 | 0.046 |
| binary | balanced acc 0.69 | 0.50 | 20 | <1e-4 |
| categorical K=3 | balanced acc 0.43 | 0.33 | 20 | <2e-3 |
| ordinal K=4 | mae_level 0.22 | 0.63 | 20 | <1e-4 |
The dev/sim_accuracy_results/multirep_v1_n20/ directory of the package source contains the per-replicate CSVs and aggregated summary used to build this table.
11 Worked example: heat tolerance of amphibians
In this example, we use BACE to impute the heat tolerance of the world’s amphibians. Pottier et al. (2025) compiled 2,661 estimates of the critical thermal maximum (CTmax), a proxy for heat tolerance limits, across 524 species of amphibians. Below, we show how to use BACE to impute the CTmax of 5,203 species, and how to cross-validate imputed values against experimental data.
11.1 The data
The dataset contains two types of rows. Experimental rows are CTmax estimates from published heat-tolerance assays. These assays differ in their methods, such as acclimation temperature and duration, heating rate, and the endpoint used to define CTmax. Because these methods affect CTmax, raw estimates are not directly comparable across species.
Standardised rows make estimates comparable across species. Each standardised row describes an assay performed under identical conditions. Adults are acclimated for 10 days and heated at 1 °C per minute in water, and CTmax is recorded at the onset of spasms. These conditions were the median or most common values in the experimental data (Pottier et al. 2025). Each species has three standardised rows, which differ only in acclimation temperature. These temperatures correspond to the 5th, 50th and 95th percentiles of the operative body temperatures that each species experiences across its range. CTmax is unknown in all 15,609 standardised rows, and BACE imputes it using the experimental data and the phylogeny.
Code
# Data and phylogeny are bundled in vignettes/amphibian_ctmax/data/
ctmax_data <- readRDS("amphibian_ctmax/data/ctmax_data.rds")
ctmax_tree <- readRDS("amphibian_ctmax/data/ctmax_tree.rds")
table(ctmax_data$row_type)
experimental standardised
2661 15609 Code
# Experimental and standardised rows for one species, Cope's grey treefrog
amph_example <- ctmax_data[ctmax_data$species == "Dryophytes chrysoscelis",
c("row_type", "CTmax", "acclimation_temp", "ramping", "endpoint", "life_stage")]
print(amph_example, digits = 3, row.names = FALSE) row_type CTmax acclimation_temp ramping endpoint life_stage
experimental 38.0 10.0 1.000000 LRR adults
experimental 40.8 10.0 1.000000 OS adults
experimental 41.5 23.0 0.000694 other adults
standardised NA 23.3 1.000000 OS adults
standardised NA 20.4 1.000000 OS adults
standardised NA 27.1 1.000000 OS adultsOne of the covariates is the medium in which temperature was recorded during CTmax trials. Studies recorded either the body temperature of the animals, the water temperature, or the ambient temperature. The medium is unknown for 36 of the 2,661 experimental rows. For demonstration purposes, we assume that these assays recorded body or water temperature, as did 2,574 of the 2,625 assays with a known medium. We do so because imputing this covariate is computationally demanding: ambient temperature is so rare that its imputation model needs very long MCMC chains. However, BACE can impute this variable if needed (see Pottier et al. 2025).
Code
ctmax_data$medium[is.na(ctmax_data$medium)] <- "body_water"11.2 Checking phylogenetic signal
Before running BACE, we checked whether the variables with missing data show phylogenetic signal (see Pre-flight: check phylogenetic signal above).
Code
exp_rows <- ctmax_data$row_type == "experimental"
signal <- phylo_signal_summary(
data = ctmax_data[exp_rows, ],
tree = ape::keep.tip(ctmax_tree, unique(ctmax_data$species[exp_rows])),
species_col = "species",
variables = c("CTmax", "ln_sd_CTmax", "ln_acclimation_time", "ramping"),
species = TRUE,
nitt = 105000, burnin = 5000, thin = 50)
signal=============================================================
Phylogenetic signal summary
=============================================================
Random-effect structure: dual (phylo + non-phylo species)
Variables : 4
--- Table ---
variable type n_obs H2_mean H2_lo H2_hi R_species_mean lambda
CTmax gaussian 2661 0.856 0.815 0.894 0.0183 0.923
ln_sd_CTmax gaussian 1960 0.424 0.258 0.608 0.2196 0.555
ln_acclimation_time gaussian 2282 0.419 0.301 0.535 0.1359 0.663
ramping gaussian 2488 0.463 0.293 0.624 0.2998 0.457
K min_ess interpretation flag
0.432 1556 high
0.141 1939 moderate
0.251 1702 moderate
0.109 1371 moderate
Interpretation bands: H2 < 0.2 = low, 0.2-0.5 = moderate, >0.5 = high.
Flag meanings:
low_ess -- effective sample size below threshold (default 1000)
geweke_fail -- Geweke diagnostic |z| > 2 on any variance component
unreliable -- H2 cell should NOT be trusted (low_ess or geweke_fail)
low_n -- n_obs < 20; posterior is data-sparse
fit_error -- MCMCglmm fit failed (check warnings)
Note: phylogenetic signal is necessary but NOT sufficient for good
imputation. A high-H2 trait can still impute poorly under MNAR missingness;
a low-H2 trait can still impute well if coupled to an observed predictor.CTmax showed high phylogenetic signal (H² = 0.86, 95% CI: 0.81 to 0.89), and the three covariates showed moderate signal (H² = 0.42 to 0.46). The phylogeny therefore carries information for every imputation model.
11.3 Running BACE
We fitted one imputation model for each variable with missing data, following Pottier et al. (2025). We modelled CTmax as a function of acclimation temperature, acclimation duration, heating rate, test medium, endpoint, acclimation status (laboratory-acclimated or field-fresh), life stage and ecotype. Acclimation duration and heating rate are missing for 379 and 173 of the 2,661 experimental rows, respectively, so each has its own imputation model. We also included the standard deviation of CTmax as an auxiliary variable. It does not predict CTmax, but it helps impute the missing assay settings. Acclimation duration and the standard deviation of CTmax were log-transformed because both were right-skewed.
All models included phylogenetic and non-phylogenetic random intercepts for species (species = TRUE), because many species have several estimates (287 of the 524 tested species). Pottier et al. (2025) also fitted random slopes for acclimation temperature, which BACE does not yet support. In this example, all species therefore share the same effect of acclimation temperature on CTmax.
Code
ctmax_formulas <- list(
"ln_sd_CTmax ~ life_stage + acclimation_temp + endpoint + acclimated + ln_acclimation_time + medium + ramping + CTmax",
"ln_acclimation_time ~ life_stage + ln_sd_CTmax + CTmax",
"ramping ~ life_stage + ln_sd_CTmax + CTmax",
"CTmax ~ acclimation_temp + ln_acclimation_time + ramping + medium + endpoint + acclimated + life_stage + ecotype"
)
ctmax_vars <- c("species", "CTmax", "acclimation_temp", "ln_acclimation_time", "ramping",
"ln_sd_CTmax", "medium", "endpoint", "acclimated", "life_stage", "ecotype")
amph_bace <- bace(fixformula = ctmax_formulas,
ran_phylo_form = "~ 1 | species",
phylo = ctmax_tree,
data = ctmax_data[, ctmax_vars],
species = TRUE,
runs = 10, n_final = 50,
nitt = 15000, burnin = 3000, thin = 24)BACE chose the model for each variable from how it is stored (see the table of variable types above). All four imputed variables are stored as numbers, so BACE fitted Gaussian models to them. The factors are complete, so they enter the models as predictors only; their listed types would apply if they had missing values.
Code
unlist(amph_bace$initial_results$types) ln_sd_CTmax life_stage acclimation_temp endpoint
"gaussian" "threshold" "gaussian" "categorical"
acclimated ln_acclimation_time medium ramping
"threshold" "gaussian" "threshold" "gaussian"
CTmax ecotype
"gaussian" "categorical"
Run time
This analysis is computationally demanding. On a desktop computer, the chained phase took 4.5 hours and the 50 final imputations took 17.7 core-hours. Because bace() keeps every fitted model in memory (about 26 GB here), we ran the same steps in chunks with bace_imp() and bace_final_imp(). The scripts are in vignettes/amphibian_ctmax/.
11.4 Checking convergence
We checked both kinds of convergence (see Two kinds of convergence above). First, we used assess_convergence() to check whether the imputed values stabilised across the 10 chained iterations.
Code
converge <- assess_convergence(amph_bace$initial_results, method = "summary")
converge
plot(converge)
=== BACE Imputation Convergence Assessment ===
Overall Convergence: YES
--- Summary Statistics Method ---
Converged: TRUE
Variables converged: 4 / 4
Use plot() to visualize convergence diagnostics
The chained loop converged (Fig. 2). Second, we checked how well each MCMC chain mixed. With verbose = TRUE, BACE prints a diagnostics summary for the models of the last chained iteration, which is also stored in the output:
Code
amph_bace$initial_results$diagnostics
========================================
BACE MCMC Diagnostics Summary
========================================
Number of models: 4
Total parameters: 44
Fixed effects: 32
Variance components: 12
Mean Effective Sample Size:
Fixed effects: 337.2
Variance components: 307.3
Convergence Assessment:
Check : 16 parameters (36.4%)
Good : 28 parameters (63.6%)
Overall: CHECK - Many parameters show convergence issues
Consider increasing nitt, adjusting priors, or checking model specification
Convergence by Response Variable:
----------------------------------------
[CH] ln_sd_CTmax : 6/13 parameters good
Check convergence: life_stagelarvae, endpointOS, endpointother, acclimatedfield-fresh, mediumbody_water, species, species2
[OK] ln_acclimation_time: 7/ 7 parameters good
[OK] ramping : 7/ 7 parameters good
[CH] CTmax : 8/17 parameters good
Check convergence: mediumbody_water, life_stagelarvae, ecotypeArboreal, ecotypeFossorial, ecotypeGround-dwelling, ecotypeSemi-aquatic, ecotypeStream-dwelling, species, species2
========================================The summary flagged 16 of the 44 parameters, all in the models for the standard deviation of CTmax and for CTmax. BACE flags a parameter when its ESS is 200 or less, or when a Geweke test suggests that its chain has not settled (p ≤ 0.05). Here, 13 of the 16 flags came from low ESS, mainly for the variance components and for the effects of ecotype on CTmax. Close relatives often share an ecotype (3,943 of 5,089 species in genera with several species, 77%, share the most common ecotype of their genus). The effects of ecotype therefore probably overlap with the phylogenetic effect, which slows the sampler.
We can check how the chain mixed with trace plots. A chain that mixes well looks like a fuzzy caterpillar without trends, whereas a chain that mixes slowly wanders. For example, Fig. 3 shows the model for the standard deviation of CTmax, from one final imputation:
Code
first_fit <- amph_bace$final_results$all_models[[1]]
plot(first_fit$ln_sd_CTmax)
ln_sd_CTmax), from one final imputation of 500 draws. The last three rows of panels are the phylogenetic (species), non-phylogenetic (species2) and residual (units) variances.
Most parameters show a fuzzy caterpillar structure, with ESS between 148 and 500 out of 500 draws. The phylogenetic variance (species) is the exception. Its chain wanders slowly between 0.1 and 1.0 (ESS = 12), and the non-phylogenetic variance (species2) also mixes slowly (ESS = 72).
Slow mixing within a chain matters less for the imputed values, because each final imputation comes from an independent chain. Pooled across the 50 chains, the ESS of every parameter of the CTmax model exceeds 2,100 (see eff.samp in the summary below). In your own analyses, increase nitt and burnin when ESS is low.
11.5 The pooled CTmax model
As in the examples above, get_pooled_model() returns the model for a variable, with the posterior draws of the 50 final imputations stacked together.
Code
ctmax_model <- get_pooled_model(amph_bace, variable = "CTmax")Code
summary(ctmax_model)
+----------------------------------------------------------------+
| BACE Pooled MCMCglmm summary (stacked per-imputation draws) |
+----------------------------------------------------------------+
Stacked from 50 imputations
Total posterior draws: 25000
(= 500 samples/imputation x 50 imputations)
Posterior summaries below are computed from stacked draws that
approximate p(theta | Y_obs) (Monte Carlo integration over imputed
values). Validity requires each per-imputation chain to have
converged on its own; check a subset via
plot(final$all_models[[1]]$<var>) and coda::effectiveSize().
----------------------------------------------------------------
Iterations = 1:25000
Thinning interval = 1
Sample size = 25000
DIC:
G-structure: ~species
post.mean l-95% CI u-95% CI eff.samp
species 1.107 0.8328 1.392 2110
~species2
post.mean l-95% CI u-95% CI eff.samp
species2 0.06734 0.03114 0.1057 2125
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.1562 0.1471 0.1658 20882
Location effects: CTmax ~ acclimation_temp + ln_acclimation_time + ramping + medium + endpoint + acclimated + life_stage + ecotype
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) -0.462003 -1.440520 0.531399 14569 0.35368
acclimation_temp 0.173695 0.163098 0.185148 17191 < 4e-05 ***
ln_acclimation_time 0.017095 0.007725 0.026726 3142 0.00040 ***
ramping -0.005003 -0.018578 0.008008 4969 0.45848
mediumbody_water -0.340305 -0.556659 -0.112725 6137 0.00264 **
endpointOS 0.390059 0.333319 0.446376 8725 < 4e-05 ***
endpointother 0.416181 0.350222 0.480290 16066 < 4e-05 ***
acclimatedfield-fresh -0.253926 -0.422659 -0.081096 10660 0.00328 **
life_stagelarvae 0.414883 0.344985 0.480877 7528 < 4e-05 ***
ecotypeArboreal -0.087319 -0.445839 0.276960 6045 0.63400
ecotypeFossorial 0.236634 -0.174448 0.637990 6668 0.25912
ecotypeGround-dwelling -0.036822 -0.367131 0.286674 6471 0.82904
ecotypeSemi-aquatic 0.019590 -0.308708 0.341528 7774 0.90016
ecotypeStream-dwelling -0.223289 -0.605744 0.130172 3985 0.23648
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BACE fits Gaussian models to standardised (z-transformed) variables, so these coefficients are expressed in standard deviations. For example, the coefficient of acclimation temperature converts to °C of CTmax per °C of acclimation temperature as follows:
Code
acc_effect <- ctmax_model$Sol[, "acclimation_temp"] *
sd(ctmax_data$CTmax, na.rm = TRUE) / sd(ctmax_data$acclimation_temp)
round(c(mean = mean(acc_effect), quantile(acc_effect, c(0.025, 0.975))), 3) mean 2.5% 97.5%
0.134 0.126 0.143 On average, CTmax increased by 0.134 °C (95% credible interval: 0.126 to 0.143) per °C of acclimation temperature. The trace and density plots show the posterior draws of the fixed effects (Fig. 4) and variance components (Fig. 5).
Code
par(mfrow = c(7, 4))
plot(ctmax_model$Sol, auto.layout = FALSE)
Code
plot(ctmax_model$VCV)
11.6 Imputed CTmax
Each final imputation is a complete dataset, so each imputed CTmax value has 50 draws. We summarised these draws with a mean and a 95% credible interval.
Code
imp <- get_imputed_data(amph_bace, format = "list")
draws <- sapply(imp, function(d) d$CTmax) # rows x imputations
ctmax_mean <- rowMeans(draws)
ctmax_pi95 <- apply(draws, 1, quantile, c(0.025, 0.975))Code
amph_imp <- amph_full$imputed[amph_full$imputed$temp_level == "p50", ] # median temperature
amph_imp$data <- factor(ifelse(amph_imp$species %in% amph_exp$species, "experimental", "imputed"),
levels = names(amph_cols))
amph_imp$width <- amph_imp$upper95 - amph_imp$lower95
amph_w <- tapply(amph_imp$width, amph_imp$data, median)
amph_s <- tapply(amph_imp$CTmax_mean, amph_imp$data, sd)
amph_n <- table(amph_imp$data)At the median acclimation temperature, imputed CTmax averaged 37.4 °C across the 5,203 species (Fig. 6). Imputations were more uncertain for species without experimental data. The median 95% credible interval spanned 6.0 °C for the 524 species with experimental data, but 8.9 °C for the 4,679 species without. The imputed means of species without data were also less variable (s.d. = 1.9 vs 2.9 °C), because BACE predicts their CTmax from their relatives only.
Code
amph_lab <- c(experimental = sprintf("Species with experimental data (n = %s)", fmt0(amph_n[["experimental"]])),
imputed = sprintf("Species without experimental data (n = %s)", fmt0(amph_n[["imputed"]])))
p_a <- ggplot(amph_imp, aes(CTmax_mean, fill = data)) +
geom_density(alpha = 0.6, colour = NA) +
scale_fill_manual(values = amph_cols, labels = amph_lab, name = NULL) +
labs(x = "Imputed CTmax (°C)", y = "Density", tag = "(a)") +
amph_theme + theme(legend.position = "top", legend.direction = "vertical")
p_b <- ggplot(amph_imp, aes(data, width, colour = data)) +
geom_jitter(width = 0.2, height = 0, alpha = 0.15, size = 0.6) +
geom_boxplot(fill = NA, colour = "black", outlier.shape = NA, width = 0.4) +
scale_colour_manual(values = amph_cols) +
scale_x_discrete(labels = c(experimental = "With data", imputed = "Without data")) +
labs(x = NULL, y = "Width of 95% credible interval (°C)", tag = "(b)") +
amph_theme
p_a + p_b + plot_layout(widths = c(1.4, 1))
11.7 Cross-validation
We assessed the accuracy of the imputations by cross-validation, following Pottier et al. (2025). We generated five datasets in which the experimental data of 16 species (~5% of the dataset) tested under conditions comparable to the imputed estimates were assigned as NA. These were estimates from adults ramped at 1 °C per minute, with CTmax recorded as the onset of spasms. We removed all entries missing CTmax, besides the observations we aimed to cross-validate, since the validation was done on experimental data. In Pottier et al. (2025), the cross-validation was done on the larger dataset (including data-deficient species), but this should not impact the accuracy of the cross-validation procedure since data-deficient species do not inform the imputation of CTmax. We then ran BACE on each cross-validation dataset, and compared the imputed values to the experimental estimates that were artificially taken as missing values.
Code
cv_folds <- readRDS("amphibian_ctmax/data/cv_folds.rds")
cv_results <- lapply(cv_folds, function(fold) {
# Experimental rows, with the CTmax of the test species hidden
d <- ctmax_data[ctmax_data$row_type == "experimental", ]
hidden <- d$row_n %in% fold$rows_masked
observed <- d$CTmax[hidden]
d$CTmax[hidden] <- NA
tree <- ape::keep.tip(ctmax_tree, unique(d$species))
fit <- bace(fixformula = ctmax_formulas, ran_phylo_form = "~ 1 | species",
phylo = tree, data = d[, ctmax_vars], species = TRUE,
runs = 10, n_final = 50, nitt = 15000, burnin = 3000, thin = 24)
# Mean and 95% credible interval of the 50 imputations of each hidden estimate
draws <- sapply(get_imputed_data(fit, format = "list"), function(x) x$CTmax[hidden])
data.frame(row_n = d$row_n[hidden], observed = observed, imputed = rowMeans(draws),
lower = apply(draws, 1, quantile, 0.025),
upper = apply(draws, 1, quantile, 0.975))
})Code
# Results of the five folds, as returned by the code above
amph_cv_df <- do.call(rbind, lapply(amph_cv, function(f) data.frame(
fold = f$fold, row_n = f$row_n, species = f$species, observed = f$truth,
imputed = rowMeans(f$draws),
lower = apply(f$draws, 1, quantile, 0.025),
upper = apply(f$draws, 1, quantile, 0.975))))
# Three species were hidden in two folds; count each hidden estimate once
twice <- amph_cv_df$row_n[duplicated(amph_cv_df$row_n)]
amph_cv_df <- amph_cv_df[!(amph_cv_df$row_n %in% twice & amph_cv_df$fold == 5), ]
amph_cv_df$error <- amph_cv_df$imputed - amph_cv_df$observed
amph_cv_tab <- data.frame(
n = nrow(amph_cv_df),
r = cor(amph_cv_df$observed, amph_cv_df$imputed),
RMSE = sqrt(mean(amph_cv_df$error^2)),
mean_error = mean(amph_cv_df$error),
coverage95 = mean(amph_cv_df$observed >= amph_cv_df$lower & amph_cv_df$observed <= amph_cv_df$upper))| n | r | RMSE | Mean error | Coverage |
|---|---|---|---|---|
| 375 | 0.86 | 1.44 | -0.45 | 93% |
BACE imputed the 375 hidden estimates of 77 species with r = 0.86 and a root-mean-square error (RMSE) of 1.44 °C (Tab. 1; Fig. 7). The 95% credible intervals contained 347 of the 375 hidden estimates (93%), close to the nominal 95%.
On average, imputed values were 0.45 °C lower than experimental values (35.74 ± 2.54 vs 36.19 ± 2.67 °C). Here, imputations primarily underestimated the heat tolerance of the most tolerant species, but there was overall good congruence with the experimental dataset.
Because BACE does not fit random slopes, all species share the same effect of acclimation temperature on CTmax. Some deviations could therefore arise because we underestimated or overestimated the plasticity of some species. Fitting random slopes for acclimation temperature could help the cross-validation, because species differ in their acclimation capacity. BACE currently does not allow such models, but see Pottier et al. (2025) for further details.
Code
lims <- range(amph_cv_df[, c("observed", "lower", "upper")])
p_a <- ggplot(amph_cv_df, aes(observed, imputed)) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
geom_linerange(aes(ymin = lower, ymax = upper), colour = amph_cols[["imputed"]], alpha = 0.25) +
geom_point(colour = amph_cols[["imputed"]], alpha = 0.7, size = 1.6) +
annotate("text", x = lims[1] + 0.5, y = lims[2] - 0.5, hjust = 0, vjust = 1, size = 3.3,
label = sprintf("n = %d\nr = %.2f\nRMSE = %.2f °C",
amph_cv_tab$n, amph_cv_tab$r, amph_cv_tab$RMSE)) +
scale_x_continuous(breaks = seq(20, 50, 5)) +
scale_y_continuous(breaks = seq(20, 50, 5)) +
coord_equal(xlim = lims, ylim = lims) +
labs(x = "Experimental CTmax (°C)", y = "Imputed CTmax (°C)", tag = "(a)") +
amph_theme
amph_long <- rbind(data.frame(source = "experimental", CTmax = amph_cv_df$observed),
data.frame(source = "imputed", CTmax = amph_cv_df$imputed))
p_b <- ggplot(amph_long, aes(CTmax, fill = source)) +
geom_density(alpha = 0.6, colour = NA) +
scale_fill_manual(values = amph_cols, name = NULL,
labels = c(experimental = "Experimental", imputed = "Imputed")) +
labs(x = "CTmax (°C)", y = "Density", tag = "(b)") +
amph_theme + theme(legend.position = c(0.02, 0.98), legend.justification = c(0, 1))
p_a + p_b + plot_layout(widths = c(1, 1.2))
12 Worked example: social behaviour of birds
The examples above impute continuous, count and categorical variables. Here we impute an ordinal variable, the social behaviour of birds, across the 10,563 species of a global bird phylogeny. Ordinal variables have ordered categories, such as low < medium < high. BACE models them on a latent (liability) scale with a threshold model: a species’ category depends on where its latent value falls relative to estimated cutpoints (Hadfield and Nakagawa 2010). Social behaviour is missing for 40% of the species, so this is a realistic imputation problem.
The example has two aims: to show how to run BACE on such data, and to show how to judge how far the imputations can be trusted. We (1) impute the social behaviour of the species without records and use the imputed datasets in a downstream analysis, (2) check BACE against known values that we hid from it, including whole families, and (3) impute an ordinal and a continuous variable together.
The checks in (2) set the expectations for the rest of the example. BACE’s imputations of social behaviour were clearly better than those of a model without phylogeny (less so for whole families), but only slightly better, and not clearly so, than a simple rule that copies the categories of each species’ closest relatives; they often missed the rarer categories and almost never identified solitary species in clades without described relatives. The example therefore illustrates the workflow, how to check the imputations, and how to carry their uncertainty into downstream analyses; it is not evidence that BACE predicts this trait better than simpler phylogenetic approaches.
12.1 The data
We took social behaviour from BIRDBASE (Şekercioğlu et al. 2025), body mass, the hand-wing index and habitat from AVONET (Tobias et al. 2022), using its eBird taxonomy, and the phylogeny from the Clements/eBird tree of McTavish et al. (2025) (version 1.6). BIRDBASE records six kinds of social behaviour for each species as separate columns: colonial; social (flocks, including mixed-species or seasonal flocks); pairs and family groups; singly and pairs; solitary; and lekking. A 1 means that the species has been described as behaving this way, and a 0 that it has not been described as doing so, not that it never does. A species can therefore have several 1s, and a species without any 1 has no description of its social behaviour.
We coded social behaviour as the largest grouping in which the species has been described: solitary < singly and pairs < pairs and family groups < flocks or colonies. These categories reflect increasing group size, so we treat them as ordinal. Colonial breeding aggregates many birds, so it joins the top category. Leks are display aggregations of males rather than social groups, so species described only as lekking (56) are treated as missing. This coding is an analytic choice: BIRDBASE mixes breeding and non-breeding behaviour, and 407 species are described both as solitary and as flocking (for example, solitary when breeding and flocking in winter); here they count as flocking. The hand-wing index measures wing pointedness, which is associated with flight efficiency and dispersal (Sheard et al. 2020); AVONET estimated it from relatives for 344 species, and we used these values as given (in the mixed-variable example below, only measured values were hidden). The bundled data in vignettes/bird_social/data/ are derived from these three sources, all licensed CC BY 4.0, and are redistributed under the same licence; please cite Şekercioğlu et al. (2025), Tobias et al. (2022) and McTavish et al. (2025) when reusing them.
Of the species with a description, 513 are solitary, 1,451 live singly and in pairs, 1,625 in pairs and family groups, and 2,709 in flocks or colonies. The categories are unbalanced: flocks or colonies account for 43% of the described species and solitary species for 8%, which matters when we judge the imputations below. Social behaviour is missing for 4,265 species: 4,205 without a description in BIRDBASE, 56 described only as lekking, and 4 that we could not match to BIRDBASE.
Code
# Data and phylogeny are bundled in vignettes/bird_social/data/
social_data <- readRDS("bird_social/data/social_data.rds")
social_tree <- readRDS("bird_social/data/social_tree.rds")
str(social_data[, c("species", "social_max", "log_mass", "hwi")])'data.frame': 10563 obs. of 4 variables:
$ species : chr "Camarhynchus_pallidus" "Camarhynchus_psittacula" "Camarhynchus_parvulus" "Camarhynchus_pauper" ...
$ social_max: Ord.factor w/ 4 levels "solitary"<"singly_pairs"<..: 2 NA NA NA NA NA NA NA 2 NA ...
$ log_mass : num 3.14 2.93 2.58 2.77 2.89 ...
$ hwi : num 14.8 15.9 15.3 17.8 16.6 17.2 14.6 19 17.4 14.8 ...Code
table(social_data$social_max, useNA = "ifany")
solitary singly_pairs pairs_family flocks_colonies <NA>
513 1451 1625 2709 4265 social_max is stored as an ordered factor. This is what tells BACE to fit a threshold model rather than an unordered (multinomial) one. Body mass was log-transformed because it is right-skewed.
12.2 Who is missing?
Multiple imputation assumes that, given the information in the imputation model, whether a value is missing does not depend on the missing value itself (missing at random). We cannot test this assumption, because the missing values are unknown, but we can look at which species are missing. In BIRDBASE, missingness follows taxonomy more closely than ecology. In logistic regressions of whether a species was missing, family (249 levels) lowered the AIC from 14,133 to 11,042, whereas log body mass, the hand-wing index and habitat together lowered it only to 13,870 (deviance explained 25% against 2%). Some orders are mostly missing, such as owls (Strigiformes, 81% of species) and hawks and eagles (Accipitriformes, 69%), whereas marine species are rarely missing (5%, against 43% of forest species).
Most species without records nevertheless have relatives with records in their own family: only 64 of them belong to the 26 families without any record. For 1,110 of them, however, at most a quarter of the species in their family have a record, so many imputations rest on sparse information.
Missingness that follows the phylogeny is what BACE is designed to use, because the phylogeny is in every imputation model. Two caveats remain. First, where most of a clade is missing, the imputations rely on more distant relatives and are less informative; we test this below by hiding whole families. Second, missingness could still depend on social behaviour itself (missing not at random), for example if solitary or inconspicuous species are less often described. The data cannot rule this out, and no imputation method can correct it without further assumptions.
12.3 Three analyses
We ran BACE three times on these data. Their results are kept in separate files and are reported separately below.
| Analysis | Social behaviour hidden from BACE | Hand-wing index hidden | Purpose | Results |
|---|---|---|---|---|
| Full | none (all 6,298 known values used) | none | the imputations to use in downstream analyses | full_results.rds |
| Validation | 1,261 known values | none | score the imputations of the hidden species | validation_results.rds |
| Mixed validation | the same 1,261 known values | 2,113 known values | score joint imputation of two variables | mixed_results.rds |
The validation analyses were run first, with the model fixed in advance, and the hidden values were scored once. Their scores are a check of this model, not a basis for choosing among models. Some summaries of the scores (such as the mean probability given to the observed category) were added afterwards to describe where the model does well or badly; they did not change the model.
12.5 How far can the imputations be trusted?
The imputations of species without records cannot be checked directly. Before fitting any model, we therefore set aside 1,261 species with known social behaviour, hid their values from BACE, imputed them together with the species without records (the validation analysis), and compared their imputations with the known values once, with the model fixed. The hidden species formed two sets, which test different things:
- Random: 630 species drawn at random among species with records. Most have close relatives with known behaviour in the data, so this set measures how well BACE fills gaps scattered across a well-sampled phylogeny.
- Whole families: all 631 species of 38 families, drawn at random among families with at most 100 species with records. Their closest relatives with known behaviour are in other families, so this set measures how well BACE extrapolates across the phylogeny.
The real gaps lie between these two cases: most species without records have some relatives with records in their family, but often few (see Who is missing?). Both sets were drawn at random from species with records, so neither represents the species without records if missingness depends on social behaviour itself.
We compared BACE with three simple baselines fitted to the species with known behaviour: the observed category frequencies (whose most likely category is “flocks or colonies” for every species); an ordinal probit regression on log body mass and the hand-wing index without phylogeny (MASS::polr()); and the category frequencies among each species’ closest relatives with known behaviour in the phylogeny. The last baseline uses the phylogeny in the simplest possible way, so it shows what BACE adds beyond “copy the relatives”. It needs two choices:
- How many relatives. We used k = 5. This value was carried over from an earlier example developed for this tutorial; it is not based on a published method, and it was not chosen from these results, but it was not pre-registered either. We report other values below.
- Ties. In an ultrametric tree, all members of a sister clade are equally distant from a species, so several relatives are often tied at the fifth-closest distance. Relatives strictly closer than that distance get weight 1, and those tied at it share the remaining weight equally, so the weights sum to five; this is the average over all ways of breaking the ties. More than five species can therefore contribute (median 13, maximum 3,279 for the hidden species, when the closest relatives with known behaviour form a whole clade). Distances that differed by less than 1e-08 million years were treated as tied: in this tree, differences due to rounding are below 10-12 and real differences above 10-6. The weighted counts of each category plus 0.5 give the probabilities.
The MCMC of the validation analysis was checked as for the full analysis: across its 20 final imputations, all \(\hat{R}\) ≤ 1.001, each chain had at least 383 effective draws for every parameter, and split \(\hat{R}\) ≤ 1.02. Geweke tests flagged more chains than expected for the second cutpoint (25%, 5 of 20 chains; binomial test against 5%, p = 0.003), and somewhat more for the phylogenetic variance and the third cutpoint (15% each). The flagged chains drifted in both directions (for the second cutpoint, 2 with positive and 3 with negative z), and the placebo test, away from the start of the chains, also flagged 10%–15% of the chains for these parameters. The excess therefore does not appear to be specific to the start of the chains; it may reflect the slow mixing of these parameters, to which the Geweke test is sensitive, but rates from 20 chains are imprecise and a short unfinished burn-in in some chains cannot be excluded. Its effect on the estimates was negligible: dropping the first 10% of every chain changed the pooled means by at most 0.006 posterior standard deviations. A longer burn-in and longer chains would be a cheap safeguard in a new analysis.
Because the imputations are probabilities over ordered categories, we scored them with the ranked probability score (RPS), which compares predicted and observed cumulative probabilities and penalises predictions more the further they are from the observed category (0 is perfect; we divide by the number of categories minus one, so the score ranges from 0 to 1). The RPS rewards predictions that are right on average, so with unbalanced categories it is dominated by the common ones. We therefore also report how often the most likely category was correct within each observed category (recall) and their mean (balanced accuracy), which weights the four categories equally. Overall accuracy is shown for completeness, but with unbalanced categories it mostly reflects the most common category: always predicting “flocks or colonies” would be right for 58% of the whole-family species.
| Hidden set | Method | Species | RPS | Balanced accuracy | Recall: solitary | Recall: singly/pairs | Recall: pairs/family | Recall: flocks/colonies | Accuracy | Mean max. probability |
|---|---|---|---|---|---|---|---|---|---|---|
| Random | BACE, probabilities averaged over posterior draws | 630 | 0.115 | 0.545 | 0.293 | 0.638 | 0.436 | 0.813 | 0.621 | 0.565 |
| Random | BACE, category frequencies across the 20 imputed datasets | 630 | 0.119 | 0.517 | 0.362 | 0.611 | 0.320 | 0.777 | 0.575 | 0.592 |
| Random | Observed frequencies | 630 | 0.181 | 0.250 | 0.000 | 0.000 | 0.000 | 1.000 | 0.398 | 0.415 |
| Random | Ordinal regression (no phylogeny) | 630 | 0.179 | 0.253 | 0.000 | 0.013 | 0.000 | 1.000 | 0.402 | 0.411 |
| Random | Closest relatives (k = 5) | 630 | 0.119 | 0.582 | 0.431 | 0.597 | 0.564 | 0.737 | 0.629 | 0.582 |
| Whole families | BACE, probabilities averaged over posterior draws | 631 | 0.137 | 0.355 | 0.018 | 0.290 | 0.198 | 0.915 | 0.612 | 0.529 |
| Whole families | BACE, category frequencies across the 20 imputed datasets | 631 | 0.137 | 0.364 | 0.089 | 0.333 | 0.233 | 0.801 | 0.564 | 0.558 |
| Whole families | Observed frequencies | 631 | 0.180 | 0.250 | 0.000 | 0.000 | 0.000 | 1.000 | 0.580 | 0.415 |
| Whole families | Ordinal regression (no phylogeny) | 631 | 0.159 | 0.250 | 0.000 | 0.000 | 0.000 | 1.000 | 0.580 | 0.457 |
| Whole families | Closest relatives (k = 5) | 631 | 0.142 | 0.363 | 0.000 | 0.312 | 0.414 | 0.727 | 0.544 | 0.556 |
| Hidden set | Observed category | Species | BACE | Closest relatives (k = 5) | Observed frequencies |
|---|---|---|---|---|---|
| Random | Solitary | 58 | 0.31 | 0.34 | 0.08 |
| Random | Singly and pairs | 149 | 0.39 | 0.43 | 0.24 |
| Random | Pairs and family groups | 172 | 0.34 | 0.43 | 0.27 |
| Random | Flocks or colonies | 251 | 0.64 | 0.55 | 0.42 |
| Whole families | Solitary | 56 | 0.14 | 0.12 | 0.08 |
| Whole families | Singly and pairs | 93 | 0.27 | 0.26 | 0.24 |
| Whole families | Pairs and family groups | 116 | 0.29 | 0.35 | 0.27 |
| Whole families | Flocks or colonies | 366 | 0.61 | 0.57 | 0.42 |
Compared with the baselines without phylogeny. As probabilistic predictions, BACE’s imputations had the lowest RPS in both sets (Tab. 6), but only slightly lower than the closest-relatives rule (below). The ordinal regression without phylogeny did hardly better than the observed frequencies in the random set (RPS 0.179 against 0.181), and BACE improved on it clearly there (RPS 0.115 against 0.179) and less in the whole families (0.137 against 0.159). Because related species are not independent, we computed intervals for the differences in RPS by resampling whole families (with replacement) rather than species: the 95% intervals for the improvement over the ordinal regression were 0.045 to 0.088 (random) and 0.0004 to 0.042 (whole families, where the lower limit nearly touches zero). Averaged over all hidden species, the imputed probabilities matched the observed frequencies of the categories, although not within every range of probabilities (see below).
Compared with the closest-relatives rule. BACE’s advantage over the closest-relatives rule was small (RPS lower by 0.005 and 0.005) and not distinguishable from zero (95% intervals -0.003 to 0.012 and -0.018 to 0.021). Its balanced accuracy was lower than that of the relatives rule in both sets (0.545 against 0.582 in the random set, 0.355 against 0.363 in the whole families). Tab. 7 shows where the difference in RPS comes from. In the random set, BACE gave more probability than the relatives rule to the observed category only for species in flocks or colonies (0.64 against 0.55), and less for the three smaller categories. In other words, BACE’s small gain in RPS comes from the most common category, at the expense of the rarer ones. In the whole families, BACE gave somewhat more probability to solitary species than the relatives rule did (0.14 against 0.12) and to species living singly and in pairs (0.27 against 0.26), but rarely enough to make “solitary” the most likely category (see below).
Other numbers of relatives. Because k = 5 is a convention rather than an estimate, we repeated the comparison with k = 1, 3, 10 and 20, with ties weighted in the same way (Tab. 8). The fixed pseudo-count of 0.5 per category smooths more strongly the fewer relatives there are: with k = 1, no category can receive a probability above 0.5. We therefore also scaled the pseudo-count with k (0.5 × k/5 per category, which leaves the main baseline at k = 5 unchanged). The smoothing does not change the most likely category, so the balanced accuracy is the same under both rules: it was higher for the relatives rule than for BACE (0.545 and 0.355) for every k in the whole families, and in the random set for k ≤ 10. For the RPS:
- With a single relative (k = 1), BACE was clearly better in both sets with the fixed pseudo-count, but in the whole families this advantage disappeared when the pseudo-count was scaled (gain 0.0081, 95% interval -0.0162 to 0.0257): it came from the heavy smoothing of a single relative.
- With three or more relatives, BACE’s RPS was lower by 0.002–0.011 under either rule.
- In the whole families, every 95% interval included zero, except at k = 1 with the fixed pseudo-count.
- In the random set, whether the interval excluded zero depended on both k and the smoothing. With the fixed pseudo-count it lay just above zero for k = 3 and 20 (lower limits 0.0001 and 0.0004); with the scaled one, for k = 10 and 20 (0.0007 and 0.0046). Both rules excluded zero at k = 1, and both included it at k = 5.
Giving every tied relative full weight instead (with the fixed pseudo-count; relatives_sensitivity in validation_results.rds) changed some details: with k = 1 in the whole families, where the closest relatives with known behaviour then often form a whole clade, the relatives rule was about as good as BACE, and in the random set BACE’s advantage was distinguishable from zero for k = 10 and 20. Across all these variants the pattern was the same: small gains in RPS for BACE, intervals that include zero for whole families, and a higher balanced accuracy for the relatives rule in most settings. The conclusions above therefore do not depend on the choice of k = 5 or of the smoothing, except that whether BACE’s small advantage in RPS in the random set is distinguishable from zero does.
| Hidden set | k | RPS gain of BACE, fixed pseudo-count (95% interval) | RPS gain of BACE, scaled pseudo-count (95% interval) | Balanced accuracy (relatives) | Relatives used: median (max) |
|---|---|---|---|---|---|
| Random | 1 | 0.0344 (0.0213 to 0.0492) | 0.0177 (0.0075 to 0.0282) | 0.589 | 1 (333) |
| Random | 3 | 0.0090 (0.0001 to 0.0186) | 0.0040 (-0.0038 to 0.0121) | 0.582 | 4 (333) |
| Random | 5 | 0.0046 (-0.0025 to 0.0115) | 0.0046 (-0.0030 to 0.0119) | 0.582 | 7 (866) |
| Random | 10 | 0.0036 (-0.0026 to 0.0094) | 0.0073 (0.0007 to 0.0138) | 0.554 | 14 (1,147) |
| Random | 20 | 0.0057 (0.0004 to 0.0114) | 0.0105 (0.0046 to 0.0165) | 0.537 | 31 (2,376) |
| Whole families | 1 | 0.0392 (0.0052 to 0.0715) | 0.0081 (-0.0162 to 0.0257) | 0.362 | 25 (3,279) |
| Whole families | 3 | 0.0110 (-0.0141 to 0.0293) | 0.0044 (-0.0182 to 0.0214) | 0.367 | 42 (3,279) |
| Whole families | 5 | 0.0047 (-0.0185 to 0.0219) | 0.0047 (-0.0189 to 0.0211) | 0.363 | 42 (3,279) |
| Whole families | 10 | 0.0025 (-0.0214 to 0.0209) | 0.0065 (-0.0167 to 0.0238) | 0.372 | 96 (4,645) |
| Whole families | 20 | 0.0026 (-0.0223 to 0.0238) | 0.0077 (-0.0162 to 0.0258) | 0.372 | 112 (4,645) |
Whole families are the hardest case. Without relatives in the same family, BACE predicted “flocks or colonies” as the most likely category for most species, and was right for 92% of the species that do live in flocks or colonies but for only 2% of the solitary species (1 of 56; Fig. 8 b); the relatives rule recovered none, and recovered more species living singly and in pairs or in pairs and family groups. Its overall accuracy in this set (0.61) was higher than that of the relatives rule (0.54) mainly because “flocks or colonies” is the most common category (58% of these species): predictions that move towards the most common category score well on overall accuracy and on the RPS, but they miss the rarer categories. Neither the relatives rule nor the ordinal regression identified any solitary species in this set, so the limit may lie partly in the information available rather than in BACE alone, although we compared only these methods.
Code
confusion_plot <- function(tab, tag) {
conf <- as.data.frame(tab)
conf$percent <- 100 * conf$Freq / ave(conf$Freq, conf$observed, FUN = sum)
ggplot(conf, aes(predicted, observed, fill = percent)) +
geom_tile(colour = "white") +
geom_text(aes(label = sprintf("%.0f%%\n(%s)", percent, fmt0(Freq))), size = 2.6) +
scale_fill_gradient(low = "white", high = "#21918c", guide = "none") +
scale_x_discrete(labels = short) + scale_y_discrete(limits = rev, labels = short) +
labs(x = "Imputed (most likely)", y = "Observed", tag = tag) + soc_theme +
theme(axis.text.x = element_text(angle = 30, hjust = 1))
}
cal <- val$calibration[val$calibration$n > 0, ]
cal$category <- factor(cal$category, levels = names(soc_cols))
p_c <- ggplot(cal, aes(mean_predicted, observed, colour = category)) +
geom_abline(linetype = "dashed") +
geom_line() + geom_point(aes(size = n)) +
scale_colour_manual(values = soc_cols, labels = short, name = NULL) +
scale_size_area(max_size = 4, guide = "none") +
coord_equal(xlim = c(0, 1), ylim = c(0, 1)) +
labs(x = "Imputed probability", y = "Observed proportion", tag = "(c)") + soc_theme
confusion_plot(val$confusion$random, "(a)") + confusion_plot(val$confusion$family_block, "(b)") + p_c
Averaged over all hidden species, the imputed probabilities matched the observed frequencies of the categories (mean imputed probabilities 0.08, 0.21, 0.25, 0.46 for the four categories, in order, against observed proportions of 0.09, 0.19, 0.23, 0.49). Within ranges of imputed probability, the agreement was looser (Fig. 8 c). For example, among species given a probability of 0.3–0.4 of living in flocks or colonies, 47% did so (193 species), and among species given 0.4–0.5 for pairs and family groups, 61% lived in pairs and family groups (87 species); at the other end, among species given a probability of at most 0.1 of being solitary, 5% were solitary, against 3% predicted. These bins are not adjusted for the dependence among related species, so their deviations are indicative only. Note that BACE’s imputed datasets do not use the most likely category; each final imputation draws a category from its predicted probabilities, so a solitary species with a probability of 0.2 of being solitary is imputed as solitary in about a fifth of the datasets. This is why the imputed category frequencies (Tab. 6) recover more solitary species than the most likely category does, at the cost of more errors elsewhere.
What this means for the imputations above. The full analysis used 1,261 more known values than the validation analysis, which may improve some imputations slightly, but it does not remove the limits found here: imputations of the rarer categories, and of species in poorly described clades, are much less reliable than those of species in flocks or colonies with well-described relatives. Downstream analyses should therefore use all the imputed datasets, which carry this uncertainty, rather than the most likely category of each species. Results that depend on the rarer categories in poorly described clades, such as the number of solitary species in mostly undescribed orders, deserve particular caution, and all results rest on the assumption that missingness does not depend on social behaviour itself (see Who is missing?).
12.6 Imputing an ordinal and a continuous variable together
BACE imputes several variables at once by chaining their imputation models. To illustrate, we started from the validation data (with the hidden species) and also hid the hand-wing index of 2,113 species (20%, drawn at random among species whose index was measured rather than inferred in AVONET). Social behaviour and the hand-wing index were then both missing for 1,037 species. Each variable with missing data gets its own model, and each model uses the other variable as a predictor:
Code
holdout <- readRDS("bird_social/data/holdout.rds")
hwi_mask <- readRDS("bird_social/data/hwi_mask.rds")
mixed_data <- social_data[, soc_vars]
mixed_data$social_max[mixed_data$species %in% holdout$species] <- NA # hidden species
mixed_data$hwi[mixed_data$species %in% hwi_mask] <- NA # hidden wing shapes
mixed_formulas <- list("social_max ~ log_mass + hwi",
"hwi ~ log_mass + social_max")
mixed_chain <- bace_imp(fixformula = mixed_formulas, ran_phylo_form = "~ 1 | species",
phylo = social_tree, data = mixed_data, runs = 25,
nitt = list(50000, 20000), burnin = list(10000, 5000), thin = list(40, 15))
mixed_final <- bace_final_imp(mixed_chain, fixformula = mixed_formulas, ran_phylo_form = "~ 1 | species",
phylo = social_tree, n_final = 20, n_cores = 20, seed = 1,
nitt = list(500000, 100000), burnin = list(50000, 10000), thin = list(450, 90))The hand-wing index is continuous, so BACE fits a Gaussian model to it, and the ordered factor social_max enters that model as a predictor (with polynomial contrasts). MCMC settings can be given per model as lists, in the order of the formulas; the Gaussian model mixes much faster, so its chains can be shorter. Because two variables are missing, each final imputation cycles through both models several times with random draws (the sweeps argument, 3 by default), so that each variable is imputed given plausible values of the other rather than fixed point estimates. The chained phase took 152 minutes and each final imputation 75 minutes.
12.6.1 Convergence
With two variables imputed from each other, both automatic checks flagged social behaviour:
=== BACE Imputation Convergence Assessment ===
Overall Convergence: NO
--- Summary Statistics Method ---
Converged: FALSE
Variables converged: 1 / 2
Use plot() to visualize convergence diagnosticsWe looked at each flag in turn.
The chained loop. assess_convergence() reported that the imputed social behaviour had not settled. Its tests use all 25 iterations. The imputed values changed most over the first four iterations (Fig. 9): the share of species imputed as living in flocks or colonies fell from 53% to 48%, the share imputed as solitary rose from 6.6% to 8.3%, and 9.3% of the imputed species changed category between the first two iterations. From iteration 5 onwards the shares fluctuated without trend, and 2.5–3.8% of species changed category between iterations, the same noise level as in the full analysis, where nothing else was imputed. Applied to iterations 5–25 only, assess_convergence() found no trend and a lag-1 autocorrelation of 0.28, below its threshold of 0.3, and reported convergence (as it did for iterations 10–25). The flag therefore reflected the initial adjustment, as the two models adapted to each other’s imputations, rather than a lack of stability at the end of the chained phase, from which the final imputations start. The mean imputed hand-wing index barely moved across iterations (25.92 to 25.95). In your own analyses, plot the series rather than relying on the flag alone, and run more iterations if the series is still moving at the end.
Code
long <- do.call(rbind, lapply(names(soc_labels), function(k)
data.frame(iteration = ms$iteration, category = k, share = ms[[k]])))
long$category <- factor(long$category, levels = names(soc_labels))
p_a <- ggplot(long, aes(iteration, 100 * share, colour = category)) +
geom_vline(xintercept = 4.5, linetype = "dashed", colour = "grey50") +
geom_line() + geom_point(size = 1) +
scale_colour_manual(values = soc_cols, labels = short, name = NULL) +
labs(x = "Chained iteration", y = "Imputed species (%)", tag = "(a)") + soc_theme
p_b <- ggplot(ms[-1, ], aes(iteration, 100 * changed)) +
geom_vline(xintercept = 4.5, linetype = "dashed", colour = "grey50") +
geom_line() + geom_point(size = 1) +
labs(x = "Chained iteration", y = "Changed category (%)", tag = "(b)") + soc_theme
p_a + p_b + plot_layout(widths = c(1.4, 1))
The models of the last chained iteration. The printed summary marked 4 parameters of the social behaviour model “Check”: the phylogenetic variance and the two cutpoints (ESS 89, 98, 77 of 1,000 draws, with Geweke p < 0.05 for the cutpoints) and the effect of body mass (Geweke p = 0.028, with ESS 1,000). The same short chains in the validation analysis, with social behaviour as the only missing variable, also had low ESS for these parameters (93, 122, 90). As in the full analysis, these chains only provide point imputations; the final chains are 10 times longer.
The final models. Within each final imputation, the social behaviour model mixed about as well as in the full analysis (Tab. 9, Fig. 10): at least 474 effective draws per chain for the phylogenetic variance and the cutpoints, and split \(\hat{R}\) ≤ 1.02. Geweke tests flagged 10%–15% of the chains for the phylogenetic variance and the cutpoints (none significantly more than 5% by a binomial test), and the placebo test away from the start of the chains flagged 0%–10%, so small early drifts in a few chains cannot be excluded; dropping the first 10% of every chain changed the pooled means by at most 0.004 posterior standard deviations.
The \(\hat{R}\) across the final imputations needs a different reading here. In this analysis, each final model of social behaviour conditions on the hand-wing indices imputed in the preceding sweep, which differ between imputations, so the 20 chains target 20 different posteriors, one per completed dataset. \(B\), the variance of the chain means, then contains the between-imputation variance of the estimates as well as Monte Carlo error, and \(\hat{R}\) above 1 is expected for parameters that depend on the imputed values. This is what happened for the effect of the hand-wing index on social behaviour (\(\hat{R}\) = 1.12, \(B/W\) = 0.25). Four checks indicate that this reflects the imputations rather than a failure to converge:
- Monte Carlo error alone would give \(B/W \approx\) 0.0011 (one over the ESS per chain); the observed value is about 230 times larger.
- In the validation analysis, the same model with the same MCMC settings but a complete hand-wing index had \(B/W\) = 0.0015 and \(\hat{R}\) = 1.000 for this parameter.
- \(B/W\) was the same in the first and second halves of the chains (0.26 and 0.24): the chains were not moving towards a common value, so longer chains would not reduce it.
- Within each chain, this parameter showed no sign of a problem (split \(\hat{R}\) ≤ 1.01, ESS ≥ 621).
Treating \(W\) as the within-imputation variance, this corresponds to a fraction of missing information of about 0.21 for this effect: \(\lambda = (1 + 1/M)B / \{W + (1 + 1/M)B\}\) with \(M\) = 20 imputations, the same to two decimals as Rubin’s fmi with the small-sample correction used by pool_mi(). Monte Carlo error, included in \(B\), contributes less than 1% of it. Parameters that depend less directly on the imputed wing shapes, such as the phylogenetic variance and the cutpoints, had \(B/W\) of 0.006–0.008. Because \(\hat{R}\) across imputations mixes the two sources of variation, use it here only together with within-chain checks, as above.
| Parameter | R-hat across imputations | ESS per chain (median) | ESS per chain (min) | Geweke flag rate | Geweke placebo rate | Split R-hat (max) | Drop first 10% (SD) |
|---|---|---|---|---|---|---|---|
| (Intercept) | 1.001 | 1000 | 911 | 0.05 | 0.15 | 1.012 | 0.000 |
| log_mass | 1.003 | 1000 | 902 | 0.05 | 0.15 | 1.012 | 0.003 |
| hwi | 1.125 | 1000 | 621 | 0.00 | 0.10 | 1.011 | -0.003 |
| phylo_var | 1.003 | 629 | 474 | 0.10 | 0.10 | 1.015 | -0.003 |
| cutpoint2 | 1.003 | 668 | 557 | 0.10 | 0.00 | 1.007 | -0.002 |
| cutpoint3 | 1.004 | 612 | 492 | 0.15 | 0.00 | 1.012 | -0.003 |
| h2_latent | 1.003 | 650 | 534 | 0.10 | 0.05 | 1.016 | -0.003 |
| hwi_phylo_var | 1.011 | 1000 | 1000 | 0.15 | 0.05 | 1.003 | -0.004 |
| hwi_residual_var | 1.010 | 1000 | 825 | 0.10 | 0.00 | 1.008 | 0.004 |
Code
tr <- mix$trace
tr$parameter <- factor(tr$parameter, levels = c("phylo_var", "cutpoint2", "cutpoint3"),
labels = c("Phylogenetic variance", "Cutpoint 2", "Cutpoint 3"))
ggplot(tr, aes(draw, value, colour = factor(chain))) +
geom_line(linewidth = 0.25, alpha = 0.8) +
facet_wrap(~ parameter, ncol = 1, scales = "free_y") +
labs(x = "Draw", y = NULL, colour = "Chain") + soc_theme
12.6.2 Accuracy of the imputations
| Method | RMSE | r | Coverage |
|---|---|---|---|
| BACE (mean of imputations) | 4.00 | 0.96 | 86% |
| Observed mean | 15.16 | ||
| Linear regression on log body mass | 14.82 | 0.21 | |
| Weighted mean of closest relatives (k = 5) | 4.29 | 0.96 |
BACE imputed the hidden hand-wing index more accurately than the baselines (Tab. 10), although the gain over the closest relatives was modest (RMSE 4.00 against 4.29), and the hidden indices were drawn at random, so this checks scattered gaps, not whole clades. BACE was also more accurate than the weighted mean of the closest relatives for every k we tried (RMSE 4.21–5.00 for k = 1–20). The improvement over the weighted mean of the closest relatives (k = 5, ties weighted as above) was clear (difference in mean squared error 2.46, 95% interval from resampling families 0.81 to 4.63), and largest for species whose social behaviour was also missing (RMSE 4.12 against 4.54, compared with 3.88 against 4.05 when it was observed). The central 90% of the imputed values (5th to 95th percentiles of the 20 imputations) contained 86% of the measured values. With 20 imputations, such intervals would contain only about 82% of the values even if the imputations were perfectly calibrated, so the imputed values were, if anything, slightly too dispersed. The hand-wing index is a continuous trait with very strong phylogenetic signal, an easier case than social behaviour.
Hiding a fifth of the hand-wing indices did not change the accuracy of the imputed social behaviour. For the hidden species, the RPS was 0.115 (random) and 0.138 (whole families), against 0.115 and 0.137 in the validation analysis, where the hand-wing index was complete. This was also true for the hidden species whose hand-wing index was hidden too (differences of 0.001 and 0.001). Social behaviour is predicted mostly by the phylogeny, so replacing measured wing shapes by imputed ones, which themselves follow the phylogeny closely, changed little. Joint imputation therefore did not make the social behaviour imputations worse here, but it did not make them better either.
How imputation uncertainty flows between variables
Within each final imputation, BACE cycles through the imputation models 3 times, drawing new values at each pass, and keeps the last pass. The social behaviour model therefore conditions on hand-wing indices that differ between imputations: at the species whose index was missing, the values it used varied across the 20 imputations (median standard deviation 0.24 on the standardised scale). The uncertainty in the imputed hand-wing index is thus carried into the imputed social behaviour, and into the pooled estimates, such as the effect of the hand-wing index whose posterior means differ between imputations (see Convergence above).
12.7 Summary: how far to trust these results
- What the example shows. How to run BACE on an ordinal variable with many missing values, alone and together with a continuous variable; how to check convergence and the imputations; and how to carry the imputation uncertainty into a downstream analysis. It does not show that BACE predicts social behaviour better than simpler phylogenetic approaches: the closest-relatives rule did almost as well. What BACE added here was a model rather than better predictions: estimates of phylogenetic signal and predictor effects, joint imputation of several variables, and imputations whose uncertainty, from the model and its parameters, can be propagated through multiple imputation. The relatives rule gives probabilities for each species but would need further assumptions to provide these.
- The MCMC. In the full and validation analyses, the independent chains of the final models agreed with each other; in all analyses, each chain had at least about 300 effective draws for the slowest parameters. Geweke tests flagged more chains than expected for one cutpoint of the validation analysis, also away from the start of the chains; dropping early draws did not change the estimates materially, and longer chains would be a cheap safeguard. In the mixed analysis, \(\hat{R}\) across imputations was above 1 for the effect of the imputed variable, as expected when the imputations differ. Other automatic flags came from the short chains of the chained phase and from the first iterations of the chained loop. We checked each flag rather than relying on the overall verdict.
- The imputations. On hidden species, BACE gave better probabilistic predictions than a model without phylogeny, with probabilities that matched the category frequencies on average but not within every range, but it was only slightly better than copying the closest relatives (for three or more relatives, by at most 0.011 in RPS, with 95% intervals that included zero in most settings and always for whole families), and it was worse at recovering the rarer categories. For whole families without described relatives, it mostly predicted the most common category and almost never identified solitary species.
- The missing species. Missingness follows taxonomy, which BACE uses, but many species without records have few described relatives, and missingness may also depend on social behaviour itself. The full analysis does not remove these limits.
- Downstream analyses. Fit the downstream model to every imputed dataset and pool the results, rather than using the most likely category, and include the variables of the downstream analysis in the imputation model. Here most of the uncertainty about the social behaviour coefficients came from the missing values (fmi about 0.7), so these results depend strongly on the imputation model and on the missing-at-random assumption; imputation does not by itself remove bias if social behaviour is missing not at random. Treat conclusions that depend on rare categories in poorly described clades with caution.
13 Citation
If you use BACE in a publication, please cite both the package (see citation("BACE")) and the underlying MCMCglmm engine (Hadfield 2010).
Bug reports, feature requests and worked examples are very welcome at the BACE issue tracker.
14 References
Barnard, John, and Donald B. Rubin. 1999. “Small-Sample Degrees of Freedom with Multiple Imputation.” Biometrika 86 (4): 948–55. https://doi.org/10.1093/biomet/86.4.948.
Blomberg, Simon P., Theodore Garland, and Anthony R. Ives. 2003. “Testing for Phylogenetic Signal in Comparative Data: Behavioral Traits Are More Labile.” Evolution 57 (4): 717–45. https://doi.org/10.1111/j.0014-3820.2003.tb00285.x.
Fritz, Susanne A., and Andy Purvis. 2010. “Selectivity in Mammalian Extinction Risk and Threat Types: A New Measure of Phylogenetic Signal Strength in Binary Traits.” Conservation Biology 24 (4): 1042–51. https://doi.org/10.1111/j.1523-1739.2010.01455.x.
Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. 3rd ed. Boca Raton, FL: Chapman; Hall/CRC. https://doi.org/10.1201/b16018.
Hadfield, Jarrod D. 2010. “MCMC Methods for Multi-Response Generalized Linear Mixed Models: The MCMCglmm r Package.” Journal of Statistical Software 33 (2): 1–22. https://doi.org/10.18637/jss.v033.i02.
Hadfield, Jarrod D., and Shinichi Nakagawa. 2010. “General Quantitative Genetic Methods for Comparative Biology: Phylogenies, Taxonomies and Multi-Trait Models for Continuous and Categorical Characters.” Journal of Evolutionary Biology 23 (3): 494–508. https://doi.org/10.1111/j.1420-9101.2009.01915.x.
Ho, Lam si Tung, and Cécile Ané. 2014. “A Linear-Time Algorithm for Gaussian and Non-Gaussian Trait Evolution Models.” Systematic Biology 63 (3): 397–408. https://doi.org/10.1093/sysbio/syu005.
McTavish, Emily Jane, Jeff A. Gerbracht, Mark T. Holder, Marshall J. Iliff, Denis Lepage, Pamela C. Rasmussen, Benjamin D. Redelings, Luna L. Sánchez Reyes, and Eliot T. Miller. 2025. “A Complete and Dynamic Tree of Birds.” Proceedings of the National Academy of Sciences 122 (18): e2409658122. https://doi.org/10.1073/pnas.2409658122.
Nakagawa, Shinichi, and Pierre de Villemereuil. 2019. “A General Method for Simultaneously Accounting for Phylogenetic and Species Sampling Uncertainty via Rubin’s Rules in Comparative Analysis.” Systematic Biology 68 (4): 632–41. https://doi.org/10.1093/sysbio/syy089.
Pagel, Mark. 1999. “Inferring the Historical Patterns of Biological Evolution.” Nature 401 (6756): 877–84. https://doi.org/10.1038/44766.
Pottier, Patrice, Michael R. Kearney, Nicholas C. Wu, Alex R. Gunderson, Julie E. Rej, A. Nayelli Rivera-Villanueva, Pietro Pollo, Samantha Burke, Szymon M. Drobniak, and Shinichi Nakagawa. 2025. “Vulnerability of Amphibians to Global Warming.” Nature 639 (8056): 954–61. https://doi.org/10.1038/s41586-025-08665-0.
Rubin, Donald B. 1987. Multiple Imputation for Nonresponse in Surveys. New York: John Wiley & Sons. https://doi.org/10.1002/9780470316696.
Şekercioğlu, Çağan H., Kyle D. Kittelberger, Flavio M. M. Mota, Amy N. Buxton, Nikolas Orton, Adara DeNiro, Evan R. Buechley, et al. 2025. “BIRDBASE: A Global Dataset of Avian Biogeography, Conservation, Ecology and Life History Traits.” Scientific Data 12: 1558. https://doi.org/10.1038/s41597-025-05615-3.
Sheard, Catherine, Montague H. C. Neate-Clegg, Nico Alioravainen, Samuel E. I. Jones, Claire Vincent, Hannah E. A. MacGregor, Tom P. Bregman, Santiago Claramunt, and Joseph A. Tobias. 2020. “Ecological Drivers of Global Gradients in Avian Dispersal Inferred from Wing Morphology.” Nature Communications 11: 2463. https://doi.org/10.1038/s41467-020-16313-6.
Tobias, Joseph A., Catherine Sheard, Alex L. Pigot, Adam J. M. Devenish, Jingyi Yang, Ferran Sayol, et al. 2022. “AVONET: Morphological, Ecological and Geographical Data for All Birds.” Ecology Letters 25 (3): 581–97. https://doi.org/10.1111/ele.13898.
van Buuren, Stef. 2018. Flexible Imputation of Missing Data. 2nd ed. Boca Raton, FL: Chapman; Hall/CRC. https://doi.org/10.1201/9780429492259.
White, Ian R., and John B. Carlin. 2010. “Bias and Efficiency of Multiple Imputation Compared with Complete-Case Analysis for Missing Covariate Values.” Statistics in Medicine 29 (28): 2920–31. https://doi.org/10.1002/sim.3944.
Zhou, Xiang, and Jerome P. Reiter. 2010. “A Note on Bayesian Inference After Multiple Imputation.” The American Statistician 64 (2): 159–63. https://doi.org/10.1198/tast.2010.09109.