Skip to contents

Getting started with pigauto

Source: vignettes/getting-started.Rmd
getting-started.Rmd

Start here

The beginner workflow reads files, checks them before fitting, and keeps the completed-data, diagnostic-prediction, uncertainty, and inference roles separate.

library(pigauto)

traits <- read_traits("traits.csv")
tree <- read_tree("tree.nwk")
check_pigauto(traits, tree)
result <- impute(traits, tree)
completed <- completed_data(result)
pigauto_report(result)

check_pigauto() is fit-free. Resolve its structured errors before fitting; inspect check$species$matched, $data_only, and $tree_only when labels do not reconcile. completed preserves observed values and contains fills only for modeled missing cells. result$prediction is an all-cell diagnostic output; its uncertainty is type-dependent, and conformal bounds are nominal held-out diagnostics. A pigauto_result does not authorize downstream inference: use multi_impute_analysis() only in its documented narrow regime.

read_traits() requires unique species identifiers. For repeated observations, use read.csv() and supply species_col = "species" to both the check and impute().

Common declarations 

# Ordered ecological states
traits$threat <- ordered(traits$threat, levels = c("LC", "NT", "VU", "EN", "CR"))

# Integer-valued continuous measurement; integer otherwise means count
result <- impute(traits, tree, trait_types = c(length_mm = "continuous"))

# Numeric 0-1 proportion and zero-inflated count require explicit declarations
result <- impute(traits, tree, trait_types = c(cover = "proportion", parasites = "zi_count"))

# A composition is complete by row and its observed components sum to one
result <- impute(traits, tree,
  multi_proportion_groups = list(diet = c("diet_insect", "diet_fruit", "diet_seed")))

What pigauto does

Many comparative datasets have missing trait values: species were not measured, specimens were unavailable, or data were lost. pigauto fills in those gaps using three sources of information:

  1. The phylogenetic tree. Closely related species tend to share similar traits, so known values of relatives are informative about missing ones.
  2. Cross-trait correlations. If body mass predicts beak length, then observed masses help impute missing beak lengths even when the species itself was never measured for that trait.
  3. Environmental covariates (optional). Climate, habitat, or geographic variables can improve imputation when trait variation has a strong environmental component beyond phylogenetic signal.

The package handles continuous measurements, counts, binary variables, ordered categories, unordered categories, bounded proportions, zero-inflated counts, and compositional multi-proportion data — all in a single call to impute().

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 threshold/liability candidates for discrete-family columns, with a per-trait Pagel’s lambda. A held-out validation split calibrates the blend and the conformal intervals.

An optional graph neural network (GNN) correction is trained when you set gnn = TRUE. It learns additional patterns from tree topology, cross-trait structure and user covariates (which only the GNN uses). A per-trait gate controls how much the GNN contributes: when the baseline is already good enough, the gate closes and the GNN stays out of the way. The GNN has been off by default since version 0.11.0.9001 because, at the current baseline defaults, it did not lower imputation error on simulated or real benchmark data and took 20 to 90 times longer to fit.

Technical note. The internal class name ResidualPhyloDAE refers to the ResNet-style skip connections inside the GNN layers, not to a statistical residual y - baseline.

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.

Installation 

# CRAN release
install.packages("pigauto")

# Development version
pak::pak("itchyshin/pigauto")

# torch backend (required; ~1 GB first-time download)
torch::install_torch()

Inspecting the internal pipeline

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

result <- impute(traits, tree)

result$completed    # data.frame: observed values preserved, NAs filled
result$imputed_mask # logical matrix: TRUE where a cell was imputed
result$prediction$se # per-cell uncertainty (original units)

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, a matching pruned Hackett MCC phylogeny from BirdTree.org.

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 traits — body mass, beak length, tarsus length, and wing length — are all log-normally distributed and highly heritable, making them ideal test cases for phylogenetic imputation.

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 trait matrix (300 × 4)
  • means, sds: per-trait normalisation parameters (for back-transformation)
  • log_transform: TRUE, recorded so predict() can invert it automatically

Step 2: Create evaluation splits

We designate 25% of all matrix cells as “missing” for evaluation. These cells are not used during training or baseline fitting; they provide an honest estimate of imputation error.

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

cat("Total held-out cells:", length(splits$val_idx) + length(splits$test_idx), "\n")
#> Total held-out cells: 1049
cat("  Validation:", length(splits$val_idx), "\n")
#>   Validation: 262
cat("  Test:      ", length(splits$test_idx), "\n")
#>   Test:       787

The missing cells are split further: validation cells guide early stopping, test cells provide the final unbiased performance estimate.

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 the type-appropriate loss on the corrupted cells plus a shrinkage penalty on delta - baseline that pulls the GNN toward the baseline. Early stopping is based on held-out validation RMSE.

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
)

Example console output during training:

Epoch  100 | train loss: 0.8431 | val RMSE: 0.8852
Epoch  200 | train loss: 0.7215 | val RMSE: 0.8601
...
Epoch  900 | train loss: 0.6108 | val RMSE: 0.8289
Early stopping at epoch 900 (no improvement for 10 evals).
print(fit)
<pigauto_fit>
  Species  : 300
  Traits   : 4 -- Mass, Beak.Length_Culmen, Tarsus.Length, Wing.Length
  Architecture: hidden_dim = 64, k_eigen = 8, dropout = 0.1
  Val RMSE : 0.828 (z-score)
  Test RMSE: 0.831 (z-score)

Step 6: Predict

predict.pigauto_fit() runs a single forward pass through the trained model, optionally repeated with MC dropout to generate multiple stochastic completions for uncertainty quantification. The result is back-transformed to the original (non-log) scale.

pred <- predict(fit, return_se = TRUE)

# pred$imputed: 300 x 4 matrix in original units
# pred$se:      300 x 4 uncertainty matrix (original units)
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)
pred$conformal_coverage           # empirical coverage on val set (target ≈ 0.95)

Expected output (values will vary by seed and convergence):

       Mass Beak.Length_Culmen Tarsus.Length Wing.Length
t1  28.4150           14.2316       18.0923     136.421
t2  41.2031           18.5034       21.3441     154.832
...

Step 7: Evaluate — and read the gate

evaluate_imputation() scores held-out cells on the scale that matches the trait type. On high-signal morphometrics (AVONET-like mass, beak, tarsus, wing) the calibrated gate often closes (r_cal = 0). Then the shipped prediction is the phylogenetic baseline; the GNN is a no-op, not a guaranteed improvement.

# Held-out scores (z-score RMSE for continuous traits)
evaluate_imputation(pred$imputed_latent, 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

Training history 

plot(fit, type = "history")

This plots the validation RMSE over epochs, showing where early stopping triggered.

Uncertainty ribbons 

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

Species are sorted by predicted mass. The blue ribbon shows the conformal prediction interval when conformal scores are available. Grey points (if truth is supplied) show known values — a useful visual check of calibration.

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 validates only the documented narrow scope; the interface remains experimental.

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, multi_impute() with its default draws_method = "auto" uses posterior draws, which can be pooled.

Phylogenetic tree uncertainty

Tree uncertainty was outside the analysis-aware validation campaign. multi_impute_trees() may be used descriptively to examine point-prediction sensitivity, but its stochastic datasets are unsupported for downstream inference and cannot currently be combined with multi_impute_analysis().

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, tree)

# 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:

ΔV(snew)=σt2∑i∈misstDik2αk \Delta V(s_{\text{new}}) = \sigma_t^2 \sum_{i \in \mathrm{miss}_t} \frac{D_{ik}^2}{\alpha_k}

where D=Rmm−RmoRoo−1RomD = R_{mm} - R_{mo} R_{oo}^{-1} R_{om} is the residual matrix at currently-missing cells, αk=Dkk\alpha_k = D_{kk} is the current relative leverage of cell kk, and σt2\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 — the cross- metric ordering 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

Solver, calibration, EM, and refinement 

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

Exact conditional prediction route 

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

Phylogenetic-signal gating

impute() runs a per-trait Pagel’s-λ check before training and routes traits with very weak phylogenetic signal to the grand-mean corner of the safety-floor simplex (i.e., it predicts the column mean for those traits). This protects against over-fitting noise on traits where the phylogeny carries little information.

It is on by default since v0.9.1.9003. The two arguments you can tune:

result <- impute(traits, tree,
                 phylo_signal_gate     = TRUE,   # default
                 phylo_signal_threshold = 0.2,   # default; min lambda to keep BM/GNN
                 phylo_signal_method   = "lambda")

If a trait’s λ falls below phylo_signal_threshold, the calibrated gate routes it to the grand-mean corner. print(fit) reports which traits triggered the gate and their λ values.

To disable (e.g., when comparing against an older snapshot or for benchmarking purposes):

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

Assembling covariates from public data

Two helpers are bundled for the common case of building a covariate matrix from species-level environmental data:

# 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)
#> 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")
#> data.frame with bio1..bio19 medians + IQRs + n_extracted per species

# Step 3: hand the whole climate frame to impute() as covariates
result <- impute(my_traits, my_tree, covariates = clim, gnn = TRUE)

Caching is on by default — both helpers write to the cache_dir you supply and re-use existing extracts on subsequent runs. GBIF requires no API key. WorldClim raster downloads are ~50–500 MB depending on resolution and are cached in worldclim_cache_dir. Both functions error with a clear message if the optional dependencies (rgbif, terra) are not installed. Both helpers return one row per requested species; alignment is enforced by name where names are present, so pass species = and keep the returned species column intact rather than reordering it.

For a worked end-to-end example with covariates, use this section alongside ?impute, ?pull_gbif_centroids, and ?pull_worldclim_per_species. The old static covariate walk-through is being refreshed and is not the current source of truth.

Tips and next steps

Use GPU if available. pigauto automatically selects CUDA > MPS > CPU. Check availability with:

torch::cuda_is_available()       # NVIDIA GPU
torch::backends_mps_is_available() # Apple Silicon GPU

Cache the phylogenetic graph. For repeated experiments, save time by writing the graph to disk:

graph <- build_phylo_graph(tree300, k_eigen = 8, cache_path = "tree300_graph.rds")

Large trees (N > 2000). build_phylo_graph() computes a dense N × N cophenetic distance matrix and a Laplacian eigendecomposition. The eigendecomposition is O(N³) in the worst case, but k_eigen = "auto" caps the number of eigenvectors at 32 (the default for all trees), which keeps computation tractable in practice — pigauto has been tested up to 10,000 species on a standard laptop. The main cost at very large N is memory for the N × N distance matrix (~800 MB at N = 10,000). If memory is tight, prune to a clade of interest or watch for a future release that switches to a sparse Lanczos eigensolver.

Supplying your own data. The workflow above works for any phylogenetic tree and trait matrix:

tree   <- ape::read.tree("my_phylogeny.nwk")
traits <- read_traits("my_traits.csv", species_col = "species")
pd     <- preprocess_traits(traits, tree)
# ...proceed as above

Trait columns can be any of eight supported types. Five are auto-detected from R column class: numeric (continuous), integer (count), factor with 2 levels (binary), factor with >2 levels (categorical), and ordered (ordinal). Three require an explicit declaration because they cannot be inferred from class alone: numeric (0–1) with trait_types = "proportion" for bounded proportions; integer with trait_types = "zi_count" for zero-inflated counts; and K numeric columns declared via the multi_proportion_groups argument for compositional (simplex) data whose rows sum to 1. pigauto reads column classes to decide how each trait is encoded and imputed. The species names in traits must overlap with tip labels in tree (partial overlap is allowed; unmatched species are dropped with a warning).

References

  • Tobias, J.A. et al. (2022). AVONET: morphological, ecological and geographical data for all birds. Ecology Letters, 25, 581–597.
  • Jetz, W. et al. (2012). The global diversity of birds in space and time. Nature, 491, 444–448.
  • Goolsby, E.W., Bruggeman, J. & Ané, C. (2017). Rphylopars: fast multivariate phylogenetic comparative methods for missing data and within-species variation. Methods in Ecology and Evolution, 8, 22–27. (pigauto’s internal BM baseline implements the same conditional-MVN approach.)
  • Cohn, D.A., Ghahramani, Z. & Jordan, M.I. (1996). Active learning with statistical models. Journal of Artificial Intelligence Research, 4, 129–145. (Variance-reduction framework that suggest_next_observation() adapts to the phylogenetic setting.)