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 Phylo-BACE
General-purpose multiple-imputation tools such as mice, Amelia and missForest assume that rows are exchangeable: 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, and throwing that information away means we impute a missing mouse trait as though the nearest source of information were the grand mean across 500 mammals rather than the handful of other mouse-like rodents in the tree.
Phylo-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 Interpreting the pooled output below).
Phylo-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 an MCMCglmm fit with a phylo random effect) and check whether the phylogenetic variance component is non-trivial.
You want posterior inference on downstream models rather than point imputations followed by ad-hoc variance adjustments.
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 cheaper.
You need sub-second imputations in a production pipeline. Each chained-equations iteration runs an MCMC fit, so runtimes are measured in minutes, not milliseconds.
For a broader discussion of the data-gap problem and imputation in comparative biology we recommend Dalia A. Conde et al. (2019) and Pottier et al. (2024); for Rubin’s-rules-based approaches to related problems see Nakagawa and de Villemereuil (2018).
3 Installation
You can install the BACE package from GitHub using the following command:
Code
# Install from github and load#install.packages("pacman") pacman::p_load(devtools)devtools::install_github("daniel1noble/BACE")
── R CMD build ─────────────────────────────────────────────────────────────────
* checking for file ‘/private/var/folders/48/07fr3p6n4nv8gq_s3kl6vj2c0000gn/T/Rtmp1OMcPa/remotes4ef0189743f6/daniel1noble-BACE-b1e07ae/DESCRIPTION’ ... OK
* preparing ‘BACE’:
* checking DESCRIPTION meta-information ... OK
* checking for LF line-endings in source and make files and shell scripts
* checking for empty or unneeded directories
Removed empty directory ‘BACE/.vscode’
Removed empty directory ‘BACE/dev’
Removed empty directory ‘BACE/ms’
Removed empty directory ‘BACE/vignettes’
* building ‘BACE_0.0.0.9000.tar.gz’
If you encounter any issues during installation, please ensure that you have the required dependencies installed.
4 Documentation
The main functions in the BACE package are bace() and bace_imp(), which performs Bayesian imputation using chained equations. You can access the documentation for these functions using the following commands:
Code
?bace
As you will see, the function allows you to specify fixed effects, random effects, 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. For most comparative datasets, this will be sufficient for modeling the phylogenetic and species-level variation in the data, but not always.
5 The BACE workflow at a glance
Fig. 1 maps the user-facing functions onto the workflow. bace() is the top-level wrapper. 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 a MCMCglmm phylogenetic mixed model, impute the missing cells, move on to the next variable, and repeat for runs iterations. 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(). Check whether the sequence of imputed datasets has stabilised. This is not within-chain MCMC convergence — see Two kinds of convergence below.
Final imputation runs — bace_final_imp(). Starting from the converged dataset, launch n_finalindependent MCMCglmm fits per variable (default n_final = 50; parallelisable via n_cores). Each run produces one complete imputed dataset and one per-variable posterior. 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 entirely.
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 Downstream analysis with Rubin’s rules below. The rest of this vignette walks through the workflow first with bace() and then shows how to run each step individually.
Figure 1— Core user-facing BACE function map. Solid arrows show the main workflow from raw comparative data through chained imputation (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()) 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. 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. Dashed boxes and arrows show optional or advanced routes: simulating benchmark data (sim_bace()), pre-flight phylogenetic-signal checks (phylo_signal_summary()), 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 reproducibilityset.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 sparsity
We 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 datahead(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).
D (binary variables only) — Fritz-Purvis D (Fritz and Purvis 2010), via caper::phylo.d().
Convergence flags — every fit is checked with coda::effectiveSize() and coda::geweke.diag(); rows with low effective sample size or failed Geweke tests are flagged unreliable and 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 honours BACE’s two-RE modes: species = FALSE fits a single phylogenetic random effect; species = TRUE fits the dual structure (phylogenetic + non-phylogenetic species) that 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)
Or, more concisely, call bace() with phylo_signal = TRUE to get the same diagnostic without running the imputation pipeline:
Code
# Signal-only preview: returns a phylo_signal object and does NOT# call bace_imp(), bace_final_imp(), or pool_posteriors().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.
Reviewer’s caveat (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.
6.3 Using BACE for Imputation
Now that we have the data and the phylogenetic tree we can run the imputation using the bace() or bace_imp() function.
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 integer whose values are all non-negative is imputed as a Poisson count — deliberately so, with no goodness-of-fit “race” against the gaussian family, because a count trait with strong phylogenetic signal is marginally overdispersed and a marginal-fit comparison would silently (and wrongly) pick gaussian exactly when the phylogenetic signal is strongest. 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). Setting ovr_categorical = TRUE instead 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 typesbace_final1$initial_results$types
Importantly, 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 users, 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 Analysis with 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: 4385.3
Variance components: 2591.7
Convergence Assessment:
Acceptable : 1 parameters (6.7%)
Check : 6 parameters (40.0%)
Fixed : 1 parameters (6.7%)
Good : 7 parameters (46.7%)
Overall: CHECK - Many parameters show convergence issues
Consider increasing nitt, adjusting priors, or checking model specification
Convergence by Response Variable:
----------------------------------------
[CH] y : 2/ 5 parameters good
Check convergence: x2.L, Species
[OK] x1 : 4/ 5 parameters good
Check convergence: (Intercept)
[CH] x2 : 1/ 4 parameters good (1 fixed)
Check convergence: x1, y, Species
========================================
Step 2: Convergence assessment
-------------------------------------------------------
Convergence not achieved. Running additional iterations...
Attempt 2 of 3
-------------------------------------------------------
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
Run 16, Imputed variable: y
Run 16, Imputed variable: x1
Run 16, Imputed variable: x2
Run 17, Imputed variable: y
Run 17, Imputed variable: x1
Run 17, Imputed variable: x2
Run 18, Imputed variable: y
Run 18, Imputed variable: x1
Run 18, Imputed variable: x2
Run 19, Imputed variable: y
Run 19, Imputed variable: x1
Run 19, Imputed variable: x2
Run 20, Imputed variable: y
Run 20, Imputed variable: x1
Run 20, Imputed variable: x2
Run 21, Imputed variable: y
Run 21, Imputed variable: x1
Run 21, Imputed variable: x2
Run 22, Imputed variable: y
Run 22, Imputed variable: x1
Run 22, 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: 4323.9
Variance components: 2915.8
Convergence Assessment:
Acceptable : 2 parameters (13.3%)
Check : 5 parameters (33.3%)
Fixed : 1 parameters (6.7%)
Good : 7 parameters (46.7%)
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 : 4/ 5 parameters good
Check convergence: x2.L
[CH] x2 : 0/ 4 parameters good (1 fixed)
Check convergence: (Intercept), x1, y, Species
========================================
Convergence not achieved. Running additional iterations...
Attempt 3 of 3
-------------------------------------------------------
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
Run 16, Imputed variable: y
Run 16, Imputed variable: x1
Run 16, Imputed variable: x2
Run 17, Imputed variable: y
Run 17, Imputed variable: x1
Run 17, Imputed variable: x2
Run 18, Imputed variable: y
Run 18, Imputed variable: x1
Run 18, Imputed variable: x2
Run 19, Imputed variable: y
Run 19, Imputed variable: x1
Run 19, Imputed variable: x2
Run 20, Imputed variable: y
Run 20, Imputed variable: x1
Run 20, Imputed variable: x2
Run 21, Imputed variable: y
Run 21, Imputed variable: x1
Run 21, Imputed variable: x2
Run 22, Imputed variable: y
Run 22, Imputed variable: x1
Run 22, Imputed variable: x2
Run 23, Imputed variable: y
Run 23, Imputed variable: x1
Run 23, Imputed variable: x2
Run 24, Imputed variable: y
Run 24, Imputed variable: x1
Run 24, Imputed variable: x2
Run 25, Imputed variable: y
Run 25, Imputed variable: x1
Run 25, Imputed variable: x2
Run 26, Imputed variable: y
Run 26, Imputed variable: x1
Run 26, Imputed variable: x2
Run 27, Imputed variable: y
Run 27, Imputed variable: x1
Run 27, Imputed variable: x2
Run 28, Imputed variable: y
Run 28, Imputed variable: x1
Run 28, Imputed variable: x2
Run 29, Imputed variable: y
Run 29, Imputed variable: x1
Run 29, Imputed variable: x2
Run 30, Imputed variable: y
Run 30, Imputed variable: x1
Run 30, Imputed variable: x2
Run 31, Imputed variable: y
Run 31, Imputed variable: x1
Run 31, Imputed variable: x2
Run 32, Imputed variable: y
Run 32, Imputed variable: x1
Run 32, Imputed variable: x2
Run 33, Imputed variable: y
Run 33, Imputed variable: x1
Run 33, Imputed variable: x2
Run 34, Imputed variable: y
Run 34, Imputed variable: x1
Run 34, 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: 4413.5
Variance components: 2779.4
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
========================================
Step 3: Final imputation runs for posterior pooling
-------------------------------------------------------
=== Running 10 final imputation iterations ===
Run 1/10 - Imputed variable: y
Run 1/10 - Imputed variable: x1
Run 1/10 - Imputed variable: x2
Run 2/10 - Imputed variable: y
Run 2/10 - Imputed variable: x1
Run 2/10 - Imputed variable: x2
Run 3/10 - Imputed variable: y
Run 3/10 - Imputed variable: x1
Run 3/10 - Imputed variable: x2
Run 4/10 - Imputed variable: y
Run 4/10 - Imputed variable: x1
Run 4/10 - Imputed variable: x2
Run 5/10 - Imputed variable: y
Run 5/10 - Imputed variable: x1
Run 5/10 - Imputed variable: x2
Run 6/10 - Imputed variable: y
Run 6/10 - Imputed variable: x1
Run 6/10 - Imputed variable: x2
Run 7/10 - Imputed variable: y
Run 7/10 - Imputed variable: x1
Run 7/10 - Imputed variable: x2
Run 8/10 - Imputed variable: y
Run 8/10 - Imputed variable: x1
Run 8/10 - Imputed variable: x2
Run 9/10 - Imputed variable: y
Run 9/10 - Imputed variable: x1
Run 9/10 - Imputed variable: x2
Run 10/10 - Imputed variable: y
Run 10/10 - Imputed variable: x1
Run 10/10 - Imputed variable: x2
=== Final imputation complete ===
Saved 10 imputed datasets with corresponding models
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: 3
Final imputations: 10
Access results via:
- $pooled_models: Pooled posterior distributions
- $imputed_datasets: List of 10 imputed datasets
- $convergence: Convergence diagnostics
- $final_results: Full final imputation results
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 that, 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 sometimes convergence fail yet convergence is likely achieved. If you think this is the case than you can do initial runs and skip the convergence assessment (see below). Note, 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_conv argument 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_size argument 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 final bace object.
Third, the species argument 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_final argument 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 embarrassingly parallel, so on multi-core machines you can set n_cores to run them concurrently.
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: 256.7651
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 0.1227 1.15e-09 0.4834 10000
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.4435 0.03414 0.869 10000
Location effects: y ~ x1 + x2
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 0.05108 -0.39972 0.47269 10000 0.7646
x1 -0.01490 -0.28788 0.27186 10000 0.9176
x2.L 0.62260 0.25336 1.02159 10663 0.0016 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
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
# You can also extract all pooled models at onceall_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)
# Or stack all imputed datasets into one data frame with an .imputation columnimputed_df <-get_imputed_data(bace_final1, format ="data.frame")table(imputed_df$.imputation)
It is worth being precise about what the pooled posterior returned by bace() actually is, because the terminology in the multiple-imputation literature can trip people up.
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 a Monte Carlo 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, and are the right tool when the per-imputation analysis produces frequentist estimates rather than full posterior draws. When you have posterior draws, stacking is both simpler and more faithful — do not apply Rubin’s formulas on top of the stacked bace() posteriors. If your downstream analysis model is frequentist (an lm, glm, nlme::gls, phylolm, …), BACE does provide a first-class Rubin’s-rules pathway — with_imputations() + pool_mi() — described in the next section.
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'first_fit <- bace_final1$final_results$all_models[[1]]$yplot(first_fit)coda::effectiveSize(first_fit$Sol)coda::autocorr.diag(first_fit$Sol)
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 .f to each of the n_final imputed datasets and returns the list of fits. It accepts the output of bace() or bace_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 matches mice::pool to 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 datasetfits <-with_imputations(bace_final1,function(d) glm(y ~ x1 + x2, family = poisson, data = d))# 2. Pool with Rubin's rulespool_mi(fits)
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 rulesfits_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 with pool_mi(); use pool_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 formulabace_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: 795.9
Variance components: 455.6
Convergence Assessment:
Acceptable : 1 parameters (6.7%)
Check : 5 parameters (33.3%)
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 : 2/ 5 parameters good
Check convergence: (Intercept), Species
[OK] x1 : 5/ 5 parameters good
[CH] x2 : 1/ 4 parameters good (1 fixed)
Check convergence: (Intercept), x1, 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 modelsummary(model_y)
Iterations = 5001:19991
Thinning interval = 10
Sample size = 1500
DIC: 256.6727
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 0.113 4.204e-10 0.4508 154.9
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.4364 0.03424 0.8518 270.2
Location effects: y ~ x1 + x2
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 0.05925 -0.35608 0.43329 668.9 0.713
x1 -0.03932 -0.32601 0.23408 761.7 0.755
x2.L 0.63752 0.28041 1.03565 837.8 <7e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# Plot trace and density plots for model parametersplot(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.
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 speciesset.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 datahead(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 speciesdim(data2)
[1] 100 6
Code
length(unique(data2$Species))
[1] 20
We 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:
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: 882.7
Variance components: 380.0
Convergence Assessment:
Acceptable : 3 parameters (16.7%)
Check : 5 parameters (27.8%)
Fixed : 1 parameters (5.6%)
Good : 9 parameters (50.0%)
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
[CH] x2 : 1/ 5 parameters good (1 fixed)
Check convergence: x1, y, Species, Species2
========================================
Code
# We can view the final model for y to check that both random effects are includedmodel_y2 <- bace_impute3$models_last_run[["y"]]summary(model_y2)
Iterations = 5001:19991
Thinning interval = 10
Sample size = 1500
DIC: 268.3542
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 0.1375 3.504e-08 0.6109 372.9
~Species2
post.mean l-95% CI u-95% CI eff.samp
Species2 0.2005 1.337e-07 0.6195 225.8
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.3378 0.01776 0.714 179.8
Location effects: y ~ x1 + x2
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 0.10581 -0.43835 0.69384 705.2 0.6360
x1 -0.02675 -0.39081 0.34771 636.6 0.8680
x2.L 0.53160 0.09310 0.96605 505.4 0.0173 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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:
Iterations = 5001:19991
Thinning interval = 10
Sample size = 1500
DIC: 75.05073
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 2.405 4.938e-07 9.064 63.57
~Species2
post.mean l-95% CI u-95% CI eff.samp
Species2 1.524 2.565e-05 6.991 67.59
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.23622 -1.63748 2.18353 683.40 0.7347
x1 -0.70334 -2.22210 0.26925 77.51 0.1520
y 0.30140 -0.01043 0.63234 582.99 0.0413 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
As 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 parametersbace_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: 4030.9
Variance components: 1405.1
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 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 MCMCglmm fit. 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 individualMCMCglmm object. bace_imp() prints a short summary of effective sample sizes and autocorrelation for the last-iteration models when verbose = TRUE, and you can always call plot() or coda::effectiveSize() on bace_impute$models_last_run[["y"]] or on any of the per-imputation fits stored in bace_final1$final_results$all_models[[m]][["y"]]. Failures here are fixed by increasing nitt, burnin, or thin, tightening priors, 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 the runs chained iterations. Failures here are fixed by increasing runs or, 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 Interpreting the pooled output above.
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 neededconverge <-assess_convergence(bace_impute2, method ="summary")plot(converge)
Code
# Now that we know that we can run BACE longer# Fit BACE with a single formulabace_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: 4524.0
Variance components: 2761.7
Convergence Assessment:
Acceptable : 1 parameters (6.7%)
Check : 6 parameters (40.0%)
Fixed : 1 parameters (6.7%)
Good : 7 parameters (46.7%)
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
Check convergence: units
[OK] x1 : 4/ 5 parameters good
Check convergence: units
[CH] x2 : 0/ 4 parameters good (1 fixed)
Check convergence: (Intercept), x1, y, Species
========================================
Code
# Check again for convergenceconverge <-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 Analysis with 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: 5158.2
Variance components: 2157.2
Convergence Assessment:
Acceptable : 8 parameters (44.4%)
Fixed : 1 parameters (5.6%)
Good : 9 parameters (50.0%)
Overall: ACCEPTABLE - Most parameters converged adequately
Convergence by Response Variable:
----------------------------------------
[~] y : 3/ 6 parameters good
[~] x1 : 4/ 6 parameters good
[~] 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 - Imputed variable: y
Run 1/10 - Imputed variable: x1
Run 1/10 - Imputed variable: x2
Run 2/10 - Imputed variable: y
Run 2/10 - Imputed variable: x1
Run 2/10 - Imputed variable: x2
Run 3/10 - Imputed variable: y
Run 3/10 - Imputed variable: x1
Run 3/10 - Imputed variable: x2
Run 4/10 - Imputed variable: y
Run 4/10 - Imputed variable: x1
Run 4/10 - Imputed variable: x2
Run 5/10 - Imputed variable: y
Run 5/10 - Imputed variable: x1
Run 5/10 - Imputed variable: x2
Run 6/10 - Imputed variable: y
Run 6/10 - Imputed variable: x1
Run 6/10 - Imputed variable: x2
Run 7/10 - Imputed variable: y
Run 7/10 - Imputed variable: x1
Run 7/10 - Imputed variable: x2
Run 8/10 - Imputed variable: y
Run 8/10 - Imputed variable: x1
Run 8/10 - Imputed variable: x2
Run 9/10 - Imputed variable: y
Run 9/10 - Imputed variable: x1
Run 9/10 - Imputed variable: x2
Run 10/10 - Imputed variable: y
Run 10/10 - Imputed variable: x1
Run 10/10 - Imputed variable: x2
=== Final imputation complete ===
Saved 10 imputed datasets with corresponding models
Pooling posteriors with sampling: 1000 samples per imputation
=======================================================
BACE Analysis Complete!
=======================================================
Convergence achieved: FALSE
Number of attempts: 1
Final imputations: 10
Access results via:
- $pooled_models: Pooled posterior distributions
- $imputed_datasets: List of 10 imputed datasets
- $convergence: Convergence diagnostics
- $final_results: Full final imputation results
Code
# Extract the pooled model for y using the accessor functiony_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: 269.8538
G-structure: ~Species
post.mean l-95% CI u-95% CI eff.samp
Species 0.1413 4.585e-10 0.6245 10000
~Species2
post.mean l-95% CI u-95% CI eff.samp
Species2 0.219 1.156e-09 0.6381 10000
R-structure: ~units
post.mean l-95% CI u-95% CI eff.samp
units 0.3188 0.01614 0.6941 10000
Location effects: y ~ x1 + x2
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) 0.11424 -0.45067 0.68750 10000 0.618
x1 -0.03130 -0.40020 0.33314 9754 0.836
x2.L 0.50093 0.06126 0.95716 10000 0.028 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
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 BACE on simulated data with known truth
For methodological work and quantitative claims about imputation accuracy, simulation-based benchmarking is essential. sim_bace() generates phylogenetically-structured comparative data with a complete pre-NA dataset returned alongside the masked one, so we can compare BACE’s imputed values against the actual held-out truth.
The pattern is:
Simulate a complete dataset with a known phylogenetic signal (sim_bace() returns both data with NAs and complete_data with the truth).
Mask a fraction of one variable as MCAR (or use the built-in missingness argument for type-aware mechanisms).
Run BACE.
Compare imputed values at the masked rows against truth using type-appropriate metrics (NRMSE / accuracy / Brier).
Compare BACE against a no-skill baseline — column-mean for continuous, modal class for categorical. BACE has to beat this floor to demonstrate it is using the phylogenetic structure.
Below we run a single replicate for one continuous (gaussian) and one categorical (multinomial3) response. For a manuscript-quality benchmark you should run many replicates per type and report mean ± SD; see dev/sim_accuracy_multirep.R in the package source for a full 20-replicate harness across all five response types.
BACE’s NRMSE is roughly half the no-skill floor on a high-signal phylogenetic trait, and the empirical 95% prediction interval covers truth ~85–95% of the time at n_final = 50 (well-calibrated; see ?bace for why smaller n_final undercovers).
The same pattern adapts to categorical responses by swapping NRMSE for balanced accuracy and the column-mean baseline for a modal-class baseline:
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), and the per-cell 95% prediction intervals are well calibrated:
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 Citation and further reading
If you use BACE in a publication, please cite both the package (see citation("BACE")) and the underlying MCMCglmm engine (Hadfield and Nakagawa 2010). Phylogenetic comparative methods and phylogenetic mixed models are reviewed in Cornwell and Nakagawa (2017), de Villemereuil, Nakagawa, and Garamszegi (2014) and Halliwell, Holland, and Yates (2025). The broader motivation for imputation in comparative biology — data gaps, trait databases, and the cost of complete-case analysis — is discussed in Dalia A. Conde et al. (2019) and Pottier et al. (2024). For a related Rubin’s-rules-based approach that handles phylogenetic uncertainty (rather than missing trait data), see Nakagawa and de Villemereuil (2018).
The Bayesian combiner used by pool_posteriors() — concatenating the per-imputation chains to approximate the marginal posterior — is discussed in Zhou and Reiter (2010) and in Chapter 18 of Gelman et al. (2013). The Rubin’s-rules combiner used by pool_mi() is from Rubin (1987) with the small-sample degrees-of-freedom correction of Barnard and Rubin (1999); van Buuren (2018) is the standard modern reference on multiple imputation practice.
Bug reports, feature requests and worked examples are very welcome at the BACE issue tracker.
12 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.
Dalia A. Conde, Johanna Staerk, Fernando Colchero, Rita Da Silva, Hugh P. Possingham, Annette Baudisch, and James W. Vaupel. 2019. “Data Gaps and Opportunities for Comparative and Conservation Biology.”Pnas 116 (10): 9658–64. https://doi.org/10.1073/pnas.1816367116.
de Villemereuil, Pierre, Shinichi Nakagawa, and L Z Garamszegi. 2014. “Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology.”Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology, 287–303. https://doi.org/10.1007/978-3-662-43550-2_11.
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, J D, and S 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.
Halliwell, Ben, Barbara R Holland, and Luke A Yates. 2025. “Multi-Response Phylogenetic Mixed Models: Concepts and Application.”Biological Reviews of the Cambridge Philosophical Society 100 (3): 1294–1316. https://doi.org/10.1111/brv.70001.
Nakagawa, Shinichi, and P de Villemereuil. 2018. “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, Daniel W A Noble, Frank Seebacher, Nicholas C Wu, Samantha Burke, Malgorzata Lagisz, Lisa E Schwanz, Szymon M Drobniak, and Shinichi Nakagawa. 2024. “New Horizons for Comparative Studies and Meta-Analyses.”Trends in Ecology & Evolution 39 (5): 435–45. https://doi.org/10.1016/j.tree.2023.12.004.
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.
Source Code
---title: "BACE: Bayesian Phylogenetic and Correlated Trait Imputation"date: "`r Sys.Date()`"author: "Daniel Noble, Szymek Drobniak, Shinichi Nakagawa"bibliography: ../bib/BACE.bibformat: html: toc: true toc-location: left toc-depth: 3 toc-title: "**Table of Contents**" output-file: "index.html" theme: cosmo embed-resources: true code-fold: show code-tools: true number-sections: true fontsize: "12" code-overflow: wrapcrossref: fig-title: Figure # (default is "Figure") tbl-title: Table # (default is "Table") title-delim: — # (default is ":") fig-prefix: Fig. # (default is "Figure") tbl-prefix: Tab. # (default is "Table")editor_options: chunk_output_type: console---# IntroductionWe 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.# When to use Phylo-BACEGeneral-purpose multiple-imputation tools such as **mice**, **Amelia** and **missForest** assume that rows are exchangeable: 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, and throwing that information away means we impute a missing mouse trait as though the nearest source of information were the grand mean across 500 mammals rather than the handful of other mouse-like rodents in the tree.**Phylo-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 *Interpreting the pooled output* below).Phylo-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 an `MCMCglmm` fit with a `phylo` random effect) and check whether the phylogenetic variance component is non-trivial.- 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).Phylo-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 cheaper.- You need sub-second imputations in a production pipeline. Each chained-equations iteration runs an MCMC fit, so runtimes are measured in minutes, not milliseconds.For a broader discussion of the data-gap problem and imputation in comparative biology we recommend @daliaa.conde2019 and @pottier2024; for Rubin's-rules-based approaches to related problems see @nakagawa2018.# InstallationYou can install the *BACE* package from GitHub using the following command:```{r}#| label: install_bace#| echo: true#| message: false#| warning: false#| eval: true# Install from github and load#install.packages("pacman") pacman::p_load(devtools)devtools::install_github("daniel1noble/BACE") pacman::p_load(BACE, ape, phytools, MCMCglmm, dplyr, magrittr)```If you encounter any issues during installation, please ensure that you have the required dependencies installed.# DocumentationThe main functions in the *BACE* package are `bace()` and `bace_imp()`, which performs Bayesian imputation using chained equations. You can access the documentation for these functions using the following commands:```{r}#| label: bace_imp_doc#| echo: true?bace```As you will see, the function allows you to specify fixed effects, random effects, 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. For most comparative datasets, this will be sufficient for modeling the phylogenetic and species-level variation in the data, but not always.# The BACE workflow at a glance@fig-function-map maps the user-facing functions onto the workflow. `bace()` is the top-level wrapper. Internally it glues together four steps that you can also run individually when you need finer control:1. **Initial chained imputation** — `bace_imp()`. For each variable with missing data, fit a `MCMCglmm` phylogenetic mixed model, impute the missing cells, move on to the next variable, and repeat for `runs` iterations. 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.2. **Convergence assessment of the chained loop** — `assess_convergence()`. Check whether the *sequence of imputed datasets* has stabilised. This is not within-chain MCMC convergence — see *Two kinds of convergence* below.3. **Final imputation runs** — `bace_final_imp()`. Starting from the converged dataset, launch `n_final` *independent* MCMCglmm fits per variable (default `n_final = 50`; parallelisable via `n_cores`). Each run produces one complete imputed dataset and one per-variable posterior. 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 entirely.4. **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-function-map 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 *Downstream analysis with Rubin's rules* below. The rest of this vignette walks through the workflow first with `bace()` and then shows how to run each step individually.{#fig-function-map}# Example Usage## Simulating DataWe'll first simulate a dataset with missing values and a phylogenetic tree. We can do this using the `sim_bace()` function.```{r}#| label: simulate_data#| echo: true# Set seed for reproducibilityset.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)```We 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.```{r}#| label: view_simulated_data#| echo: true# Object contains both the simulated data and phylogeny data <- sim_data$data tree <- sim_data$tree# View first few rows of the datahead(data)```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.## Pre-flight: check phylogenetic signalBefore 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 [@pagel1999] and Blomberg's [@blomberg2003] classical statistics, via `phytools::phylosig()`.- **D** (binary variables only) — Fritz-Purvis D [@fritz2010], via `caper::phylo.d()`.- **Convergence flags** — every fit is checked with `coda::effectiveSize()` and `coda::geweke.diag()`; rows with low effective sample size or failed Geweke tests are flagged `unreliable` and 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 honours BACE's two-RE modes: `species = FALSE` fits a single phylogenetic random effect; `species = TRUE` fits the dual structure (phylogenetic + non-phylogenetic species) that BACE uses when there are within-species replicates.```{r}#| label: phylo_signal_preflight#| echo: true#| eval: false# 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)```Or, more concisely, call `bace()` with `phylo_signal = TRUE` to get the same diagnostic without running the imputation pipeline:```{r}#| label: phylo_signal_via_bace#| echo: true#| eval: false# Signal-only preview: returns a phylo_signal object and does NOT# call bace_imp(), bace_final_imp(), or pool_posteriors().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.**Reviewer's caveat (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.## Using *BACE* for ImputationNow that we have the data and the phylogenetic tree we can run the imputation using the `bace()` or `bace_imp()` function.:::::: {.callout-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 `integer` whose values are all non-negative is imputed as a Poisson count — deliberately so, with no goodness-of-fit "race" against the gaussian family, because a count trait with strong phylogenetic signal is marginally overdispersed and a marginal-fit comparison would silently (and wrongly) pick gaussian exactly when the phylogenetic signal is strongest. 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`). Setting `ovr_categorical = TRUE` instead 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:```{r}#| label: check_variable_types#| echo: true#| eval: false# After running bace() (see next section), inspect the inferred variable typesbace_final1$initial_results$types```Importantly, 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 users, but later we will show how different components of the workflow can be run separately for more control over the imputation process.## Simple BACE WorkflowThe `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.```{r}#| label: simple_bace_workflow1#| echo: true#| eval: true#| warning: false#| message: false# 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)```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 that, 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 sometimes convergence fail yet convergence is likely achieved. If you think this is the case than you can do initial runs and skip the convergence assessment (see below). Note, 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_conv` argument 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_size` argument 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 final `bace` object.- Third, the `species` argument 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_final` argument 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 embarrassingly parallel, so on multi-core machines you can set `n_cores` to run them concurrently.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:```{r}#| label: simple_bace_workflow2#| echo: true#| eval: true#| warning: false# Extract the pooled model for the response variable using get_pooled_model()y_model <-get_pooled_model(bace_final1, variable ="y")summary(y_model)plot(y_model) # 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")# You can also extract all pooled models at onceall_models <-get_pooled_model(bace_final1)names(all_models)```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:```{r}#| label: access_imputed_datasets#| echo: true#| eval: true#| warning: false# 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)head(imputed_datasets[[1]])# Or stack all imputed datasets into one data frame with an .imputation columnimputed_df <-get_imputed_data(bace_final1, format ="data.frame")table(imputed_df$.imputation)```## Interpreting the pooled outputIt is worth being precise about what the pooled posterior returned by `bace()` actually is, because the terminology in the multiple-imputation literature can trip people up.`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 a Monte Carlo 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 @zhou2010 and Chapter 18 of @gelman2013.Two practical notes follow from this:1. **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, and are the right tool when the per-imputation analysis produces frequentist estimates rather than full posterior draws. When you have posterior draws, stacking is both simpler and more faithful — do not apply Rubin's formulas *on top of* the stacked `bace()` posteriors. If your downstream analysis model is frequentist (an `lm`, `glm`, `nlme::gls`, `phylolm`, ...), *BACE* does provide a first-class Rubin's-rules pathway — `with_imputations()` + `pool_mi()` — described in the next section.2. **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:```{r}#| label: check_per_imputation_mcmc#| echo: true#| eval: false# Inspect MCMC diagnostics on the first per-imputation fit for 'y'first_fit <- bace_final1$final_results$all_models[[1]]$yplot(first_fit)coda::effectiveSize(first_fit$Sol)coda::autocorr.diag(first_fit$Sol)```If these look poor, increase `nitt` and/or `burnin` and re-run `bace()` rather than simply trusting the pooled banner.# Downstream analysis with Rubin's rulesThe 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-function-map):1. `with_imputations(object, .f)` applies your model-fitting function `.f` to **each** of the `n_final` imputed datasets and returns the list of fits. It accepts the output of `bace()` or `bace_final_imp()` (or a plain list of completed data frames), shows progress, and captures per-fit errors rather than aborting the whole loop.2. `pool_mi(fits)` combines the per-imputation coefficient estimates with **Rubin's rules** [@rubin1987]: the pooled variance is the within-imputation variance plus (1 + 1/M) times the between-imputation variance, degrees of freedom use the @barnard1999 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 matches `mice::pool` to machine precision on identical inputs.Continuing with the `bace_final1` object fitted above (recall `y` is a count):```{r}#| label: rubin_pathway#| echo: true#| eval: true#| warning: false# 1. Fit a downstream GLM to each imputed datasetfits <-with_imputations(bace_final1,function(d) glm(y ~ x1 + x2, family = poisson, data = d))# 2. Pool with Rubin's rulespool_mi(fits)```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:```{r}#| label: rubin_pathway_gls#| echo: true#| eval: false# Phylogenetic GLS on each imputed dataset, pooled with Rubin's rulesfits_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 with `pool_mi()`; use `pool_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 @rubin1987 and Chapter 2 of @vanbuuren2018.# Advanced: running the workflow step by stepFor 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.## Controlling the imputation models with a list of formulasThe 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.```{r}#| label: bace_imputation#| echo: true #| eval: false# 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):```{r}#| label: bace_imputation_simple#| echo: true# Fit BACE with a single formulabace_impute2 <-bace_imp(fixformula ="y ~ x1 + x2",ran_phylo_form ="~ 1 | Species",phylo = tree,data = data,runs =5, nitt =20000, burnin =5000, thin =10) ```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.## 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:```{r}#| label: evaluate_models#| echo: true# 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 modelsummary(model_y)# Plot trace and density plots for model parametersplot(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.# Phylogenetic and Non-Phylogenetic Random EffectsBACE 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.```{r}#| label: non_phylo_random_effects#| echo: true#| eval: true#| warning: false#| message: false# Simulate data with multiple replicates per speciesset.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 datahead(data2) # We can also make a table of the number of cases per species to confirm that there are multiple replicates per speciesdim(data2)length(unique(data2$Species))```We 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:```{r}#| label: bace_non_phylo_random_effects#| echo: true#| eval: truebace_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) # We can view the final model for y to check that both random effects are includedmodel_y2 <- bace_impute3$models_last_run[["y"]]summary(model_y2)```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:```{r}#| label: bace_non_phylo_random_effects_error#| echo: true#| eval: true#| error: true#| warning: false# 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)```## Adjusting MCMC Parameters For ConvergenceAs 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:```{r}#| label: check_ess#| echo: true#| eval: truemodel_x1 <- bace_impute3$models_last_run[["x1"]]summary(model_x1)model_x2 <- bace_impute3$models_last_run[["x2"]]summary(model_x2)```As 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:```{r}#| label: mcmc_parameters#| echo: true#| eval: true#| warning: false#| message: false# Define MCMC parameters for each model nitt_list <-list(60000, 60000, 105000) thin_list <-list(10, 10, 10)# Run BACE with model-specific MCMC parametersbace_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)```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 to improve mixing, however, it's still good to look at the final models to check the model visually.## Two kinds of convergenceWhen 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.1. **Within-chain MCMC convergence** of each per-variable `MCMCglmm` fit. 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 *individual* `MCMCglmm` object. `bace_imp()` prints a short summary of effective sample sizes and autocorrelation for the last-iteration models when `verbose = TRUE`, and you can always call `plot()` or `coda::effectiveSize()` on `bace_impute$models_last_run[["y"]]` or on any of the per-imputation fits stored in `bace_final1$final_results$all_models[[m]][["y"]]`. Failures here are fixed by increasing `nitt`, `burnin`, or `thin`, tightening priors, or (less often) simplifying the model.2. **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 the `runs` chained iterations. Failures here are fixed by increasing `runs` or, 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 *Interpreting the pooled output* above.## Assessing chained-equations convergenceWe 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:```{r}#| label: check_convergence#| echo: true#| eval: true#| warning: false#| message: false# Assess whether number of runs was enough by looking at the stability of imputed points across runs. Adjust runs if neededconverge <-assess_convergence(bace_impute2, method ="summary")plot(converge)# Now that we know that we can run BACE longer# Fit BACE with a single formulabace_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) # Check again for convergenceconverge <-assess_convergence(bace_impute2_2, method ="summary")plot(converge)```## Worked example: bace() with species random effectsNow 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.```{r}#| label: simple_bace_workflow#| echo: true#| eval: true#| warning: false#| message: false# 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)# Extract the pooled model for y using the accessor functiony_model <-get_pooled_model(bace_final, variable ="y")summary(y_model)plot(y_model) # 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")```# Validating BACE on simulated data with known truthFor methodological work and quantitative claims about imputation accuracy, simulation-based benchmarking is essential. `sim_bace()` generates phylogenetically-structured comparative data with a complete pre-NA dataset returned alongside the masked one, so we can compare BACE's imputed values against the actual held-out truth.The pattern is:1. Simulate a complete dataset with a known phylogenetic signal (`sim_bace()` returns both `data` with NAs and `complete_data` with the truth).2. Mask a fraction of one variable as MCAR (or use the built-in `missingness` argument for type-aware mechanisms).3. Run BACE.4. Compare imputed values at the masked rows against truth using type-appropriate metrics (NRMSE / accuracy / Brier).5. Compare BACE against a **no-skill baseline** — column-mean for continuous, modal class for categorical. BACE has to beat this floor to demonstrate it is using the phylogenetic structure.Below we run a single replicate for one continuous (`gaussian`) and one categorical (`multinomial3`) response. For a manuscript-quality benchmark you should run many replicates per type and report mean ± SD; see `dev/sim_accuracy_multirep.R` in the package source for a full 20-replicate harness across all five response types.```{r}#| label: simulated_benchmark#| echo: true#| eval: falselibrary(BACE)set.seed(2026)# 1. Simulate a complete dataset with phylogenetic signal 0.9sim <-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) # no NAs yet — we mask manually below)truth_data <- sim$complete_datatree <- sim$tree# 2. Mask 20% of y as MCARmiss_idx <-sample.int(nrow(truth_data),size =floor(nrow(truth_data) *0.20))masked <- truth_datamasked$y[miss_idx] <-NA# 3. Run BACEres <-bace(fixformula ="y ~ x1 + x2",ran_phylo_form ="~ 1 | Species",phylo = tree,data = masked,runs =4,n_final =50,nitt =2500, thin =10, burnin =500,skip_conv =TRUE,verbose =FALSE)# 4. Compute BACE imputed-vs-truth accuracytruth_y <- truth_data$y[miss_idx]imp_mat <-sapply(res$imputed_datasets,function(d) as.numeric(d$y[miss_idx])) # missing × n_finalpost_mean <-rowMeans(imp_mat)sd_full <-sd(truth_data$y, na.rm =TRUE)nrmse_bace <-sqrt(mean((post_mean - truth_y)^2)) / sd_full# 95% prediction interval coveragelo <-apply(imp_mat, 1, quantile, probs =0.025)hi <-apply(imp_mat, 1, quantile, probs =0.975)coverage95 <-mean(truth_y >= lo & truth_y <= hi)# 5. Compare against a column-mean baselinebase_pred <-mean(masked$y, na.rm =TRUE)nrmse_base <-sqrt(mean((base_pred - truth_y)^2)) / sd_fullcat(sprintf("BACE NRMSE = %.3f (95%% PI coverage %.2f)\n", nrmse_bace, coverage95))cat(sprintf("baseline NRMSE = %.3f (column-mean imputer)\n", nrmse_base))```Expected output (numbers vary slightly with the seed):```BACE NRMSE = 0.45 (95% PI coverage 0.85)baseline NRMSE = 1.00 (column-mean imputer)```BACE's NRMSE is roughly half the no-skill floor on a high-signal phylogenetic trait, and the empirical 95% prediction interval covers truth ~85–95% of the time at `n_final = 50` (well-calibrated; see `?bace` for why smaller `n_final` undercovers).The same pattern adapts to categorical responses by swapping NRMSE for balanced accuracy and the column-mean baseline for a modal-class baseline:```{r}#| label: simulated_benchmark_cat#| echo: true#| eval: falseset.seed(2026)sim <-sim_bace(response_type ="multinomial3",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))truth_data <- sim$complete_datatree <- sim$treemiss_idx <-sample.int(nrow(truth_data),size =floor(nrow(truth_data) *0.20))masked <- truth_datamasked$y[miss_idx] <-NAres <-bace(fixformula ="y ~ x1 + x2",ran_phylo_form ="~ 1 | Species",phylo = tree, data = masked,runs =4, n_final =50,nitt =2500, thin =10, burnin =500,skip_conv =TRUE, verbose =FALSE)truth_y <-as.character(truth_data$y[miss_idx])imp_mat <-sapply(res$imputed_datasets,function(d) as.character(d$y[miss_idx]))# Modal class per missing cell (BACE's prediction)bace_pred <-apply(imp_mat, 1, function(r)names(sort(table(r), decreasing =TRUE))[1])acc_bace <-mean(bace_pred == truth_y)# Modal class baselinemodal <-names(sort(table(masked$y), decreasing =TRUE))[1]acc_base <-mean(modal == truth_y)cat(sprintf("BACE accuracy = %.2f\n", acc_bace))cat(sprintf("baseline accuracy = %.2f (modal-class imputer)\n", acc_base))```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), and the per-cell 95% prediction intervals are well calibrated:| 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.# Citation and further readingIf you use *BACE* in a publication, please cite both the package (see `citation("BACE")`) and the underlying `MCMCglmm` engine [@hadfield2010]. Phylogenetic comparative methods and phylogenetic mixed models are reviewed in @cornwell2017, @devillemereuil2014a and @halliwell2025. The broader motivation for imputation in comparative biology — data gaps, trait databases, and the cost of complete-case analysis — is discussed in @daliaa.conde2019 and @pottier2024. For a related Rubin's-rules-based approach that handles *phylogenetic* uncertainty (rather than missing trait data), see @nakagawa2018.The Bayesian combiner used by `pool_posteriors()` — concatenating the per-imputation chains to approximate the marginal posterior — is discussed in @zhou2010 and in Chapter 18 of @gelman2013. The Rubin's-rules combiner used by `pool_mi()` is from @rubin1987 with the small-sample degrees-of-freedom correction of @barnard1999; @vanbuuren2018 is the standard modern reference on multiple imputation practice.Bug reports, feature requests and worked examples are very welcome at the [BACE issue tracker](https://github.com/daniel1noble/BACE/issues).# References::: {#refs}:::
