Getting started with pigauto

Start here

These examples use functions and defaults from pigauto 0.11.0. Check packageVersion("pigauto"): if CRAN has an earlier version, use the GitHub installation described below before following this workflow.

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() and completed_data() were added in 0.11.0. If your installed version is 0.10.0, call impute(traits, tree) directly; that version does not have them.

check_pigauto() is fit-free. Its default gnn = FALSE matches impute() and does not probe torch; pass gnn = TRUE when planning a GNN fit to include runtime availability. Resolve structured input errors before fitting, and 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().

# 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().

Installation

The CRAN release can lag behind the GitHub version described on this site. Check your installed version before running the examples above.

# CRAN release
install.packages("pigauto")
packageVersion("pigauto")

# GitHub version, if CRAN has an earlier version
# install.packages("pak") # once, if needed
pak::pak("itchyshin/pigauto")

# Install the torch runtime only if you plan to use gnn = TRUE
torch::install_torch()

library(pigauto)
library(ape)

Inspecting the internal pipeline

The workflow 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, posterior member 69 from the BirdTree Hackett-backbone sample, pruned to those 300 species. It is included as an example tree, not as an MCC or consensus tree or as a recommended tree for substantive analyses. For your own analysis, obtain a tree suited to your taxa and cite Jetz et al. (2012). 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 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)

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() 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-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. 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

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 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 and no covariates, multi_impute() with its default draws_method = "auto" uses posterior draws. Downstream pooling is supported for linear analyses of the imputed traits on their imputation scale; see the multiple imputation article for the conditions.

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:

\[ \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 — 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

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

Tips and next steps

Optional GNN: use a GPU if available. When you set gnn = TRUE, pigauto selects CUDA when available, then Apple MPS, and otherwise CPU. The default gnn = FALSE fits the phylogenetic baseline only and does not call torch. To check runtime availability before opting into GNN training:

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. build_phylo_graph() creates dense N × N distance, adjacency, Laplacian, and phylogenetic-correlation matrices. Each double- precision matrix needs about 800 MB at N = 10,000, and several are live during graph construction, so memory use can be substantial. For trees with more than 7,500 tips, the spectral step uses the optional RSpectra Lanczos solver when it is installed; if it is unavailable or does not converge, the function falls back to dense eigendecomposition. k_eigen = "auto" limits the number of returned spectral features to 32, but this does not reduce the storage needed for the dense matrices. For very large trees or limited memory, consider pruning to the clade relevant to your analysis.

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 should overlap with tip labels in tree. Partial overlap is allowed: unmatched input species remain in completed but are not modeled or filled. Tree tips with no input row can still be used during model fitting.

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.
• Hackett, S.J. et al. (2008). A phylogenomic study of birds reveals their evolutionary history. Science, 320(5884), 1763–1768. https://doi.org/10.1126/science.1157704
• 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.)