Mixed-Type Trait Imputation

Overview

Real comparative datasets contain many kinds of traits: body mass (continuous), clutch size (count), migratory status (yes/no), diet type (carnivore/herbivore/omnivore), threat status (LC < VU < EN). pigauto handles all of them in a single model: you do not need separate imputation runs for different column types.

Choose a workflow

Fill missing cells in one trait table: start with impute(). The GNN is off by default. Review the returned predictions and the limits for each trait type.
Inspect preprocessing, fit stages, or validation scores: follow this article’s lower-level workflow. Keep validation and test cells separate from predictions, and check each metric’s scale and sample size.
Carry missing-data uncertainty into a downstream model: start with the multiple imputation guide. Choose a route that matches your data and analysis, then check its eligibility and validation limits. Random variation among completed tables alone does not establish valid inference.
Compare predictions across posterior trees: use the tree sensitivity guide. This is a prediction-sensitivity check; tree-aware downstream inference is not supported by this workflow.

The mixed-type example below teaches the lower-level steps. For a first analysis with your own data, check the species names and trait declarations before fitting; the checklist below shows what to verify.

Type	R class	Example	Notes
Continuous	`numeric`	Body mass, wing length	Auto-detected
Count	`integer`	Clutch size, litter size	Auto-detected
Binary	`factor` (2 levels)	Migratory yes/no	Auto-detected
Categorical	`factor` (>2 levels)	Diet type, lifestyle	Auto-detected
Ordinal	`ordered`	Threat status (LC < VU < EN)	Auto-detected
Proportion	`numeric` (0–1)	Habitat cover, diet fraction	Requires a named override, e.g. `trait_types = c(cover = "proportion")`
ZI count	`integer` (zero-inflated)	Parasite load, rare behaviour	Requires a named override, e.g. `trait_types = c(parasite_load = "zi_count")`; experimental
Multi-proportion	K `numeric` columns summing to 1	Diet composition, plumage-colour fractions, microbiome relative abundances	Requires `multi_proportion_groups = list(<name> = c("col1", ..., "colK"))`; encoded via centred log-ratio (CLR) + per-component z-score

The first five rows in this table are generally inferred from R column class, subject to the exceptions below. proportion, zi_count, and multi_proportion must be declared explicitly (trait_types or multi_proportion_groups). Two-level ordered factors are detected as binary; integer columns with negative observed values are treated as continuous, while other integer columns are treated as counts. Use a named trait_types override when that inference does not match the measurement. All eight types share a latent representation, but type also determines baseline eligibility and dispatch, masking rules, loss, calibration, uncertainty summaries, and back-transformation. For examples, see the baseline section and prediction and decoding. Names in trait_types must match the corresponding trait columns.

Before fitting a mixed table, use check_pigauto(traits, tree). Its trait table records the detected type, declaration source, usable matched rows, and copy-ready corrections; its species report distinguishes matched, data-only, and tree-only identities.

For a real dataset, first confirm that row names (or the declared species column) match the tree tips. Set factor levels in the intended order, explicitly declare proportions and zero-inflated counts, and declare compositional groups with multi_proportion_groups. Then resolve any issues reported by check_pigauto() before fitting.

Synthetic example

This teaching example uses 60 simulated species. Validation evidence is described separately. The validation and test splits can contain only a few cells for some traits, so their scores and calibration are imprecise. For a real analysis, use enough observed values to support the validation split and inspect the reported n before interpreting a score.

library(pigauto)
library(ape)

set.seed(42)
n <- 60
tree <- rtree(n)

traits <- data.frame(
  row.names = tree$tip.label,
  mass      = exp(rnorm(n, 3, 0.5)),
  clutch    = as.integer(rpois(n, 3) + 1L),
  migr      = factor(sample(c("no", "yes"), n, replace = TRUE)),
  diet      = factor(sample(c("herb", "carn", "omni"), n, replace = TRUE)),
  threat    = ordered(sample(c("LC", "VU", "EN"), n, replace = TRUE),
                      levels = c("LC", "VU", "EN"))
)

Preprocessing

preprocess_traits() auto-detects column types from R classes:

pd <- preprocess_traits(traits, tree)
#> Auto-detected trait type "count" (log1p-encoded) for integer column(s): clutch. If any is a continuous measurement stored as whole numbers (common after read.csv), override it: trait_types = c(clutch = "continuous").
print(pd)
#> pigauto_data
#>   Species: 60 
#>   Traits:  5 
#>   Types:   binary=1, categorical=1, continuous=1, count=1, ordinal=1 
#>   Latent columns: 7

The trait_map records each trait’s type, levels, latent column indices, and normalisation parameters:

str(pd$trait_map, max.level = 1)
#> List of 5
#>  $ mass  :List of 9
#>  $ clutch:List of 9
#>  $ migr  :List of 8
#>  $ diet  :List of 8
#>  $ threat:List of 8
pd$trait_map$diet
#> $name
#> [1] "diet"
#> 
#> $type
#> [1] "categorical"
#> 
#> $n_latent
#> [1] 3
#> 
#> $latent_cols
#> [1] 4 5 6
#> 
#> $levels
#> [1] "carn" "herb" "omni"
#> 
#> $log_transform
#> [1] FALSE
#> 
#> $mean
#> [1] NA
#> 
#> $sd
#> [1] NA

Creating splits

When trait_map is supplied, make_missing_splits() operates at the original-trait level. For categorical traits, all K one-hot columns are held out together. A zero-inflated count uses two latent columns: a recorded zero has an observed zero/non-zero gate and no magnitude, while a genuinely missing count has no observed gate. The split uses the gate to identify observed cells; the missing magnitude for a recorded zero does not mean the raw count was missing.

spl <- make_missing_splits(pd$X_scaled, missing_frac = 0.20,
                           seed = 1, trait_map = pd$trait_map)
cat("Val cells (latent):", length(spl$val_idx), "\n")
#> Val cells (latent): 21
cat("Test cells (latent):", length(spl$test_idx), "\n")
#> Test cells (latent): 57

Baseline fitting

The baseline depends on trait type. Continuous-family latent columns use phylogenetic conditional multivariate-normal imputation. Discrete traits use phylogenetic label propagation or, for eligible combinations of traits, a threshold model. The architecture article describes the dispatch and its fallbacks:

bl <- fit_baseline(pd, tree, splits = spl)
dim(bl$mu)
#> [1] 60  7

Fitting

fit_pigauto() calibrates the phylogenetic baseline for each trait type on validation cells. It computes conformal scores for continuous, count, ordinal, and proportion traits, and for the positive-count magnitude of a zero-inflated count. That last score describes magnitude conditional on a non-zero count; it does not describe uncertainty about whether the count is zero. Binary, categorical, and multi-proportion traits do not receive these conformal scores. The lower-level fitter defaults to gnn = FALSE, so this example fits the baseline without constructing a network. Set gnn = TRUE to train the GNN correction as well. The example fixes both continuous and discrete Pagel’s lambda at 1 to keep this small illustration quick; it does not compare lambda settings. The impute() convenience function also defaults to the baseline-only route.

fit <- fit_pigauto(
  pd, tree,
  splits = spl,
  gnn = FALSE,
  lambda_mode = "fixed_1",
  discrete_lambda = "fixed_1",
  verbose = FALSE,
  seed = 1
)
#> Warning: Small validation set for 5 trait(s): mass (n=2), clutch (n=4), migr
#> (n=3), diet (n=3), threat (n=3). Calibrated gate and conformal scores will be
#> noisy for these trait(s). 95% split-conformal coverage is NOT achievable for 3
#> trait(s) with fewer than 19 validation cells (mass (n=2), clutch (n=4), threat
#> (n=3)): the achievable ceiling is n_val / (n_val + 1), which only reaches 0.95
#> at n_val >= 19. See `?fit_pigauto` under 'Calibration at small n' for smoothing
#> options.
print(fit)
#> pigauto_fit
#>   GNN     : off (baseline only)
#>   Species : 60 
#>   Traits  : 5 -- mass, clutch, migr, diet, threat 
#>   Types   : binary=1, categorical=1, continuous=1, count=1, ordinal=1 
#>   Architecture: hidden_dim = 64 | k_eigen = 4 
#>   Best val loss : 1.4259 
#>   Test loss     : 1.3562 
#>   Gate calibration: yes
#>   Conformal scores: 3 traits

Prediction and decoding

predict() decodes latent predictions back to original types:

pred <- predict(fit, return_se = TRUE)
head(pred$imputed)
#>          mass clutch migr diet threat
#> t58  9.518026      4  yes carn     LC
#> t8   9.629026      5  yes carn     LC
#> t36 21.377761      3  yes herb     VU
#> t4  12.202983      8  yes carn     LC
#> t22 20.067241      4   no herb     VU
#> t18 16.213927      4   no omni     EN

For binary and categorical traits, class probabilities are available:

# Binary: probability of "yes"
head(pred$probabilities$migr)
#>       t58        t8       t36        t4       t22       t18 
#> 0.7875313 0.7875313 0.7875313 0.7875313 0.2124687 0.2124687

# Categorical: probability of each diet class
head(pred$probabilities$diet)
#>          carn      herb      omni
#> t58 0.6495271 0.1752365 0.1752365
#> t8  0.6495271 0.1752365 0.1752365
#> t36 0.1752365 0.6495271 0.1752365
#> t4  0.6495271 0.1752365 0.1752365
#> t22 0.1752365 0.6495271 0.1752365
#> t18 0.1752365 0.1752365 0.6495271

The uncertainty matrix is type-dependent. Entries for binary and categorical traits are uncertainty scores, not Gaussian standard errors or measurement units. In this example, the ordinal threat column uses the threshold-joint baseline, which currently returns zero in pred$se; that zero is a placeholder, not evidence of certainty or an ordinal standard error. Continuous-family columns have scale-specific standard errors; see the uncertainty guide before using them:

Class probabilities (pred$probabilities): predicted probabilities for binary and categorical classes. Compare classes or report their probabilities.
Uncertainty scores (pred$se for binary or categorical traits): type-specific scores, not standard errors. Use them to rank predictions by uncertainty; do not form a normal interval from them.
Model-based standard errors (pred$se for supported continuous-family traits): values on a trait-specific scale. Check the scale and assumptions in the uncertainty guide; these are not conformal intervals.
Conformal bounds (pred$conformal_lower and pred$conformal_upper): held-out residual-calibrated prediction bounds for supported traits. They are not Rubin-pooled uncertainty for a downstream model.

head(pred$se)
#>     mass clutch      migr      diet threat
#> t58    0      0 0.2124687 0.3504729      0
#> t8     0      0 0.2124687 0.3504729      0
#> t36    0      0 0.2124687 0.3504729      0
#> t4     0      0 0.2124687 0.3504729      0
#> t22    0      0 0.2124687 0.3504729      0
#> t18    0      0 0.2124687 0.3504729      0

Evaluation

evaluate_imputation() dispatches type-specific metrics. Use evaluation = TRUE to keep validation and test values out of 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.

pred_eval <- predict(fit, evaluation = TRUE)
ev <- evaluate_imputation(pred_eval, pd$X_scaled, spl)
print(ev)
#>    split  trait        type  n      rmse   pearson_r coverage_95       mae
#> 1    val   mass  continuous  2 1.2367092 -1.00000000   0.5000000        NA
#> 2    val clutch       count  4 1.6710009 -0.44216026          NA 1.6452599
#> 3    val   migr      binary  3        NA          NA          NA        NA
#> 4    val   diet categorical  3        NA          NA          NA        NA
#> 5    val threat     ordinal  3 0.9688590          NA          NA        NA
#> 6   test   mass  continuous 13 1.2786925 -0.03154819   0.8461538        NA
#> 7   test clutch       count 10 0.9346293 -0.14660792          NA 0.8382015
#> 8   test   migr      binary 10        NA          NA          NA        NA
#> 9   test   diet categorical  6        NA          NA          NA        NA
#> 10  test threat     ordinal  6 1.4604038          NA          NA        NA
#>    spearman_rho accuracy     brier zero_accuracy aitchison rmse_clr simplex_mae
#> 1            NA       NA        NA            NA        NA       NA          NA
#> 2            NA       NA        NA            NA        NA       NA          NA
#> 3            NA      1.0 0.1481444            NA        NA       NA          NA
#> 4            NA      0.0        NA            NA        NA       NA          NA
#> 5     0.8660254       NA        NA            NA        NA       NA          NA
#> 6            NA       NA        NA            NA        NA       NA          NA
#> 7            NA       NA        NA            NA        NA       NA          NA
#> 8            NA      0.5 0.2481624            NA        NA       NA          NA
#> 9            NA      0.0        NA            NA        NA       NA          NA
#> 10   -0.8332381       NA        NA            NA        NA       NA          NA

Read split to distinguish validation from test results and n for the number of scored cells. For continuous and count traits, RMSE is on a transformed latent scale; for zero-inflated counts, RMSE and MAE are on the original count scale. Lower error scores and Brier scores indicate closer predictions; higher accuracy and correlation indicate better agreement. coverage_95 here uses a symmetric standard-error interval on a latent scale, not the separate conformal-bound coverage reported by evaluate(). NA means that metric is not available for that trait and split.

Experimental analysis-aware multiple imputation

multi_impute_analysis() requires the downstream formula before imputation. The initial backend supports exactly one incomplete continuous covariate under MAR. It does not yet impute the mixed discrete and compositional outcomes shown earlier in this vignette for downstream inference.

This code is a template, not a standalone example. analysis_data must contain a binary outcome y, one incomplete continuous predictor x, and a fully observed continuous predictor z; here model = "glm" specifies a binomial-logit model. For a runnable minimal example, see the multi_impute_analysis() function reference. The multiple-imputation guide explains eligibility and the current validation limits.

# Precompute nonlinear auxiliaries explicitly when required.
analysis_data$z_sq <- analysis_data$z^2
mi <- multi_impute_analysis(
  data = analysis_data, formula = y ~ x + z,
  missing = "x", model = "glm", m = 50L,
  auxiliary = "z_sq", seed = 1L
)

# Fit a downstream model to each
fits <- with_imputations(mi, function(d) {
  glm(y ~ x + z, data = d, family = binomial)
})

# Pool with Rubin's rules
pool_mi(fits)

pool_mi() returns a tidy data.frame with estimate, std.error, p.value, df (classical Rubin degrees of freedom by default, or Barnard-Rubin when df_fun supplies complete-data residual degrees of freedom), fmi (fraction of missing information), and riv (relative increase in variance) per coefficient.

The "lm" branch uses proper Bayesian normal-regression MI, "glm" uses smcfcs for a binomial-logit model, and "lmer" uses jomo.smc for one Gaussian random intercept. Fixed effects only: multiple incomplete columns, missing discrete covariates, MNAR, interactions/nonlinear terms involving the incomplete covariate, random slopes, variance components, correlations, BLUPs/conditional modes, and latent loadings are unsupported.

For prediction diagnostics, the lower-level predict(fit, n_imputations = 5L) interface returns stochastic datasets. These are not supported for downstream inference. The analysis-aware multi_impute_analysis() route admits its own narrow class of objects to with_imputations() and pool_mi(); the posterior route from multi_impute() supports the separately documented continuous-trait analyses.

The analysis-aware backend is implemented, but implementation alone does not establish validated inferential performance. In the registered tests, 23 of 24 method-by-scenario results were indeterminate and one binomial-logit GLM result failed; none was classified as supported, and 12 controls used complete data. See the registered inferential results and the multiple-imputation guide before choosing an inferential route.