Package {pigauto}

AVONET morphological and ecological trait data for 300 bird species

Description

A data frame with 300 rows and 8 columns: 4 continuous morphometric traits and 3 ecological traits (2 categorical, 1 ordinal) for demonstrating mixed-type imputation. Species names are in the Species_Key column. Set rownames(df) <- df$Species_Key; df$Species_Key <- NULL before passing to preprocess_traits.

Usage

avonet300

Format

A data frame with 300 rows and 8 variables:

Species_Key

Character. Species name in BirdTree format (spaces replaced by underscores).

Mass

Numeric. Body mass (grams).

Beak.Length_Culmen

Numeric. Beak length from culmen (mm).

Tarsus.Length

Numeric. Tarsus length (mm).

Wing.Length

Numeric. Wing length (mm).

Trophic.Level

Factor. Dietary category: Carnivore, Herbivore, Omnivore, or Scavenger.

Primary.Lifestyle

Factor. Lifestyle category: Aerial, Aquatic, Generalist, Insessorial, or Terrestrial.

Migration

Ordered factor. Migration strategy: Resident < Partial < Full.

Source

Tobias et al. (2022) AVONET: morphological, ecological and geographical data for all birds. Ecology Letters, 25, 581-597. Species labels are aligned to the bundled example phylogenies.

Full AVONET morphological and ecological trait data for 9,993 bird species

Description

The full-scale counterpart to avonet300: every bird species for which AVONET3 and the bundled example phylogeny agree on both a species label and a complete set of continuous morphometric measurements. The schema is identical to avonet300 (same trait columns, same factor encodings, same Species_Key column) so any code that runs on avonet300 runs on avonet_full with no modification.

Usage

avonet_full

Format

A data frame with 9,993 rows and 8 variables. Columns match avonet300: Species_Key, Mass, Beak.Length_Culmen, Tarsus.Length, Wing.Length, Trophic.Level, Primary.Lifestyle, Migration.

Details

Unlike the 300-species subset, native AVONET missingness is PRESERVED in the two ecological columns: 4 NAs in Trophic.Level and 20 NAs in Migration. The four continuous columns are complete by construction.

Use avonet_full + tree_full to exercise pigauto at real-world scale (the AVONET missingness sweep benchmark uses this dataset). For quick examples or unit tests, prefer the 300-species avonet300 subset – it is ~30x smaller and loads instantly.

See avonet300 for the full column schema (same traits and encodings).

Source

Tobias et al. (2022) AVONET: morphological, ecological and geographical data for all birds. Ecology Letters, 25, 581-597. Species labels are aligned to the bundled example phylogeny.

See Also

avonet300, tree_full

Build a phylogenetic graph representation from a tree

Description

Computes a symmetric-normalised Gaussian kernel adjacency matrix and Laplacian spectral node features from the cophenetic distance matrix. Results are optionally cached to an .rds file.

Usage

build_phylo_graph(tree, k_eigen = "auto", sigma_mult = 0.5, cache_path = NULL)

Arguments

Details

The graph is built in three steps: (1) cophenetic distances between all pairs of tips; (2) Gaussian kernel A_{ij} = \exp(-d_{ij}^2 / (2\sigma^2)) with \sigma = \mathrm{median}(D) \times s (where s is the user-supplied sigma_mult); (3) symmetric normalisation \tilde{A} = D^{-1/2} A D^{-1/2} with self-loops added before normalisation.

Spectral node features are the k_eigen smallest non-zero eigenvectors of the unnormalised Laplacian, which encode the broad cluster structure of the phylogeny.

Value

A list with:

adj

Numeric matrix (n x n). Symmetric-normalised adjacency.

coords

Numeric matrix (n x k_eigen). Spectral node features.

n

Integer. Number of tips.

sigma

Numeric. Bandwidth used.

D

Numeric matrix (n x n). Cophenetic (patristic) distances between tips, row/column-ordered by tree$tip.label. Returned so downstream functions can reuse it instead of calling ape::cophenetic.phylo() again.

D_sq

Numeric matrix (n x n). Element-wise squared cophenetic distances (D * D). Used by GraphTransformerBlock for the learnable per-head Gaussian bandwidth (B2 rate-aware attention). Not freed during impute() training — only D is freed.

R_phy

Numeric matrix (n x n). Phylogenetic correlation matrix cov2cor(ape::vcv(tree)) (diagonal = 1). Used by fit_baseline for BM conditional imputation.

Scaling

For trees with more than 7500 tips, build_phylo_graph() uses RSpectra's sparse Lanczos eigensolver (eigs_sym) instead of the dense base::eigen(). This reduces the spectral step from O(n^3) (hours at n = 10{,}000) to O(n \cdot k \cdot \mathrm{iters}) (seconds to a minute). The threshold is conservative because the dense eigensolver is competitive on small and mid-size trees, and on pathological ultrametric simulations (ape::rcoal) the Laplacian spectrum has tight degenerate clusters that slow down Lanczos until the tree is fairly large. On realistic asymmetric phylogenies with variable branch lengths (e.g. the AVONET BirdTree Stage2 Hackett tree) the sparse path wins much earlier, and at n = 9{,}993 it delivers a ~10x speedup over dense. The dense path is kept as the default for smaller trees and as a fallback when RSpectra is not installed, or if Lanczos fails to converge. The cophenetic distance matrix is computed once per call and reused for both the adjacency and the spectral features; it is also returned in the result so that downstream consumers (notably fit_baseline) can skip recomputation.

Examples

set.seed(1)
tree <- ape::rtree(30)
g <- build_phylo_graph(tree, k_eigen = 4)
dim(g$adj)    # 30 x 30
dim(g$coords) # 30 x 4
dim(g$D)      # 30 x 30

Compute calibration data for probability predictions

Description

Bins predicted probabilities and computes observed frequencies within each bin. Useful for calibration plots of binary classifiers.

Usage

calibration_df(truth, prob, n_bins = 10L)

Arguments

Value

A data.frame with columns bin_mid, obs_freq, mean_pred, n.

Examples

calibration_df(
  truth = c(0, 0, 1, 1),
  prob = c(0.1, 0.3, 0.7, 0.9),
  n_bins = 2
)

Check whether traits and a tree are ready for pigauto

Description

Inspect declared inputs without constructing a phylogenetic graph or fitting a model. The check is intended to make input problems and likely compute scale visible before calling impute(). Its fingerprint identifies the declared input object for comparison within a workflow; it is not a security hash.

Usage

check_pigauto(
  traits,
  tree,
  species_col = NULL,
  trait_types = NULL,
  multi_proportion_groups = NULL,
  covariates = NULL,
  gnn = FALSE
)

Arguments

Value

A list of class pigauto_check containing a schema version, status, stable input fingerprint, input summaries, and structured messages.

Examples

tree <- ape::read.tree(text = "((a:1,b:1):1,c:2);")
traits <- data.frame(x = c(1, NA, 3), row.names = tree$tip.label)
check_pigauto(traits, tree)

Compare BM baseline and pigauto methods across replicates

Description

Runs a full comparison of the BM baseline and pigauto (with attention + calibration) over multiple random seeds. Returns a tidy data.frame suitable for plotting or downstream analysis.

Usage

compare_methods(
  data,
  tree,
  splits = NULL,
  seeds = NULL,
  epochs = 500L,
  verbose = TRUE,
  gnn = TRUE,
  ...
)

Arguments

Details

For each seed the function:

1. Creates train/val/test splits.

2. Fits the phylogenetic BM baseline.

3. Fits pigauto (with attention and calibration enabled by default).

4. Evaluates both methods on the test split.

5. Collects results into a single data.frame.

The default comparison trains the optional torch GNN. The example checks that torch can create a tensor before fitting; it reports a skip when the Lantern runtime is unavailable.

Value

A data.frame with columns: method, trait, type, metric, value, rep.

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
torch_ready <- requireNamespace("torch", quietly = TRUE) &&
  isTRUE(tryCatch(torch::torch_is_installed(), error = function(e) FALSE)) &&
  isTRUE(tryCatch({
    probe <- torch::torch_tensor(1)
    as.numeric(probe$item()) == 1
  }, error = function(e) FALSE))
if (torch_ready) {
  cmp <- compare_methods(preprocess_traits(traits, tree), tree, seeds = 1L,
                         epochs = 5L, verbose = FALSE, gnn = TRUE)
  aggregate(value ~ method + trait + metric, data = cmp, FUN = mean)
} else {
  message("Skipping GNN comparison: install torch and its Lantern runtime.")
}

Extract completed trait data

Description

Extract completed trait data

Usage

completed_data(result)

Arguments

Value

Exactly the stored result$completed data.frame.

Examples

## Not run: completed_data(result)

Compute a confusion matrix for categorical or binary predictions

Description

Compute a confusion matrix for categorical or binary predictions

Usage

confusion_matrix(truth, predicted, levels = NULL)

Arguments

Value

A list with:

table

Confusion matrix (rows = truth, columns = predicted).

accuracy

Overall accuracy.

per_class

Data.frame with per-class precision, recall, F1.

Examples

confusion_matrix(
  truth = c("a", "a", "b", "b"),
  predicted = c("a", "b", "b", "b")
)

k-fold cross-validation for pigauto trait imputation

Description

Performs stratified k-fold cross-validation by rotating which cells serve as the test set. Returns per-fold and aggregated metrics.

Usage

cross_validate(
  data,
  tree,
  k = 5L,
  seeds = NULL,
  epochs = 500L,
  verbose = TRUE,
  ...
)

Arguments

Value

A list of class "pigauto_cv" with:

results

Data.frame with columns: fold, rep, trait, type, metric, value.

summary

Data.frame with mean and sd across folds/reps for each trait + metric.

conformal_coverage

Data.frame of coverage per trait across folds (if available).

k

Number of folds.

n_reps

Number of replicates.

Examples

data(avonet300, tree300, package = "pigauto")
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(12L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
cv <- cross_validate(preprocess_traits(traits, tree), tree, k = 3L,
                     seeds = 1L, epochs = 1L, verbose = FALSE,
                     k_eigen = 4L, hidden_dim = 8L,
                     use_transformer_blocks = FALSE)
print(cv)
summary(cv)

Simulated multi-observation-per-species CTmax data

Description

A simulated dataset mimicking a thermal tolerance study where each species has multiple measurements of critical thermal maximum (CTmax) taken at different acclimation temperatures. The data-generating process is

Usage

ctmax_sim

Format

A data frame with 1,464 rows and 3 variables:

species

Character. Species name matching tree300 tips.

acclim_temp

Numeric. Acclimation temperature (degrees C). This is an observation-level covariate that varies within species.

CTmax

Numeric. Critical thermal maximum (degrees C). Contains NA values for unobserved species and MCAR missingness.

Details

CTmax_{ij} = 38 + phylo_i + 0.50 \times acclim\_temp_{ij} + \epsilon_{ij}

where phylo_i follows Brownian motion on tree300, the within-species acclimation response ratio is 0.50, and \epsilon_{ij} \sim N(0, 1.5). Thirty percent of species are entirely unobserved (all CTmax values are NA), and an additional 15\ observations are missing at random.

This dataset demonstrates pigauto's multi-observation and observation-level covariate support. Use with tree300 and set species_col = "species" in impute.

Source

Simulated. See data-raw/make_ctmax_sim.R for the generation script.

See Also

tree300, impute

Evaluate a fitted pigauto model on its test set

Description

Computes type-specific performance metrics for each trait. By default the model is evaluated on the test split stored in the fit object, but an alternative data / splits can be supplied.

Usage

evaluate(fit, data = NULL, splits = NULL)

Arguments

Details

Metrics by trait type:

continuous

RMSE, Pearson r, MAE

count

RMSE, MAE, Pearson r

binary

Accuracy, Brier score

categorical

Accuracy (overall)

ordinal

RMSE (on integer scale), Spearman rho

When conformal scores are present in the fit, conformal coverage at the 95\ ordinal traits.

When the fit includes a baseline, baseline metrics are appended with method = "baseline" for direct comparison.

Value

A data.frame with columns: method, trait, type, metric, value, n_test.

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
pd <- preprocess_traits(traits, tree)
splits <- make_missing_splits(pd$X_scaled, trait_map = pd$trait_map)
fit <- fit_pigauto(pd, tree, splits = splits, epochs = 5L,
                   verbose = FALSE)
eval_df <- evaluate(fit)
eval_df[eval_df$metric == "rmse", ]

Evaluate imputation performance against known values

Description

Computes type-specific metrics for each trait on the validation and test splits. When a trait_map is supplied, metrics are dispatched per trait type; otherwise the function falls back to continuous-only metrics (RMSE, Pearson r, 95% coverage).

Usage

evaluate_imputation(pred, truth, splits, pred_se = NULL, trait_map = NULL)

Arguments

Details

Metrics by trait type:

continuous

RMSE, Pearson r, 95\ supplied)

proportion

RMSE, Pearson r, 95\ supplied)

count

RMSE, MAE, Pearson r

ordinal

RMSE, Spearman rho

binary

Accuracy, Brier score

categorical

Accuracy

zi_count

RMSE, MAE, Pearson r, zero-accuracy, Brier score (on gate)

For binary and categorical traits the function accepts either a pigauto_pred object (preferred, gives access to probabilities) or raw matrices (latent scale).

Value

A data.frame with columns split, trait, type, n, and type-specific metric columns.

Examples

truth <- matrix(c(1, 2, 3, 4), nrow = 2L)
splits <- list(test_idx = 1L, val_idx = 2L, n = 2L, p = 2L)
evaluate_imputation(truth, truth, splits)

Fit the phylogenetic baseline

Description

Dispatches to pigauto's phylogenetic baseline machinery and returns imputed latent-scale means and standard errors for every species.

Usage

fit_baseline(
  data,
  tree,
  splits = NULL,
  model = "BM",
  graph = NULL,
  multi_obs_aggregation = c("hard", "soft"),
  lambda_mode = c("estimate", "fixed_1", "cv", "bayes"),
  lambda_fixed = NULL,
  discrete_lambda = c("estimate", "fixed_1"),
  em_iterations = 0L,
  em_tol = 0.001,
  em_offdiag = FALSE,
  joint_solver = c("inhouse", "rphylopars"),
  predict_method = c("auto", "exact", "per_column"),
  joint_refine_iter = 0L,
  predict_route = NULL,
  seed = NULL
)

Arguments

Details

When splits is supplied the val and test cells are masked to NA before fitting, so the baseline is evaluated under the same conditions as fit_pigauto.

Continuous-family columns use Brownian-motion conditional MVN baselines on the phylogenetic correlation matrix, either independently or through the joint MVN path when the data and optional dependencies support it. Binary, ordinal, categorical, and zero-inflated gate columns use the appropriate label-propagation or threshold/liability baseline candidates, with per-column fallbacks when a joint path is not available.

Covariates and the joint baseline (P1-8): data$covariates is only used by the per-column BM path (bm_impute_col_with_cov()). The joint MVN and threshold-joint (Rphylopars) baselines do not accept a covariate design matrix, so when a joint path is selected (BM-eligible columns >= 2, or binary/ordinal cols present, with Rphylopars available) any supplied covariates are ignored for the BASELINE and a warning is emitted; covariates still reach the GNN correction via fit_pigauto regardless of which baseline path fires.

Per-type lambda dispatch (arc/lambda-per-type; corrected in S4/S6, feat/joint-lambda-default): lambda_mode only ever governs the baseline for CONTINUOUS-FAMILY columns (continuous, count, proportion, zi_count magnitude) – NOT ordinal. Binary, ordinal, categorical, and zero-inflated gate columns always stay at lambda = 1 in EVERY path (threshold-joint, OVR-categorical, label propagation, and the per-trait ordinal path-selection fallback below), regardless of lambda_mode – there is no discrete-trait analogue of Pagel's lambda, and previously forcing these columns onto label propagation any time lambda_mode != "fixed_1" cost 19pp of Trophic.Level accuracy on AVONET (0.789 -> 0.600; see docs/dev-log/2026-08-16-external-comparison-results.md). When lambda_mode != "fixed_1" and the threshold-joint baseline fires for a dataset with binary/ordinal AND continuous-family columns, the joint liability fit's own continuous-column OUTPUT is now USED (each continuous-family column gets its own lambda_k via the joint solver's lambda_cols), rather than being discarded for a separate lambda-aware per-column BM re-fit as in the pre-S4 design. Ordinal liability columns inside that same joint fit, and the BM-via-MVN / K-class-OVR-LP alternatives the "Per-trait ordinal path selection" block below compares against it, all stay at lambda = 1.

Value

A list with:

mu

Numeric matrix (n_species x p_latent), baseline means in latent scale.

se

Numeric matrix (n_species x p_latent), standard errors.

path

Named character vector, one entry per trait_map entry (name = trait/group name), recording which dispatch branch produced that trait's baseline: "joint_mvn" (Phase 2 continuous-only joint), "threshold_joint" (Phase 3/5 binary + ordinal + continuous-family liability joint), "ovr_categorical" (Phase 6 one-vs-rest K-fits), "per_column_bm" (per-column Brownian-motion conditional MVN, including the lambda-aware path and the per-trait ordinal fallback that beats threshold-joint on held-out val MSE), "multi_proportion_bm" (per-component per_column_bm on the CLR columns of a multi_proportion group – same kernel as "per_column_bm", labelled separately because multi_proportion is architecturally distinct), "label_propagation" (phylogenetic label propagation for binary/categorical/ordinal, and the per-trait ordinal fallback that beats the other two candidates), "zi_gate_lp" (label propagation for a zero-inflated count's gate column), or "zi_mag_constant" (global mean/sd fallback for a zero-inflated count's magnitude column when fewer than 5 non-zero observations exist). For zi_count traits (2 latent columns: gate then magnitude), path reports the GATE column's dispatch – the magnitude column can independently land on "joint_mvn", "threshold_joint", "per_column_bm", or "zi_mag_constant" and is not separately surfaced here.

lambda_per_trait

Named numeric vector (length p_latent, names = colnames(data$X_scaled)), the lambda actually used for each latent column: 1 for every column not eligible for lambda estimation (all discrete columns; continuous-family columns under lambda_mode %in% c("cv", "bayes")), else the estimated or fixed value.

lambda_block

Numeric scalar, the shared lambda used internally by whichever joint fit ran for its Sigma M-step / opt-in exact conditional / opt-in EM refine (see R/joint_mvn_solver.R's $lambda_block); NA when no joint fit ran.

lambda_mode

Character, echoes the resolved lambda_mode argument.

predict_method_used

Character scalar, "exact", "per_column", or "auto". Under a concrete predict_method ("exact" / "per_column"), aggregated across every joint fit that ran: "exact" only if every joint fit achieved exact; "per_column" if any fell back, if none ran (a pure per-column baseline), or if every joint fit used joint_solver = "rphylopars" (no exact/per_column dichotomy there). Under predict_method = "auto" (or a forced predict_route), always "auto"; see $predict_method_by_trait for the resolved per-trait routes.

predict_method_by_trait

Named character vector (names = trait names), one entry per trait_map entry, giving the route ("exact" or "per_column") actually used for that trait's cells. Under a concrete predict_method, every entry equals predict_method_used. Under "auto" (or a forced predict_route), the per-trait validation choice (see predict_method's "auto" entry above).

Examples

data(avonet300, tree300, package = "pigauto")
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
pd     <- preprocess_traits(traits, tree)
splits <- make_missing_splits(pd$X_scaled, trait_map = pd$trait_map)
bl     <- fit_baseline(pd, tree, splits)

Fit a pigauto model for trait imputation

Description

Trains a pigauto model: a gated ensemble of a phylogenetic baseline and an attention-based graph neural network correction, implemented as an internal torch module (ResidualPhyloDAE; "Residual" here refers to the ResNet-style skip connections in the GNN layers, not to a statistical residual). For continuous, count, and ordinal traits the baseline is Brownian motion (phylogenetic correlation matrix); for binary and categorical traits it is phylogenetic label propagation. Supports all eight trait types via a unified latent space.

Usage

fit_pigauto(
  data,
  tree,
  splits = NULL,
  graph = NULL,
  baseline = NULL,
  gnn = FALSE,
  baseline_full = NULL,
  hidden_dim = 64L,
  k_eigen = "auto",
  n_gnn_layers = 2L,
  gate_cap = 0.8,
  use_attention = TRUE,
  use_transformer_blocks = TRUE,
  n_heads = 4L,
  ffn_mult = 4L,
  use_trait_attention = FALSE,
  n_trait_heads = 2L,
  trait_embed_dim = 32L,
  dropout = 0.1,
  lr = 0.003,
  weight_decay = 1e-04,
  epochs = 3000L,
  corruption_rate = 0.55,
  corruption_start = 0.2,
  corruption_ramp = 500L,
  refine_steps = 8L,
  lambda_shrink = 0.03,
  lambda_gate = 0.01,
  warmup_epochs = 200L,
  edge_dropout = 0.1,
  eval_every = 100L,
  patience = 10L,
  clip_norm = 1,
  conformal_method = c("split", "bootstrap", "mondrian"),
  conformal_bootstrap_B = 500L,
  conformal_split_val = FALSE,
  gate_method = c("cv_folds", "median_splits", "single_split"),
  gate_splits_B = 31L,
  gate_cv_folds = 5L,
  safety_floor = FALSE,
  phylo_signal_gate = FALSE,
  phylo_signal_threshold = 0.2,
  phylo_signal_method = c("lambda", "blomberg_k"),
  min_val_cells = 20L,
  lambda_mode = c("estimate", "fixed_1", "cv", "bayes"),
  discrete_lambda = c("estimate", "fixed_1"),
  joint_solver = c("inhouse", "rphylopars"),
  predict_method = c("auto", "exact", "per_column"),
  joint_refine_iter = 0L,
  verbose = TRUE,
  seed = NULL
)

Arguments

Details

Blend formulation: The prediction is \hat{x} = (1-r)\mu + r\delta, where \mu is the BM baseline, \delta is the model's direct prediction, and r = \sigma(\rho) \times \mathrm{cap} is a per-column learnable gate bounded in (0, \mathrm{gate\_cap}). When r = 0, the prediction collapses to the baseline. The gate is regularised toward zero via the shrinkage penalty on \delta - \mu, so the model defaults to the baseline unless the GNN's correction demonstrably helps on the validation set.

Training objective (per epoch):

1. A random subset of observed cells is corrupted with a learnable mask token.

2. The model predicts \delta from graph context.

3. Loss = type-specific reconstruction on corrupted cells + lambda_shrink * MSE(\delta - \mu) + lambda_gate * MSE(r).

The gate penalty on r is necessary because when \delta = \mu (the BM-optimal solution for observed cells), the reconstruction and shrinkage losses both equal zero regardless of r, leaving no gradient to close the gate. The explicit penalty ensures gates default toward zero when the GNN correction provides no benefit.

Type-specific losses:

continuous/count/ordinal

MSE

binary

BCE with logits

categorical

cross-entropy over K latent columns

Value

An object of class "pigauto_fit".

Calibration at small n

Conformal bounds use the empirical (1 - alpha) quantile of |y - \hat y| on held-out validation cells. They are nominal held-out diagnostics and do not certify package-wide coverage.

Small validation sets make both the residual quantile and gate selection unstable. The min_val_cells warning identifies traits where this is especially relevant. conformal_split_val = TRUE can separate gate selection from residual estimation for eligible columns, but it does not remove the small-sample limitation. Use cross_validate() on the intended data and report the resulting held-out diagnostics rather than treating an interval target as certified.

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
data <- preprocess_traits(traits, tree)
splits <- make_missing_splits(data$X_scaled, trait_map = data$trait_map,
                              seed = 1)
fit <- fit_pigauto(data, tree, splits = splits, epochs = 5L,
                   verbose = FALSE, seed = 1)

Impute missing phylogenetic traits (convenience wrapper)

Description

One-call interface for preprocessing, phylogenetic baseline fitting, and prediction. The default gnn = FALSE uses the baseline alone. Set gnn = TRUE to also train a GNN and calibrate its contribution. For fine-grained control, use the individual functions (preprocess_traits, fit_baseline, fit_pigauto, etc.) directly.

Usage

impute(
  traits,
  tree,
  species_col = NULL,
  trait_types = NULL,
  multi_proportion_groups = NULL,
  k_eigen = "auto",
  log_transform = TRUE,
  missing_frac = 0.25,
  n_imputations = 1L,
  covariates = NULL,
  epochs = 2000L,
  verbose = TRUE,
  seed = NULL,
  multi_obs_aggregation = c("hard", "soft"),
  lambda_mode = c("estimate", "fixed_1", "cv", "bayes"),
  discrete_lambda = c("estimate", "fixed_1"),
  joint_solver = c("inhouse", "rphylopars"),
  predict_method = c("auto", "exact", "per_column"),
  joint_refine_iter = 0L,
  em_iterations = 0L,
  em_tol = 0.001,
  em_offdiag = FALSE,
  pool_method = c("median", "mean", "mode"),
  clamp_outliers = FALSE,
  clamp_factor = 5,
  match_observed = c("none", "pmm"),
  pmm_K = 5L,
  safety_floor = FALSE,
  phylo_signal_gate = FALSE,
  phylo_signal_threshold = 0.2,
  phylo_signal_method = "lambda",
  conformal_split_val = FALSE,
  gnn = FALSE,
  ...
)

Arguments

Value

An object of class "pigauto_result" with components:

completed

The input traits data.frame with observed values preserved and only missing cells filled in. This is the primary output – typically what users want.

imputed_mask

Logical matrix (same shape as completed) that is TRUE for cells that were imputed (originally NA) and FALSE for observed cells.

prediction

A pigauto_pred object from predict.pigauto_fit containing raw model predictions for every cell (observed + missing), per-cell uncertainty ($se – see predict.pigauto_fit's "What $se means per trait type" section; it is not uniformly a standard error), class probabilities, and conformal intervals.

fit

The trained pigauto_fit object.

baseline

The phylogenetic baseline.

data

The preprocessed pigauto_data object.

splits

The val/test splits (or NULL if missing_frac = 0).

evaluation

Evaluation metrics on test set (or NULL).

check

A compact pigauto_check preflight summary.

Trait type auto-detection

pigauto infers each trait's type from its R class — no trait_types argument is needed for most data:

Two types cannot be inferred from class alone and must be declared via trait_types:

"proportion"

A numeric bounded 0–1, e.g. survival rate: trait_types = c(Survival = "proportion").

"zi_count"

An integer with excess zeros, e.g. parasite count: trait_types = c(Parasites = "zi_count").

Use the trait_types argument directly (it is an explicit parameter, not a ... pass-through).

Traits vs covariates

The distinction is functional, not ontological: a trait is something you want to impute (NA values allowed in traits); a covariate is something you use to sharpen imputation accuracy (must be fully observed, passed via covariates). The same variable can be either depending on the scientific question.

Examples:

• IUCN status with Data Deficient species → put it in traits as ordered(c("LC","NT","VU","EN","CR")) so pigauto predicts the unknown categories.

• IUCN status fully known for all species → pass as a covariate to inform imputation of other traits (e.g. body mass, range size).

• Realm / biome (factor) → pass as a covariate; pigauto one-hot encodes factor columns automatically (v0.6.1+).

Variables that belong in traits: anything with missing values you care about predicting. Variables that belong in covariates: fully observed, exogenous to the trait space (geography, climate, habitat, experimental treatment).

Safety floor (v0.9.1.9002+)

With safety_floor = TRUE (opt-in since 0.11.0.9002; it was the default from v0.9.1.9002), the post-training calibration grid searches a 3-way convex combination of the Brownian-motion baseline, the GNN delta, and the per-trait grand mean. The simplex is sampled at step 0.05 (231 candidates per latent column). Because the corner (0, 0, 1) — pure grand mean — is always in the grid, the selected candidate cannot be worse than the grand-mean corner on the validation cells under the calibration metric. The fit object gains four new slots: r_cal_bm, r_cal_gnn, r_cal_mean (each a named numeric of length p_latent), and mean_baseline_per_col.

Set safety_floor = FALSE to reproduce the pre-v0.9.1.9002 1-D calibration bit-identically (no mean term; r_cal_mean = 0; r_cal_bm = 1 - r_cal_gnn). See specs/2026-04-23-safety-floor-mean-gate-design.md for the design rationale and plans/2026-04-23-safety-floor-mean-gate.md for the implementation plan.

What gets imputed (read this first)

pigauto only imputes cells that are NA in the input. Observed cells are preserved as-is in result$completed. The slot result$prediction$imputed contains the model's prediction for every cell – observed and missing alike – and is intended for diagnostics (e.g. checking calibration on training cells), not as the imputed-values output. The imputed values themselves are result$completed[result$imputed_mask].

Common pitfall. If you call impute() on a fully observed trait matrix (no NAs among matched target rows), there is nothing to impute. impute() stops before graph construction or training and directs you to cross_validate() for predictive assessment. To exercise imputation, supply a real missing target cell.

Imbalanced K-class traits. At default settings (n_imputations = 1L, pool_method = "median"), a small ordinal / categorical trait whose marginal distribution is heavily skewed (e.g. AVONET Migration is ~78\ ~14\ collapse onto a corner that predicts the majority class everywhere. Compare candidate settings, including pool_method = "mode", with cross_validate() on the intended trait table rather than applying a generic setting.

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
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_
result <- impute(df, tree, epochs = 5L, verbose = FALSE)
head(result$completed)

Load a saved pigauto model

Description

Reads a pigauto_fit object previously saved with save_pigauto.

Usage

load_pigauto(path)

Arguments

Value

An object of class "pigauto_fit".

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
traits$Mass[seq_len(3L)] <- NA_real_
fit <- impute(traits, tree, epochs = 5L, verbose = FALSE)$fit
path <- tempfile(fileext = ".pigauto")
save_pigauto(fit, path)
pred <- predict(load_pigauto(path))

Split cells into train/val/test for imputation evaluation

Description

Randomly designates a fraction of cells as "missing" and splits them into validation and test sets. When a trait_map is supplied, masking operates at the original trait level – all latent columns belonging to one trait are held out together (important for categorical traits).

Usage

make_missing_splits(
  X,
  missing_frac = 0.25,
  val_frac = 0.25,
  seed = NULL,
  trait_map = NULL,
  mechanism = c("MCAR", "MAR_trait", "MAR_phylo", "MNAR"),
  mechanism_args = list(),
  tree = NULL
)

Arguments

Details

The returned index vectors use linear (column-major) indexing. Both original-trait-space and latent-space indices are returned when a trait_map is present.

Value

A list with:

val_idx

Integer vector of linear indices (latent space).

test_idx

Integer vector of linear indices (latent space).

val_idx_trait

Integer vector in original-trait space (if trait_map supplied).

test_idx_trait

Integer vector in original-trait space (if trait_map supplied).

n

Number of species (rows).

p

Number of latent columns.

n_traits

Number of original traits.

mask

Logical matrix (n x p_latent). TRUE = observed.

mechanism

Character string of the mechanism used.

Examples

X <- matrix(rnorm(100), nrow = 20)
splits <- make_missing_splits(X, missing_frac = 0.25, seed = 1)
length(splits$val_idx)

# MAR: missingness depends on another trait
splits_mar <- make_missing_splits(X, mechanism = "MAR_trait",
  mechanism_args = list(driver_col = 1, beta = 2))

Create an observed/missing mask matrix

Description

Returns a logical matrix of the same dimensions as X, with TRUE where values are observed (not NA).

Usage

mask_missing(X)

Arguments

Value

Logical matrix, TRUE = observed.

Examples

X <- matrix(c(1, NA, 3, 4), nrow = 2)
mask_missing(X)

Generate experimental stochastic completion datasets

Description

Run pigauto's full imputation pipeline and return M stochastic completions of the trait matrix instead of a single point estimate. The conformal-width and Brownian/MC-dropout draws returned here are experimental prediction-diagnostic draws. Do not use these datasets for downstream inference or Rubin pooling. For continuous traits, draws_method = "posterior" instead returns proper Bayesian posterior imputations for the analyses described in "Posterior draws" below: fit the analysis to them with with_imputations() and pool the resulting fits with pool_mi(). The separate analysis-aware backend, multi_impute_analysis(), covers its own documented narrow regime.

Usage

multi_impute(
  traits,
  tree,
  m = 100L,
  draws_method = c("auto", "conformal", "mc_dropout", "posterior"),
  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,
  gnn = FALSE,
  lambda_mode = c("estimate", "fixed_1", "cv", "bayes"),
  posterior_control = list(),
  ...
)

Arguments

Details

With the default draws_method = "auto", eligible continuous-trait data use posterior draws; other inputs use conformal prediction-diagnostic draws. Check result$draws_method before downstream fitting.

The conformal and MC-dropout draws do not condition on a declared substantive analysis model. Consequently, stochastic variation alone does not make them proper or congenial multiple imputations. The analysis-aware backend requires the analysis model before generating draws and dispatches only across its documented supported model classes.

draws_method = "conformal": Run the model once; missing cells are sampled from x_{ij}^{(k)} \sim \mathrm{N}(\hat\mu_{ij},\; (q_{j}/1.96)^2) where q_j is the trait-level split-conformal residual quantile. Dividing this quantile by 1.96 is a pragmatic Normal-scale construction, not an inference consequence of a nominal held-out conformal diagnostic. For discrete traits (binary, categorical) it uses Bernoulli / categorical draws from the estimated probability vector. For zi_count the gate is Bernoulli from P(nonzero) and the conditional magnitude is drawn on the log1p-z latent using the conformal score (or the magnitude latent SE when scores are unavailable) — not the reported SE of expected count. For multi_proportion groups it draws the K CLR latent columns with their BM latent SEs, projects back to sum-zero CLR space, and decodes to the simplex.

draws_method = "mc_dropout": Run M GNN forward passes in training mode (dropout active) on top of stochastic BM baseline draws. When r_cal = 0, the GNN-dropout term disappears but the BM draw still contributes between-draw variance.

draws_method = "posterior": see the "Posterior draws" section.

Nakagawa & Freckleton (2008, 2011) review the consequences of ignoring missing data in ecological and comparative analyses and argue for multiple imputation as the default.

Value

An object of class "pigauto_mi" with components:

datasets

A list of length m. Each element is a data.frame with the same shape and column types as the input traits; observed cells are preserved and missing cells are filled with the corresponding stochastic draw. These datasets are for prediction diagnostics, not downstream inference, except with draws_method = "posterior" (unless posterior_control$param_uncertainty = "none"), whose datasets are posterior imputations: fit the analysis to them with with_imputations() and pool the fits with pool_mi() (see posterior below and "Posterior draws").

m

Number of stochastic completion datasets.

pooled_point

A single data.frame whose missing cells are replaced by the MC-averaged point estimate. Convenient for reporting but does not provide a valid downstream MI analysis.

se

Matrix of per-cell uncertainty summaries combining the baseline SE and the between-draw standard deviation. A descriptive spread for reporting and ranking only — it is not a Rubin's-rules pooled standard error (it contains no within-imputation variance component and no small-sample df correction) and must not be used for downstream pooling.

mi_workflow

"pigauto_diagnostic_mi", recording that these are prediction-diagnostic completions and cannot be passed to with_imputations() or pool_mi().

imputed_mask

Logical matrix; TRUE where a cell was originally missing.

fit

The underlying pigauto_fit object, retained for diagnostics and for calls to predict() on new data.

data

The pigauto_data object.

tree

The input phylogeny.

species_col

Passed-through species-column name or NULL.

posterior

draws_method = "posterior" only. A list with cell_interval (data.frame: row (row of traits), trait, lower, upper, median; 95% posterior predictive interval on the original scale), diagnostics (data.frame: parameter, rhat, ess_bulk, with attributes "converged", "n_extensions" (the number of automatic extensions made) and "sweeps_per_chain" (the sweeps after burn-in per chain that the diagnostics use)), params (list: Sigma_P and Sigma_E as ⁠K x K x draws⁠ arrays, lambda and mu as ⁠draws x K⁠ matrices, on the latent scale), converged, control, hyper (priors), start (REML lambda and chain-1 start lambda, before clamping to [0.02, 0.98]), draw_index (kept sweeps used for the m datasets), wall_s and sweeps (all sweeps, burn-in included, summed over chains). For this method fit is NULL, se holds the posterior predictive SD of each imputed cell, the class is c("pigauto_posterior_mi", "pigauto_mi", "list") and mi_workflow is "pigauto_posterior_mi_v1" ("pigauto_posterior_plugin_diagnostic" with posterior_control$param_uncertainty = "none", which with_imputations() and pool_mi() refuse).

posterior_improper

Only with posterior_control$param_uncertainty = "both" (validation only). A list with datasets (m completed data.frames drawn with the covariance matrices fixed at their posterior means, from the same chain run as datasets), cell_interval (same columns as posterior$cell_interval), and the fixed Sigma_P and Sigma_E. Not for downstream inference: it ignores parameter uncertainty.

When to use this

With draws_method = "conformal" or "mc_dropout", this function is useful for comparing stochastic prediction behavior from one tree; those draws are not for downstream inference. With draws_method = "posterior" (continuous traits), it returns proper posterior imputations for the analyses described in "Posterior draws": with_imputations() fits the analysis to them and pool_mi() pools those fits. multi_impute_analysis() is the separate analysis-aware backend for its own documented regime.

multi_impute_trees() provides an experimental posterior-tree sensitivity path, but tree uncertainty is not supported by multi_impute_analysis().

Posterior draws

With draws_method = "posterior", the latent (z-scored, optionally log-transformed) traits follow the multivariate phylogenetic mixed model

\mathrm{vec}(Y) \sim \mathrm{N}(1\mu^\top,\; \Sigma_P \otimes R + \Sigma_E \otimes I_n),

with R the phylogenetic correlation matrix of the tree, full K \times K phylogenetic (\Sigma_P) and residual (\Sigma_E) covariance matrices, and a flat prior on \mu. The implied Pagel's lambda of trait k is \Sigma_P[k,k] / (\Sigma_P[k,k] + \Sigma_E[k,k]).

The phylogenetic covariance uses a parameter-expanded inverse-Wishart prior on its standardized covariance matrix and Normal priors on the expansion parameters (Hadfield 2010; Gelman 2006). The default residual covariance prior uses a separation strategy: residual standard deviations have half-Cauchy priors scaled by the observed latent-trait standard deviations, and the residual correlation has an LKJ(1) prior. The legacy posterior_control$residual_prior = "iw" option instead uses an inverse-Wishart prior on the residual covariance. Prior settings are returned in mi$posterior$hyper.

The sampler draws the phylogenetic effects, the trait means and the missing cells as one block with a sparse Cholesky factor of the Hadfield and Nakagawa (2010) node precision, and updates the covariance matrices by parameter-expanded Gibbs steps plus Metropolis moves with the phylogenetic effects integrated out. Convergence is checked with rank-normalised split R-hat and bulk effective sample size (Vehtari et al. 2021); mi$posterior$diagnostics reports them. When any R-hat is at least 1.05 or any ESS is at most 400, the chains are extended automatically (posterior_control$auto_extend and max_extend), and a warning is raised if the rule still fails after the last extension.

The m datasets are posterior predictive draws spaced evenly across the kept sweeps of all chains. Per-cell 95% intervals in mi$posterior$cell_interval come from all kept sweeps (keep_draws of them by default), never from the m datasets.

These draws come from a linear-Gaussian model of the imputed traits on the scale where they were imputed (the log scale for traits log-transformed by log_transform). They are proper imputations for analyses that are linear in the imputed traits on that scale and whose variables are all among the imputed traits, for example a phylogenetic regression of one imputed trait on others (on the log scale for log-transformed traits). Not covered: covariates from outside the imputed traits, because the imputation model does not contain them; nonlinear terms (such as squares) or interactions among imputed traits; and analysing a log-transformed trait on its raw scale. The GNN is not used, and gnn, epochs, missing_frac, lambda_mode and other fitting arguments are ignored (a message lists any that were supplied). Non-continuous traits, multiple observations per species (or any species_col), covariates, and input with no missing cells in the rows of traits (tree tips absent from traits do not count) are errors.

Optional safety floor

When fit_pigauto() is called with safety_floor = TRUE (opt-in; the default is FALSE), the 3-way blend r_BM * BM + r_GNN * GNN + r_MEAN * MEAN propagates through prediction-diagnostic draws through predict.pigauto_fit(). The posterior sampler uses its own model and does not use this blend. For draws_method = "mc_dropout" the mean term contributes no between-draw variance (it is a deterministic scalar per column); between-draw variance comes from the BM-draw and GNN-dropout terms only. For draws_method = "conformal" the blend centre is the 3-way prediction and conformal scores remain calibrated on the blended residuals.

References

Rubin DB (1987). Multiple Imputation for Nonresponse in Surveys. Wiley.

Nakagawa S, Freckleton RP (2008). "Missing inaction: the dangers of ignoring missing data." Trends in Ecology & Evolution 23(11): 592-596.

Gelman A (2006). "Prior distributions for variance parameters in hierarchical models." Bayesian Analysis 1(3): 515-534.

Hadfield JD (2010). "MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package." Journal of Statistical Software 33(2): 1-22.

Hadfield JD, Nakagawa S (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.

Vehtari A, Gelman A, Simpson D, Carpenter B, Buerkner P-C (2021). "Rank-normalization, folding, and localization: an improved R-hat for assessing convergence of MCMC." Bayesian Analysis 16(2): 667-718.

Nakagawa S, Freckleton RP (2011). "Model averaging, missing data and multiple imputation: a case study for behavioural ecology." Behavioral Ecology and Sociobiology 65(1): 103-116.

See Also

impute() for point imputation and multi_impute_analysis() for the narrow analysis-aware inferential backend.

Examples

library(pigauto)
data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
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(df, tree, m = 2L, draws_method = "conformal",
                   epochs = 5L, verbose = FALSE)
print(mi)
lapply(mi$datasets, head)

Analysis-aware multiple imputation for narrow regression models

Description

Generate multiple imputations conditional on a specified downstream regression model. Validity requires MAR, a correctly specified conditional imputation model, and adequate SMCFCS or jomo convergence. This experimental interface is deliberately narrow: exactly one continuous covariate may contain missing values, while the outcome, all other predictors, auxiliary variables, and any grouping variable must be fully observed.

Usage

multi_impute_analysis(
  data,
  formula,
  missing,
  model = c("lm", "glm", "lmer"),
  m = 50L,
  auxiliary = character(),
  seed = NULL,
  control = list()
)

Arguments

Value

An object with classes "pigauto_analysis_mi" and "pigauto_mi". Its datasets component is a list of m completed data frames compatible with with_imputations() and pool_mi(). The object also records the analysis model, formula, engine, controls, engine version, seed, runtime provenance, and imputed-variable metadata. Any engine warning aborts the call instead of returning an inference-ready object. mi_workflow is a stable analysis-aware provenance marker propagated by with_imputations() and checked by pool_mi().

Experimental scope

This function does not fit pigauto's phylogenetic baseline or graph neural network and does not propagate posterior-tree uncertainty. The "lm" route uses proper Bayesian Normal-regression draws. The "glm" route requires the optional smcfcs package, and the "lmer" route requires the optional jomo and lme4 packages. Interactions, transformed terms, random slopes, nested or crossed random effects, multiple incomplete variables, missing outcomes, and pooling of variance components or BLUPs are unsupported.

See Also

multi_impute(), with_imputations(), pool_mi()

Examples

set.seed(11)
dat <- data.frame(
  y = stats::rnorm(80),
  x = stats::rnorm(80),
  z = stats::rnorm(80)
)
dat$x[seq(4, 80, by = 5)] <- NA_real_
mi <- multi_impute_analysis(
  dat, y ~ x + z, missing = "x", model = "lm", m = 5, seed = 12
)
fits <- with_imputations(
  mi, function(d) stats::lm(y ~ x + z, data = d), .progress = FALSE
)
pool_mi(fits)

Posterior-tree prediction sensitivity

Description

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

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.

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.

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.

Share-GNN (tree-sharing) mode

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.

Safety floor + share_gnn (v0.9.1.9002+)

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.

See Also

multi_impute() for single-tree prediction diagnostics or, with draws_method = "posterior", single-tree inference on continuous traits; multi_impute_analysis() for analysis-aware MI; and trees300

Examples

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)
print(mi)
mass_by_tree <- vapply(mi$datasets, function(dat) dat$Mass,
                       numeric(nrow(df)))
apply(mass_by_tree, 1L, stats::sd)

Generate an HTML imputation report from a pigauto fit

Description

Produces an HTML file with interactive charts (loading Chart.js from its public CDN). For fits created with the default gnn = FALSE, the report describes baseline-only imputation. For fits created with gnn = TRUE, it also compares the GNN against the phylogenetic baseline and includes training history. Per-trait metrics, gate values, and held-out diagnostic conformal coverage are included when available.

Usage

pigauto_report(
  fit,
  data = NULL,
  splits = NULL,
  output_path = "pigauto_report.html",
  title = "pigauto Imputation Report",
  open = TRUE
)

Arguments

Value

The output path (invisibly).

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
traits$Mass[seq_len(3L)] <- NA_real_
result <- impute(traits, tree, epochs = 5L, verbose = FALSE)
pigauto_report(result, output_path = tempfile(fileext = ".html"),
               open = FALSE)

Plot a pigauto benchmark

Description

Creates a multi-panel comparison plot showing BM baseline vs pigauto performance across all simulation scenarios.

Usage

## S3 method for class 'pigauto_benchmark'
plot(x, metric = "rmse", ...)

Arguments

Value

Invisible NULL.

Examples

benchmark <- structure(list(summary = data.frame(
  scenario = rep("BM", 2), trait = rep("trait1", 2),
  method = c("baseline", "pigauto"), metric = rep("rmse", 2),
  mean = c(1, 0.8)
)), class = "pigauto_benchmark")
plot(benchmark)

Plot diagnostics for a fitted pigauto model

Description

S3 plot method for pigauto_fit objects, using base R graphics. Three plot types are available: training history, calibrated gate values, and conformal prediction scores.

Usage

## S3 method for class 'pigauto_fit'
plot(x, type = "history", ...)

Arguments

Value

Invisible NULL. Called for its side effect (plotting).

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
data <- preprocess_traits(traits, tree)
splits <- make_missing_splits(data$X_scaled, trait_map = data$trait_map)
# The package default disables the optional GNN, so this example works
# without torch or its Lantern runtime.
fit <- fit_pigauto(data, tree, splits = splits, epochs = 5L,
                   eval_every = 1L, verbose = FALSE, gnn = FALSE)
plot(fit)
plot(fit, type = "gates")
plot(fit, type = "conformal")

Plot predictions from a pigauto model

Description

S3 plot method for pigauto_pred objects, using base R graphics. Three plot types are available: scatter plots of observed vs predicted values, interval plots per species, and probability charts for discrete traits.

Usage

## S3 method for class 'pigauto_pred'
plot(x, data = NULL, splits = NULL, type = "scatter", trait = NULL, ...)

Arguments

Value

Invisible NULL. Called for its side effect (plotting).

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
traits$Mass[seq_len(3L)] <- NA_real_
pred <- impute(traits, tree, epochs = 5L, verbose = FALSE)$prediction
plot(pred)
plot(pred, type = "scatter", data = traits)
plot(pred, type = "intervals", trait = "Mass")

Plot diagnostics for a pigauto result

Description

S3 plot method for pigauto_result objects. Delegates to plot.pigauto_fit on x$fit.

Usage

## S3 method for class 'pigauto_result'
plot(x, ...)

Arguments

Value

Invisible NULL. Called for its side effect (plotting).

Forest-plot style comparison of benchmark results

Description

Takes a data.frame of benchmark results and produces a forest-plot style figure comparing baseline and GNN performance per trait. Each trait appears on the y-axis; points mark RMSE (continuous/count) or accuracy (binary/categorical) for the two methods, connected by a horizontal line colour-coded by whether the GNN improves over the baseline.

Usage

plot_comparison(
  results,
  metric = NULL,
  methods = c("baseline", "pigauto"),
  ...
)

Arguments

Details

When both RMSE and accuracy metrics are present and metric is NULL, a two-panel figure is produced automatically.

Value

Invisible NULL. Called for its side effect (plotting).

Examples

results <- data.frame(
  trait = rep("Mass", 2L), type = rep("continuous", 2L),
  metric = rep("rmse", 2L), method = c("baseline", "pigauto"),
  value = c(1, 0.8)
)
plot_comparison(results)
plot_comparison(results, metric = "rmse")

Plot training history (ggplot2, deprecated)

Description

Deprecated. Use the base-R S3 method instead: plot(fit, type = "history").

Produces a ggplot2 validation loss training curve.

Usage

plot_history_gg(x, ...)

Arguments

Value

A ggplot object.

Examples

fit <- structure(
  list(history = data.frame(epoch = 1:3, val_loss = c(1, 0.7, 0.5))),
  class = "pigauto_fit"
)
plot_history_gg(fit)

Plot uncertainty ribbons for imputed trait values

Description

Plots conformal prediction intervals when present, otherwise an approximate prediction +/- 1.96 * SE ribbon, sorted by predicted value. If truth is supplied, observed values are overlaid as points. Works for continuous, count, and ordinal traits. For binary traits, plots predicted probabilities with an uncertainty ribbon.

Usage

plot_uncertainty(pred_result, truth = NULL, trait_name)

Arguments

Value

A ggplot object.

Examples

pred <- list(
  imputed = data.frame(mass = c(10, 12, 15)),
  se = matrix(c(1, 1.5, 2), ncol = 1,
              dimnames = list(NULL, "mass"))
)
plot_uncertainty(pred, trait_name = "mass")

Pool downstream model fits across multiple imputations (Rubin's rules)

Description

Combine regression coefficients from M model fits – one per imputed dataset – into a single pooled table using Rubin's rules. The pooled standard errors combine within-imputation sampling variance and between-imputation variance. Correct Rubin arithmetic does not make an incompatible imputation model inferentially valid. For pigauto inference, fits should come from the documented multi_impute_analysis() workflow, or, for continuous traits, from multi_impute() with draws_method = "posterior" (see with_imputations() for the analyses those draws cover), via with_imputations(). The backend is experimental and supports fixed-effect coefficients only. Variance components, correlations, random-effect predictions, BLUPs/conditional modes, latent loadings, and other structured parameters are unsupported.

Usage

pool_mi(
  fits,
  conf.level = 0.95,
  coef_fun = NULL,
  vcov_fun = NULL,
  df_fun = NULL,
  tidy_fun = NULL
)

Arguments

Details

Let \hat\theta_i be the coefficient vector from fit i and U_i = \mathrm{vcov}(\mathrm{fit}_i), for i = 1, \ldots, M. Rubin's rules (Rubin 1987) give

\bar\theta = M^{-1} \sum_i \hat\theta_i

W = M^{-1} \sum_i \mathrm{diag}(U_i)

B = (M-1)^{-1} \sum_i (\hat\theta_i - \bar\theta)^2

T = W + (1 + 1/M) B

with pooled standard error \sqrt{T}. The relative increase in variance is r = (1 + 1/M) B / W, the classical pooled df is \nu_{\text{old}} = (M-1)(1 + 1/r)^2, and the fraction of missing information is

\mathrm{fmi} = (r + 2/(\nu + 3)) / (r + 1).

When df_fun returns finite complete-data df nu_com, the Barnard-Rubin (1999) correction combines \nu_{\text{obs}} = ((\nu_{\text{com}}+1)/(\nu_{\text{com}}+3)) \nu_{\text{com}} (1 - \lambda) with nu_old via \nu_{\text{BR}} = 1/(1/\nu_{\text{old}} + 1/\nu_{\text{obs}}). With no between-imputation variation (B = 0), the classical limit is df = Inf, riv = 0, and fmi = 0. If finite complete-data df are supplied, Barnard–Rubin instead retains finite observed-data df and its small-sample FMI adjustment. A completely deterministic quantity (B = W = 0) always has zero FMI and infinite df.

Only fixed-effect coefficients are pooled. Random-effect variances and correlations, BLUPs/conditional modes, latent loadings, and other structured parameters require parameter-specific transformations and are not supported by the automatic pool_mi() adapters. Custom extractors are an expert escape hatch and cannot be inspected by pigauto; using them to select unsupported structured parameters is outside the documented scope. The glmmTMB adapter selects conditional fixed effects only. The drmTMB adapter includes named distributional fixed-effect blocks such as regression coefficients for mu and sigma; those are fixed coefficients, not random-effect variance components. drmTMB and gllvmTMB need to be installed and loadable only when their automatic adapters are used. For either recognized backend class, pooling also requires the fit to report scalar convergence code 0 and pdHess = TRUE; missing or invalid status is rejected before coefficient extraction, including with custom extractors.

Supported classes. lm, glm, gls, lme, merMod (lme4), glmmTMB (conditional component), drmTMB, and gllvmTMB_multi are handled automatically. Other classes are pooled through their coef() / vcov() methods, or through coef_fun / vcov_fun / tidy_fun.

Bayesian fits (brms, MCMCglmm) are rejected because Rubin's rules are not the right tool for posterior samples: they are pooled by concatenating posterior draws across imputations. For brms, use brms::brm_multiple(formula, data = mi$datasets), or brms::combine_models() on a list of existing fits. For MCMCglmm, concatenate the Sol samples across fits. No Bayesian downstream workflow is supported by the analysis-aware backend.

Value

A data.frame with one row per coefficient and columns:

term

Coefficient name.

estimate

Pooled point estimate (mean across fits).

std.error

Pooled standard error sqrt(T) where T = W + (1 + 1/M) * B.

df

Pooled degrees of freedom (Barnard-Rubin if df_fun supplied, else classical Rubin).

statistic

estimate / std.error.

p.value

Two-sided p-value from a t distribution on df.

conf.low, conf.high

Pooled conf.level interval.

fmi

Fraction of missing information.

riv

Relative increase in variance due to non-response.

References

Rubin DB (1987). Multiple Imputation for Nonresponse in Surveys. Wiley.

Barnard J, Rubin DB (1999). "Small-sample degrees of freedom with multiple imputation." Biometrika 86(4): 948-955.

Nakagawa S, Freckleton RP (2008). "Missing inaction: the dangers of ignoring missing data." Trends in Ecology & Evolution 23(11): 592-596.

Nakagawa S, Freckleton RP (2011). "Model averaging, missing data and multiple imputation: a case study for behavioural ecology." Behavioral Ecology and Sociobiology 65(1): 103-116.

See Also

multi_impute_analysis(), with_imputations()

Examples

dat <- data.frame(y = stats::rnorm(30L), x = stats::rnorm(30L),
                  z = stats::rnorm(30L))
dat$x[seq(3L, 30L, by = 5L)] <- NA_real_
mi <- multi_impute_analysis(
  data = dat, formula = y ~ x + z, missing = "x",
  model = "lm", m = 2L
)
fits <- with_imputations(mi, function(d) lm(y ~ x + z, data = d))
pool_mi(fits)

Impute missing traits using a fitted pigauto model

Description

Runs a single forward pass through the fitted model and returns imputed trait values back-transformed to the original scale. Supports all trait types (continuous, binary, categorical, ordinal, count, proportion, zero-inflated count, multi-proportion) and MC dropout for stochastic prediction diagnostics (when n_imputations > 1). The fitted model is a gated ensemble of a phylogenetic baseline and a graph neural network correction; prediction is the per-trait blend (1 - r_cal) * baseline + r_cal * delta_GNN.

Usage

## S3 method for class 'pigauto_fit'
predict(
  object,
  newdata = NULL,
  return_se = TRUE,
  n_imputations = 1L,
  baseline_override = NULL,
  pool_method = c("median", "mean", "mode"),
  clamp_outliers = FALSE,
  clamp_factor = 5,
  match_observed = c("none", "pmm"),
  pmm_K = 5L,
  ...
)

Arguments

Details

When n_imputations > 1, each imputation m draws a BM posterior sample t_BM_draw ~ N(BM_mu, BM_se) on the latent scale for originally-missing cells (BM_se = 0 for observed cells so they are never perturbed). The model runs in train mode (GNN dropout active) using t_BM_draw as input. The final blend is (1 - r_cal) * t_BM_draw + r_cal * GNN_delta(t_BM_draw): when the calibrated gate is zero the imputation is a pure BM posterior draw; when r_cal > 0 both BM draws and GNN dropout contribute variance. Point estimates are the mean (continuous, count) or mode (binary, categorical, ordinal) across passes. The M complete datasets are returned in imputed_datasets for prediction diagnostics. These datasets are not supported for downstream inference. Use multi_impute_analysis() for the narrow analysis-aware backend after consulting its supported-model and lifecycle documentation.

gnn = FALSE fits (traditional-stats mode): when the fit was produced with gnn = FALSE, no GNN was ever trained, so this function makes no torch:: calls at all – prediction is the calibrated baseline blend r_cal_bm * baseline + r_cal_mean * mean_baseline_per_col (r_cal_gnn is always zero) computed on plain R matrices. Two baselines are available: in production mode (the default, no .mask_observed_idx) predictions use baseline_full – fit with splits = NULL, so every observed cell can inform every other cell; in evaluation mode (.mask_observed_idx supplied, as used internally for val/test scoring) predictions use the held-out baseline instead, so scored cells never see their own value. n_imputations > 1 draws are BM posterior samples MU_m ~ N(baseline_mu, baseline_se) at originally-missing cells (observed cells are always restored to their true value), with no GNN dropout contribution.

Decoding per type:

continuous

reverse z-score, then exp() if log-transformed

binary

sigmoid(latent) to probability, round to 0/1

count

reverse z-score of log1p, expm1(), round, clip >= 0

ordinal

reverse z-score, round to nearest valid integer level

categorical

softmax() over K latent columns, argmax

Conformal intervals: conformal_lower / conformal_upper use whatever conformal_method the fit was trained with (fit_pigauto). For "mondrian", the half-width is not constant across rows: each target cell's phylogenetic locality (mean cophenetic distance to its 5 nearest training-observed species for that trait) picks the near- or far-stratum score stored on the fit, so predictions in undersampled clades get wider intervals. See fit_pigauto's conformal_method argument for details.

zi_count conformal intervals are conditional on non-zero. conformal_lower[, nm] / conformal_upper[, nm] for a zi_count trait back-transform the magnitude latent (log1p-z of the count given it is non-zero) plus/minus the conformal half-width that was calibrated on that same magnitude column (fit_pigauto's compute_conformal_scores() scores zi_count on latent_cols[2], never on the decoded E[X]). They are not an interval for the expected value E[X] = P(nonzero) * E[count | nonzero] – constructing a naive \pmquantile interval directly on E[X] would conflate gate and magnitude uncertainty into one number with no valid coverage interpretation. The gate probability P(nonzero) is available separately in pred$probabilities[[nm]]; combine it with the conditional interval yourself if you need an E[X]-scale bound.

Value

A list of class "pigauto_pred" with:

imputed

data.frame of imputed values in original scale with proper R types (numeric, integer, factor, ordered).

imputed_latent

Numeric matrix (n x p_latent) of predictions in latent scale.

se

Numeric matrix (n x n_original_traits) of per-cell uncertainty. See the "What $se means per trait type" section below for what each trait type's column actually contains. NULL if return_se = FALSE.

probabilities

Named list. Binary traits: numeric probability vector. Categorical traits: n x K probability matrix. Other types: not present.

imputed_datasets

List of M stochastic prediction data.frames when n_imputations > 1; NULL otherwise. Not supported for downstream inference.

trait_map

Trait map from the fitted model.

species_names

Character vector.

trait_names

Character vector.

n_imputations

Integer, number of imputations performed.

What $se means per trait type

pred$se is a single n x n_original_traits matrix, but its per-column meaning differs by trait type — do not treat every column as a Gaussian standard error. This table is the definitive reference; other pigauto documentation (e.g. impute's return value) points here instead of repeating it.

These families are incompatible objects sharing one column name: only the continuous-family entries are standard errors in the Gaussian sense, and none of the binary / categorical / zi_count entries may be used in Rubin's-rules arithmetic. The supported downstream multiple-imputation routes start with multi_impute_analysis or, for continuous traits, multi_impute with draws_method = "posterior", then use with_imputations and pool_mi.

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
data <- preprocess_traits(traits, tree)
splits <- make_missing_splits(data$X_scaled, trait_map = data$trait_map)
fit <- fit_pigauto(data, tree, splits = splits, epochs = 5L,
                   verbose = FALSE)
pred <- predict(fit, return_se = TRUE)
pred$imputed
pred$se
pred$probabilities
pred2 <- predict(fit, n_imputations = 2L)
pred2$imputed_datasets

Preprocess trait data: align to tree, encode into latent space

Description

Aligns species in the trait data frame to the tree, detects or accepts trait types (continuous, binary, categorical, ordinal, count, proportion, zi_count), and encodes each trait into a continuous latent matrix.

Usage

preprocess_traits(
  traits,
  tree,
  species_col = NULL,
  trait_types = NULL,
  multi_proportion_groups = NULL,
  log_cols = NULL,
  log_transform = TRUE,
  center = TRUE,
  scale = TRUE,
  covariates = NULL
)

Arguments

Details

When each species has one observation (the default), output rows match tree$tip.label order. When species_col is supplied, multiple observations per species are supported: the output matrix has one row per observation, plus an obs_to_species mapping for the GNN (which operates at species level).

Automatic type detection (when trait_types = NULL) follows the R class of each column — no user input is required for most data:

Two types require an explicit override because they cannot be distinguished from their R class alone:

"proportion"

A numeric column bounded 0–1 looks identical to continuous. Declare it explicitly: trait_types = c(SurvivalRate = "proportion"). Encoded via qlogis(clamp(x, 0.001, 0.999)).

"zi_count"

An integer column with excess zeros looks identical to count. Declare it explicitly: trait_types = c(Parasites = "zi_count"). Encoded as a binary zero/non-zero gate plus log1p-z magnitude.

Practical examples of type assignment:

• Body mass (numeric, all positive) → continuous, auto-log-transformed.

• Clutch size (integer) → count.

• Migratory (factor with levels "Yes"/"No") → binary.

• Diet (factor with >2 levels) → categorical.

• IUCN status (ordered factor, LC < NT < VU < EN < CR) → ordinal. If left as an unordered factor, it becomes categorical — both are valid depending on the question.

• Parasite load (integer with many zeros) → needs trait_types = c(Parasites = "zi_count").

• Survival rate (numeric, values in 0 to 1) → needs trait_types = c(Survival = "proportion").

Latent encoding per type:

continuous

optional log(), then z-score (1 latent column)

binary

0/1 encoding (1 latent column)

count

log1p(), then z-score (1 latent column)

ordinal

integer coding (0 to K-1), then z-score (1 latent column)

categorical

one-hot encoding (K latent columns)

proportion

qlogis(clamp(x, 0.001, 0.999)), then z-score (1 latent column)

zi_count

gate (0/1) + log1p-z of non-zeros (2 latent columns)

multi_proportion

centred log-ratio (CLR) + per-component z-score (K latent columns per group). Rows must sum to 1. Declared via the multi_proportion_groups argument, not trait_types.

Value

A list of class "pigauto_data" with components:

X_scaled

Numeric matrix (n_obs x p_latent), latent encoding. When species_col is NULL, n_obs = n_species.

X_raw

Numeric matrix of continuous traits after optional log but before z-scoring (for backward compatibility).

X_original

Original data.frame (aligned to tree, before encoding).

means

Named numeric vector of column means used for z-scoring (continuous/count/ordinal traits only).

sds

Named numeric vector of column SDs.

species_names

Character vector of unique species matching tree$tip.label order (length = n_species).

obs_species

Character vector of species labels per observation (length = n_obs). When multi-obs, can have duplicates. NULL when species_col is NULL.

obs_to_species

Integer vector (length = n_obs) mapping each observation to its species index in species_names. NULL when species_col is NULL.

n_species

Integer, number of unique species.

n_obs

Integer, number of observations (= n_species when single-obs).

multi_obs

Logical, TRUE when multiple observations per species are present.

trait_names

Character vector of original trait names.

latent_names

Character vector of latent column names.

trait_map

List of trait descriptors (see Details).

p_latent

Integer, total number of latent columns.

log_transform

Logical, legacy field (TRUE if any continuous trait was log-transformed).

Examples

# Single-obs per species (backward compatible)
data(avonet300, tree300, package = "pigauto")
traits <- avonet300
rownames(traits) <- traits$Species_Key
traits$Species_Key <- NULL
pd <- preprocess_traits(traits, tree300)
dim(pd$X_scaled)   # 300 x p_latent

# Multi-obs per species (via species_col)
pd2 <- preprocess_traits(avonet300, tree300, species_col = "Species_Key")
pd2$n_obs      # number of observations
pd2$n_species  # number of unique species

Print a pigauto input check

Description

Print a pigauto input check

Usage

## S3 method for class 'pigauto_check'
print(x, ...)

Arguments

Value

x, invisibly.

Print a pigauto result summary

Description

Print a pigauto result summary

Usage

## S3 method for class 'summary_pigauto_result'
print(x, ...)

Arguments

Value

x, invisibly.

Fetch species range-centroid covariates from GBIF

Description

Pulls occurrence records from the Global Biodiversity Information Facility (GBIF) for each species in species, aggregates them to a median lat / lon centroid, and returns a data.frame ready to be passed to impute via its covariates argument.

Usage

pull_gbif_centroids(
  species,
  cache_dir = NULL,
  occurrence_limit = 500L,
  sleep_ms = 100L,
  verbose = TRUE,
  refresh_cache = FALSE,
  store_points = FALSE
)

Arguments

Details

Caching is strongly recommended — GBIF rate-limits anonymous calls and the per-species fetch is expensive. With cache_dir set, each species gets one RDS file; subsequent calls skip the API.

For each species: resolve taxon via name_backbone, fetch up to occurrence_limit records via occ_search (paginated at 300 per GBIF call), filter out records with hasGeospatialIssues = TRUE and basisOfRecord in c("FOSSIL_SPECIMEN", "LIVING_SPECIMEN"), drop out-of-range coordinates, then take the median latitude and longitude as the species centroid.

Species with zero post-filter records receive NA centroids; their rows are still included in the returned data.frame so it aligns with the input species list.

Requires the optional rgbif package (in DESCRIPTION Suggests). If rgbif is not installed the function errors with an installation message.

Value

A data.frame with columns species, centroid_lat, centroid_lon, n_occurrences. Rownames are set to species. NA centroids are returned for species with no GBIF hits; the caller should decide how to handle them (typical: drop or impute).

See Also

impute (pass the return value as covariates).

Examples

sp <- c("Quercus alba", "Pinus taeda", "Acer saccharum")
cov <- pull_gbif_centroids(sp, cache_dir = "script/data-cache/gbif")
cov_num <- cov[, c("centroid_lat", "centroid_lon"), drop = FALSE]

Fetch per-species bioclim covariates from WorldClim v2.1

Description

Extends pull_gbif_centroids by extracting 19 WorldClim v2.1 bioclim variables at each species' GBIF occurrence points, aggregating per species (median + IQR), and returning a data.frame ready for impute via its covariates argument.

On first call, downloads the WorldClim 10-arc-minute raster stack (~130 MB compressed, ~500 MB unzipped) to worldclim_cache_dir. A sentinel file makes subsequent calls fully offline.

Per-species extracts are also cached (one RDS per species in worldclim_cache_dir/extracts/), so repeated calls for the same species list are instantaneous.

Since v0.9.1.9006 (2026-04-30) extraction operates on the raw per-occurrence point list when pull_gbif_centroids was called with store_points = TRUE. For legacy caches written by older pull_gbif_centroids (centroid-only) calls, extraction silently falls back to the species centroid.

Requires the optional terra package (in DESCRIPTION Suggests). If terra is not installed the function errors with an installation message.

Usage

pull_worldclim_per_species(
  species,
  gbif_cache_dir,
  worldclim_cache_dir,
  resolution = "10m",
  verbose = TRUE,
  refresh_cache = FALSE
)

Arguments

Value

A data.frame with length(species) rows and columns:

• species

• bio1_median, bio1_iqr, ..., bio19_median, bio19_iqr

• n_extracted (integer, number of valid occurrence points that contributed)

Rownames are set to species. Species with no GBIF hits get all-NA bio values and n_extracted = 0.

References

Fick SE, Hijmans RJ (2017). WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37, 4302–4315.

See Also

pull_gbif_centroids (B.1, required input), impute (consume as covariates).

Examples

sp <- c("Quercus alba", "Pinus taeda", "Acer saccharum")
gbif_df <- pull_gbif_centroids(sp,
  cache_dir = "script/data-cache/gbif")
wc_df   <- pull_worldclim_per_species(
  species             = sp,
  gbif_cache_dir      = "script/data-cache/gbif",
  worldclim_cache_dir = "script/data-cache/worldclim")
# Combine for impute()
cov <- cbind(gbif_df[, c("centroid_lat", "centroid_lon")],
             wc_df[, grep("^bio", colnames(wc_df))])

Read trait data from a CSV file or data frame

Description

Loads trait data and sets species names as row names. If a CSV path is supplied, it is read with read.csv. All species identifiers must be present and unique because row names identify one trait row per species. For repeated observations, use read.csv() and pass the resulting data frame with species_col to check_pigauto() and impute(). By default, all non-species columns are returned (numeric, factor, integer, etc.). Use trait_cols to select a subset.

Usage

read_traits(
  x,
  species_col = "species",
  trait_cols = NULL,
  factor_cols = NULL,
  ordered_cols = NULL
)

Arguments

Value

A data.frame with unique species names as row names.

Examples

df <- data.frame(species = c("Sp_a", "Sp_b"), mass = c(10, 20))
traits <- read_traits(df, species_col = "species")

Read a phylogenetic tree from a file

Description

Reads a Newick or NEXUS tree file using ape. Tries ape::read.tree first; falls back to ape::read.nexus if that fails.

Usage

read_tree(path)

Arguments

Details

For bird phylogenies, BirdTree provides a species-subset download tool for up to 2,500 species and full-tree distributions for larger sets. Download and unzip the trees, then read a Newick tree file here. BirdTree documents its distributions as Newick trees; this function also accepts NEXUS files. A file with one tree returns a "phylo" object; a file with several trees returns a "multiPhylo" object.

pigauto does not download BirdTree data automatically. Pass a single-tree result to impute(traits, tree). A multi-tree result can be passed to multi_impute_trees(traits, trees) for experimental prediction- sensitivity checks only; this path has not been validated for downstream inference.

BirdTree's subset tool defaults to at least 100 draws and recommends a reasonable sample (more than 100) for analyses using its full-tree distributions. pigauto's multi_impute_trees() remains experimental: it supports descriptive checks of prediction sensitivity only, and its stochastic datasets are not validated for downstream inference. See vignette("tree-uncertainty").

BirdTree asks researchers using full or partial tree data to cite Jetz et al. (2012). If you use the BirdTree web tool, also credit BirdTree.org.

Value

An object of class "phylo" when the file contains one tree, or "multiPhylo" when it contains multiple trees.

References

Jetz W, Thomas GH, Joy JB, Hartmann K, Mooers AO (2012). The global diversity of birds in space and time. Nature, 491, 444-448. doi:10.1038/nature11631. BirdTree.org. Downloads, phylogeny subsets, and FAQ. https://birdtree.org/downloads/ https://birdtree.org/subsets/ https://birdtree.org/faq/.

See Also

multi_impute_trees

Examples

path <- tempfile(fileext = ".tre")
ape::write.tree(ape::rtree(10L), path)
tree <- read_tree(path)

Save a fitted pigauto model

Description

Serialises a pigauto_fit object to disk, including the torch model state. The saved file can be loaded on any machine with the same pigauto version.

Usage

save_pigauto(fit, path, compress = TRUE)

Arguments

Details

Standard saveRDS() does not work for pigauto fits because torch tensor states cannot be serialised by R's native mechanism. This function converts the model state to a portable raw format before saving.

Value

Invisible path.

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
traits$Mass[seq_len(3L)] <- NA_real_
fit <- impute(traits, tree, epochs = 5L, verbose = FALSE)$fit
path <- tempfile(fileext = ".pigauto")
save_pigauto(fit, path)
fit2 <- load_pigauto(path)

Run a simulation benchmark for pigauto

Description

Generates trait data under various evolutionary models, introduces missing data, fits both the Brownian motion baseline and the full pigauto GNN, and compares performance. This is the recommended way to assess pigauto on data with known properties before applying it to real data.

Usage

simulate_benchmark(
  n_species = 100L,
  n_traits = 4L,
  scenarios = c("BM", "OU", "regime_shift", "nonlinear", "mixed"),
  missing_frac = 0.25,
  n_reps = 3L,
  epochs = 500L,
  verbose = TRUE,
  seed = NULL,
  gnn = TRUE,
  ...
)

Arguments

Details

Available scenarios:

"BM"

Pure Brownian motion – the baseline is exact, so the GNN should tie or slightly improve via inter-trait correlations.

"OU"

Ornstein-Uhlenbeck – stabilising selection constrains variation. BM over-estimates evolutionary variance.

"regime_shift"

Two-regime BM – clade-specific optima create bimodal distributions that BM cannot capture.

"nonlinear"

Non-linear inter-trait relationships – the GNN's multi-layer message passing can capture quadratic and interaction effects that BM's linear covariance misses.

"mixed"

Mixed trait types: 2 continuous + 1 binary + 1 categorical (3 levels). Tests the full type pipeline.

Value

An object of class "pigauto_benchmark" with:

results

data.frame with columns: scenario, rep, method, trait, type, metric, value, n_test.

summary

data.frame averaged across replicates.

scenarios

character vector of scenarios run.

n_reps

integer.

n_species

integer.

Examples

bench <- simulate_benchmark(n_species = 20L, n_traits = 2L,
                            scenarios = "BM", epochs = 5L, n_reps = 1L,
                            verbose = FALSE, gnn = FALSE)
bench$summary
plot(bench)

Simulate non-BM trait data for benchmarking

Description

Generates species trait data under models that deviate from Brownian Motion. Supports OU (stabilising selection), regime shifts (clade-specific optima), and non-linear inter-trait correlations.

Usage

simulate_non_bm(
  tree,
  n_traits = 4,
  scenario = c("OU", "regime_shift", "nonlinear"),
  alpha = 2,
  theta = 0,
  sigma = 1,
  shift_magnitude = 2,
  seed = NULL
)

Arguments

Value

A data.frame with species as rownames.

Examples

tree <- ape::rtree(12)
traits <- simulate_non_bm(tree, n_traits = 2, scenario = "OU", seed = 1)
head(traits)

Suggest model-based candidate observations

Description

For a fitted pigauto_result (returned by impute), compute model-based BM/label-propagation proxy reductions in predictive variance or entropy across currently-missing cells if each candidate cell were observed next. Use the returned ranking as sampling-design guidance under those model assumptions; it is not an optimal-design guarantee or a demonstrated field gain.

Usage

suggest_next_observation(
  result,
  top_n = 10L,
  by = c("cell", "species"),
  types = c("continuous", "count", "ordinal", "proportion", "binary", "categorical",
    "zi_count", "multi_proportion")
)

Arguments

Details

Two metrics are supported, dispatched by trait type:

Continuous-family traits (continuous, count, ordinal, proportion) use the BM variance-reduction formula. The variance- reduction formula is derived from a Sherman-Morrison rank-1 inverse update on the BM conditional MVN: adding species s to the observed set updates the inverse correlation matrix by a known closed form. For each candidate cell (s, t),

\Delta V(s, t) = \sigma_t^2 \sum_{i \in \mathrm{miss}_t} \frac{D_{ik}^2}{\alpha_k}

where D = R_{mm} - R_{mo} R_{oo}^{-1} R_{om} is the residual matrix at currently-missing cells, \alpha_k = D_{kk} is the current relative leverage of cell k, and \sigma_t^2 is the REML BM variance for trait t.

Discrete traits (binary, categorical) use a label- propagation expected-entropy-reduction formula. The current LP probability at species i is p_i = \mathrm{sim}[i, \mathrm{obs}] y_{\mathrm{obs}} / \sum \mathrm{sim}[i, \mathrm{obs}], with entropy H(p_i) = -\sum_k p_{i,k} \log p_{i,k}. After observing s_{\mathrm{new}} with unknown class y_{\mathrm{new}}, the new LP probability has a closed form, and the expected entropy is averaged over P(y_{\mathrm{new}}) = current LP estimate at s_{\mathrm{new}}. Total expected entropy reduction sums across all currently-missing cells (the entropy at s_{\mathrm{new}} itself drops to 0).

Variance and entropy are NOT directly comparable. The output sorts within each metric and the cross-metric ordering by delta is approximate. When you want a strict ranking, filter by metric first.

Reductions are summed across the included traits for each species when by = "species", supporting the typical use case where measuring a species observes all of its currently-missing traits at once. At by = "species", the per-trait variance and entropy reductions are summed separately into delta_var_total and delta_entropy_total columns; the delta column is whichever is non-NA (or delta_var_total when both are populated). Cross-type species-level ranking is approximate – see the variance-vs- entropy caveat above.

The ranking is conditional on the fitted BM or label-propagation proxy and its missing-data assumptions. It does not establish that collecting a selected observation improves a future field study or downstream analysis.

zi_count (v2): observing a missing zi_count cell reveals the gate value (entropy reduction at the gate column, computed via the LP binary formula) AND, with probability p_{\mathrm{gate}} (current LP estimate at s_{\mathrm{new}}), reveals a magnitude (variance reduction at the magnitude column, computed via the BM Sherman-Morrison formula on the gate=1 subset). Output rows for zi_count populate BOTH delta_var_total (= expected magnitude variance reduction = p_{\mathrm{gate}} \times \Delta V_{\mathrm{mag}}) AND delta_entropy_total (= gate entropy reduction). metric is set to "variance" so the row sorts on the magnitude scale; delta_entropy_total is available for users who care about gate-uncertainty separately.

multi_proportion (v2): observing a row reveals all K simplex components simultaneously. Per-component variance reductions are computed via BM Sherman-Morrison on each CLR-z latent column, summed across components. metric is "variance"; delta_var_total is the K-component sum.

Value

A data.frame of class "pigauto_active". Columns when by = "cell": species, trait, type, metric ("variance" or "entropy"), delta, delta_var_total (NA for discrete rows), and delta_entropy_total (NA for continuous rows), sorted by delta descending. When by = "species": species, delta_var_total, delta_entropy_total, n_traits_missing, sorted by the SUM of available metrics.

See Also

impute

Examples

data(avonet300, tree300, package = "pigauto")
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
traits$Mass[seq_len(3L)] <- NA_real_
res <- impute(traits, tree, epochs = 5L, verbose = FALSE)
suggest_next_observation(res, top_n = 5L)
suggest_next_observation(res, top_n = 5L, by = "species")

Summary method for pigauto_fit objects

Description

Prints a formatted evaluation table including per-trait metrics, gate calibration, and conformal coverage. Requires the original pigauto_data object to compute test-set performance.

Usage

## S3 method for class 'pigauto_fit'
summary(object, ..., data = NULL)

Arguments

Value

Invisibly returns the evaluation data.frame (or NULL if data is not supplied).

Examples

data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
pd <- preprocess_traits(traits, tree)
splits <- make_missing_splits(pd$X_scaled, trait_map = pd$trait_map)
fit <- fit_pigauto(pd, tree, splits = splits, epochs = 5L,
                   verbose = FALSE)
summary(fit, data = pd)

Summarise a pigauto result

Description

Summarise a pigauto result

Usage

## S3 method for class 'pigauto_result'
summary(object, ...)

Arguments

Value

A small summary_pigauto_result object containing input-check, completion, output-role, and evaluation summaries.

Example bird phylogeny for the 300 species in avonet300

Description

An object of class 'phylo' from the ape package. It is posterior sample member 69 from the Hackett-backbone half of BirdTree's 100-tree bird collection, retrieved through the megatrees package and pruned to the 300 species in avonet300. The MIT software licence for megatrees does not establish rights for the underlying BirdTree data; see the installed ‘NOTICE’. This object is included for examples, not as an MCC or consensus tree or a recommended tree for substantive analysis.

Usage

tree300

Format

An object of class phylo with 300 tips.

Source

Li (2026), megatrees 1.0.0; member 69 of megatrees::get_tree_bird_n100() is from the Hackett-backbone set. The BirdTree distribution is described by Jetz et al. (2012).

References

Hackett SJ, Kimball RT, Reddy S, et al. (2008). A phylogenomic study of birds reveals their evolutionary history. Science, 320(5884), 1763-1768. doi:10.1126/science.1157704.

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

Example bird phylogeny for the species in avonet_full

Description

An object of class 'phylo' from the ape package. It is posterior sample member 69 from the Hackett-backbone half of BirdTree's 100-tree bird collection, pruned to the species present in avonet_full rather than the 300-species example subset. The MIT software licence for megatrees does not establish rights for the underlying BirdTree data; see the installed ‘NOTICE’. This object is included for examples, not as an MCC or consensus tree or a recommended tree for substantive analysis.

Usage

tree_full

Format

An object of class phylo with 9,993 tips.

Details

Row order in avonet_full matches tip order in tree_full: all(avonet_full$Species_Key == tree_full$tip.label) returns TRUE.

Source

Li (2026), megatrees 1.0.0; member 69 of megatrees::get_tree_bird_n100() is from the Hackett-backbone set. The BirdTree distribution is described by Jetz et al. (2012).

References

Hackett SJ, Kimball RT, Reddy S, et al. (2008). A phylogenomic study of birds reveals their evolutionary history. Science, 320(5884), 1763-1768. doi:10.1126/science.1157704.

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

See Also

tree300, avonet_full

50 posterior phylogenies for the 300 species in avonet300

Description

A multiPhylo list of 50 phylogenetic trees sampled from the 100-tree BirdTree bird collection retrieved through megatrees. The collection's first 50 trees use the Ericson backbone and its next 50 use the Hackett backbone. The fixed sample generated with seed 42 contains 28 Ericson and 22 Hackett trees, each pruned to the 300 species in avonet300. Topologies and branch lengths vary across this posterior sample. The MIT software licence for megatrees does not establish rights for the underlying BirdTree data; see the installed ‘NOTICE’.

Usage

trees300

Format

An object of class multiPhylo containing 50 phylo objects, each with 300 tips.

Details

Use with multi_impute_trees for experimental sensitivity of point imputations to posterior-tree choice. Tree uncertainty is not supported by the analysis-aware MI backend, and these stochastic datasets are not validated for downstream inference.

Source

Li (2026), megatrees 1.0.0; trees pruned from megatrees::get_tree_bird_n100(). The source collection and its Ericson and Hackett backbones are described by Ericson et al. (2006), Hackett et al. (2008), and Jetz et al. (2012).

References

Ericson PGP, Anderson CL, Britton T, et al. (2006). Diversification of Neoaves: integration of molecular sequence data and fossils. Biology Letters, 2(4), 543-547. doi:10.1098/rsbl.2006.0523.

Hackett SJ, Kimball RT, Reddy S, et al. (2008). A phylogenomic study of birds reveals their evolutionary history. Science, 320(5884), 1763-1768. doi:10.1126/science.1157704.

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

See Also

tree300, avonet300, multi_impute_trees

Fit a downstream model on every imputed dataset

Description

Apply a user-supplied model-fitting function .f to each of the M complete datasets stored in an analysis-aware MI object and return the list of fits. This is the supported fitting entry for multi_impute_analysis() and for multi_impute() with draws_method = "posterior". Prediction-diagnostic ("conformal", "mc_dropout") and posterior-tree sensitivity completions stop with an actionable error before any model is fitted. The workflow is experimental and limited to fixed-effect coefficients and their covariance matrices.

Usage

with_imputations(
  mi,
  .f,
  ...,
  .progress = interactive(),
  .on_error = c("continue", "stop")
)

Arguments

Value

A list of length M with class "pigauto_mi_fits". Each element is either a model fit or, if .f errored on that imputation and .on_error = "continue", an object of class "pigauto_mi_error" containing the captured condition. pool_mi() filters error elements automatically.

Posterior draws and congeniality

The draws_method = "posterior" completions come from a linear-Gaussian phylogenetic mixed model of the imputed traits only, on the scale where they were imputed (the log scale for traits log-transformed by log_transform). They are proper imputations for analyses that are linear in the imputed traits on that scale and whose variables are all among the imputed traits, for example a phylogenetic regression of one imputed trait on others (on the log scale for log-transformed traits). Not covered: covariates from outside the imputed traits, nonlinear terms (such as squares) or interactions among imputed traits, and analysing a log-transformed trait on its raw scale. The formula passed to .f is not checked; this is the analyst's responsibility. Plug-in draws from posterior_control$param_uncertainty = "none" (validation only) are refused.

See Also

multi_impute_analysis(), pool_mi()

Examples

dat <- data.frame(y = stats::rnorm(30L), x = stats::rnorm(30L),
                  z = stats::rnorm(30L))
dat$x[seq(3L, 30L, by = 5L)] <- NA_real_
mi <- multi_impute_analysis(
  data = dat, formula = y ~ x + z, missing = "x",
  model = "lm", m = 2L
)
fits <- with_imputations(mi, function(d) stats::lm(y ~ x + z, data = d))
pool_mi(fits)