Getting started with pigauto

Common declarations for your own data

These are separate templates, not additional steps in the bundled bird example. Choose only the pattern that matches your data, and replace the example column names with columns in your own trait table. Each example assumes that traits and tree are already loaded.

Ordered categories

Declare the scientific order of the levels; pigauto never guesses it.

result <- impute(traits, tree, ordinal_levels = list(
  threat = c("LC", "NT", "VU", "EN", "CR")
))

# Equivalent: construct the ordered factor yourself.
traits$threat <- ordered(
  traits$threat,
  levels = c(
    "LC", "NT", "VU", "EN", "CR"
  )
)
result <- impute(traits, tree)

An integer-valued measurement treated as continuous

result <- impute(
  traits, tree,
  trait_types = c(
    length_mm = "continuous"
  )
)

Proportions and zero-inflated counts

Numeric proportions and zero-inflated counts need explicit type declarations.

result <- impute(
  traits, tree,
  trait_types = c(
    cover = "proportion",
    parasites = "zi_count"
  )
)

Compositional data

These grouped columns must be observed or missing together by row and sum to one when observed.

result <- impute(
  traits, tree,
  multi_proportion_groups = list(
    diet = c(
      "diet_insect",
      "diet_fruit",
      "diet_seed"
    )
  )
)

How it works (briefly)

By default pigauto fits a phylogenetic baseline: Brownian-motion conditional imputation for continuous-family latent columns and phylogenetic label propagation or a threshold model for discrete traits, with a per-trait Pagel’s lambda. The architecture article explains how pigauto selects among these baselines. A held-out validation split calibrates the blend and the conformal intervals. Split conformal assumes the calibration cells resemble the scored cells. Masking a whole clade can violate that assumption: residuals from other clades may not represent the held-out clade, so nominal coverage need not hold there. For the current implementation, treat these intervals as nominal diagnostics, not a general 95% coverage guarantee: validation data also inform model selection and gate calibration. See the interval limitations.

An optional graph neural network (GNN) is trained when you set gnn = TRUE. It can learn additional patterns from tree topology, cross-trait structure, and user covariates. Covariates can also enter a covariate-aware per-column numeric baseline when that fallback is selected. Joint MVN and threshold-joint routes omit them and emit a separate baseline warning. One-versus-rest categorical routes also omit covariates but do not trigger that warning themselves. In a mixed fit, another selected joint route can still trigger it. Validation selects how much the GNN contributes for each trait. The default is gnn = FALSE. In the tested single-observation benchmarks, adding the GNN did not improve held-out prediction and took longer to fit. A later repeated-observation study with covariates found gains over the same fitted model without its neural contribution in two prespecified scenarios, but not in a no-signal control; separately fitted GLS had lower error in the interactive scenario. These are different settings and comparisons, not evidence of general GNN superiority. See the validation results and their limits for details. For the implementation name and code map, see the optional architecture article.

The phylogenetic structure is encoded as:

• A symmetric-normalised Gaussian-kernel adjacency matrix derived from cophenetic (patristic) distances.
• k Laplacian eigenvectors as node positional features: these encode phylogenetic clusters at multiple scales.

library(pigauto)
library(ape)

How the imputation pipeline works

The workflow above is the primary interface. The following is the fine-grained pipeline for inspection, caching, or benchmarking.

This template assumes you already have a trait data frame named traits and a matching phylogenetic tree named tree. If you are following the bundled bird example, run the AVONET setup in the next section first and use avonet300 and tree300 in place of those placeholders.

result <- impute(traits, tree)

result$completed    # observed values preserved; some missing values may remain NA
result$imputed_mask # TRUE marks missing targets in fitted rows, not successful fills
result$prediction$se # per-cell, type-dependent uncertainty summaries

impute() runs the full pipeline internally: preprocessing, baseline fitting, graph construction, calibration, and prediction (plus GNN training when gnn = TRUE). The rest of this vignette walks through each step individually for users who need fine-grained control, want to cache intermediate objects, or are benchmarking components.

The AVONET bird data

pigauto ships with avonet300, a 300-species subset of the AVONET v3 bird morphology database (Tobias et al. 2022, Ecology Letters), and tree300. The bundled tree is member 69 of the 100-tree BirdTree collection, from the Hackett-backbone set, and pruned to those 300 species; see ?tree300 for the package-specific extraction record. Hackett et al. (2008) describes the backbone. The tree is an example, not an MCC or consensus tree or a recommendation for substantive analyses. For your own analysis, obtain a tree suited to your taxa and cite Jetz et al. (2012) for the BirdTree distribution. If you use the BirdTree website, also cite BirdTree.org. You can download a tree file and read it with read_tree(); see BirdTree downloads and the tree-uncertainty article for options.

data(avonet300, tree300, package = "pigauto")

head(avonet300)
#>                Species_Key    Mass Beak.Length_Culmen Tarsus.Length Wing.Length
#> 2   Nothoprocta_pentlandii   300.7               27.9          37.0       148.2
#> 1         Eudromia_formosa   772.9               32.5          49.5       223.9
#> 3           Rhea_americana 23000.0               86.5         308.0       604.5
#> 4     Ptilopachus_petrosus   193.0               16.9          28.0       117.2
#> 9      Bambusicola_fytchii   313.7               21.6          45.2       142.8
#> 10 Francolinus_psilolaemus   478.3               30.2          43.1       167.3
#>    Trophic.Level Primary.Lifestyle Migration
#> 2      Herbivore       Terrestrial  Resident
#> 1      Herbivore       Terrestrial  Resident
#> 3       Omnivore       Terrestrial  Resident
#> 4      Herbivore       Terrestrial  Resident
#> 9      Herbivore       Terrestrial  Resident
#> 10     Herbivore       Terrestrial  Resident
ape::Ntip(tree300)
#> [1] 300

The four morphological columns are example measurements. The dataset also includes two categorical traits and one ordinal trait, so this walk-through uses seven original traits in total. Their distributions and phylogenetic signal vary; this tutorial does not use them to establish a general performance claim.

Step 1: Preprocess

preprocess_traits() log-transforms and z-scores the traits, and reorders rows to match the phylogeny’s tip order (required by the graph operations).

traits <- avonet300[, -1]              # drop Species_Key column
rownames(traits) <- avonet300$Species_Key

pd <- preprocess_traits(traits, tree300, log_transform = TRUE)
print(pd)
#> pigauto_data
#>   Species: 300 
#>   Traits:  7 
#>   Types:   categorical=2, continuous=4, ordinal=1 
#>   Latent columns: 14 
#>   Missing values: 0 %

The object pd is a pigauto_data list containing:

• X_scaled: the z-scored, log-transformed latent matrix (300 × 14). The four continuous traits and the ordinal trait each use one latent column; the two categorical traits expand to one column per level.
• means, sds: per-trait normalisation parameters (for back-transformation)
• log_transform: TRUE, recorded so predict() can invert it automatically

log_transform = TRUE (the default) logs every continuous trait whose observed values are all positive. If a trait is already on a log scale, such as log body mass, set log_transform = FALSE (and log any raw traits yourself), or it will be imputed on a log(log) scale; see ?impute.

Step 2: Create evaluation splits

We hold out 25% of eligible observed species-by-trait entries for evaluation. All latent columns that encode a selected trait are held out together. These values are not used during training or baseline fitting; they let you assess prediction error under this masking design. The scores may differ when real-world missingness follows another pattern.

splits <- make_missing_splits(pd$X_scaled, missing_frac = 0.25, seed = 42,
                              trait_map = pd$trait_map)

cat("Held-out latent cells:", length(splits$val_idx) + length(splits$test_idx), "\n")
#> Held-out latent cells: 1040
cat("  Validation:", length(splits$val_idx), "\n")
#>   Validation: 260
cat("  Test:      ", length(splits$test_idx), "\n")
#>   Test:       780

The missing cells are split further: validation cells guide model selection, while test cells assess predictions under this masking design. Test scores can differ from performance on naturally missing values if the masking design does not resemble the real pattern of missingness.

Step 3: Fit the BM baseline

fit_baseline() fits the phylogenetic BM baseline after masking the held-out cells. This gives the benchmark against which an optional GNN would be compared.

baseline <- fit_baseline(pd, tree300, splits = splits, model = "BM")

cat("Baseline object contains:\n")
#> Baseline object contains:
cat("  mu matrix:", nrow(baseline$mu), "x", ncol(baseline$mu), "\n")
#>   mu matrix: 300 x 14
cat("  se matrix:", nrow(baseline$se), "x", ncol(baseline$se), "\n")
#>   se matrix: 300 x 14

Step 4: Build the phylogenetic graph

build_phylo_graph() computes the adjacency matrix and Laplacian spectral features from the tree. The result can be cached to disk for reuse.

graph <- build_phylo_graph(
  tree300,
  k_eigen    = 8,
  sigma_mult = 0.5
)

cat("Graph: n =", graph$n, "species\n")
#> Graph: n = 300 species
cat("Adjacency: [", nrow(graph$adj), "x", ncol(graph$adj), "]\n")
#> Adjacency: [ 300 x 300 ]
cat("Spectral coords: [", nrow(graph$coords), "x", ncol(graph$coords), "]\n")
#> Spectral coords: [ 300 x 8 ]
cat("Kernel bandwidth sigma:", round(graph$sigma, 3), "\n")
#> Kernel bandwidth sigma: 80.293

For repeated runs, supply cache_path = "graph_cache.rds" to skip recomputation.

Step 5: Train the model

With the default gnn = FALSE, fit_pigauto() fits no network and makes no torch calls; it calibrates the baseline and the conformal intervals. The chunk below sets gnn = TRUE to show the optional GNN. Its training loop corrupts a random 55% of observed cells each epoch with a learnable mask token, then minimises a type-appropriate loss on those known values. A shrinkage penalty keeps GNN predictions closer to the phylogenetic baseline. For typed data, early stopping monitors a composite held-out validation loss: each trait contributes a loss suited to its type. Lower is better. This value is not an RMSE in the original trait units.

fit <- fit_pigauto(
  data            = pd,
  tree            = tree300,
  splits          = splits,
  graph           = graph,
  baseline        = baseline,
  hidden_dim      = 64,
  k_eigen         = 8,
  dropout         = 0.10,
  lr              = 3e-3,
  epochs          = 2000,
  corruption_rate = 0.55,
  lambda_shrink   = 0.03,
  eval_every      = 100,
  patience        = 10,
  verbose         = TRUE,
  seed            = 1,
  gnn             = TRUE
)

With verbose = TRUE, each evaluation reports training-loss components (rec, shrink, and gate), the learning rate, and, when validation data are available, the current and best validation loss and remaining patience. The exact values depend on your data and settings. For this mixed-type example, val is a composite loss, not an RMSE in Mass units.

print(fit)

The fit contains seven original traits and 14 latent columns. Its summary reports the original trait names and the mix of continuous, categorical, and ordinal types. Do not interpret the latent-column count as the number of input traits.

Step 6: Predict

predict.pigauto_fit() returns predictions on the original (non-log) scale. If you request n_imputations > 1, it also returns stochastic prediction-diagnostic draws. These draws are not validated multiple imputations for downstream inference. For that question, first check whether your analysis meets a supported route in the multiple-imputation guide.

pred <- predict(fit, return_se = TRUE)

# pred$imputed: 300-row data.frame with the 7 original trait columns
# pred$imputed_latent: 300 x 14 matrix on the latent scale
# pred$se: 300 x 7 per-trait uncertainty matrix
head(pred$imputed)

# Conformal prediction intervals (nominal held-out diagnostic)
pred$conformal_lower[, "Mass"]    # lower bound per species (original units)
pred$conformal_upper[, "Mass"]    # upper bound per species (original units)
# Coverage is not stored on the prediction object. Score it with evaluate()
# (or evaluate_imputation() on held-out cells).

The categorical columns in pred$imputed are returned as factors. For those traits, pred$se contains uncertainty scores rather than Gaussian standard errors. The latent prediction has one column per encoded category, while the decoded data frame has one column per original trait.

Step 7: Evaluate and read the gate

evaluate_imputation() scores held-out cells on the scale that matches the trait type. For a single latent column in this mixed-trait fit, no direct covariates are supplied and the mean safety floor is off. On high-signal morphometrics (AVONET-like mass, beak, tarsus, wing), the calibrated neural weight for a latent column can be zero (fit$r_cal_gnn = 0), leaving the phylogenetic baseline in this example. In other routes, a zero neural weight closes only the GNN component; inspect the BM, GNN, and mean weights and any direct covariate terms to understand the full prediction.

CRAN 0.11.0: use evaluate(fit) for held-out scores. That release silently ignores evaluation = TRUE in predict(), so the call below would score production predictions and can give misleading results. The GitHub version supports the argument. evaluate(fit) returns a long table (method, trait, type, metric, value, n_test), whereas evaluate_imputation() returns a wide table. Do not compare production predict(fit) values with held-out truth as an evaluation shortcut.

# Held-out scores (z-score RMSE for continuous traits)
# Mask validation/test cells before scoring; plain predict() is for imputation.
pred_eval <- predict(fit, evaluation = TRUE)
evaluate_imputation(pred_eval, pd$X_scaled, splits)

# Per-trait calibrated gate. 0 means baseline-only.
fit$r_cal

Do not treat a GNN-versus-BM table as the default story. If the gate is closed, a lower GNN training loss on masked cells is not a claim that the ensemble beat BM on the test cells you ship. How much (if anything) the GNN adds depends on the trait, clade, missingness pattern, and whether calibration left the gate open.

Step 8: Visualise

plot(fit, type = "history")

plot(pred, type = "intervals", trait = "Mass")

Experimental analysis-aware multiple imputation

Downstream inference requires imputations compatible with the substantive analysis. multi_impute_analysis() therefore receives the formula and model class before drawing one incomplete continuous covariate under MAR.

This interface implements a narrow, experimental workflow. In the registered analysis-aware tests, 23 of 24 method-by-scenario results were indeterminate and one binomial-logit GLM result failed; none was classified as supported. These 24 results include 12 complete-data controls. See the registered inferential results.

This is a template, not a standalone example: analysis_data must contain a continuous outcome y, one incomplete continuous predictor x, and a fully observed continuous predictor z. This example selects the Gaussian lm route; change both the outcome and model setting if your analysis uses another supported outcome type.

analysis_data$z_sq <- analysis_data$z^2

mi <- multi_impute_analysis(
  data = analysis_data,
  formula = y ~ x + z,
  missing = "x",
  model = "lm",
  m = 50L,
  auxiliary = "z_sq",
  seed = 1L
)

fits <- with_imputations(mi, function(d) lm(y ~ x + z, data = d))
pool_mi(fits)

The initial dispatch is deliberately narrow: proper Bayesian normal-regression MI for Gaussian lm, smcfcs for binomial-logit glm, and jomo.smc for a Gaussian lmer with one random intercept. Users must precompute auxiliary terms such as z_sq. Fixed effects are the only supported pooling targets.

Conformal-width, Brownian/MC-dropout, and PMM draws remain experimental prediction diagnostics. They failed the downstream fixed-effect gate and must not be passed to with_imputations()/pool_mi() for inference. For continuous traits with one row per species and no covariates, multi_impute() with its default draws_method = "auto" selects posterior draws when the data meet its eligibility conditions. This is a software capability. The documented downstream route fits analyses linear in the imputed continuous traits on their imputation scale and uses variables among those traits; it does not establish calibrated inference. In the registered validation, five of eight posterior slope regimes failed and three were indeterminate; no inferential scenario was classified as supported. See the registered inferential results and multiple imputation article for scope and conditions.

Phylogenetic tree uncertainty

Use multi_impute_trees() to see whether point predictions change across a posterior sample of trees. This is a prediction-sensitivity check only; its outputs do not provide inferential uncertainty for downstream estimates and cannot be pooled for inference. See the tree-sensitivity guide for a step-by-step example.

Active imputation: where to measure next

suggest_next_observation() returns model-based BM/label-propagation proxy rankings under its stated assumptions. It is not an optimal sampling guarantee or a demonstrated field gain.

res <- impute(traits, tree300)

# Top-10 individual (species, trait) cells, ordered by expected
# total variance reduction:
suggest_next_observation(res, top_n = 10, by = "cell")

# Or aggregated by species (sum of variance reductions across the
# species' currently-missing continuous-family traits):
suggest_next_observation(res, top_n = 10, by = "species")

The delta_var_total column is the model-based expected reduction in total predictive variance across all currently-missing cells if you observed that candidate next:

\[ \Delta V(s_{\text{new}}) = \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. Treat it as a proxy ranking conditional on the fitted model, rather than a claim about a future sampling campaign.

Scope. suggest_next_observation() v2 (2026-05-01) supports all eight pigauto trait types:

• Variance reduction for continuous-family traits (continuous, count, ordinal, proportion) and the magnitude column of zi_count, computed by closed-form Sherman-Morrison rank-1 update on the BM phylogenetic-correlation matrix.
• Entropy reduction for binary, categorical, and the gate column of zi_count, computed by closed-form label-propagation posteriors.
• Multi-component variance reduction (per-component BM, summed across the K CLR-z latent columns) for multi_proportion.

For zi_count rows the output populates BOTH delta_var_total (probability-weighted magnitude variance reduction) AND delta_entropy_total (gate entropy reduction); the metric column is set to "variance" so rows sort on the magnitude scale. See ?suggest_next_observation for full details.

Variance vs entropy are not directly comparable: ordering across metrics by delta_var_total is approximate. When you want a strict ranking on a single metric, filter by metric first or restrict types to continuous-family-only or discrete-only.

Multi-obs (multiple observations per species) inputs are not supported and error with a clear message.

Advanced controls

These templates assume your own trait table is named traits and its matching phylogeny is named tree. In the bundled bird walkthrough above, use tree300 in place of tree.

result <- impute(
  traits, tree,
  joint_solver = "rphylopars",
  em_iterations = 2L,
  joint_refine_iter = 1L,
  conformal_split_val = TRUE
)

result <- impute(traits, tree, predict_method = "exact")

result <- impute(traits, tree,
                 phylo_signal_gate     = TRUE,   # opt-in (default FALSE)
                 safety_floor          = TRUE,   # enables the grand-mean fallback
                 phylo_signal_threshold = 0.2,   # default; min lambda to keep BM/GNN
                 phylo_signal_method   = "lambda")

result <- impute(traits, tree, phylo_signal_gate = FALSE)

# Step 1: GBIF occurrence centroids (median lat/lon across cleaned points)
gbif <- pull_gbif_centroids(species          = rownames(my_traits),
                            cache_dir        = "cache/gbif",
                            occurrence_limit = 500L,
                            store_points     = TRUE,
                            refresh_cache    = TRUE)
#> data.frame with columns species, centroid_lat, centroid_lon, n_occurrences

# Step 2: WorldClim bioclim summaries (median + IQR of all 19 bio variables
# across each species' GBIF occurrences)
clim <- pull_worldclim_per_species(species             = rownames(my_traits),
                                    gbif_cache_dir      = "cache/gbif",
                                    worldclim_cache_dir = "cache/worldclim",
                                    resolution          = "10m",
                                    refresh_cache       = TRUE)
#> data.frame with bio1..bio19 medians + IQRs + n_extracted per species

# Check the species order before fitting: covariate rows must match trait rows.
stopifnot(identical(clim$species, rownames(my_traits)))

# Step 3: keep the intended numeric climate summaries, not species or count
climate_cols <- grep("^bio[0-9]+_(median|iqr)$", names(clim), value = TRUE)
clim_cov <- clim[, climate_cols, drop = FALSE]
if (anyNA(clim_cov)) stop("Resolve missing climate values before fitting")
result <- impute(my_traits, my_tree, covariates = clim_cov, gnn = TRUE)