The plot.pigauto_fit() example now follows the package default and fits without
the optional GNN runtime, while still demonstrating the gate and conformal plots.
The compare_methods() example now checks that torch and its Lantern runtime
can create a tensor before fitting the GNN. It reports a skip when the runtime
is unavailable, while compare_methods() keeps its documented
gnn = TRUE default.
The README and tutorials now identify tree300 as a posterior sample tree and
describe how to obtain and cite BirdTree trees. The README clarifies that
BirdTree distributes Newick files and read_tree() also accepts NEXUS files.
The installed NOTICE
distinguishes the optional megatrees package licence from rights to the
underlying data, and the bundled trees300 help documents its Ericson/Hackett
mixture. The retired simulation study remains in the repository's history but
is excluded from the public site, sitemap, and search index.
The package description no longer makes an unqualified claim that the current release has been tested at 10,000 species. The historical scaling note below now distinguishes the completed 9,993-species missingness study from the v0.9.0 timing run, which stopped at 5,000 tips.
Two defaults change, and results change for anyone using the defaults.
First, binary, ordinal and categorical traits now get their own Pagel's
lambda in the baseline (on the liability scale; categorical traits through
their one-vs-rest fits). Before, categorical traits were held at lambda = 1 and
binary and ordinal traits at 1, or at the joint lambda shared with the
continuous traits on the "exact" route. The new argument
discrete_lambda = c("estimate", "fixed_1") (on impute(), fit_pigauto()
and fit_baseline(), and through ... on multi_impute() and
multi_impute_trees()) defaults to "estimate" and is recorded in the fit's
model_config. Binary and ordinal lambdas appear in $lambda_per_trait; the
categorical one-vs-rest lambdas do not. The zero-inflated-count gate column
also gets its own lambda.
Second, safety_floor and phylo_signal_gate now default to FALSE in
impute() and fit_pigauto(). With lambda estimated, the gate and the floor
did not improve accuracy in any setting tested: at Pagel's lambda of 0 to 0.1
the old defaults were worse than predicting the mean (continuous zRMSE 1.01
against 0.89 to 0.91 with them off), and on the BIEN plant data (4,745
species) that the floor was built for, turning them off changed error by
-0.002 (SE 0.001). Fits
are also 3 to 6 times faster at n = 1,000.
To restore the previous behaviour, pass
discrete_lambda = "fixed_1", safety_floor = TRUE, phylo_signal_gate = TRUE.
Evidence (docs/dev-log/discrete-lambda/README.md; simulated Brownian-motion
and OU data, n = 100 to 1,000, plus AVONET300 and BIEN): estimating discrete
lambda alone added about 0.01 to 0.04 discrete accuracy at Pagel's lambda of 0.3 or
below; together with the gate and floor off the gain was 0.04 to 0.12, and on
the four-arm study's datasets binary and categorical accuracy matched or beat
BACE's. Costs: with lambda near 1 and n of 1,000 or more, accuracy -0.002 to
-0.005 and Brier score +0.010 to +0.016; on AVONET300's two categorical
traits, Brier +0.014 to +0.020 with accuracy unchanged or better. Ordinal
accuracy still trails BACE's at Pagel's lambda of 0.3 (0.38 against 0.44 at
n = 300).
?impute, ?multi_impute, ?multi_impute_trees, ?preprocess_traits and the
getting-started vignette now say that the default log_transform = TRUE logs
every all-positive continuous trait, so traits that are already on a log scale
are logged twice. Set log_transform = FALSE for such data. In a real-data
check on PanTHERIA (4,027 species, four already-logged traits, posterior
draws), the default shifted Rubin-pooled regression slopes by up to 3.9
standard errors of the complete-data slope, in a direction that depended on
which values were missing; with FALSE the largest shift was 1.8
(docs/dev-log/mi-posterior/pantheria_logtf/). No behaviour changes.
plot_comparison() now defaults to methods = c("baseline", "pigauto"),
the labels that evaluate() and compare_methods() return. The previous
default (BM_baseline, pigauto_GNN) made
plot_comparison(compare_methods(...)) error with no matching rows.
multi_impute(draws_method = "posterior") now puts a separation-strategy
prior on the residual covariance: half-Cauchy residual standard deviations
(scaled by each trait's observed SD) and a uniform (LKJ(1)) residual
correlation. Posterior draws change. To reproduce earlier results, pass
posterior_control = list(residual_prior = "iw").
Why. The inverse-Wishart prior kept the residual covariance away from 0
and its correlation near 0 when traits are close to Brownian motion, which
attenuated a downstream relationship between correlated traits. In the
Rubin study simulation (continuous traits, 30% MCAR, 200 datasets per cell;
docs/dev-log/arc/2026-10-02-rubin-pigauto-campaign.md on branch
arc/rubin-freq-bace) at lambda = 1 and trait correlation 0.5, PGLS slope
coverage was:
| n | 100 | 300 | 1000 |
|---|---|---|---|
| "iw" (old default) | 0.955 | 0.875 | 0.750 |
| "sep" (new default) | 0.955 | 0.940 | 0.965 |
| frequentist MI | 0.955 | 0.920 | 0.935 |
Over all 18 cells (n 100 to 1000, lambda 0.3 to 1, correlation 0 or 0.5)
"sep" had 0 sampler errors and 0 unconverged fits in 3,600 ("iw": 12
and 107), downstream coverage 0.925 to 0.975 and per-cell coverage 0.942 to
0.959. The #189 acceptance simulation (8,000 cells) passes both gates under
"sep" (docs/dev-log/mi-posterior/sep_validation/). Costs: at n = 100
with correlated traits the pooled slope is biased by about -0.02 (coverage
0.96), and fits take about 5% more compute.
multi_impute() now defaults to draws_method = "auto". It uses
"posterior" when every trait is continuous, there is one row per species,
and no covariates or multi_proportion_groups are given; these are proper
multiple imputations that with_imputations() and pool_mi() accept.
Otherwise it uses "conformal" and prints one line saying why and that
those draws cannot be pooled. result$draws_method records the choice. To
get the previous default, set draws_method = "conformal".
Why. With the old default, a user who called multi_impute() and followed
the documented with_imputations() / pool_mi() route was stopped at the
last step, because conformal draws are refused for inference: in a
16-regime simulation of a phylogenetic regression slope they biased the
pooled slope towards zero by 0.20 to 0.46 with 0 to 17% interval coverage
(see the draws_method documentation). Note that the posterior sampler is
slower than conformal draws: about ten minutes per fit for 300 species and
four traits on one core.
impute(), fit_pigauto(), multi_impute() and multi_impute_trees() now
default to gnn = FALSE: they fit the calibrated phylogenetic baseline,
with its conformal intervals, and train no graph neural network. Results
from code that relied on the old default will change. To restore the
previous behaviour, pass gnn = TRUE.
Why. A rerun of the with/without-GNN campaign at the current baseline
defaults (docs/dev-log/arc/2026-10-02-campaign-gnn-rerun.md; Totoro,
200 cells, 20 seeds, single-observation data, 30% missing at random) found:
What to know.
covariates are used only by the GNN. With gnn = FALSE they are
ignored, and the warning now says to set gnn = TRUE to use them.conformal_method = "mondrian" needs the GNN; its error now says so.multi_impute(draws_method = "mc_dropout") and the
multi_impute_trees() default draw from the baseline posterior when
gnn = FALSE, with a message.compare_methods() and simulate_benchmark() gain a gnn argument
that defaults to TRUE, so they keep comparing the baseline with the
GNN-corrected model.gnn = FALSE; torch is still in Imports.vignettes/multiple-imputation.Rmd): which draws_method can be pooled,
and worked multi_impute(draws_method = "posterior") ->
with_imputations() -> pool_mi() examples for nlme, lme4, glmmTMB,
drmTMB and gllvmTMB, plus brms and MCMCglmm by posterior concatenation.pool_mi() now refuses brmsfit objects, as it already did for
MCMCglmm, and points to brms::brm_multiple(). Previously a brmsfit
fell through to coef() and vcov(), and coef() on a brmsfit does not
return its fixed effects.script/cran-0.11-integration/, outside the built package.pool_mi() now requires recognized drmTMB and gllvmTMB fits to report
convergence code zero and a positive-definite Hessian, including when
custom extractors are supplied. Workflow markers are caller-controlled
metadata and do not authenticate a fit list.fit_pigauto()'s gate_method default is
"cv_folds" and min_val_cells defaults to 20; the help pages said
"single_split" and 10.docs/dev-log/defaults-audit/2026-10-02-defaults-audit.md.multi_impute() gains draws_method = "posterior", which returns proper
Bayesian multiple imputations for continuous traits. It fits the
multivariate phylogenetic mixed model
vec(Y) ~ N(1 mu', Sigma_P %x% R + Sigma_E %x% I_n) on pigauto's latent
scale, with full phylogenetic and residual covariance matrices (so a separate
Pagel's lambda per trait) and MCMCglmm-style parameter-expanded priors. Each
of the m completed datasets is a draw from the posterior predictive
distribution of the missing cells, so imputation uncertainty is correlated
across tips and traits and includes parameter uncertainty. The sampler draws
the phylogenetic effects, trait means and missing cells as one block with a
sparse Cholesky factor of the Hadfield and Nakagawa (2010) node precision, so
cost grows roughly linearly in the number of species. Metropolis moves with
the effects integrated out keep it mixing near lambda = 0 and lambda = 1.
posterior_control (chains, burn-in, iterations, thinning,
keep_draws, param_uncertainty, seed, auto_extend, max_extend).mi$posterior holds per-cell 95% posterior predictive intervals
(cell_interval, computed from all kept draws, not from the m
datasets), convergence diagnostics (diagnostics: rank-normalised split
R-hat and bulk ESS; attribute "converged"), and the parameter draws
(params). A warning is raised if the chains have not converged after
any automatic extensions (next item).n_iter sweeps, up to posterior_control$max_extend times (default 3, so
at most 4 times the default length), and the diagnostics are recomputed on
all sweeps after burn-in. Burn-in is not repeated and the Metropolis step
sizes stay as tuned in burn-in, so an extended chain is the same chain as
one run longer from the start. The number of kept draws does not change:
after k extensions they are every thin * (k + 1)-th sweep.
posterior_control$auto_extend = FALSE turns this off. A fit that meets
the rule first time is not extended; its draws and diagnostics are the
same as with auto_extend = FALSE.
mi$posterior$diagnostics records the extensions made (attribute
"n_extensions") and the sweeps after burn-in per chain
("sweeps_per_chain"); print() and the non-convergence warning report
the extensions.with_imputations(), and the fits it returns
to pool_mi(); both now accept the provenance marker
"pigauto_posterior_mi_v1" alongside "pigauto_analysis_mi_v1". The
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 for
analyses that are linear in the imputed traits on that scale and whose
variables are all among the imputed traits. Not covered: external
covariates, because the imputation model does not contain them; nonlinear
terms or interactions among imputed traits; and analysing a
log-transformed trait on its raw scale.posterior_control$param_uncertainty = "none" (plug-in draws with the
covariance matrices fixed at their posterior means, for validation only)
returns the marker "pigauto_posterior_plugin_diagnostic", which
with_imputations() and pool_mi() refuse.species_col), covariates, and input with no missing cells in the rows
of traits (tree tips absent from traits do not count). The GNN is
not used, and fitting arguments such as gnn and epochs are ignored
with a message.When this posterior method was introduced, the default draws method remained
"conformal". The 0.11.0 default is now "auto", which selects posterior
draws for eligible all-continuous data and otherwise falls back to conformal
prediction-diagnostic draws.
Caveat for downstream inference with the other draws methods. The
"conformal" and "mc_dropout" draws perturb missing cells around a point
prediction rather than drawing them jointly from their conditional
distribution given the observed data. In a 16-regime simulation of a
downstream phylogenetic regression (PGLS slope; branch
arc/mi-gls-attenuation, docs/dev-log/mi-gls/results.md), conformal
draws biased the pooled slope by -0.20 to -0.46, and the pooled 95%
intervals covered the truth in 0 to 17% of replicates, in all 16 regimes;
MC-dropout draws biased it by -0.03 to -0.38. For downstream inference on
continuous traits, use draws_method = "posterior".
impute(), fit_pigauto() and fit_baseline() now default predict_method to "auto" (previously
"per_column"). "auto" fits the phylogenetic baseline two ways and keeps, for each trait, the one with the
lower error on that trait's validation cells:
"exact": the full cross-trait conditional of vec(L) ~ MVN(0, Sigma x R(lambda_block)), solved in sparse
precision form (Hadfield & Nakagawa, 2010), with each trait centred at its phylogenetic GLS mean. It uses
the other traits' observed values to predict a missing one."per_column": each trait predicted from its own observed values only (the previous default).Validation error is squared error on the latent scale for continuous, count, proportion, ordinal and
zero-inflated magnitude traits, and log-loss for binary, zero-inflated gate and categorical traits. The unit of
validation is the held-out trait row (the species, in data with several observations per species). A trait
with at least 38 held-out rows, whose two candidate fits differ, splits its rows in two: one half chooses the
route and the other calibrates the gate and the conformal intervals, so no row does both. With fewer rows the
same rows do both, because a 95% conformal interval needs at least 19 calibration rows; when the two fits
agree (for example a single trait) no choice is made and every row calibrates. A trait with fewer than 5
validation cells, or a fit with no validation split, uses "exact". The choice is recorded in
fit$model_config$predict_method_by_trait, and the final refit on all data reuses it (new predict_route
argument of fit_baseline()). "exact" and "per_column" remain available as explicit choices. In two
timings the "auto" baseline took about 1.4 to 1.5 times as long as a single baseline fit.
Measured against the previous default with the GNN off (18 core simulation cells, 200 seeds each, results
averaged over two cross-trait correlations): mean z-RMSE of the continuous-family traits falls by 3.1 to 6.6%
at true lambda 0.3 and by 6.1 to 9.4% at lambda 0.7. At lambda 1 it falls by 1.2% (100 species) and 0.3%
(300) and rises by 0.4% at 1,000 species. Mean discrete accuracy rises by 0.009 to 0.032 at lambda 0.3 and
0.7 and is unchanged (within 0.001) at lambda 1; no single scenario falls by more than 0.004. Interval
coverage of the continuous traits rises by 0.003 to 0.010 at lambda 0.3 and 0.7. At lambda 1 with 100
species, where the previous default already covered only 0.86 to 0.92 of cells, coverage falls by a further
0.006 to 0.018 (largest for the count and proportion traits); with 300 species it falls by at most 0.008.
The opt-in GNN-on arm of impute() was not re-benchmarked for this change.
On 13 real-data cases (continuous traits, 5 seeds each, GNN off, no Monte Carlo error computed), z-RMSE falls
by 17 to 51% on AVONET and PanTHERIA and by 2.5% on AmphiBIO (2,000 species); the other seven cases change by
-0.9 to +0.8%, the largest rise on the full GlobTherm set (1,969 species). Evidence:
docs/dev-log/exact-default/.
Details:
"exact", discrete traits (binary, zero-inflated gate, ordinal) share the block lambda used for the
whole joint fit instead of staying at lambda = 1; under "per_column" they still stay at 1."exact" falls back to "per_column" above about 20,000 unknown cells (roughly 4,000 species with 5
traits), for a numerically unusable Sigma, or with fewer than two joint traits. Under the default this is
reported once per session with message(); an explicit predict_method = "exact" still warns.
model_config$predict_method_used records the route that ran.em_iterations >= 1 silently ignored predict_method inside the threshold-joint and categorical EM
loops; both now pass it through.fit_baseline(..., lambda_fixed = )) now reproduces the original
fit under "exact", because the block lambda is stored with the per-trait values.The route choice scores ordinal traits by squared error on the latent scale. Two alternatives, class error
rate and the log-probability of the true class, were tested on the 18 core simulation cells and both lowered
ordinal accuracy (on average by 0.002 and 0.019), so neither was adopted (docs/dev-log/ordinal-route/).
For each ordinal trait, the threshold-joint baseline picks between the threshold-joint fit and a per-column BM
fit on validation error. A third candidate, label propagation (added in May 2026), decoded the 0 to K-1 class
coding as 1 to K, so it ignored every observation in the lowest class. It has been removed. On the 18 core
simulation cells (GNN off, 200 seeds) ordinal accuracy changes by -0.002 to +0.007 (mean +0.001) and macro F1
by -0.009 to +0.002 (mean -0.004); no other trait changes. A version with the decoding fixed lowered ordinal
accuracy (mean -0.003) and was not kept. Evidence: docs/dev-log/ordinal-route/.
The profile likelihood used the REML variance divisor but dropped the REML
determinant term 0.5 * log(1' R(lambda)^-1 1) (0.5 * log|X' R^-1 X|
with covariates), which biased the lambda estimate low at weak signal.
Measured at n = 300, 30% missing, 200 seeds: bias -0.061 -> -0.022 at true
lambda 0.3, and -0.038 -> -0.027 at 0.7; unchanged at lambda = 1. Prediction
error changes little (knowing the true lambda improves it by about 0.3%),
so this mainly makes the reported lambda_per_trait trustworthy.
lambda_mode = "fixed_1" is unaffected.
impute(), fit_pigauto(), and multi_impute() now default lambda_mode
to "estimate" (previously "fixed_1"); "fixed_1", "cv", and "bayes"
remain available. multi_impute() gains an explicit lambda_mode argument
(previously only reachable via ...). multi_impute_trees()'s per-tree
baseline refit now falls back to "estimate" (not "fixed_1") when the
reference fit's model_config$lambda_mode is unset, so tree-specific
baselines estimate their own per-tree Pagel's lambda by default.
Under lambda_mode = "estimate", each continuous-family (BM-eligible)
latent column gets its own per-trait Pagel's lambda via the existing
per-column ML estimator (bm_impute_col(lambda = "estimate")); discrete
traits (binary, categorical, zi gate, OVR synthetic columns) AND ordinal
stay at lambda = 1, unaffected by lambda_mode (see fit_baseline()'s
"Per-type lambda dispatch" details).
This is a real change to default baseline numbers, not just a relabelling.
Measured max abs difference in baseline mu of ~0.38 (z-scored latent
scale) between lambda_mode = "fixed_1" and "estimate" on a 2-trait,
40-tip smoke fixture, with the joint MVN baseline active on both sides
(bl$path reports "joint_mvn" for both) -- the joint solver itself now
estimates lambda, rather than falling back to the per-column path. See
the joint-solver lane reports under docs/dev-log/lambda-default/ for
the estimation mechanism. Binary/ordinal/categorical/zi-gate columns are
unaffected: still threshold-joint / OVR, still lambda = 1 in every path.
The opt-in predict_method = "exact" / joint_refine_iter > 0 paths use
the shared lambda_block only for cross-trait computations that need one
common phylogenetic correlation matrix; they never override an individual
discrete or ordinal column's own lambda = 1.
fit$model_config gains lambda_per_trait (named numeric per latent
column) and lambda_block (a single shared lambda used by the joint/exact
prediction paths, i.e. predict_method = "exact" or joint_refine_iter > 0, which need one common phylogenetic correlation matrix across traits
rather than per-trait values), sourced from fit_baseline()'s own
$lambda_per_trait / $lambda_block. fit_baseline() also gains
lambda_fixed (a named numeric vector), letting a caller rebuild a
baseline at previously-estimated per-trait lambda values instead of
re-running ML estimation -- the mechanism a future predict-time rebuild
would use. Columns not named in lambda_fixed default to lambda = 1, on
both the per-column and the joint MVN / threshold-joint paths.
Covariate-aware baselines now estimate lambda too. When
data$covariates is supplied, the per-column covariate path
(bm_impute_col_with_cov()) accepts a numeric lambda or "estimate", so
lambda_mode %in% c("estimate", "fixed_1") reaches it and each
covariate-aware BM-eligible column gets its own estimated lambda.
"cv" / "bayes" are not supported there and fall back to lambda = 1
for BM-eligible columns, with a warning.
joint_solver = "rphylopars" now runs model = "lambda" under
lambda_mode = "estimate" (previously always model = "BM"), letting
Rphylopars estimate its own shared lambda instead of pigauto's in-house
per-trait estimator. Measured about 26x slower than model = "BM" on the
lane's benchmark datasets, with occasional explosive tip predictions on
small or weak-signal data; a plausibility guard (tip predictions more
than 10x the observed range) falls back to the in-house solver with a
warning when that happens. The same guard also now fires under
lambda_mode = "fixed_1" (model = "BM"), not just under "estimate".
cross_validate(), compare_methods(), simulate_benchmark(), and
multi_impute_trees() all inherit the new "estimate" default, because
none of them pass an explicit lambda_mode to fit_pigauto() /
fit_baseline(). Known gap: simulate_benchmark() fits its own baseline
directly (fit_baseline(pd, tree, splits = spl), always "estimate")
and then passes it to fit_pigauto(..., baseline = bl, ...); a
lambda_mode supplied through simulate_benchmark(...)'s ... reaches
fit_pigauto() and is recorded in fit$model_config$lambda_mode, but
never reaches the already-fitted bl, so the recorded lambda_mode can
disagree with the baseline that was actually used.
Known limitation. suggest_next_observation() still scores candidate
observations' expected variance/entropy reduction using the tree's raw
phylogenetic correlation matrix (fit$graph$R_phy, lambda = 1); it does
not read the fit's own model_config$lambda_per_trait. Under the new
default the imputation baseline itself estimates lambda, so this
sampling-design helper and the baseline it is meant to complement can now
disagree on the phylogenetic correlation structure. Not changed in this
release.
Local candidate after CRAN pigauto 0.10.0 (Date/Publication 2026-07-30).
No CRAN submission or public release is implied. The installed BACE bridge has
been removed for 0.11: BACE remains an in-tree comparator only, with no
Suggests dependency, export, or installed help page.
The pre-registered real-data confirmation of conformal_method = "mondrian"
(docs/dev-log/mondrian-realdata/) ran on PanTHERIA (n = 4,027; random and
phylogenetically structured masks, 3 each), AVONET (n = 1,500; random mask, 3)
and FishBase (n = 10,484; structured mask, 1).
Rule verdict for the mondrian default: KEEP_SPLIT. The default remains "split".
What the data showed, in every dataset and arm: Mondrian raises coverage in the
far (undersampled) stratum from about 0.92-0.93 to 0.94-0.96 by widening those
intervals, and narrows near-stratum intervals by about 20%, lowering near
coverage from 0.97-0.98 to 0.95-0.97 (pooled over traits). The pre-registered
rule failed on its near-stratum non-inferiority condition (Mondrian near
coverage no more than 2 points below split's), which failed for AVONET (pooled
near coverage 2.3 points below split's) and could not be demonstrated for
FishBase from its single mask (1.6 points below, p = 0.22); PanTHERIA passed.
"mondrian" remains available and opt-in. In the far stratum it raised
coverage in all 18 trait-by-arm rows at the cost of wider intervals there, and
it lowered near-stratum coverage by 1 to 2 points, in six trait rows to below
0.95. Users whose missing cells sit in undersampled clades may prefer that
trade; the study measured it on masked observed cells, not on real missing
cells.
compute_conformal_scores() now records n_val, n_near and n_far for each
Mondrian trait, and a verbose fit names the trait and its stratum sizes when
Mondrian falls back to the split score.
multi_impute(draws_method = "conformal") now draws each missing cell with its
own Mondrian stratum score when the fit used conformal_method = "mondrian"
(the per-cell scores come from the new internal mondrian_cell_scores(), also
used by predict(), whose output is unchanged). Split fits are unchanged.
Caveat, measured on simulated bivariate Brownian traits (500 replicates,
n = 1000): conformal draws, split or Mondrian, halved a phylogenetic GLS slope
(0.35 and 0.32 against a true 0.70) while OLS on the same draws was nearly
unbiased. A one-tree diagnostic points to two causes: the draw spread is
calibrated on point-prediction error, which was 1.3 times the proper
conditional spread there, and the default prediction route does not use the
other observed traits. Treat conformal MI followed by a phylogenetic
GLS with caution until this is resolved; an investigation is under way on
arc/mi-gls-attenuation. See useful/mondrian-mi-se-justification.md.
impute(), fit_pigauto(), multi_impute() and multi_impute_trees() gain
gnn = TRUE. With gnn = FALSE the GNN is never constructed or trained and
the fit, prediction and every downstream call make no torch:: call at all
(no device probe, no tensors, no GPU): the latent multivariate phylogenetic
baseline (joint MVN / threshold-joint / OVR, or per-column BM and label
propagation) is what runs through preprocessing, calibration, conformal
scores, prediction, completed data, multiple imputation, cross_validate(),
simulate_benchmark() and evaluate(). On the 40-tip test fixture
impute(gnn = FALSE) completes in about one second.
Semantics, chosen so that the GNN-on and GNN-off arms differ only in the GNN term:
r_cal_gnn is then folded into r_cal_bm.
safety_floor and phylo_signal_gate therefore behave identically in both
arms. The pure traditional-statistics arm is
gnn = FALSE, safety_floor = FALSE, phylo_signal_gate = FALSE.fit$baseline is the val/test-masked fit and is
what evaluate(), cross_validate() and every other scorer read, so
held-out metrics stay leakage-free. fit$baseline_full is a refit with no
held-out cells (splits = NULL); predict() uses it only in production
mode (no baseline_override, no .mask_observed_idx), so completed data
no longer pay the held-out-cell cost that gate calibration imposes. Under
gnn = TRUE predictions are unchanged; the fit object gains three
inert slots (baseline_full = NULL, model_config$gnn = TRUE,
baseline$path).baseline_full, so it is conservative
when the BM model is right (inference, not a guarantee) and should be
measured on user-masked cells rather than through evaluate().draws_method = "mc_dropout" under gnn = FALSE draws from the BM
posterior (no dropout exists); one message says so, no error.multi_impute_trees(gnn = FALSE) fits each per-tree baseline without
held-out cells, and the per-tree baseline replay now also carries
joint_solver, predict_method and joint_refine_iter (pre-existing
omission).fit_baseline() now returns path, a named character vector recording
which dispatch produced each trait's baseline (joint_mvn,
threshold_joint, ovr_categorical, per_column_bm,
multi_proportion_bm, label_propagation, zi_gate_lp,
zi_mag_constant).summary()/print() show GNN: off (baseline only) (and no architecture
line); plot(fit, "history") explains that no history exists for a
gnn = FALSE fit; pigauto_report() labels the non-BM column
Blend (GNN off) and shows no BM-vs-GNN verdict.missing_frac = 0, or fit_pigauto(splits = NULL)) the GNN-off fit is pure baseline (r_cal_bm = 1) with no
conformal scores, and still predicts.covariates enter pigauto only through the GNN; under gnn = FALSE
they are ignored, with a warning.Scope note: for continuous-only data at lambda_mode = "fixed_1" the
GNN-off arm is a phylopars-style joint BM fit by construction, so a
comparison against raw Rphylopars there measures the pipeline around the
solver, not a different method. The arm differs from Rphylopars for mixed
types, other lambda_mode values, multi-proportion and multi-observation
data, and in its intervals. This entry makes no performance claim; the
with/without-GNN comparison against Rphylopars and BACE is a separate,
pre-registered campaign.
check_pigauto() is now the documented fit-free preflight before
impute(), and completed_data() is the documented completed-data
extractor. The beginner journey reads traits and a tree, checks inputs,
imputes, extracts completed data, and writes a report.gnn = FALSE, matching impute() and
skipping the torch runtime probe; pass gnn = TRUE to check runtime
availability when planning a GNN fit.All chronological entries below are historical unless this 0.11 section says
otherwise. In particular, older instructions to pass multi_impute() or
multi_impute_trees() completions to with_imputations() or pool_mi(), and
older claims of Rubin-valid downstream standard errors for those diagnostic or
tree-sensitivity completions, are superseded. The current supported downstream
route begins with multi_impute_analysis() in its documented narrow regime.
Older chronological statements about broad interval coverage, active-selection
optimality or novelty, and solver performance comparisons are historical
context, not current package-wide claims; use the current README and vignettes
for user guidance.
fit_baseline() previously set force_per_column <- TRUE for
lambda_mode %in% c("estimate", "cv", "bayes"), which disabled the joint
MVN, threshold-joint, AND OVR-categorical baselines for ALL traits --
including binary, ordinal, categorical, and zero-inflated gate columns
that have no lambda concept and were always fit at lambda = 1 regardless.
This was a dispatch-condition defect: affected discrete columns could take an
unintended per-column fallback when a continuous-family lambda mode was
selected. The correction restores the intended type-specific baseline routing;
it does not make a comparative performance claim.
lambda_mode now governs ONLY the baseline for continuous-family columns
(continuous, count, ordinal, proportion, zi_count magnitude); the
threshold-joint and OVR-categorical baselines for binary/ordinal/
categorical/zi_gate columns fire regardless of lambda_mode, at
lambda = 1 as they always have. When the threshold-joint baseline runs
for lambda_mode != "fixed_1", its joint liability fit still uses the
continuous-family columns internally (they inform the joint Sigma the
binary/ordinal posteriors condition on) -- only the fit's continuous-column
baseline output is discarded in favour of the lambda-aware per-column
bm_impute_col(..., lambda = bm_lambda) fit. lambda_mode = "fixed_1"
output is unchanged (see tests/testthat/test-lambda-per-type.R).
fit_baseline(), fit_pigauto(), and impute() gain a new joint_solver
argument, c("inhouse", "rphylopars"). The default, "inhouse", is
byte-identical to prior releases. "rphylopars" routes the joint MVN,
threshold-joint, and OVR categorical baselines through Rphylopars::phylopars()
instead, whose output contract ($anc_recon / $anc_var / $pars$phylocov)
the in-house solver was built to match. If the phylopars call errors or
produces a non-finite tip prediction, pigauto warns and falls back to the
in-house solver automatically -- phylopars' own solve(): system is singular
warnings are not themselves treated as failure, only a hard error or
non-finite output is.
Also fixed: the fit_joint_threshold_baseline() / fit_joint_mvn_baseline()
roxygen previously said the joint baseline ran "via Rphylopars::phylopars()"
while the code actually called the in-house solver; that stale claim (and
the matching one in impute()'s em_iterations docs) is corrected.
fit_mvn_bm_inhouse() and fit_joint_solver() (R/joint_mvn_solver.R)
gain a sigma_method = c("single_pass", "fisher_ml") argument. The
default, "single_pass", is byte-identical to prior releases.
"fisher_ml" ports (adapted, not cherry-picked) the never-merged
commit e7ca41c: it estimates Sigma by maximising the observed-data
marginal MVN log-likelihood over a Cholesky-parameterised Sigma
(species treated as iid; phylogeny re-enters through the unchanged
per-column-BM init and cross-trait EM refinement), as an alternative
to the default's R-credited Kronecker M-step. Falls back to the
single_pass Sigma estimate, with a warning, if optim() errors or
fails to converge. Not yet threaded through fit_baseline() /
impute() -- exposure is gated on a recovery-simulation check
(dev/sigma_recovery_sim.R).
fit_pigauto() / impute() gain a new conformal_method option,
"mondrian", alongside the existing "split" (default) and "bootstrap".
It targets the exchangeability failure diagnosed in the mechanism-coverage
campaign (docs/dev-log/2026-08-16-mechanism-coverage-results.md): a single
global conformal quantile per trait is calibrated on the observed
complement, which sits systematically closer (in cophenetic distance) to
other observed cells than genuinely-missing cells do -- worst under
clade-structured missingness (MAR_phylo, -3.4pp coverage vs MCAR at
n=300/1000).
"mondrian" conditions calibration on phylogenetic sampling locality: each
validation cell's locality is the mean cophenetic distance to its 5 nearest
species with an observed value for that trait, cells are split into 2
strata at the median locality, and a split-conformal quantile is computed
per stratum. Cells in the far (undersampled) stratum get the wider
quantile at prediction time -- intervals widen in poorly-sampled clades,
which is the point. Traits with fewer than 19 validation residuals in
either stratum fall back to the global "split" score (19 is the
smallest stratum whose achievable coverage ceiling n/(n+1) reaches 0.95;
smaller strata are arithmetically guaranteed to undercover). Single-obs only:
fit_pigauto() errors immediately for multi-observation data, since
locality is defined at species level and multi-obs validation cells are
per-observation.
Four places where the pipeline silently did less than the user asked are now loud, plus one documentation debt:
lambda_mode = "estimate"/"cv"/"bayes" + covariates: the covariate-aware
GLS path (bm_impute_col_with_cov()) has no lambda argument, so those
modes silently ran at lambda = 1. fit_baseline() now warns once.fit_baseline() now warns once when a joint path fires with covariates
supplied. Threading covariates through the joint solver is future work.impute() now exposes conformal_split_val explicitly (previously
reachable only via ..., i.e. undocumented at the primary entry point).Documentation debt: lambda_mode values "cv" (per-column CV-selected
lambda) and "bayes" (Bayesian model-averaged prediction over a lambda
grid, with posterior-mean lambda as a diagnostic; the SE combines
within- and between-lambda variance) shipped in PR #109 but were never
recorded in NEWS. The 2026-08-16 lambda-attribution campaign
(docs/dev-log/2026-08-16-lambda-attribution-results.md) measured
"bayes" as the mode that eliminates the ML boundary collapse the
v0.10.0 "Known limitation" section describes, superseding that section's
"v0.11 CV-lambda" remedy plan.
compute_phylo_signal_per_trait() indexed species_names (one entry per
species) with an observation-length logical drawn from X_scaled (one row
per observation). In multi-obs mode this produced NA tip names,
ape::keep.tip() errored, the caller's tryCatch swallowed it to NA, and
phylo_signal_gate therefore never fired for any multi-obs trait —
silently, with no warning. Signal is now computed on species means, matching
how fit_baseline() aggregates multi-obs input. Verified against the
pre-fix code: a BM-simulated trait returned NA before and a finite
lambda > 0.5 after. Single-obs behaviour is unchanged, and with one
observation per species the two paths agree exactly.
Multi-obs users with phylo_signal_gate = TRUE may see gate decisions that
previously could not occur. That is the documented behaviour finally taking
effect, not a change of policy.
make_missing_splits() treated a (species, trait) cell as observed only when
all its latent columns were non-NA. zi_count encodes an observed zero as
gate = 0 with the log1p-z magnitude column NA — the NA is data ("this count
is zero"), not missingness — so every observed zero failed that test and could
never be drawn into validation or test. Because zero-inflation is the defining
feature of the type, this silently restricted all zi_count gate calibration,
conformal scoring and evaluation to the non-zero subset. Observation status for
zi_count is now decided by the gate column alone.
Measured on a 60-species fixture with 23 observed zeros: 0 of them could be
held out before the fix; they are eligible after. Any previously reported
zi_count validation or test metric was computed on non-zeros only and should
be re-measured.
Gate calibration accepted blends from half-set comparisons (or half-B
fallback corners) that could be materially worse than the pure BM
baseline on the full validation set — for continuous-family traits no
full-val BM comparison existed at all. calibrate_gates() now checks
both pure corners on the full validation set for every family; the
continuous-family BM comparison is margin-based (cal_min_rel_gain / 2,
default 1%) so near-ties survive (the 2026-04-29 AVONET Mass case)
while material losses close the gate. A second, identical floor runs in
fit_pigauto() on the post-refine calibrated surface — the one
predict() delivers — because a gate that wins on the pre-refine
surface can lose to pure BM after the calibrated refine pass (+2.0%
measured on the #157 fixture); a gate closed there costs nothing since
the delta contribution is zeroed. The strict val-floor test now
evaluates on the calibration-matched surface
(.mask_observed_idx = c(val_idx, test_idx)), matching evaluate() /
cross_validate() / impute() / report() — a plain predict() pins
each held-out cell's own truth as input context, a surface no genuinely
missing cell has in production.
This development note recorded an intended fit_baseline_bace() bridge,
including an experimental bace_final_imp() path. The installed bridge is
removed in 0.11; BACE is retained only as an in-tree comparator. Neither the
historical bridge nor its final-draw path has current pigauto inference
validation, so it must not be treated as a supported MI or pooling route.
The current supported downstream route begins with multi_impute_analysis()
in its documented narrow regime.
fit_ovr_categorical_fits() no longer describes its one-vs-rest
decomposition as "the OVR strategy BACE uses". BACE's
ovr_categorical default became FALSE (true multinomial) in
2026-08. pigauto's OVR is a different estimator — K threshold-joint
Rphylopars fits, used so that Rphylopars stays well-conditioned — and
that motivation is unchanged. Comment/roxygen only; no estimator
change.
preprocess_traits() / impute() reordered trait rows to tree tip
order but only checked that covariates had the same nrow. Shuffled
trait tables with matching-nrow environmental data silently mispaired
species and climate in the GNN covariate tensor. Covariates are now
aligned to the same species (single-obs rownames) or observation
(multi-obs species_col) order as X_scaled. Unmatched names error.
Unnamed covariates still follow the input-row permutation of traits.
GNN math is unchanged.
multi_impute(draws_method = "conformal") now scores zi_count traits on the log1p-z magnitude latent (held-out non-zeros) and draws count|nonzero from that score, or from the magnitude latent SE when scores are unavailable. The zero-inflation gate stays Bernoulli from P(nonzero). Previously scores were skipped and draws used pred$se of E[X], mixing gate uncertainty into the count SD. draws_method = "mc_dropout" is unchanged. Reported $se of expected count is not rewritten (P1).
Under max_iter = 0 (shipped default; cross-trait EM is not restored),
the Henderson init used for K≥2 joint MVN and mixed threshold-joint
baselines recycled one empirical SD for every missing tip. Missing-cell
SE is now BM-style σ√(1−h_i) from the same sparse Q already used for μ.
Binary threshold decode (pnorm(μ / sqrt(1+se²))) therefore sees
per-tip uncertainty. Conformal intervals are unchanged (they do not
consume baseline se).
Three related wiring bugs in the ResidualPhyloDAE fit path (review 2026-08-07). Defaults and the public API are unchanged; behaviour is more honest about what the model sees at train and cal time.
mask_heldout_with_baseline();
model_config$train_mask_heldout = TRUE).phylo_signal_gate triggers and safety_floor = FALSE (so no mean
baseline is available), the override falls back to pure BM
(r_cal_bm = 1) with a warning, instead of the mean corner that
would become latent 0 at predict time.refine_steps refine loop that predict() uses
(pigauto_refine_forward(); model_config$cal_refine_steps).The new multi_impute_analysis() backend requires the substantive model
before imputation and is deliberately narrow: one incomplete continuous
covariate under MAR, with proper Bayesian normal-regression MI for Gaussian
lm, smcfcs for binomial-logit glm, and jomo.smc for a Gaussian lmer
with one random intercept. It pools fixed effects only.
A warning-free 6,000-task package-level campaign at clean source SHA
2e3809d passed all 24 predeclared fixed-effect cells. Added bias was
0.00086-0.05042 oracle SD, pooled-SE/empirical-SD ratios were 0.942-1.030,
coverage was 93.9%-96.3%, and all results were finite. Engine warnings are
fatal, so no inference-ready object is returned after a rejection or
convergence warning. This validates the backend only within its documented
narrow scope.
The earlier conformal-width and Brownian/MC-dropout methods each passed 0/12 downstream fixed-effect cells in the 3,000-fit package campaign. PMM also failed its redesign pilot. Those methods remain experimental prediction diagnostics and are unsupported for downstream inference. Tree uncertainty was outside the analysis-aware campaign and remains unsupported by the new backend.
Adds first-pass support for Pagel's lambda < 1 in the BM baseline.
Spec: specs/2026-05-18-pagel-lambda-baseline-design.md (APPROVED).
New parameters:
fit_pigauto(..., lambda_mode = c("fixed_1", "estimate"))fit_baseline(..., lambda_mode = c("fixed_1", "estimate"))bm_impute_col(..., lambda) accepts numeric in [0, 1]
OR the string "estimate" (in-house ML via REML profile + 11-point
grid scan + local optim)Default lambda_mode = "fixed_1" is bit-identical to v0.9.x — no
silent behaviour change for existing callers.
What lambda < 1 does mathematically:
R(lambda) = lambda * R + (1 - lambda) * I
shrinks the off-diagonal of the phylogenetic correlation matrix toward zero. lambda = 1 is full BM kriging; lambda = 0 reduces the GLS predictor to the grand mean. Intermediate lambda allows partial phylogenetic shrinkage for traits with weak phylo signal.
No phylopars dependency on the lambda path. Both transform_tree_pagel
(internal scaler with Pagel terminal-edge compensation) and
ml_lambda_for_col (per-column ML) are in-house. When
lambda_mode = "estimate", the dispatcher forces per-column BM
(joint MVN / threshold-joint / OVR are gated off because phylopars
model = "BM" is silently lambda = 1).
On the v0.9.4 BIEN h2h CI bench (n=2000, 30% MCAR, m=20):
The sla case is the well-known weak-phylo-signal pathology of ML for Pagel's lambda: the profile likelihood becomes flat for small lambda when phylogenetic signal is weak, finite-sample noise dominates, and the point estimator can collapse to the boundary. The math is correct; the predictor is just suboptimal.
BACE's MCMCglmm achieves z-RMSE ≈ 0.69 on sla because Bayesian model averaging integrates over the variance-component posterior rather than relying on a point ML estimate.
Practical guidance: keep lambda_mode = "fixed_1" (the default)
unless you've verified on your data that ML lambda doesn't pick
degenerate boundary values. The spec's hard target ("sla z-RMSE
≤ 0.75") is not achieved by this release; it's a v0.11 target
via cross-validation lambda selection (spec to follow).
tests/testthat/test-pagel-lambda.R (31 expectations)bm-internal tests pass — back-compat verifiedA two-part Monte Carlo study (no code change) testing whether pigauto's
gated GNN holds up against the plain BM-kriging baseline under realistic
non-random missingness. Drivers:
script/sim_bench/overnight_2026_05_19_mechanisms.R (240 reps) and
script/sim_bench/pigauto_full_sweep.R (600 reps, Dan-comparable n=500,
50 sims/cell). Methodology writeup:
script/sim_bench/mechanism_sweep_methodology.md.
Design: 4 process scenarios (bm_strong, ou_strong, nonlinear_cov,
weak_signal) × 3 missingness conditions (phylo_MAR, trait_MAR,
trait_MNAR, DEP_STRENGTH = 1.5 per Penone et al. 2014) × 30% missing.
Caveat: in this single-trait DGP, trait_MAR depends on the env
covariate, which exists only for nonlinear_cov; for the other three
scenarios trait_MAR has no variable to depend on and degenerates to
MCAR. The conditions genuinely exercised are therefore phylo-MAR, MCAR,
and value-dependent MNAR for all four scenarios, plus a genuine
covariate-MAR for nonlinear_cov only.
Result on the 600-rep n=500 sweep — across all 12 cells: 6 material wins for pigauto, 6 ties, 0 material losses vs BM kriging (z-RMSE, material = |Δ| > 0.05). The win mechanism is scenario-specific and only partly attributable to the GNN:
weak_signal: the calibrated gate closes to ≈0; pigauto falls back to
column-mean while BM kriging mis-extrapolates 90%-noise as phylo signal
(Δ up to −0.22). This is the safety floor, not GNN learning.ou_strong: gate stays low (≈0.07–0.09) but pigauto beats both
baselines by avoiding BM's process-misspecification (Δ ≈ −0.13).nonlinear_cov: gate most open (≈0.18–0.25); the GNN + environmental
covariate genuinely beat column-mean and BM (Δ ≈ −0.09).bm_strong: gate harmless; pigauto ties BM (the correct answer).The missingness condition scales difficulty (trait_MNAR is uniformly
hardest) but never flips a verdict: pigauto's per-scenario standing
relative to BM kriging is stable across phylo-MAR, MCAR, and MNAR. The
2 apparent losses seen in the smaller 10-sim pilot were both
n=200/trait_MNAR and did not replicate at n=500.
External 60-replicate ablation (b1805, #106, code in #116)
using a corrected multi-trait DGP that fires the joint-MVN baseline.
N=2000, 5 seeds per cell, three arms — including the "lazy-optimiser
disarmed" pigauto_ON_L0 configuration (baseline_mu masked, lambda_shrink = 0).
Result (z-RMSE, averaged across all scenarios):
| arm | z_rmse | coverage | width |
|-----------------|--------|----------|-------|
| column_mean | 1.291 | — | — |
| bm_kriging | 1.163 | — | — |
| pigauto_OFF | 1.038 | 0.887 | 2.65 |
| pigauto_ON | 1.056 | 0.884 | 2.67 |
| pigauto_ON_L0 | 1.057 | 0.887 | 2.68 |
(Measured before the 2026-08 GNN train/cal-symmetry fix, which removed held-out truth from the training-time input context; not yet re-verified on the corrected pipeline.)
At BIEN scale (N=2000) the joint-MVN baseline already encodes the
cross-trait Σ structure; the within-row attention mechanism injects
noise and slightly widens conformal intervals. Even with the lazy-
optimiser trap disarmed (pigauto_ON_L0), the network regresses
against the baseline rather than finding new structure.
use_trait_attention stays default-FALSE. Kept as opt-in for
small-N regimes where the baseline hasn't converged (b1805's local
smoke test at N=300 showed it dropping validation loss before the
joint-MVN catches up at larger N).
Wires the existing WorldClim bioclim covariate pipeline (shipped in v0.9.1.9005) into the BACE-compatible head-to-head CI bench for the BIEN plant dataset. Closes the gap surfaced by the v0.9.3 autoresearch sweep — at that scale, GNN architecture knobs were exhausted; the remaining improvement had to come from environmental covariates, not capacity.
BIEN pooled per-trait log-RMSE on the h2h CI bench (n=2000 species, 30% MCAR, m=20 imputations):
baseline (full 19,109-pool, no cov, gate ON): 8.286
control (filtered 3,450-pool, no cov, gate OFF): 7.455
treatment (filtered pool, WC covariates, gate OFF): 7.236
Pool-filter effect = −0.83. Net covariate effect = −0.22. Combined gain vs v0.9.3 baseline: −1.05 (−12.7%).
Per-trait gates with covariates active:
height_m = 0.68 (raw |r| vs bioclim = 0.64)
leaf_area = 1.00 (raw |r| vs bioclim = 0.41)
wood_density = 0.40 (raw |r| vs bioclim = 0.16)
sla = 0 (raw |r| vs bioclim = 0.27; gate-closed)
seed_mass = 0 (raw |r| vs bioclim = 0.43; gate-closed)
Two of five traits (sla, seed_mass) remain gate-closed even with
covariates active. sla is the expected case (low raw correlation).
seed_mass is anomalous — its raw bioclim correlation is comparable
to height_m's, but the val-set gate calibration finds no benefit.
Diagnosis to follow.
script/gha/_bace_compat.R: BIEN config now sets
covariate_cols to the 38 bioclim columns (19 vars × median + IQR),
phylo_signal_gate = FALSE, and external_cov_rds pointing to the
new bioclim cache. .run_bench_bace_compat() learned an
external_cov_rds join that cbinds covariates onto traits_df
before subsetting and filters out species without complete bioclim
coverage. cov_cols AND external-cov column names are removed
from the default trait-subset so they don't get masked + imputed
as traits.useful/bace_data_snapshot/data/bien_worldclim_covariates.rds
(19,109 species × 38 cols; 3,450 with complete coverage).The W-series sweep tested six tweaks on top of the W0 (treatment) baseline; all regressed by 0.09–0.12, except W6 which was tied within noise:
W6's tie shows the phylo_signal_gate = FALSE recommendation
from specs/2026-04-24-worldclim-covariates-design.md §5.4 is
unnecessary when WC covariates are active — the calibrator handles
the gating either way. The config keeps FALSE for stability but
the constraint is softer than the spec implied.
Small opt-in feature release from an autoresearch sweep on the BIEN
plant bench. New: within-row cross-trait self-attention in the GNN
(use_trait_attention, default FALSE). Eighteen knob-tweak
experiments + three structural variants on BIEN n=2000 were run
between v0.9.2 and v0.9.3; none improved BIEN pooled RMSE
meaningfully over the default architecture. The cross-trait attention
path is shipped as opt-in because it is structurally distinct and may
help on datasets with strong within-row functional coupling that the
joint MVN / threshold-joint baseline cannot capture; it did NOT help
on BIEN (the baseline already encodes Σ for the dominant trait types).
fit_pigauto() and multi_impute() accept use_trait_attention = FALSE (default), n_trait_heads = 2L, trait_embed_dim = 32L.
When TRUE, the GNN's encoder receives an additional pooled
trait-context feature built by per-trait token embedding +
one multi-head self-attention block over the trait sequence.
Backward-compatible: saved fits without the field default to
FALSE and reconstruct identically.vignettes/gnn-architecture.Rmd). Documents the
end-to-end data flow, the baseline dispatcher (BM / LP / joint MVN /
threshold-joint / OVR), the transformer-block GNN with B2 rate-aware
attention and the new B3 within-row cross-trait attention, the gate-
and-calibrate safety mechanism, conformal prediction intervals, the
three uncertainty quantification mechanisms, and the two-step
Nakagawa & de Villemereuil (2019) tree-uncertainty workflow.Negative results, documented for the trail:
n_gnn_layers 2→4, hidden_dim 64→128,
ffn_mult 4→8, n_heads 4→16, and all combos regress (+0.02 to
+0.05 on pooled RMSE).n_heads 4→2, hidden_dim 64→32,
legacy attention path also regress.dropout 0.10→0.20,
lambda_gate 0.01→0.001, gate_cap 0.8→1.0.refine_steps 8→16,
corruption_rate 0.55→0.30.Marginal positive: n_heads 4→8 in the CI bench config (-0.002,
at-noise; bench-only).
Strategic implication: at BIEN's scale and signal regime, further RMSE lift needs environmental covariates (WorldClim) or larger posterior tree ensembles, not GNN-architecture changes.
Historical comparator snapshot — withdrawn as current evidence. This release recorded an in-house-versus-BACE comparison snapshot while the in-house Σ solver and reporting alignment were changing. It is not a current ranking or release claim; the external comparator gate remains open. The underlying imputation algorithm was unchanged from v0.9.1.
fit_mvn_bm_inhouse() had three independent bugs that surfaced
together on real phylogenetic data:
Henderson R⁻¹ scale mismatch. henderson_R_inv_apply() was
documented to return R⁻¹b but actually returned A⁻¹b (where
A = vcv(tree), R = cov2cor(vcv(tree))). For variable-depth
trees the mismatch is a per-tip rescale; for ultrametric trees
it is a global scalar. Fixed by adding cor_scale = TRUE
option that pre-/post-multiplies by sqrt(diag(A)).
henderson_R_inv_apply(..., cor_scale = TRUE) now matches
dense solve(R) to ~1e-6 relative on non-ultrametric trees;
the existing cor_scale = FALSE (A-scale) tests continue to
pass unchanged.
Σ init via species-iid sample covariance. The previous init
used cov(L, use = "pairwise.complete.obs") which treats
species as iid and double-counts sister tips as independent
observations. On strongly-phylo-conserved liability matrices
this produced near-singular Σ̂ initialisations whose off-diagonals
were biased toward zero. Replaced with the closed-form
matrix-normal MLE Σ̂ = (1/n) L̂ᵀR⁻¹L̂ on per-column-BM-completed
L̂ — properly credits R for the species correlation and recovers
Σ̂ to the right order of magnitude.
EM E-step diverged at near-sister tips. build_conditional_prior()
computes a within-row cross-trait conditional posterior but ignores
R-mediated cross-row correlation. For near-sister tips, the proper
joint posterior would strongly pull a missing tip's value toward its
observed sister; the within-row approximation cannot see that, and
combined with henderson_bm_predict()'s constant conservative
variance for missing cells, the inverse-variance pool over-weights
the wrong cross-trait prior. R⁻¹ at near-sister pairs (values of
4000–7500 observed at sister pairs on real AVONET) then amplifies
the discrepancy and Σ̂ inflated 10–50× per EM iteration. Fixed by
changing the default max_iter to 0L so the solver returns the
per-column BM L̂ + closed-form Σ̂ MLE without the broken cross-trait
EM. Callers can still opt in to the EM via max_iter > 0.
Historical comparator observation — withdrawn as current evidence. This section recorded a changing in-house/BACE snapshot while the solver was being repaired. It is not a current cross-package ranking; the external comparator gate remains open.
Major changes to script/gha/_bace_compat.R and
script/gha/make_headtohead_report.R so pigauto's CI numbers
are directly comparable to BACE's snapshot.
mi$pooled_point to match BACE's
reporting policy. BACE uses majority-vote-across-M-imputations
for categorical / ordinal accuracy and rowMeans(imp_mat) for
continuous RMSE. pigauto's mi$pooled_point carries the
analogous central tendency (argmax-of-average-probability for
categorical / binary, mode-or-mean for ordinal, median-or-mean
for continuous depending on log-transform). Before this fix,
.bace_eval_per_imputation scored every stochastic categorical
sample individually — bounded by E[p_max] across cells —
badly underestimating the model's accuracy. The previous
per-draw default also inflated continuous RMSE by Jensen's
inequality. A new pool_method = "auto" argument is the default;
legacy "per_draw" behaviour stays available for diagnostics.coverage_95, interval_width, and brier to the
canonical CI schema. Coverage and interval width are
pigauto-only (BACE has no posterior-interval snapshot column).
Brier is now computed on both sides for discrete traits so the
two methods can be compared on calibration.ordinal as a continuous-family type and tried to pull RMSE
from both sides, returning NA for every ordinal trait. Removed
ordinal from CONTINUOUS_TYPES so the report uses accuracy
for ordinal too. AVONET migration, PanTHERIA diet_breadth /
habitat_breadth now show up properly.load_pigauto_results() and
load_bace_snapshot() now route every loaded part through
.normalize_eval() so rbind across schema-version mismatches
no longer crashes the aggregator. The BACE snapshot, pinned to
a pre-Tier-1.5 schema, is automatically padded with NA in the
new pigauto-only columns.n_imputations in PIGAUTO_CI_CONFIG raised from 10
to 20. Mild improvement on BIEN total RMSE (8.312 → 8.288 in
an autoresearch sweep of 17 settings; m=30/50/100 all
regressed, so 20 is the empirical sweet spot at the current
CI scale).cov_encoder /
obs_refine paths) was tried and reverted: CI runs showed
AVONET categorical accuracy dropped sharply when covariates
were passed (trophic_level 0.825 → 0.731, primary_lifestyle
0.828 → 0.657) without lifting continuous traits. The
underlying covariate infrastructure stays in place — when
covariate_cols is empty the path is a no-op. Re-enabling
for AVONET (or BIEN) needs a categorical-safe covariate
routing first; planned for a future minor release.tests/testthat/test-bace-compat-eval.R covers the new
pool semantics + coverage / brier computation paths.cor_scale = TRUE
assertions on non-ultrametric trees and per-cell variance
checks for henderson_bm_predict().predict.pigauto_fit() now seeds DAE refinement with stored observed latent cells and uses the phylogenetic baseline only for originally missing cells. Previously prediction initialized every cell from the baseline, so the GNN did not receive the observed trait context it was trained to use. Evaluation, cross-validation, and report metrics explicitly mask validation/test cells from that context so held-out scores remain leakage-free.cov_linear fixed-effect path is now included in training, validation, gate calibration, conformal scoring, and prediction whenever user covariates are present. Previously that path was constructed and returned by the model, but production callers ignored it unless using the optional low-level baseline_mu forward path.r_cal_gnn alone.multi_impute_trees() completions. Tree completions are descriptive only;
the supported inference route begins with multi_impute_analysis() in its
documented narrow regime.multi_impute_trees(share_gnn = FALSE) now computes pooled_point by averaging all completed T * m_per_tree datasets, matching the documented return contract. Previously this path averaged one per-tree pooled result, which could ignore within-tree draw-level variation when m_per_tree > 1.multi_impute_trees(share_gnn = TRUE) now computes pooled_point from the completed imputation datasets, preserving species_col, user input row order, observed cells, and all m_per_tree draws. Previously this path pooled directly from raw prediction output, which could drop non-trait columns and return the point estimate in internal prediction order.multi_impute_trees(share_gnn = TRUE) now recomputes each posterior
tree's baseline from that tree's own cophenetic distances and
phylogenetic covariance. Previously the shared-GNN path accidentally
reused the reference tree's cached graph when fitting per-tree baselines,
which could collapse the intended tree-uncertainty signal.predict.pigauto_fit() now returns conformal_lower and conformal_upper
columns for explicit proportion traits when conformal scores are
available. Cross-validation/report coverage and prediction interval plots
now treat proportions as numeric back-transformable traits on the original
0-1 scale.multi_impute(draws_method = "conformal") now generates stochastic
draws for multi_proportion groups by sampling their CLR latent columns
with BM latent uncertainty and decoding back to the simplex. Previously
the default conformal MI path skipped these groups because their decoded
output columns are the component names rather than the group name.User-facing docs reorganisation following the issues by @b1805 (#67, #68)
that surfaced confusion the README and ?impute should have prevented.
PR #69 fixed the function-level docstring; this release reorganises the
pkgdown site itself:
useful/pkgdown_page_rethink.md.docs/dev/ on the current build are documented as
pending regeneration in the YAML and will re-enter the navbar once
their script/make_*_html.R driver has been re-run; the Phase-2
math vignette will land in a separate PR alongside this list refresh.vignettes/common-pitfalls.Rmd: five sections covering
prediction$imputed semantics, K=3 ordinal majority-class collapse,
the safety-floor closed-gate behaviour, phylogenetic signal diagnosis,
and clamp_outliers tail safety. Each follows a fixed
Symptom / Why / Diagnose / Fix / See-also template.DOCS.md retired; README + live pkgdown site are the single source of
truth for documentation navigation. README gains a one-line
"Live documentation" link near the top.No benchmark reruns or test changes.
Two silent-degradation paths surfaced by the pre-shipping coverage audit are now loud errors instead of "warn + degraded result":
rownames(traits) and tree$tip.label
used to drop every trait row, emit a "tree tips have no trait
data" message, and return a pigauto_data with all-NA X_scaled.
Downstream fit_pigauto() would then fail in a less-obvious place.
Now preprocess_traits() errors immediately with "No overlap between trait row names and tree tip labels: 0 species in common.".pigauto_data with no phylogenetic signal. Now errors with
"tree has fewer than 2 tips. Phylogenetic imputation requires at least 2 species.".No effect on any non-pathological input.
Two production bugs surfaced by the pre-shipping coverage audit and fixed in PR #65:
suggest_next_observation(by = "species") no longer dies with
"no rows to aggregate" on continuous-only traits. The species
aggregator now uses stats::aggregate()'s data.frame interface
(which tolerates all-NA columns) instead of the formula interface
(which routes through model.frame() and na.omits every row
when the lhs is all NA). See R/active_impute.R::692.pool_imputations() for count and zi_count traits now
honours pool_method = "mean" instead of silently falling back
to median. A leftover-from-refactor vals <- row_median_decoded(nm)
line had been overriding the pool_method-aware assignment. The
zi_count branch additionally had a duplicate (and therefore
unreachable) else if (tm$type == "zi_count") that has now been
collapsed into a single branch which (a) honours pool_method on
the magnitude column and (b) always probability-pools the gate
column for P(non-zero). See R/predict_pigauto.R::809+865.tests/testthat/test-shipping-coverage.R (+17 test_that blocks,
+68 expectations) covers:
confusion_matrix, calibration_df, simulate_non_bm,
read_tree, plot_history_gg, plot_uncertainty.multi_impute(draws_method = "mc_dropout"),
fit_pigauto(use_attention = FALSE),
fit_pigauto(conformal_method = "bootstrap").The bench evidence collected after PR #61 (Phase G' PMM) and the unsuccessful Phase G'' attempts (per-draw and post-pool conditional PMM) showed that PMM is not a tail-safety tool for pigauto. PMM in mice works because linear-regression predictions are unbiased estimators of truth across the prediction range. pigauto's GNN can be systematically biased on rare-clade species (the AVONET Casuarius case) AND simultaneously accurate on legitimate-but-out- of-observed-range species (the AVONET seed-2032 case). PMM cannot distinguish these two cases by donor-distance alone, so it can hurt RMSE just as easily as it can help.
This historical inference framing is withdrawn. PMM is an experimental
stochastic prediction mechanism, not a supported pooling or downstream
inference route. The current supported downstream route begins with
multi_impute_analysis() in its documented narrow regime.
This release updates the ?impute and ?predict.pigauto_fit
documentation + this NEWS entry to reflect the corrected framing.
No code changes.
For tail safety on single-imputation point estimates, the
recommended tool remains clamp_outliers = TRUE (Phase G PR #60).
Adds Predictive Mean Matching (PMM; Little 1988; Buuren & Groothuis-Oudshoorn 2011 mice) as an opt-in stochastic prediction mechanism. PMM replaces each missing cell's back-transformed prediction with an actual observed value drawn from a donor pool ranked by predicted-value proximity.
\strong{Originally framed as "tail safety", but the bench result showed mixed prediction behavior.} This historical implementation does not establish a downstream inference route; see the superseding note above.
Two new arguments on impute() and predict.pigauto_fit():
match_observed (character, default "none"). Set to
"pmm" to enable PMM.pmm_K (integer, default 5L -- mice convention). Donor
pool size.When match_observed = "pmm", for each missing cell of an at-risk
trait type the function:
pmm_K observed cells whose own predictions are
closest to the missing cell's prediction (the "donor pool").By construction, PMM-imputed values are always in the observed data range. No extrapolation is possible.
PMM is type-restricted to at-risk trait types where the back-transform amplifies tail overshoots:
expm1 decode)Discrete-class types (binary / categorical / ordinal / multi_proportion) are no-op: their decoders already return values from the observed level set by construction. Untransformed continuous: no-op (no expm1 risk).
match_observed = "none" remains the default for backward
compatibility (PR #60 added clamp_outliers opt-in; this PR adds
PMM opt-in alongside it). Cross-dataset evidence in v0.9.2 may
justify flipping the default to "pmm" for at-risk types.
Recommended use:
multi_impute_analysis() in its documented narrow regime.clamp_outliers = TRUE, clamp_factor = 5 for users who want
a deterministic point estimate without PMM stochasticity.R/pmm.R (NEW) -- pmm_impute_one_trait(), apply_pmm_to_decoded(),
recover_X_orig(), pmm_is_eligible(). RNG state is saved and
restored around per-call set.seed() so PMM does not leak the
user's RNG state.R/predict_pigauto.R -- predict.pigauto_fit() and
decode_from_latent() extended with match_observed and pmm_K
arguments. PMM is applied between decode and pool, per draw.R/impute.R -- user-facing wrapper extended with same args + roxygen.R/fit_pigauto.R -- the fit object now retains X_scaled so
predict() can recover original-units observed values for the
donor pool. Adds ~8 * n_obs * p_latent bytes (~600 KB at AVONET
full).tests/testthat/test-pmm.R -- 5 layers of tests, 65
assertions, all pass:
script/bench_phase_g_prime_pmm.R -- AVONET 4 log-cont traits +
synthetic heavy-tail bench across 3 configs (none / clamp / pmm).script/bench_phase_g_prime_pmm.R reports RMSE, max prediction, and
the fraction of imputations falling in the observed range, for
each (dataset, seed, trait, config) cell. See
useful/MEMO_2026-05-01_phase_g_prime_results.md for the verdict.
Closes the AVONET Mass tail-extrapolation mode documented in
useful/MEMO_2026-05-01_avonet_mass_diag.md: at N_IMP=20 on seed
2030, Casuarius bennetti was predicted at 538 kg vs truth 35 kg
because all 20 MC-dropout draws produced a +3-4 SD latent overshoot
that expm1() amplified to 250-1700 kg.
Two new arguments on impute() and predict.pigauto_fit():
clamp_outliers (logical, default FALSE). When TRUE,
post-back-transform predictions for log-transformed continuous,
count, and zi_count magnitude traits are capped at
tm$obs_max * clamp_factor (where tm$obs_max is the observed
maximum on the original scale, recorded at preprocess time).clamp_factor (numeric, default 5). Multiplicative factor
on the observed maximum. Anything > obs_max * clamp_factor is
considered implausible and capped.For non-at-risk types (untransformed continuous, ordinal, binary, categorical, proportion, multi_proportion) the clamp is a no-op.
clamp_outliers = FALSE remains the default for backward
compatibility. Users imputing log-transformed traits with
n_imputations > 1 and concerned about tail blow-ups should pass
clamp_outliers = TRUE explicitly. See
useful/MEMO_2026-05-01_phase_g_results.md for the AVONET
acceptance bench:
| seed | clamp OFF | clamp ON | reduction | |------|-----------|----------|-----------| | 2030 | 24,330 | 6,273 | -74 % | | 2031 | 312 | 318 | +1.7 % (MC noise) | | 2032 | 591 | 433 | -27 % |
Seed-2030 Casuarius prediction: 551,050 g (OFF) -> 167,847 g (ON, capped at 35,000 * 5 = 175,000).
R/preprocess_traits.R: records entry$obs_max and
entry$obs_min (original-scale) for log-transformed continuous,
count, and zi_count types. Untransformed continuous: NA (not at
risk of expm1 amplification).R/predict_pigauto.R: extends predict.pigauto_fit() and
decode_from_latent() with clamp_outliers / clamp_factor
arguments. New .clamp_high() helper applies the cap at the
three back-transform sites (continuous-log, count, zi_count
magnitude).R/impute.R: extends the user-facing wrapper with the same
arguments + roxygen.tests/testthat/test-clamp-outliers.R: 5 test_that blocks
(24 assertions), all pass. Covers obs_max / obs_min recording,
default-OFF no-op, ON behaviour on synthetic tail, clamp_factor
scaling, and input validation.
script/bench_phase_g_clamp.R runs the AVONET seed-2030/2031/2032
N_IMP=20 sweep with clamp ON and OFF. Output:
script/bench_phase_g_clamp.{rds,md}.
Closes the AVONET Migration N_IMP-dependent regression flagged in
useful/MEMO_2026-05-01_phase_f_smoke_results.md.
predict.pigauto_fit() and impute() now accept a third value for
pool_method:
"median" (default, unchanged) — per-cell median for count /
proportion / zi_count magnitude; per-cell mean(integer-class) -> round for ordinal."mean" (unchanged) — pre-v0.9.2 arithmetic-mean pooling."mode" (new) — per-cell majority vote across the M
draws for ordinal traits. For continuous-family traits ("mode"
falls back to "median" since mode does not apply to continuous
decoders). Binary, categorical, and multi_proportion are
unaffected (already probability-averaged + argmax).The previous ordinal pooling (mean of integer class indices, then round) biases predictions toward middle classes when MC-dropout
spreads the M draws across adjacent classes. Concrete K=3 failure:
True class 1 (sedentary), 5 dropout draws produce {1,1,1,2,3}. Mean = 1.6, round = 2 -> "partial migrant" (wrong). Mode = 1 -> "sedentary" (correct).
Empirical impact on AVONET Migration (3 seeds, n=1500, miss=0.30):
| pool_method | N_IMP=1 | N_IMP=20 | |-------------|-------------|-------------| | median (default) | 0.767 +/- 0.011 | 0.713 +/- 0.032 | | mode (new) | 0.767 +/- 0.011 | 0.779 +/- 0.011 |
Mean baseline (per-class mode): 0.800 +/- 0.013.
Switching to "mode" improves Migration accuracy at N_IMP=20 by
+6.6 pp and reduces the regression vs baseline from -10.9 pp
to -2.1 pp. Pre-registered acceptance criterion (Migration acc at
N_IMP=20 >= N_IMP=1) passes 3 / 3 seeds.
"median" remains the default for backward compatibility.
Users running ordinal traits with n_imputations > 1 should pass
pool_method = "mode" explicitly until cross-dataset evidence
justifies a default flip. See
useful/MEMO_2026-05-01_phase_h_results.md.
R/predict_pigauto.R: extended pool_method enum; new ordinal
branch using tabulate(class_indices, nbins = K) -> which.max().
Existing median-vs-mean guards on count / proportion / zi_count
branches extended to treat "mode" as "median" (fall back).R/impute.R: extended user-facing pool_method enum + roxygen.tests/testthat/test-mi-pool-robust.R: 5 new assertions covering
the mode-pool ordinal path and the median-fallback for continuous.script/bench_phase_h_pool.R runs the 3-seed AVONET-1500 Migration
sweep at N_IMP in {1, 20} with pool_method = "mode". Output:
script/bench_phase_h_pool.{rds,md} and the verdict memo at
useful/MEMO_2026-05-01_phase_h_results.md.
Closes the v1.5 scope gap: suggest_next_observation() now supports
all 8 pigauto trait types.
Hybrid metric: 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 \eqn{p_{\mathrm{gate}}}
(current LP estimate at \eqn{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 (= probability-
weighted magnitude variance reduction) 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.
Per-component variance reduction summed across K CLR-z-scored
latent columns. Per CLAUDE.md, multi_proportion rows are
fully observed OR fully missing (no partial K-component
observations), so delta_var_total is the K-component sum at each
fully-missing row.
R/active_impute.R:
zi_count and multi_proportion
in suggest_next_observation().bm_variance_reduction() and
lp_entropy_reduction_binary() helpers.needs_R_phy / needs_sim_phylo derived flags clarify
which kernels each trait type needs (zi_count needs both;
multi_proportion needs R_phy only).4 new test_that blocks (tests/testthat/test-active-impute.R):
[active v2] zi_count: hybrid variance + entropy populated[active v2] zi_count: types argument can exclude zi_count[active v2] multi_proportion: per-component BM variance summed[active v2] multi_proportion: delta equals K-component sum
(verifies via independent recomputation against
bm_variance_reduction() per component)Total active-impute tests: 17 / 60 expects, all pass.
For a fitted pigauto_result (from impute()), compute the
closed-form expected reduction in TOTAL predictive variance across
all currently-missing cells if each candidate cell were observed
next. Useful for sampling-design guidance when measurement budget
is limited: which species, if observed next, would most reduce
imputation uncertainty across the rest of your tree?
suggest_next_observation(result, top_n = 10L,
by = c("cell", "species"))
Returns a pigauto_active data.frame (descending by
delta_var_total). Columns:
by = "cell": species, trait, type, delta_var_total.by = "species": species, delta_var_total,
n_traits_missing.Built on a Sherman-Morrison rank-1 inverse update of the BM
conditional MVN: adding species s_new to the observed set updates
R_oo^{-1} by a known closed form, and the variance reduction at
every other missing cell is a single scalar per candidate.
Closed form per candidate cell:
delta_V_total(s_new) = sigma2 * sum_{i in miss} D[i, k]^2 / alpha_k
where D = R_mm - R_mo R_oo^{-1} R_om is the residual matrix and
alpha_k = D[k, k] is the current relative leverage of cell k.
The formula assumes sigma2 (REML variance) is held fixed -- the
standard active-learning convention. Variance reduction depends
only on the geometry of which cells are observed, not on the
unknown value at s_new.
continuous, count, ordinal,
proportion) -- variance reduction via the BM Sherman-Morrison
closed form above.zi_count and multi_proportion not yet supported (silently
skipped).Variance and entropy are not directly comparable. The
default output mixes them by delta value, which is approximate
across metrics. When precise ranking matters, filter by
metric ("variance" or "entropy") first.
This historical note described why the active-imputation helper was explored. The helper is now documented only as a model-based BM/label-propagation proxy ranking under its stated assumptions. It makes no uniqueness, novelty, optimal-design, or field-gain claim.
R/active_impute.R (new) -- 304 LOC, hosts
bm_variance_reduction() (internal) and
suggest_next_observation() (exported), plus
print.pigauto_active S3 method.R/impute.R -- adds tree to the pigauto_result so
suggest_next_observation() can recompute R_phy if it was
freed in the post-fit cleanup block.tests/testthat/test-active-impute.R (new) -- 10 tests, 29
expects. Verifies closed-form against an independent fixed-
sigma2 brute-force refit on a 20-tip fixture; descending order;
by-species aggregation; multi-obs error path; discrete-skip.On AVONET 300 with 90 cells masked across 4 continuous traits:
top-1 cell is Chloephaga_rubidiceps Beak.Length_Culmen,
delta_var_total = 2.99; top-1 species is Chloephaga_rubidiceps
with 2 currently-missing traits combining to delta = 3.23.
Closes the longest-standing open NEWS item: "A future improvement (cross-validated gate selection, ...) could close this further; not in scope for this release."
calibrate_gates() and fit_pigauto() gain gate_method = "cv_folds" (opt-in alongside the existing "single_split" default
and "median_splits" option). When selected, val cells are
partitioned into gate_cv_folds = 5L (default) deterministic
non-overlapping folds; the existing grid + half-B-verify procedure
runs once per fold (half_a = K-1 training folds, half_b =
held-out fold), and the componentwise median of K winning weight
vectors becomes w_final. Post-cal full-val invariant fires
unchanged.
"median_splits"| Property | median_splits | cv_folds | |---|---|---| | half_a per split | 1/2 of n_val | (K-1)/K of n_val | | Number of splits | B=31 random | K=5 deterministic | | Overlap | Random (each cell in many splits) | None | | Statistical interpretation | Median of B noisy gate selections | Standard k-fold CV |
Larger training sets per split → less noisy w_k. Standard CV gives a cleaner story for users and reviewers.
The strict val-floor enforces pigauto_val <= baseline_val, but
val→test drift survives in 4/32 binary cells (per
useful/MEMO_2026-04-29_discrete_bench_reruns.md). Per-trait gate
calibration overfits to the specific val cells because the same
cells score every candidate weight. Larger training sets per
split reduce overfitting; non-overlapping folds give a true held-
out evaluation.
gate_cv_folds < 2 errors with a clear message.n_val < gate_cv_folds caps K_eff = n_val so each fold has at
least 1 cell.n_val < 2 (CV ill-defined) falls back to a single split.gate_method = "single_split", and the "median_splits" path is bit-identical
to pre-fix.5 new tests in tests/testthat/test-safety-floor.R (under
"# CV-fold gate calibration"):
n_val < K graceful fallbackgate_cv_folds < 2 errorAll pass. Existing 24 safety-floor tests + B4 edge cases unchanged.
| gate_method | wall (s) | r_cal_gnn[col_3] | |---|---|---| | single_split | 82.6 | 0.000 | | median_splits | 77.2 | 0.050 | | cv_folds | 75.9 | 0.150 |
cv_folds is slightly faster (smaller per-split work) and selects
a non-zero GNN gate on the col where the others snap to BM.
script/bench_cv_vs_median.R runs a focused 3-way comparison on
a synthetic continuous + binary + categorical fixture (BM, lambda=1,
n=300, 30 % mask, 5 reps). Test-set means across reps:
| Method | bin (acc) | cat3 (acc) | cont (RMSE) | |---|---|---|---| | single_split | 0.8222 | 0.7356 | 1.0778 | | median_splits | 0.8222 | 0.7356 | 1.0720 | | cv_folds | 0.8222 | 0.7356 | 1.0611 |
useful/MEMO_2026-04-29_discrete_bench_reruns.md) showed:
the val→test drift on 4/32 binary cells is val→test
extrapolation noise, not gate-selection noise; cv_folds does
not close it on this regime. The hypothesis that larger
training sets per split would close MORE of the val→test
drift than median_splits is not supported by this
fixture. cv_folds is preferable for continuous traits but
not transformative for discrete.Full per-type bench reruns (bench_binary.R, bench_categorical.R,
bench_zi_count.R at multiple signal levels) deferred to a
later session.
A type-restriction bug in R/fit_helpers.R::calibrate_gates was
silently exempting binary, categorical, and zi_count traits from the
post-calibration "strict val-floor" invariant. The buggy behaviour
let pigauto underperform its own baseline on the package's own
simulator benches (bench_binary.md, bench_categorical.md,
bench_zi_count.md):
signal_0.6 mean accuracy: −4.8 pp below baselineK_3 cat2: −12 pp below baselinezf_0.2 zi1: +23 % RMSE above baselineThese violate the safety-floor's "pigauto cannot underperform its
own baseline" promise. The bug was introduced 2026-04-23 with the
initial safety-floor commit (28d8e45); the type filter mirrored an
earlier continuous-only MSE evaluator and was not extended when the
simplex grid landed.
The fix splits the post-calibration invariant into two type-specific blocks:
Empirical evidence (CORRECTED 2026-04-30 per the Opus adversarial review; the original tabulation here mistakenly averaged val + test splits in the bench RDS, overstating the closure). The user-facing held-out test-set mean acc gap pigauto − baseline (lower-magnitude = closer to baseline-parity floor):
| bench | Apr 17 (pre-bug) | Apr 28 (buggy) | Apr 29 (post-fix v3) | |---|---|---|---| | binary | −0.029 | −0.027 | −0.014 | | categorical | −0.031 | −0.026 | −0.016 | | zi_count (RMSE ratio) | 1.058 | 1.043 | 1.012 |
The strict val-floor closes ~50 % of the test-side regression on binary and categorical (from −2.7 / −2.6 pp to −1.4 / −1.6 pp); the remaining ~1.4 pp is val→test sampling drift. zi_count test ratio drops from 1.06 to 1.01 -- most of the regression closes, with a small residual. 3 of 32 binary cells still drift ≥ 5 pp on test (val→test extrapolation that no val-only floor can prevent).
The strict val-floor's invariant is enforced at calibration time
on the calibration-time loss surface
(compute_corner_loss(g, val_row_idx, ...)). At predict time the
model runs additional refine_steps and mask_token substitution that
introduce a small (~5 %) drift; the
tests/testthat/test-safety-floor.R strict-floor invariant test
documents this with a * 1.05 slack. So the user-facing guarantee
is "pigauto val-loss ≤ baseline val-loss + 1e-12 at calibration time,
modulo ~5 % refine_step drift at predict time" -- not the bit-exact
"≤ 1e-12 everywhere" reading.
Continuous benchmarks (AVONET n=1500 single seed × N_IMP=5): the originally-reported "v3 within MC-dropout noise of pre-fix" landed on a particularly close pair of draws. Re-running the same seed × n × N_IMP three times produced Mass test RMSE = 377, 430, 504 -- the comparison is non-deterministic at this evidence level (pooling order over the 5 MC dropout draws varies between R sessions). A multi-seed × N_IMP=20 verification is queued; the single-seed single-run evidence does NOT support strong claims about whether v3 changed continuous performance, only that the change (if any) is within run-to-run pooling noise on this dataset.
Per-cell variance survives in multiple binary cells (val→test
sampling drift). The strict val-floor only guarantees
pigauto_val_loss ≤ baseline_val_loss + 1e-12 (calibration time);
it does not guarantee pigauto_test_loss ≤ baseline_test_loss. A
future improvement (cross-validated gate selection, epsilon-margin
tightening, or LP added as a corner option for ordinal) could
close this further; not in scope for this release.
R CMD check after the fix: 0 errors / 0 warnings / 0 notes.
Full devtools::test(): FAIL 0 / SKIP 2 / PASS 1081.
New regression test in tests/testthat/test-safety-floor.R:
"strict val-floor: pigauto val-loss <= baseline val-loss per trait, all types" — fits on a mixed continuous + binary +
categorical + count fixture, recomputes val-loss for the
calibrated blend and the pure baseline, asserts blend ≤ baseline ×
1.05 for every trait.
Calibration evidence preserved at
script/_strict_floor_v3_calibration/ (4-way comparison: pre-fix /
v1 over-strict / v2 intermediate / v3 current). Memos:
useful/MEMO_2026-04-29_*.md.
fit_baseline() now computes BOTH the threshold-joint baseline (B3,
commit a541dbd) AND a per-column BM-via-MVN alternative for each
ordinal trait, then picks the lower-val-MSE path per trait against
splits$val_idx. The chosen path is exposed as
fit_baseline()$ordinal_path_chosen (named character: latent col →
"threshold_joint" or "bm_mvn").
Motivation: the AVONET-300 Phase 6 bench's Migration RMSE regressed
from 0.879 (Apr 16) to 0.975 (Apr 29) when B3 swapped the ordinal
baseline path from a label-propagation-equivalent BM-via-MVN to the
truncated-Gaussian threshold-joint approach. At K=3 ordinals the
K−1 thresholds are pinned to a narrow band by phylopars EM,
producing systematically worse predictions than BM-via-MVN on
z-scored integer class. See
useful/MEMO_2026-04-29_phase6_migration_bisect.md.
Per Opus's adversarial review #6, we did not ship a K ≤ 3 → LP
heuristic. Instead the per-trait selection runs at every K and
the gate's safety-floor simplex no longer has to route around
threshold-joint output for low-K ordinals -- the lower-val-MSE
baseline is selected before the gate calibration sees it.
Scope: single-obs only (multi-obs ordinal selection is out of scope
for this fix; it would require species-level aggregation of the BM
alternative). The code path is gated by !multi_obs && !is.null(splits)
and only activates when threshold-joint actually populated the
ordinal column.
Two regression tests in
tests/testthat/test-joint-threshold-baseline.R:
"fit_baseline picks per-trait ordinal path and reports it" --
smoke test on K=3 ordinal + continuous mixed data; verifies finite
output and that $ordinal_path_chosen is one of the two valid
values."fit_baseline ordinal selection picks BM when threshold-joint loses on val" -- recomputes both paths' val MSE independently
and asserts the selection picked the lower one.Four correctness fixes surfaced by the Opus adversarial pass on top of PR #49 (the safety-floor / mean-gate work). All four are silent-failure modes that produced plausible-looking but invalid output; none of them threw a useful error before this release.
HIGH — multi_impute(draws_method = "conformal") no longer
collapses to zero between-draw variance on shuffled input.
R/multi_impute.R::.sample_conformal_draw now threads
input_row_order so it correctly maps user-input mask indices
to internal (tree-tip) row indices when perturbing
pred$imputed. Previously, on any data frame whose rows are
not already in tree-tip order, the sampler perturbed the
wrong cells and build_completed re-aligned the unperturbed
values back into the user's missing cells, producing M
identical datasets and var = 0 after pooling. Any prior
Rubin's-rules pooling on shuffled-input multi-imputation was
unreliable. Caught by a new shuffled-input variance regression
test; pre-fix per_cell_sd ≡ 0, post-fix > 1e-8 at every
masked cell.
MEDIUM (opt-in) — conformal_split_val argument on
fit_pigauto() lets users avoid the validation-set
double-dipping that breaks split-conformal exchangeability.
When TRUE, the validation cells are split per-column into a
calibration half (used by calibrate_gates() to pick the
blend weights) and a conformal half (used to estimate the
residual quantile). Reusing the same val cells for both
selects the gate to minimise residual MSE on the very cells
whose residuals drive the conformal quantile, producing
systematic undercoverage that's most visible at small n_val
(matches the coverage_investigation memo).
Default is FALSE (pre-fix single-set behaviour) because
forcing the split regresses the AVONET300 / OVR-categorical /
BIEN safety-floor smoke benches by 2-26% on small-val
datasets. Use conformal_split_val = TRUE when accurate 95%
coverage matters more than bench-grade RMSE.
MEDIUM — calibrate_gates half-A loop is now NA-safe.
Mirrors the half-B fix from commit 573decd: when
loss_a_pure_bm is non-finite (e.g. multi_proportion CLR
latents whose half-A subset is too small for some component),
best_la is now seeded with Inf and the legacy
pure-BM-relative-gain filter is skipped rather than throwing
missing value where TRUE/FALSE needed.
LOW — predict.pigauto_fit no longer produces a cryptic
torch shape error when calibrated_gates is overridden with a
vector of the wrong length. Length 0 is now treated as "no
calibration" (degrades gracefully to the learned-gate path);
any other non-matching length errors with an explicit
"calibrated_gates has length k but the model has p latent column(s)" message. Unblocks the AE-attribution ablation
workflow used by script/phase1_gnn_ablation.R.
Each fix has a dedicated regression test:
tests/testthat/test-multi-impute.R::"multi_impute(draws_method=\"conformal\") produces non-zero between-draw variance on shuffled input",
tests/testthat/test-fit-predict.R::"conformal_split_val toggle changes conformal scores (no double-dipping)",
and
tests/testthat/test-fit-predict.R::"predict.pigauto_fit catches wrong-length calibrated_gates override".
Fixes the core limitation of PR #47 (centroid-only bioclim extraction): bioclim is now extracted at every GBIF occurrence point per species (not just the centroid), so the IQR columns carry real range-breadth signal instead of being constant zero. This is what actually lets the safety-floor calibrator see a meaningful covariate input and open the GNN gate on plants.
Validated by the honest-sim diagnostic
(experiment/covariate-honest-sim) which showed pigauto's architecture
can produce 32 % RMSE lift on strong-environment traits when
covariates carry real signal — something centroid-only bioclim
failed to deliver on plants.
pull_gbif_centroids(..., store_points = FALSE) — new arg,
default FALSE preserves pre-v0.9.1.9006 cache format. When
TRUE, persists the raw filtered lat/lon points in each species'
cache RDS, enabling per-occurrence bioclim extraction.pull_worldclim_per_species() — no signature change. Internally,
.wc_extract_one() now prefers cached per-occurrence points,
falling back to centroid only when points field is absent
(legacy caches).Pre-v1.1 cache RDS files remain fully readable. Users who want per-occurrence extraction upgrade their cache with one call:
pull_gbif_centroids(sp_list,
cache_dir = "script/data-cache/gbif",
store_points = TRUE,
refresh_cache = TRUE) # one-time refresh
After that, subsequent calls are instantaneous (cache hit).
tests/testthat/fixtures/worldclim_plants_300.rds re-generated with
per-occurrence aggregation. IQR columns are now informative.
script/bench_bien_worldclim.R. Expected: SLA r >= 0.35,
leaf_area r >= 0.30.New exported helper pull_worldclim_per_species(species, gbif_cache_dir, worldclim_cache_dir, ...) extracts 19 WorldClim v2.1
bioclim variables at each species' GBIF occurrence points, aggregates
per species (median + IQR), and returns a 38-column numeric
data.frame ready for impute(..., covariates = ...).
pull_gbif_centroids() (B.1).worldclim_cache_dir on first call; sentinel file makes subsequent
calls offline.terra::extract().worldclim_cache_dir/extracts/.The phylo_signal_gate introduced in v0.9.1.9003 fires on raw trait
values, so traits with weak phylogenetic signal (e.g. plants SLA,
leaf_area, height_m) get routed directly to the grand-mean corner
BEFORE the GNN sees any covariates. To see bioclim lift those
traits, users must set phylo_signal_gate = FALSE:
cov_bio <- pull_worldclim_per_species(
rownames(traits),
gbif_cache_dir = "script/data-cache/gbif",
worldclim_cache_dir = "script/data-cache/worldclim")
res <- impute(traits, tree, covariates = cov_bio,
phylo_signal_gate = FALSE) # <-- required
A follow-up spec ("covariate-aware phylo-signal gating") will compute lambda on residuals-after-covariates so the gate composes cleanly with covariates without the manual opt-out.
terra added to Suggests. pull_worldclim_per_species() errors
with an install hint when terra is absent.
tests/testthat/test-worldclim-covariates.R: 20 offline
expectations (cache key, aggregation, download sentinel, extract
cache, NA handling).tests/testthat/fixtures/worldclim_plants_300.rds: pre-extracted
fixture for 300 BIEN species.phylo_signal_gate = FALSE.New exported helper pull_gbif_centroids(species, cache_dir, ...)
fetches species occurrence records from the Global Biodiversity
Information Facility (GBIF) and returns a data.frame of per-species
median lat/lon centroids ready for impute(..., covariates = ...).
For each species, pull_gbif_centroids():
Plant traits like SLA and leaf area are environment-driven. Even the simplest geographic covariate -- a species range centroid from GBIF -- carries biogeographic signal that phylogeny alone doesn't. This is the B.1 proof-of-pipe for the pigauto covariate path; the full B.2 spec (WorldClim bioclim + SoilGrids raster extraction) follows.
Adding GBIF centroids as covariates keeps SLA and leaf_area RMSE within +10% of the no-cov baseline (ratios ~1.09 and ~1.01 at n=200). At this sample size the val-to-test sampling noise dominates any covariate lift signal; production-grade positive lift reserved for the B.2 bioclim run at larger n.
rgbif added to Suggests (not Imports). pull_gbif_centroids() fails with a helpful install hint when rgbif is absent.
B.2: full WorldClim bioclim + SoilGrids extraction at each GBIF occurrence, aggregated per species. Expected to lift plants RMSE significantly more than centroids alone.
Default change since v0.9.1.9003: impute() / fit_pigauto() now
compute per-trait Pagel's lambda on training-observed cells and route
traits with lambda < 0.2 to the grand-mean corner (0, 0, 1) of the
safety-floor simplex, skipping BM fitting and GNN training effort on
those traits. Complements the v0.9.1.9002 safety floor with an
explicit per-trait decision: "pigauto knows when phylogeny helps".
phylo_signal_gate = FALSE reproduces v0.9.1.9002 safety-floor-only behaviour bit-identically.
print.pigauto_fit() now surfaces a "Phylogenetic signal" section listing gated vs passed traits when any trait was gated.
phytools added to Suggests. When absent, phylo_signal_gate degrades gracefully to FALSE (safety-floor-only) with a one-time warning.
Default change: impute() / fit_pigauto() / multi_impute* now run
with safety_floor = TRUE. The calibration grid post-training searches
a three-way convex combination of (BM baseline, GNN delta, grand mean)
on a 231-point simplex at step 0.05. Because the pure-mean corner
(0, 0, 1) is always a candidate in the grid, validation RMSE is
guaranteed by construction to satisfy pigauto_val_RMSE <= mean_val_RMSE
on every trait. On held-out test, the invariant holds up to
sampling-noise slack (tolerance +2% in the production canary,
verified on five taxa by script/regress.R).
height_m, leaf_area, seed_mass,
sla) the calibrator opens r_MEAN > 0 and the prediction falls
toward the grand mean, fixing the boundary-case regression where
pigauto was previously 15-101% worse than mean; wood_density lift
preserved.impute(..., safety_floor = TRUE) — new arg, defaults to TRUE.fit_pigauto(..., safety_floor = TRUE) — same.r_cal_bm, r_cal_gnn,
r_cal_mean (each a named numeric of length p_latent), and
mean_baseline_per_col.multi_impute() / multi_impute_trees() propagate the mean term
automatically via predict.pigauto_fit(); no signature change.%||% fallback in
predict.pigauto_fit(): when the new slots are absent, the 2-way
blend (1 - r_cal) * BM + r_cal * GNN is reconstructed from the
scalar r_cal.safety_floor = FALSE reproduces the v0.9.1 1-D calibration
bit-identically on the legacy path.
tests/testthat/test-safety-floor.R: smoke canary, ~3 min on MPS.
23 test_that blocks, 75 expectations covering: simplex grid,
mean-baseline scalar, calibrate_gates 3-way path, fit_pigauto
integration, predict 3-way blend + legacy fallback, multi_impute
propagation, multi_impute_trees share_gnn propagation, AVONET300
regression (continuous RMSE within +2%, discrete accuracy within
-2pp), plants safety smoke (cached BIEN n=1000, RMSE within +15%
of grand mean).script/regress.R: full canary, ~20 min on MPS at n=1000 per taxon.
Writes script/regress_result.json with per-bench pass/fail.R/calibration_three_way.R: new — simplex_grid(),
mean_baseline_scalar() helpers.R/fit_helpers.R: calibrate_gates() now returns list of 3 weight
vectors and optionally searches the simplex grid.R/fit_pigauto.R: computes mean_baseline_per_col from training
cells; passes safety_floor through; stores the new slots.R/predict_pigauto.R: 3-way blend with %||% backward compat.R/impute.R: pass-through arg + roxygen.R/multi_impute.R, R/multi_impute_trees.R: roxygen only.inst/extdata/legacy_fit_v091.rds: fixture for backward-compat
test.data-raw/make_legacy_fit_v091.R: one-shot fixture builder.specs/2026-04-23-safety-floor-mean-gate-design.md: design spec.plans/2026-04-23-safety-floor-mean-gate.md: implementation plan.pool_method = c("median", "mean") on impute() and
predict.pigauto_fit(). Default "median" takes the per-cell
median of the M decoded draws when n_imputations > 1 for count,
proportion, and zi_count magnitude traits — robust to
dropout-noisy latents that get amplified by expm1() / plogis()
decoders. "mean" restores pre-v0.9.2 arithmetic-mean pooling.litter_size at
n_imputations = 20 (default RMSE 7.7, em5 RMSE 60.0) — a single
dropout-noisy latent draw at ~+5 SD decoded to thousands and
contaminated the mean. With median pooling the RMSE returns to
n_imputations = 1 scale. Coverage was unaffected (conformal ≥ 0.95,
MC-dropout ≥ 0.92); this is purely a point-estimate fix.pool_method has no effect on them.tests/testthat/test-count-mi-pooling.R guards the
robustness: 5-draw synthetic with one injected outlier of 1000 on a
count trait returns median = 2 (expected) but mean ≈ 202 (bug).Second real-data benchmark for pigauto after AVONET (birds). Tests the "general-purpose phylogenetic trait imputation" claim on a different taxon with a different trait mix.
script/fetch_pantheria_and_tree.R — one-shot downloader. Fetches
PanTHERIA (Jones et al. 2009, Ecology) from the ESA archive and
builds a taxonomy-derived cladogram from the Order/Family/Genus/
Species columns via ape::as.phylo.formula() + Grafen branch lengths.
Not a time-calibrated phylogeny, but encodes the hierarchical
taxonomic structure pigauto's phylo prior uses. Self-contained; no
external tree-database dependency.script/bench_pantheria_full.R — full-scale bench on ~4,629 mammals
× 8 canonical traits (body mass, head-body length, gestation length,
litter size, max longevity, diet breadth, habitat breadth,
terrestriality), MCAR 30% per trait. pigauto_default +
pigauto_em5 vs mean_baseline.script/bench_pantheria_bace_head_to_head.R — stratified 500-
species subset with optional BACE::bace() for cross-package parity
(graceful requireNamespace skip when BACE isn't available).script/make_pantheria_summary_html.R — aggregate landing page
linked from the validation suite alongside the Phase 8 summary.Five reproducible bench scripts + one aggregate report, all linked from
the pkgdown validation suite. No R/ code changes; bench infrastructure
only.
script/bench_signal_sweep.R — Pagel's λ ∈ {0.1, 0.3, 0.5, 0.7,
0.9, 1.0} × mixed-type traits (2 cont + 1 bin + 1 cat K=3) × 4 methods
× 3 reps. Headline: at λ = 1.0 pigauto reaches 94.4% binary accuracy
vs mean_baseline 76.3%; at λ = 0.1 all methods collapse near chance
(expected — no signal).script/bench_bace_avonet_head_to_head.R — pigauto_default vs
pigauto_em5 vs optional BACE::bace() on avonet300/tree300 with
identical splits. Graceful requireNamespace guard when BACE is
unavailable. Honest finding: Phase 6 EM does NOT universally beat the
plug-in at n = 300 (default 80.5% vs em5 65.9% on Trophic.Level),
consistent with the "Calibration at small n" caveat from v0.9.1 — EM
is opt-in for a reason.script/bench_correlation_sweep.R (Phase 8.1) — cross-trait
correlation ρ ∈ {0, 0.2, 0.4, 0.6, 0.8} × 4 continuous traits × 3 reps.
pigauto delivers consistent 6–10× RMSE reduction over mean_baseline
across all ρ (0.07–0.14 vs 0.64–0.80). Pearson r averages 0.97 at all ρ.script/bench_evo_model_sweep.R (Phase 8.2) — 4 evolutionary
models × 3 reps. BM RMSE 0.13 (6× lift over mean); OU 0.32 (3×);
regime_shift 0.07 (14×); nonlinear 0.63 (1.6× — pigauto struggles
on genuinely non-BM correlations, Pearson r drops to 0.69). Honest
graceful-degradation story.script/bench_clade_missingness.R (Phase 8.3) — realistic MAR
pattern via random-clade masking on focal trait. At every target_frac
tested pigauto maintains Pearson r ≥ 0.96 and RMSE 5–8× lower
than mean_baseline: 0.10→0.18 vs 0.99 (r 0.97); 0.25→0.14 vs 1.11
(r 0.97); 0.40→0.16 vs 1.16 (r 0.96). Strong validation of pigauto's
phylo prior in the realistic missingness regime (real databases
undersample clades, not random species). Clade-picker caps at
target_frac to avoid degenerate 100%-masked runs.script/phase8_summary.html — aggregate TL;DR + sub-report links
sessionInfo()).pkgdown/assets/validation_suite.html) gains five
new rows, one per bench, with graceful pending fallback on missing
RDS files.Real-data benches now cover four vertebrate classes plus a plants run that establishes the scope of phylogenetic imputation. All use 30% MCAR held-out, seed 2026, and report both conformal and MC-dropout 95% coverage alongside RMSE / accuracy. The vertebrate benches consistently lift point estimates by 27-75% RMSE / 10-40 pp accuracy; the plants bench is the first honest boundary-case result, with a clean wood_density lift but r ≤ 0.21 on four other traits, attributable to weak phylogenetic signal in pooled BIEN species means and random polytomy resolution in V.PhyloMaker2 (details below). Conformal coverage still lands at 0.90-0.98 on plants, so the uncertainty story generalises across the kingdom jump even where point estimates do not.
script/bench_avonet9993_bace_n3000.rds, Vulcan L40S):
AVONET n=3,000 subset x 7 traits. Mass RMSE -45%, Beak -55%,
Tarsus -46%, Wing -57%; Trophic.Level accuracy +16 pp;
Primary.Lifestyle +13 pp. Conformal 0.94-0.96.
New driver script/bench_avonet_full_local.R runs the full
n=9,993 locally on Mac MPS (Vulcan hit a predict-stage OOM at
this scale); overnight job, HTML generator
script/make_bench_avonet_full_local_html.R ready.feature/bench-pantheria): PanTHERIA n=4,027 x
8 traits. terrestriality accuracy 0.55 -> 0.95 (+40 pp);
body_mass_g RMSE -27%, head_body_length -75%; diet/habitat
breadth +10 pp each. Conformal 0.93-0.96.script/bench_fishbase.rds, Mac MPS):
FishBase x fishtree n=10,654 x 6 traits.
Vulnerability RMSE -45%, Length -38% (r=0.81), BodyShapeI
accuracy +31 pp; Troph -23%, Depth -10%; Weight improved from
r=0.04 (noise) to r=0.42 under the new median-pool MI fix.script/bench_amphibio.rds, Mac MPS): AmphiBIO
n=5,237 x 2 continuous traits. Body_size_mm RMSE 89 -> 44
(-51%), Body_mass_g -27% (r=0.98). Continuous-only v1 --
AmphiBIO binary / categorical columns hit a character/double
type mismatch in the threshold-joint baseline, parked for
v0.9.2. Calibrated gates = 0 on both traits: pigauto's GNN
correctly recognises that the Level-C joint MVN baseline is
already optimal, and still beats the mean baseline by 27-51%
because the BM baseline captures cross-trait correlation +
phylogenetic signal.script/bench_bien.R, honest boundary-case finding):
BIEN x V.PhyloMaker2 at n=4,745 species, n_imputations=20.
Only wood_density lifts (-6% RMSE, r=0.43); height_m, sla,
seed_mass, leaf_area show r ≤ 0.21 and are 15-101% worse
than grand mean on RMSE. Root cause is in the data, not pigauto:
(i) BIEN species means are pooled from heterogeneous individual
observations, diluting phylogenetic signal, and (ii) V.PhyloMaker2
scenario S3 resolves within-genus polytomies randomly, so most
species sit behind random branches relative to the Smith & Brown
2018 backbone. pigauto's gate safety contains the damage: on
weak-signal traits the calibrated gate closes to 0 and the
pipeline degrades gracefully to the (also weak) BM prior rather
than blowing up. Conformal coverage still lands at 0.90-0.98
on all 5 traits -- the distribution-free uncertainty story
holds even when point estimates are weak, separating the UQ
narrative (always robust) from the point-estimate narrative
(taxon-dependent). A first pass at n_imputations = 1 showed
height_m RMSE blowing up to 47.7 via Jensen exp()-decode bias;
rerunning at n_imputations = 20 activated the dc8cffa
median-pool MI correction and dropped it to 15.15. The weak-
signal result is what remains after that fix -- it is a property
of BIEN x V.PhyloMaker2 at this species pool, not a pigauto
regression. See script/bench_bien.html for the full honest
reading and the in-page Jensen note.Validation suite (pkgdown/assets/validation_suite.html) and the
paper-draft summary (useful/results_summary.md) collect all
numbers in one place.
script/bench_fishbase.R benchmarks pigauto on the intersection
of fishtree::fishtree_phylogeny() (Rabosky et al. 2018) and
rfishbase's load_taxa() + species() + ecology() -- 10,654 species
with 6 mixed-type traits (Length, Weight, BodyShapeI, DepthRangeDeep,
Vulnerability, Troph). Completes the vertebrate breadth triad:
birds (AVONET), mammals (PanTHERIA), fish (FishBase + fishtree).script/make_bench_fishbase_html.R + pkgdown
mirror at pkgdown/assets/dev/bench_fishbase.html.pool_imputations() now uses median pooling across decoded
values for all trait types where the decode step amplifies
dropout-noisy draws: log-transformed continuous (new), count
(extended), zi_count (extended), and proportion (new). Untransformed
continuous keeps mean pooling (backward-compat default).n_imputations > 1, a few dropout-noisy latent
draws with high-tail values decode via expm1() / plogis() to
absurd magnitudes that dominate rowMeans(). Median is the
outlier-robust pooling choice that matches PR #41's approach for
count / proportion, now extended across the full set of
decode-amplifying trait types.test-mi-pool-robust.R cover the four
trait types that switched to median, verifying that a single
extreme draw does not dominate the pooled output.fit_pigauto() now moves its returned model_state to CPU and
calls torch::cuda_empty_cache() before returning, so the training
run's activations, optimizer state, and per-epoch intermediate
tensors can be reclaimed by torch's CUDA caching allocator. Prior
to this fix, the fit object held GPU tensor refs that kept
~40 GB of training-time memory resident, and the downstream
predict() path OOM'd on any card with <= 80 GB at n >= 5000.impute() now calls gc() + torch::cuda_empty_cache() between
fit_pigauto() and predict() as belt-and-braces.pigauto_fit shape,
predict() output identical; just no longer OOM at scale.em_offdiag = FALSE on impute() and fit_baseline().
When TRUE AND em_iterations >= 2L, each binary/ordinal liability
cell's prior at iteration k + 1 is the full conditional-MVN
(mu, sd) given the posterior liability of other traits at iteration
k. This means observing one discrete trait shifts (not just
tightens) the prior on correlated other traits — the off-diagonal
entries of Σ now enter the E-step, instead of being ignored as in
Phase 6.em_offdiag = FALSE default preserves Phase 6 behaviour byte-for-byte.em_offdiag = TRUE with em_iterations < 2L is silently clamped
(iter 0 has no EM; iter 1 is plug-in with no previous Σ to condition
on).em_state gains em_offdiag field.em_iterations = 0L on impute() and fit_baseline().
When >= 1, the threshold-joint path (binary + ordinal + OVR
categorical) iterates: the BM rate Σ learned by
Rphylopars::phylopars() at iteration k is fed back as the
per-trait prior SD at iteration k + 1, up to em_iterations times
(default 5 when enabled) or until
||Σ_k - Σ_{k-1}||_F / ||Σ_{k-1}||_F < em_tol (default 1e-3). Closes
the final open item from the v0.9.0 "Deferred to future releases" list.em_iterations = 0L preserves v0.9.1 output byte-for-byte.Σ is deferred to a future Phase 7.
Phase 6 uses sqrt(diag(Σ)) only.pigauto 0.9.1 is a consolidation release on top of 0.9.0. Two headline
improvements: (a) tree-sharing GNN is now the default in
multi_impute_trees() for historical descriptive sensitivity work; and (b) opt-in
calibration smoothers (bootstrap conformal, median-over-splits gate,
low-val-cells warning) for small-n regimes. This release also ships the
three post-0.9.0 feature branches (B1 soft-liability, B2 rate-aware
attention, B3 full-threshold ordinal baseline), a 9,993-species AVONET
missingness sweep, a separate single-replicate timing curve through 5,000
tips, and a clean R CMD check (0 errors / 0 warnings / 1 note, down from
1 / 4 / 3 on v0.9.0).
fit_pigauto() argument min_val_cells = 10L. Warns at fit time
if any trait had fewer than this many validation cells available for
gate calibration and conformal-score estimation. Default of 10 is the
floor of pathological territory (operational target is 20-30).gate_method = "median_splits" runs the half-A / half-B
grid search on gate_splits_B = 31L random splits and takes the
median best-gate. Small empirical benefit at small n_val (gate SD
0.406 → 0.360 across 10 seeds on the calibration sim). Default stays
"single_split" for backward compat.conformal_method = "bootstrap" averages the conformal
quantile across conformal_bootstrap_B = 500L bootstrap resamples of
the val residuals. Empirically reduces conformal-score variance ~30%
but does not reliably improve 95%-interval coverage stability on
the evaluation sim. Shipped as experimental; default stays
"split".fit_pigauto() gains a @section Calibration at small n:
explaining when the 95% intervals should be treated as approximate.
Short version: at n_val < 20 the split-conformal guarantee degrades;
increase missing_frac or collect more species rather than relying
on estimator smoothing.share_gnn = TRUE in multi_impute_trees(). The GNN is
trained once on a reference tree (MCC via phangorn::maxCladeCred,
fallback trees[[1]]) and reused across posterior trees, with only
the BM baseline recomputed per tree. This historical feature is now
descriptive sensitivity only: do not send its completions to downstream
fitting or pooling, and do not interpret them as MI, Rubin results, SEs,
or FMI.m_per_tree setting remains part of the API but does not
create an inference-supported tree workflow.share_gnn = FALSE restores the pre-v0.9.1 per-tree fit
path for users who need exact per-tree model independence.reference_tree lets users override the MCC choice.predict.pigauto_fit() gains a baseline_override argument (mostly
internal — used by the shared-GNN path).phangorn (for maxCladeCred()).decode_ordinal_liability() converts the
Rphylopars liability posterior back to z-scored integer class for
downstream GNN compatibility.fit_baseline() threshold-joint dispatch now fires when ordinal
cols are present (in addition to binary / categorical). Ordinal
cols that can't be populated (e.g., <2 observations) fall back to
per-column BM unchanged.GraphTransformerBlock attention bias now uses learnable per-head
Gaussian bandwidth computed from raw squared cophenetic distances:
bias_h = -D_sq / (2 * softplus(log_bw_h)^2). Each head can learn
its own phylogenetic scale — tight for fast-evolving traits, broad
for conserved traits. Replaces the fixed adj_bias_scale * log(adj)
formulation.build_phylo_graph() returns D_sq (squared cophenetic distances)
alongside adj. Persists through training (not freed like D).D_sq = NULL (pre-B2 saved fits), the old
log(adj) path fires unchanged. Legacy attention path (use_transformer_blocks = FALSE) is completely unaffected.multi_obs_aggregation = c("hard", "soft") on impute()
and fit_baseline(). Default "hard" preserves v0.9.0 behaviour.
"soft" uses Rao-Blackwell convex-combination E-step for binary and
categorical traits: a species with 6/10 class-1 observations now
produces a less-extreme posterior liability than 10/10, fixing the
threshold-at-0.5 / argmax-one-hot information loss from Phase 10.estep_liability_binary_soft(p, mu_prior, sd_prior) and
estep_liability_categorical_soft(p_vec, mu_prior, sd_prior) in
R/liability.R. Both reduce to their hard counterparts at boundary
(one-hot) proportions; both return the prior mean at maximum ambiguity
(p = 0.5 for binary, uniform for categorical).aggregate_to_species() gains soft_aggregate flag: when TRUE,
preserves species-level proportions instead of thresholding/argmax.script/bench_avonet_missingness.R).
Mean / mode vs BM baseline vs pigauto at 20% / 50% / 80% missing, run
from Compute Canada (Narval; one run per missingness setting). At 80%
missing, continuous RMSE matches BM for all four traits. Across the two
categorical traits and three settings, pigauto accuracy differs from BM
by −4.1 to +4.2 percentage points; the direction varies by trait and
setting. These are descriptive results from this run, with no uncertainty
estimates from repeated runs. Output: bench_avonet_missingness.rds +
bench_avonet_missingness.md, rendered into the pkgdown validation
suite.script/bench_scaling_v090.R).
This v0.9.0 run measured per-stage wall time and R heap use at
n ∈ {100, 300, 1,000, 3,000, 5,000}, with four traits and 25% MCAR
missingness. It completed the 5,000-tip case in 23.1 minutes overall;
graph construction used about 1,971 MB of R heap at that size. The
recorded run did not include 7,500 or 10,000 tips, so it does not support
a current 10,000-tip runtime claim. Results are in
script/bench_scaling_v090.md and script/bench_scaling_v090.log.bench_covariate_sim.R rerun on v0.9.0 (previously pinned to v0.6.0).
Now the reference entry for the environmental-covariates path in the
validation suite.bench_multi_obs_mixed.R output regenerated with the B1 soft E-step
path enabled; captures the +12–28pp categorical lift under soft
aggregation vs hard.vignettes/tree-uncertainty.Rmd for
the current descriptive boundary.zi_count and multi_proportion
are now listed everywhere alongside the other six trait types
(Getting Started, Mixed Types, the validation suite, and the README).R CMD check clean: 0 errors / 0 warnings / 1 note, down from
1 / 4 / 3 on v0.9.0. The remaining note is the BACE in-tree
"Package sources not in a standard location" flag, which is a
deliberate .Rbuildignore choice so BACE ships as a comparison
reference without being part of the installed pigauto package.submit_v090_cloud/ (Narval;
renamed clusters after the 2026 Alliance Canada infrastructure
renewal — Cedar → Fir etc.) ships one-shot scripts for the AVONET
missingness sweep and the scaling curve. submit_v090_vulcan/ (PAICE
account aip-snakagaw) ships a SLURM array driver for the 60-cell
calibration-coverage grid.R/, vignettes/, pkgdown/,
README.md, NEWS.md, DESCRIPTION). Saves ~50% of the
GitHub Actions minutes previously consumed by pigauto PR checks.bench_multi_obs_mixed.R with a portable relative path; refreshed
the bench_multi_obs_mixed .rds.pigauto 0.9.0 is a large release that ships two complementary upgrades: a Level-C phylogenetic baseline (captures cross-trait correlation and threshold structure), and a Graph Transformer GNN backbone (replaces the 2-layer single-head attention stack with a proper multi-head transformer). Both are gate-safe — pigauto still closes its per-trait gate when the baseline is already optimal, so high-phylo-signal datasets retain identical behaviour to v0.7.0.
A new family of joint multivariate baselines replaces per-column Brownian motion + label propagation on datasets where cross-trait correlation carries signal.
Liability encoding (R/liability.R). Every trait type now has an
explicit liability interpretation: continuous types are their own
liability; binary uses a 1D threshold at 0; ordinal uses K-1
thresholds; categorical uses K liabilities with an argmax constraint.
A single dispatcher estep_liability(tm, observed, mu_prior, sd_prior)
returns posterior mean/variance of the underlying liability given an
observation. This is the foundation the joint paths below build on.
Joint multivariate BM baseline (R/joint_mvn_baseline.R). When
Rphylopars is installed and the dataset has ≥2 BM-eligible latent
columns, the baseline now fits a single joint multivariate BM across
continuous, count, ordinal, proportion, and ZI-magnitude cols via
Rphylopars::phylopars(), capturing phylogenetic correlation across
traits. Benchmark on correlated simulated BM data:
33.7% RMSE lift vs per-column BM (bench_joint_baseline.R).
Threshold-joint binary baseline (R/joint_threshold_baseline.R).
Binary traits get a truncated-Gaussian E-step posterior mean via
estep_liability_binary and join the joint MVN alongside continuous
traits. Decoded back to logit(P(y=1)) so BCE and gate math are
unchanged. Benchmark: 3/3 scenarios with ≥2pp accuracy lift over the
LP baseline, peaking at +16pp at rho = 0.6
(bench_binary_joint.R).
OVR categorical baseline (R/ovr_categorical.R). Each K-class
categorical trait is decomposed into K independent "is_class_k vs
rest" binary threshold fits, then normalised into a row-stochastic
distribution. Historical comparator figures for this path are withdrawn as
current evidence; the external comparator gate remains open. The K-independent-fits path
sidesteps the rank-(K-1) numerical instability that made earlier
single-fit approaches unstable; each individual fit has only one
categorical-related column.
ZI-count gate routing. The binary gate col of zero-inflated count traits now joins the threshold-joint path alongside binaries (the magnitude col already rode the Phase-2 joint MVN).
Graceful fallback throughout. Any Level-C path that can't fit (Rphylopars missing, too few observations per column, multi_proportion trait present) falls through to the legacy per-column BM + LP paths without user-visible breakage.
Graph Transformer backbone (R/graph_transformer_block.R). The
2-layer single-head attention stack is replaced by n_gnn_layers
instances of a GraphTransformerBlock — multi-head attention (default
4 heads) with learnable per-head log(adj + eps) bias (so the phylo
prior is respected by each head independently), feed-forward network
(default 4× width), pre-norm, two residual skips. FFN output
projection initialises to zero so each block starts near-identity,
preserving pigauto's gate-closed-at-init safety.
New hyperparameters on ResidualPhyloDAE and fit_pigauto():
use_transformer_blocks = TRUE (new default), n_heads = 4L,
ffn_mult = 4L. The legacy architecture is fully preserved behind
use_transformer_blocks = FALSE and reproduces pre-0.9.0 numbers
exactly. Pre-0.9.0 saved pigauto_fit objects reconstruct via the
legacy path automatically.
Empirical characterisation. On the 4-scenario discriminative benchmark the transformer and legacy GNN produce identical predictions in high-signal regimes (gate closes to zero — architecture below the gate is irrelevant) and the transformer shows small RMSE lift in moderate/low signal scenarios where the gate stays partly open. Runtime is ~25–30% slower per fit. Larger transformer gains likely need self-supervised pretraining on a large tree corpus (future work).
Multi-obs data now uses Level-C paths. Previously the joint MVN,
threshold-joint, and OVR categorical paths refused multi-obs input via
four !multi_obs guards in fit_baseline(). Those guards are all
removed; multi-obs data is aggregated to species level inside each
Level-C helper via a new internal aggregate_to_species() helper.
Aggregation rules: mean for continuous-family, threshold-at-0.5 for
binary/ZI-gate, argmax-one-hot for categorical. Splits are
mask-then-aggregated so val/test leakage is impossible; single-obs
data passes through unchanged.
Benchmark (bench_multi_obs_mixed.R). On a synthetic multi-obs
mixed-type simulation (150 species, 2 continuous + 1 binary + 1 K=4
categorical, Poisson-5 obs per species), the Level-C path beats
the LP path by +12–28pp categorical accuracy and +5–7.5pp
binary accuracy across all three phylogenetic-signal regimes.
Aggregation caveat. Binary/categorical aggregation is lossy — a species with split observations (e.g., 6/10 class-1) becomes a single class-1 observation at species level. For datasets with strong within-species variability on discrete traits, prefer single-obs mode with a user-computed species-level summary.
Nine new or reran benchmark reports quantifying the release's impact:
bench_joint_baseline.R — Phase-2 joint MVN A/B on correlated BM data.bench_binary_joint.R — Phase-3 threshold-joint A/B on binary traits.bench_discriminative.R + bench_discriminative_phase9.R — 4-scenario
mixed-signal sweep, including the transformer vs legacy GNN A/B.bench_avonet_phase6.R — post-Phase-6 AVONET 300 vs LP baseline.bench_categorical_joint.R — OVR K-fits A/B on synthetic categorical
data.bench_covariate_sim.R — rerun on 0.9.0. Shows −10.6% RMSE on the
low-phylo-strong-env scenario vs pre-Phase-6. Pre-0.9.0 snapshot in
bench_covariate_sim_preP6.md.bench_multi_obs.R — rerun on 0.9.0. Numbers unchanged from pre-0.9.0
(bench is single-trait so it can't exercise Level-C dispatch
conditions; bench_multi_obs_mixed.R is the correct validator).bench_multi_obs_mixed.R — new synthetic mixed-type multi-obs
benchmark specifically to validate Phase 10's multi-obs unlock.joint_mvn_available() gates the Level-C dispatch — returns TRUE
iff Rphylopars is installed. All Level-C paths fall back gracefully
when it returns FALSE.aggregate_to_species() is the sole multi-obs → species-level helper
and is the right place to hook future soft-liability aggregation
(currently threshold-at-0.5 for binary, argmax for categorical).test-liability.R,
test-joint-threshold-baseline.R, test-joint-baseline.R,
test-ovr-categorical.R, test-graph-transformer-block.R,
test-phase9-integration.R, and test-multiobs-levelc.R. Full suite
is 778 PASS, 0 FAIL.script/).multi_proportion: a new trait type for compositional data — K columns
per row that sum to 1 (e.g. plumage-colour proportions, diet composition,
microbiome relative abundances, allele frequencies). Declared via a new
multi_proportion_groups = list(colour = c("black", ..., "yellow"))
argument to impute(), preprocess_traits(), and multi_impute().
Encoding: centred log-ratio (CLR) + per-component z-score. Small
epsilon (1e-6) is applied to zero components before log to keep the
transform safe; rows are re-normalised so the composition is preserved.
Baseline: K independent Brownian-motion imputations on the CLR columns. BM is well-defined on CLR space (Euclidean), unlike the raw simplex.
GNN output: K logits projected to sum-zero, then softmax at decode time — guaranteed to lie on the simplex.
Loss: MSE on CLR columns (averaged across K), applied group-wise so the whole composition is masked together during DAE corruption.
Metrics: Aitchison distance (the natural compositional-data metric), RMSE on CLR space (comparable to continuous traits), and simplex MAE (mean absolute error on decoded proportions).
Benchmark: script/bench_multi_proportion.R with 5-scenario
signal sweep and K ∈ {3, 5, 8, 12} secondary sweep. Report rendered
by script/make_bench_multi_proportion_html.R.
Tests: tests/testthat/test-multi-proportion.R (5 test blocks,
20 expectations) — encoding, validation, baseline, end-to-end impute().
evaluate_imputation() output adds three columns: aitchison,
rmse_clr, simplex_mae. Existing metrics columns are unchanged.
impute() gains explicit trait_types and multi_proportion_groups
arguments (previously trait_types was documented as passed through
... but was not actually forwarded; this fixes that).
multi_impute_trees() now
supports descriptive prediction sensitivity only; its completions must not
be passed to downstream fitting or pooling and do not provide MI, Rubin,
standard-error, or FMI output.n = 300 / 5,000 / 10,000 species with T = 10 / 50
trees, plus guidance on reducing T for large trees and
parallelising across machines.MC dropout zero variance fixed (R/predict_pigauto.R): when
n_imputations > 1, the blend equation now uses a BM posterior draw
t_BM_draw ~ N(BM_mu, BM_se) instead of the deterministic t_MU.
When the calibrated gate is zero (BM dominates), each imputation now
samples from the correct BM posterior rather than returning the same
point estimate, giving non-zero between-imputation variance.
Conformal draw scale bug fixed (R/multi_impute.R): for
log-transformed traits (Mass, Wing length, etc.), draws are now made
on the latent z-score scale and back-transformed, avoiding near-zero
variance caused by dividing a log-scale SE by an original-scale value.
The same fix applies to log1p-scaled counts and logit-scaled
proportions.
draws_method = "conformal" is now the default for multi_impute().
Historical wording about downstream pooling is superseded: these draws are
prediction diagnostics and are not supported for downstream inference.draws_method = "mc_dropout" remains available and is now correct
(see bug fix above).man/figures/logo.png: pink pig driving a car with a 7-tip
ultrametric pectinate cladogram above the pig's head. Forest green /
mint palette (v3) selected as the official sticker.When datasets have multiple observations per species measured under
different conditions (e.g. CTmax at different acclimation temperatures),
the GNN now produces covariate-conditional predictions within species.
A new refinement MLP (obs_refine) re-injects user-supplied covariates
after species-level phylogenetic message passing, so that different
observations of the same species receive different predicted values
based on their covariate context. This is controlled by a residual
connection that preserves the phylogenetic signal from message passing.
R/model_residual_dae.R: new n_user_cov parameter and obs_refine
MLP (activated only when species_col + covariates are both supplied)R/fit_pigauto.R: n_user_cov stored in model configR/predict_pigauto.R: backward-compatible model reconstructionSimulated multi-observation-per-species CTmax data (1,464 observations
across 300 species, with acclimation temperature as an observation-level
covariate). 30% of species are entirely unobserved. Uses tree300.
See data-raw/make_ctmax_sim.R for the generation script.
script/bench_multi_obs.R simulates CTmax-like datasets under varying
phylogenetic signal and within-species acclimation response ratios,
comparing species_mean, pigauto (no covariates), and pigauto (with
covariates). The observation-level refinement gives measurable RMSE
improvement when within-species covariate effects are strong.
pigauto now explicitly combines three sources of information for
imputation: (1) the phylogenetic tree, (2) cross-trait correlations,
and (3) optional environmental covariates. The covariates argument
is accepted by impute(), multi_impute(), multi_impute_trees(),
and the lower-level pipeline functions. Covariates are threaded
through the GNN with the same gated safety that protects against
GNN degradation — they only contribute when they demonstrably
improve accuracy.
The Brownian-motion baseline for continuous, count, ordinal, and
proportion traits is now computed internally using conditional
multivariate normal imputation on the phylogenetic correlation matrix
R = cov2cor(vcv(tree)) (Goolsby et al. 2017). This removes the
hard dependency on Rphylopars, which is now in Suggests only.
The internal implementation (R/bm_internal.R) uses univariate
per-column imputation with Cholesky decomposition and nugget
regularisation for near-singular submatrices. Validation against
Rphylopars shows r = 0.97–0.98 (expected given univariate vs
multivariate difference).
This historical feature description is superseded. Posterior-tree completions may be compared as descriptive prediction sensitivity, but they are not multiple-imputation datasets and must not be used for Rubin pooling, standard errors, or FMI.
trees300: 10 posterior trees from the BirdTree Hackett backbone,
pruned to the same 300 tips as tree300delhey5809: plumage lightness and 4 environmental covariates
(absolute latitude, mean temperature, mean precipitation,
tree cover) for 5,809 bird species (Delhey 2019)tree_delhey: matching BirdTree phylogeny for delhey5809Thirteen benchmark reports are now accessible from the pkgdown site under Articles > Benchmarks:
Each report is a self-contained HTML page generated from
script/make_bench_*_html.R and backed by reproducible
script/bench_*.R drivers.
DOCS.md updated with new benchmarks, datasets, and functionsfit_baseline_tabpfn() and setup_tabpfn() are gone, along with the
reticulate suggested dependency, the script/bench_tabpfn.* files,
script/run_tabpfn.py, tests/testthat/test-tabpfn.R, and the
dev/bench_tabpfn.html benchmark page. The TabPFN wrapper was
originally added as a curiosity — a test of whether a general-purpose
tabular transformer could compete with pigauto as a drop-in baseline.
In practice it is a poor fit for what pigauto is trying to do:
n = 9,993 species (the full
AVONET dataset) this is already near the edge of what a single
GPU can handle, and most users of pigauto care precisely about
scales at which a baseline must be cheap.torch and Python torch share libtorch and cannot coexist in
the same process, so TabPFN had to run via system2() on a
separately-installed Python venv. Every time the user's Python,
tabimpute, or torch versions drifted, the wrapper broke. This
is not maintenance the project can afford.Users who still want to benchmark pigauto against TabPFN can install
tabimpute separately and call it directly from Python. The removal
simplifies pigauto's installation, CI, and documentation.
DOCS.md: rewritten benchmark and developer-reports sections to
reflect the v0.3.2+ benchmark set (scaling to 10,000 tips and the
AVONET missingness sweep). Dropped the obsolete benchmark_report
and benchmark_final_report links that v0.3.2 had already pruned
from _pkgdown.yml but had not yet removed from DOCS.md.CLAUDE.md: removed the TabPFN Python-venv setup recipe and the
"R torch ↔ Python torch cannot share a process" gotcha (no longer
applicable). Version line bumped to 0.4.0.The multiple-imputation example in README.md, the two HTML
tutorials (pigauto_intro.html, pigauto_workflow_mixed.html), and
tests/testthat/test-multi-impute.R all constructed the phylogenetic
random-effect structure with
Vphy <- ape::vcv(tree) # WRONG: covariance, diag = tree height
glmmTMB(... + propto(0 + species | dummy, Vphy), ...)
propto() estimates sigma^2 freely, so passing the raw covariance
matrix (whose diagonal is the tree height) silently rescales the
phylogenetic variance component by tree height. The correct usage
passes a correlation matrix:
Vphy <- cov2cor(ape::vcv(tree)) # diag = 1
All four sites are now fixed. Thanks to the user who flagged this.
The user-facing prose throughout the package historically described pigauto as a Residual Phylogenetic Denoising Autoencoder (ResidualPhyloDAE) that learns a residual correction on top of a BM baseline. This wording was misleading for binary and categorical traits in two ways:
y - baseline. So the GNN's output delta is a full per-cell
prediction, not a statistical residual.What pigauto actually is, in honest terms, is a gated ensemble of two predictors:
delta from spectral node
features and an attention-biased message passing over the
phylogeny).Prediction is the per-trait blend
(1 - r_cal) * baseline + r_cal * delta, where r_cal is calibrated
on a held-out validation split (with a split-validation cross-check
for discrete traits). When the baseline is already optimal the gate
closes and the GNN becomes a no-op.
The internal torch class name ResidualPhyloDAE is preserved because
its GNN layers use ResNet-style residual skip connections — that is a
correct, well-defined ML term. But user-facing prose no longer
describes the GNN as "learning a residual on top of the baseline".
Files updated: README.md, DESCRIPTION, CLAUDE.md, _pkgdown.yml,
DOCS.md, vignettes/getting-started.Rmd, vignettes/mixed-types.Rmd,
R/fit_pigauto.R (roxygen), R/predict_pigauto.R (roxygen),
R/model_residual_dae.R (clarifying header comment),
R/simulate_traits.R, R/plot.R (plot titles), and
script/make_intro_html.R.
The honest fix above is a documentation correction, not a model
change. A genuinely residual baseline for binary and ordinal traits
would require a liability-scale model (threshold BM, phylolm::phyloglm,
or MCMCglmm with family = "threshold") — a design decision that
should be made deliberately rather than rolled into a same-day patch.
A v0.4.0 tracking issue lays out the trade-offs.
The full AVONET3 + BirdTree Stage2 Hackett MCC phylogeny, aligned to
9,993 species with 7 mixed-type traits (4 continuous morphometric +
2 categorical + 1 ordinal) and native missingness preserved.
Complements the 300-species avonet300 / tree300 bundled data:
avonet300 / tree300 remain the right choice for quick examples,
unit tests, and the getting-started vignette.avonet_full / tree_full are the right choice for scale
benchmarks, realistic-missingness experiments, and anything that
wants to exercise the v0.3.1 sparse Lanczos / cophenetic caching
code paths on a real-world phylogeny.The schemas are identical, so any code that runs on avonet300 runs
on avonet_full with no modification. Combined tarball increment is
~328 KB (xz-9).
Articles -> Benchmarks -> AVONET missingness sweep is a new pkgdown
page comparing three methods on the full avonet_full dataset at
three MCAR missingness levels:
Per-trait RMSE, Pearson r, Spearman rho, and accuracy on held-out
test cells. Single seed per cell, trait-level masking so categorical
traits are held out as whole groups. Hyperparameters match the
validate_avonet_full.R scaling run for comparability.
pigauto_intro and pigauto_workflow_mixed HTML
tutorials against the v0.3.2 package state.benchmark_report.html, benchmark_final_report.html,
test_report.html) and from pkgdown/assets/dev/.script/ (trial artefacts from
pre-v0.3.1 benchmarking phases: bench_v2/v3/v4.*,
benchmark_avonet2000.*, benchmark_final.*,
benchmark_simulation.*, scaling duplicates, and Apr-6 demo
artefacts).README.md citation now points at 0.3.2.DESCRIPTION updated to mention both bundled datasets.pigauto now runs end-to-end on real 10,000-tip phylogenies. Before this
release the pipeline was tested and tuned around the 300-species bundled
AVONET data, and at 10,000 tips it wedged on a single dense eigendecomposition
of the graph Laplacian (hours of CPU time). Two targeted fixes land in this
release; no rewrite was needed.
build_phylo_graph()build_phylo_graph() now uses RSpectra::eigs_sym() sparse Lanczos to
compute the k + 1 smallest Laplacian eigenvectors when the tree has more
than 7,500 tips and RSpectra is installed. This replaces
base::eigen(L, symmetric = TRUE), which was O(n^3) in time and discarded
99% of its work (we only need the bottom ~32 eigenvectors).
ape::rcoal) have degenerate spectra that need
large Krylov subspaces to resolve.Suggests. Installing it is strongly recommended
for any tree larger than ~7,500 tips.The cophenetic distance matrix D = ape::cophenetic.phylo(tree) was previously
computed four separate times in a single impute() call (three times inside
build_phylo_graph() and once inside fit_baseline()). Each call allocates a
dense n*n double matrix and runs an O(n^2) tree traversal -- wasted work at
scale. In 0.3.1:
build_phylo_graph() computes D exactly once and returns it as part of
the result list (graph$D).fit_baseline() gains an optional graph argument; when supplied, it
reuses graph$D instead of recomputing the cophenetic matrix.impute() and fit_pigauto() pass the graph through automatically, so
users get the speedup without any code changes.impute() also explicitly releases graph$D after fit_baseline()
finishes, so the ~800 MB cophenetic matrix does not sit in R memory
during training (a subtle side-effect that caused a large per-epoch
regression at n = 10,000 until it was dropped).See pkgdown/assets/dev/scaling.html (linked from the site as
Articles -> Benchmarks -> Scaling to 10,000 tips) for a side-by-side
comparison of v0.3.0 vs v0.3.1 at n in {300, 1k, 2k, 3k, 5k, 7.5k, 10k},
plus an end-to-end validation run on the full AVONET3 + BirdTree dataset
(9,993 species, 7 mixed-type traits). The pkgdown article also walks
through both a quick 300-species example with the bundled avonet300
data and a full 10k-species example using the BirdTree Stage2 Hackett
MCC phylogeny.
build_phylo_graph(): the N > 2000 memory warning was out of date (the
real blocker above that size was wall-clock time, not memory). The
warning is now issued at N > 10000 and its text points at the true
bottleneck.The GNN now uses learned attention weights (with phylogenetic adjacency
as positional bias) instead of fixed adjacency matrix multiplication.
The model learns which phylogenetic neighbours are most informative for
each species, while the log-adjacency bias ensures initialisation
respects evolutionary distances. Controlled via use_attention = TRUE
(now the default).
After training, per-trait gate values are optimised on the validation set via grid search. This prevents the GNN from adding noise when the Brownian motion baseline is already optimal. On the AVONET benchmark, this eliminates the RMSE regression seen in earlier versions.
This historical statement about guaranteed coverage is superseded. Current
conformal bounds are nominal held-out diagnostics, not package-wide certified
coverage. They remain available for continuous, count, and ordinal traits via
pred$conformal_lower and pred$conformal_upper.
Binary and categorical trait baselines now use phylogenetic similarity- weighted label propagation instead of flat population-level frequencies. Each species gets a personalised baseline reflecting its phylogenetic neighbourhood.
k_eigen now defaults to "auto", scaling with tree size:
min(max(ceil(n/20), 4), 32). This gives 4 features for small trees,
15 for 300 tips, and 32 for 640+ tips.
pigauto_report() produces self-contained HTML reports with interactive
Chart.js visualisations: per-trait metrics, calibrated gate values,
training history, and conformal coverage.
evaluate(): per-trait metrics on held-out datacompare_methods(): automated BM vs GNN comparisoncross_validate(): proper k-fold cross-validationsummary.pigauto_fit(): formatted console summaryplot.pigauto_fit(): training history, gate values, conformal scoresplot.pigauto_pred(): scatter plots, conformal intervals, probability
distributionsplot_comparison(): forest-plot style method comparisonsimulate_benchmark() runs a complete simulation study with configurable
scenarios (BM, OU, regime shift, non-linear, mixed types), returning tidy
results with print/summary/plot methods. Useful for methods papers and
for verifying pigauto's behaviour on data with known properties.
build_phylo_graph() default k_eigen changed from 8L to "auto"fit_pigauto() default k_eigen changed from 8L to "auto"all_grads_finite() now handles non-standard tensor types gracefullyimpute() uses adaptive k_eigen by defaultR CMD check with 0 errors and 0 notesspecies_col)