Posterior-tree prediction sensitivity

Run pigauto's full imputation pipeline on each of T posterior phylogenies, generating m_per_tree stochastic completions per tree for a total of T * m_per_tree completed datasets. Each completed dataset is conditional on a specific posterior tree (recorded in mi$tree_index).

Usage

multi_impute_trees(
  traits,
  trees,
  m_per_tree = 1L,
  species_col = NULL,
  trait_types = NULL,
  multi_proportion_groups = NULL,
  log_transform = TRUE,
  missing_frac = 0.25,
  covariates = NULL,
  epochs = 2000L,
  verbose = TRUE,
  seed = NULL,
  share_gnn = TRUE,
  reference_tree = NULL,
  draws_method = c("mc_dropout", "conformal"),
  gnn = FALSE,
  ...
)

Arguments

traits: data.frame. Same format as multi_impute() and impute().
trees: list of phylo objects (class multiPhylo or plain list). Each tree must contain the species in traits as tips. Posterior samples from BirdTree.org (Jetz et al. 2012) are ideal; the bundled trees300 dataset provides 50 posterior trees for avonet300.
m_per_tree: integer. Number of stochastic completion draws per tree (default 1). Total datasets = length(trees) * m_per_tree.
species_col: character or NULL. See impute().
trait_types: named character vector overriding auto-detected trait types. Required for "proportion" and "zi_count". See impute() / preprocess_traits(). Default NULL (auto-detect).
multi_proportion_groups: named list declaring compositional trait groups (rows summing to 1), forwarded to impute() / preprocess_traits(). Default NULL.
log_transform: logical. Log-transform every continuous trait whose observed values are all positive (default TRUE), then back-transform the imputations. Set it to FALSE when traits are already on a log scale (for example log body mass): otherwise they are logged a second time and imputed on a log(log) scale, which can bias analyses done on the log scale. With a mix of raw and already-logged traits, log the raw ones yourself and set FALSE. In a real-data check (PanTHERIA, 4,027 mammal species, four already-logged traits, posterior draws), the default shifted Rubin-pooled regression slopes by up to 3.9 standard errors of the complete-data slope, in a direction that depended on which values were missing; with FALSE the largest shift was 1.8.
missing_frac: numeric. Fraction held out for validation/test during training (default 0.25).
covariates: data.frame or matrix of environmental covariates (fully observed, numeric). Passed through to impute(). Default NULL (no covariates).
epochs: integer. Maximum GNN training epochs per tree (default 2000).
verbose: logical. Print progress (default TRUE).
seed: optional integer. Base random seed; when supplied, each tree uses seed + t - 1 for reproducible results. The default NULL uses the current RNG stream.
share_gnn: logical. With gnn = TRUE, TRUE (the default) trains one GNN on a reference tree and reuses it while recomputing the baseline for each posterior tree. Set FALSE to refit the GNN on each tree. With the default gnn = FALSE, no GNN is trained in either mode. See the "Share-GNN" section below.
reference_tree: optional phylo used for the reference fit when share_gnn = TRUE. Default NULL selects the maximum-clade-credibility tree via phangorn::maxCladeCred(trees). If phangorn is not installed, falls back to trees[[1]] with a warning.
draws_method: character. How the m_per_tree stochastic draws are generated for each tree (P1-11). One of "mc_dropout" (default, back compatible with pre-P1-11 behaviour) or "conformal". See the "Which draw mechanism this uses" section below for the trade-off between the two.
gnn: logical. Passed through to impute() / fit_pigauto() for every tree. When TRUE, the usual GNN correction is trained (once, on the reference tree, when share_gnn = TRUE). When FALSE, no GNN is used anywhere in the pipeline: baseline-only fits throughout (the default since 0.11.0.9001; see impute()). With share_gnn = TRUE, each tree's baseline is fit without held-out splits (splits = NULL), matching the production baseline_full semantics fit_pigauto() uses.
...: additional arguments forwarded to fit_pigauto() via impute().

Value

An object of class "pigauto_mi_trees", inheriting from "pigauto_mi", with components:

datasets: List of T * m_per_tree completed data.frames. Observed cells are preserved; missing cells are filled with stochastic draws for prediction-sensitivity diagnostics.
m: Total number of datasets (T * m_per_tree).
n_trees: Number of posterior trees used.
m_per_tree: Imputations per tree.
tree_index: Integer vector of length m; element i gives the tree index (1..T) for dataset i.
pooled_point: Single data.frame averaging across all T * m_per_tree datasets. For reporting, not inference.
se: Matrix of per-cell pooled SEs (NA if not available).
imputed_mask: Logical matrix; TRUE where a cell was originally missing.
share_gnn: Logical; TRUE if the shared-GNN path was used.
draws_method: Character; "mc_dropout" or "conformal", echoing the argument this call used (P1-11).
mi_workflow: "pigauto_tree_sensitivity_diagnostic", recording that these are prediction-sensitivity completions and cannot be passed to with_imputations() or pool_mi().
fit: Single pigauto_fit fitted on the reference tree when share_gnn = TRUE; NULL otherwise.
fits: List of T pigauto_fit objects (one per tree) when share_gnn = FALSE; NULL when share_gnn = TRUE.
reference_tree: The reference phylo used for the shared fit when share_gnn = TRUE; NULL otherwise.
trees: The input posterior trees.
species_col: Passed-through species column name.

Details

This is an experimental prediction-sensitivity workflow. It compares stochastic completions across a posterior tree sample; it does not provide a validated downstream inferential workflow. See vignette("tree-uncertainty").

Every completed dataset carries a different tree's signal so that between-tree variation enters the returned stochastic completions. This path is not supported by multi_impute_analysis() and has not passed the analysis-aware inferential gate. Do not use it for downstream inference.

For each tree the function runs preprocessing, baseline fitting and prediction. A GNN is also fitted only when gnn = TRUE; with share_gnn = TRUE, that GNN is trained once on a reference tree. Topologies and branch lengths vary across trees, so the phylogenetic baseline covariance differs for each tree.

The returned datasets may be compared descriptively to assess sensitivity of point imputations to the tree sample. No calibrated downstream standard error or Rubin-pooling claim is made for this path.

Computation time. The following historical planning estimates apply to gnn = TRUE, with one shared GNN fit plus T baseline passes when share_gnn = TRUE. They are not timings for the current GNN-off default; estimate the intended tree and sample setting before a campaign.

Species n	1 fit	T = 50 share_gnn=TRUE	T = 50 share_gnn=FALSE
300	~30-60 s	~3-5 min	25-50 min
5,000	~5-10 min	~10-20 min	4-8 hr
10,000	~20-40 min	~30-60 min	17-33 hr

When to use this

Choose the prediction-diagnostic function based on how many trees you have:

One tree (single published phylogeny, single time-calibrated tree): use multi_impute() and select a draw method explicitly when needed.
Multiple posterior trees (BirdTree samples, BEAST posterior, etc.): use multi_impute_trees() only for prediction sensitivity. Tree uncertainty is unsupported by the analysis-aware backend.

When gnn = TRUE, share_gnn = TRUE, the GNN weights and spectral features are trained once on the reference tree (MCC by default). For each posterior tree the BM / joint-MVN baseline is recomputed, and the prediction is the blend (1 - r_cal) * baseline_t + r_cal * gnn_shared. Because r_cal is calibrated once on held-out data at the reference tree and applied uniformly, the tree-uncertainty contribution is:

Fully preserved when the gate is closed (r_cal near 0): the GNN contributes nothing, and the baseline varies per tree.
Partially preserved when the gate is open: the baseline portion still varies, but the GNN portion is a tree-invariant constant — this slightly under-estimates tree variance in the GNN channel.
Lost in the GNN channel when the gate is fully open (rare on real data; the baseline channel still carries tree variation).

On every real dataset benchmarked in the v0.9.0 campaign the gate closed partially or fully. This evidence is specific to those benchmark regimes and does not guarantee tree-variance calibration elsewhere. Set share_gnn = FALSE with gnn = TRUE if you need a separate GNN fit for each tree. With the default gnn = FALSE, there is no GNN channel.

Which draw mechanism this uses (and its limitation)

The default draws_method = "mc_dropout" passes n_imputations = m_per_tree to impute(). With gnn = TRUE, each draw combines a fresh baseline draw with a stochastic GNN pass. With the default gnn = FALSE, no GNN or dropout runs; only the baseline draw varies within each tree.

This differs from the conformal option in multi_impute(). That function now defaults to draws_method = "auto", selecting posterior draws for eligible continuous-trait inputs and conformal draws otherwise. Baseline and optional GNN draws reflect model and prior uncertainty, whereas conformal draws use held-out prediction residuals. Neither mechanism validates downstream pooling or inference for posterior-tree completions.

multi_impute_trees() now exposes a draws_method argument (P1-11). draws_method = "conformal" runs impute with n_imputations = 1 on each tree and draws the m_per_tree completions from the calibrated conformal scores instead, via the same internal sampling helper multi_impute() uses. Default remains "mc_dropout" for back compatibility. Within-tree calibration changes the within-tree stochastic completion mechanism; both settings are prediction-sensitivity diagnostics. Between-tree variation is present under either setting because the baseline (and, when share_gnn = FALSE, the GNN) still varies per tree.

When share_gnn = TRUE with safety_floor = TRUE, the grand-mean baseline mean_baseline_per_col and the three calibrated weights (r_cal_bm, r_cal_gnn, r_cal_mean) are computed ONCE on the reference tree and reused across all posterior trees. They are properties of the observed training traits, not of the tree topology. Each posterior tree only recomputes the BM baseline; the GNN delta and the three weights stay fixed. This preserves prediction sensitivity to the per-tree baseline without re-calibrating the safety floor; it does not establish downstream inferential validity.

References

Nakagawa S, de Villemereuil P (2019). "A general method for simultaneously accounting for phylogenetic and species sampling uncertainty via Rubin's rules in comparative analysis." Systematic Biology 68(4): 632-641.

Jetz W, Thomas GH, Joy JB, Hartmann K, Mooers AO (2012). "The global diversity of birds in space and time." Nature 491(7424): 444-448.

Examples

# \donttest{
data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
trees <- structure(list(tree, tree), class = "multiPhylo")
df <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                c("Mass", "Wing.Length"), drop = FALSE]
rownames(df) <- tree$tip.label
df$Mass[seq_len(3L)] <- NA_real_
mi <- multi_impute_trees(df, trees, m_per_tree = 1L, epochs = 5L,
                         verbose = FALSE)
#> Warning: Few tree-sensitivity draws (T * m_per_tree = 2). Descriptive sensitivity summaries may be unstable.
#> Warning: phangorn required for MCC tree selection; falling back to trees[[1]]. Install with install.packages('phangorn') for MCC-based selection.
#> Warning: Small validation set for 2 trait(s): Mass (n=3), Wing.Length (n=4). Calibrated gate and conformal scores will be noisy for these trait(s). 95% split-conformal coverage is NOT achievable for 2 trait(s) with fewer than 19 validation cells (Mass (n=3), Wing.Length (n=4)): the achievable ceiling is n_val / (n_val + 1), which only reaches 0.95 at n_val >= 19. See `?fit_pigauto` under 'Calibration at small n' for smoothing options.
#> Warning: phangorn required for MCC tree selection; falling back to trees[[1]]. Install with install.packages('phangorn') for MCC-based selection.
print(mi)
#> pigauto experimental posterior-tree prediction sensitivity
#>   Trees     : 2 posterior phylogenies
#>   Per tree  : 1 completion draws (MC dropout)
#>   Total     : 2 completed datasets
#>   Species   : 30
#>   Traits    : 2 -- Mass, Wing.Length
#>   Cells     : 3 imputed / 60 total (5.0%)
#> 
#>   Access diagnostic draws:  mi$datasets[[i]]
#>   Tree index for draw i:    mi$tree_index[i]
#>   Downstream inference:     unsupported for this path
mass_by_tree <- vapply(mi$datasets, function(dat) dat$Mass,
                       numeric(nrow(df)))
apply(mass_by_tree, 1L, stats::sd)
#>  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# }