Fit the phylogenetic baseline — fit

Dispatches to pigauto's phylogenetic baseline machinery and returns imputed latent-scale means and standard errors for every species.

Usage

fit_baseline(
  data,
  tree,
  splits = NULL,
  model = "BM",
  graph = NULL,
  multi_obs_aggregation = c("hard", "soft"),
  lambda_mode = c("estimate", "fixed_1", "cv", "bayes"),
  lambda_fixed = NULL,
  em_iterations = 0L,
  em_tol = 0.001,
  em_offdiag = FALSE,
  joint_solver = c("inhouse", "rphylopars"),
  predict_method = c("auto", "exact", "per_column"),
  joint_refine_iter = 0L,
  predict_route = NULL,
  seed = NULL
)

Arguments

data: object of class "pigauto_data".
tree: object of class "phylo".
splits: list (output of make_missing_splits) or NULL.
model: character. Evolutionary model: "BM" (default) or "OU".
graph: optional list returned by build_phylo_graph. When supplied, graph$D (cophenetic distances) is reused for label propagation and graph$R_phy (phylogenetic correlation matrix) is reused for BM imputation, avoiding duplicate $O(n^2)$ allocations. When NULL (default), both matrices are computed here.
multi_obs_aggregation: character. How to aggregate multiple observations per species before the Level-C joint baseline: "hard" (default) thresholds binary proportions at 0.5 and uses argmax for categorical, matching Phase 10 behaviour. "soft" preserves species-level proportions and dispatches the truncated-Gaussian soft E-step (estep_liability_binary_soft) so that intermediate class frequencies contribute fractional liability evidence. Only relevant for multi-obs data with binary or categorical traits when the Level-C joint baseline is active.
lambda_mode: character. Pagel-lambda mode for the CONTINUOUS-FAMILY baseline (continuous, count, proportion, zi_count magnitude columns; NOT ordinal, which the threshold-joint path treats as a liability column via estep_liability_ordinal() – see “Per-type lambda dispatch” in Details). "estimate" (default, S4) fits each continuous-family column's own Pagel's lambda via profile REML; when the joint MVN or threshold-joint baseline fires, this now runs INSIDE that joint fit (R/joint_mvn_solver.R's per-column lambda_cols) rather than being discarded in favour of a separate per-column re-fit. "fixed_1" preserves the classic Brownian correlation matrix (lambda = 1 everywhere). "cv" and "bayes" have no joint analogue and force the per-column BM path for continuous-family columns (as before); binary/ordinal/categorical/ zi_gate columns are unaffected by lambda_mode in every case and keep the threshold-joint / OVR-categorical baseline at lambda = 1. Covariate caveat: when data$covariates is supplied, the per-column path switches to bm_impute_col_with_cov(), which accepts a numeric lambda or "estimate" but not "cv" / "bayes"; those two modes are silently ignored (fit at lambda = 1) for BM-eligible columns and a warning is emitted.
lambda_fixed: optional named numeric vector (names = latent column names, i.e. colnames(data$X_scaled)) giving a FIXED lambda per continuous-family column, overriding lambda_mode entirely for those columns (spec 4.5: reproduce a previous fit's exact per-trait lambda at predict time without re-estimating – typically supply that fit's own $lambda_per_trait). Columns not present in lambda_fixed keep their lambda = 1 default. NULL (default) means every continuous-family column follows lambda_mode normally.
em_iterations: integer. Number of Phase 6 EM iterations for the threshold-joint baseline (binary + ordinal + OVR categorical). Default 0L disables the EM loop and preserves v0.9.1 output byte-for-byte. When >= 1, the BM rate $\Sigma$ learned by the in-house joint solver (R/joint_mvn_solver.R) at iteration $k$ is fed back as the per-trait prior SD at iteration $k+1$, up to em_iterations times or until em_tol convergence. em_iterations = 1L is a degenerate single-pass run and produces the same baseline output as 0L; >= 2L is needed for actual iteration. Only affects the threshold-joint path (continuous-only traits pass through the existing joint MVN path unchanged).
em_tol: numeric. Relative-Frobenius convergence tolerance for the Phase 6 / 7 EM loop. Early-stops when $||\Sigma_k - \Sigma_{k-1}||_F / ||\Sigma_{k-1}||_F < $ em_tol. Default 1e-3.
em_offdiag: logical. Phase 7 opt-in: when TRUE AND em_iterations >= 2L, each liability cell's prior at iteration $k+1$ is the conditional-MVN $(\mu, sd)$ given the posterior liability of other traits at iteration $k$, using the full off- diagonal entries of $\Sigma$. Binary + ordinal only (OVR categorical stays on Phase 6 diagonal). Default FALSE preserves Phase 6 behaviour.
joint_solver: character. Which solver estimates the joint Sigma / posterior for the joint MVN, threshold-joint, and OVR categorical baselines. "inhouse" (default) uses the single-pass in-house solver (R/joint_mvn_solver.R); under lambda_mode = "fixed_1" it is byte-identical to prior releases, but lambda_mode = "estimate" (the current default) is new behaviour, not a byte-compatibility guarantee. "rphylopars" delegates to Rphylopars::phylopars(); under "estimate" this calls it with model = "lambda" (substantially slower than model = "BM"), with automatic fallback to "inhouse" on failure or implausible output (a plausibility guard that also fires under lambda_mode = "fixed_1" – see NEWS). Only affects the joint MVN / threshold-joint / OVR categorical paths above; ignored when those paths don't fire. lambda_mode does NOT disable the continuous-only joint MVN path: both fit_joint_mvn_baseline() and fit_joint_threshold_baseline() accept a lambda_mode / lambda_fixed argument and estimate lambda inside the joint fit via lambda_cols – see “Per-type lambda dispatch” in Details. Only "cv" / "bayes" force continuous-family columns off the joint path entirely (no joint analogue for those two modes).
predict_method: character. Prediction route for the in-house joint solver. "auto" (default, S5b/S5c) fits the baseline once with the "exact" route and once with the "per_column" route, then for each TRAIT (each trait_map entry; a categorical trait's K latent columns and a zi_count trait's gate + magnitude columns are chosen together) picks whichever route has the lower loss on that trait's cells: mean squared error on the z-scored latent scale for continuous/count/proportion/ordinal/zi_count magnitude, mean log-loss of plogis(mu) against the observed 0/1 truth for binary/zi_count gate, and mean multinomial log-loss for categorical. Which cells the choice sees: a trait's held-out validation rows (species in multi-observation data) are split into two halves (seeded via seed, or the ambient RNG state when seed = NULL) when the trait has at least 38 of them and the two candidate fits predict it differently. The route choice sees only half A; half B is reserved so fit_pigauto()/impute() calibrate the GNN gate and compute conformal scores on rows that never informed the choice. Below 38 rows each half would keep fewer than the 19 a 95\ split-conformal interval needs, so the same rows both choose the route and calibrate; when both fits agree, no choice is made and all rows calibrate. With fewer than 5 validation cells, or on a genuine tie, a trait keeps "exact", as it does when splits is NULL (no validation cells at all makes every trait "exact"). The chosen route per trait is returned as $predict_method_by_trait. "exact" uses the full cross-trait conditional mean and variance of vec(L) ~ MVN(0, Sigma %x% R(lambda_block)) in the sparse precision form (Hadfield & Nakagawa, 2010), with each column centred at its own GLS phylogenetic mean at lambda_block before the solve (mean-model consistency; see docs/dev-log/exact-default/). Discrete liability columns (binary, zi gate, ordinal-via-OVR synthetic columns) share lambda_block under "exact" rather than staying fixed at lambda = 1, so the whole joint fit uses one internally-consistent R(lambda_block). Falls back to "per_column" above roughly 20000 unknown cells (roughly 4000 species at 5 traits), on a singular/unusable Sigma, or when fewer than 2 joint columns or no Henderson sparse precision are available; the fallback prints a one-time message() (not a warning) per R session when predict_method was left at its default, or a warning() every time when "exact" was requested explicitly. "per_column" retains the original per-column conditional prediction route (each column's own posterior, no cross-trait borrowing in the prediction step; se is that column's own conditional SE). Neither option changes covariance estimation or the "rphylopars" solver.
joint_refine_iter: integer, default 0L. Enables cross-trait refinement of the joint baseline's cell imputations using the estimated Sigma (the in-house solver's max_iter EM cell-refinement; R/joint_mvn_solver.R). 0L preserves current behaviour byte-for-byte. The refinement is guarded: the Sigma step must shrink each iteration, or the loop rolls back to the last good iterate and sets $diverged. Has no effect on any trait predicted by the "exact" route: the exact conditional solve returns its analytic posterior directly and never enters this EM loop. Under the default predict_method = "auto", joint_refine_iter therefore only ever applies to traits that "auto" routes to "per_column" (see $predict_method_by_trait below).
predict_route: optional named character vector (names = trait names, i.e. vapply(data$trait_map, function(tm) tm$name, character(1)); values "exact" or "per_column"). Forces the prediction route for each named trait, bypassing the "auto" validation-loss comparison entirely (and overriding predict_method, if a value other than "auto" was also supplied). Traits not named default to "exact". A name that matches no trait in data$trait_map is ignored, with a warning (once per call). Used internally to replay a validation-split "auto" choice onto a production refit that has no validation cells (e.g. fit_pigauto()'s baseline_full, fit with splits = NULL) so that refit reuses the SAME per-trait route rather than re-deciding with zero validation evidence. NULL (default) disables forcing; predict_method governs normally.
seed: optional integer. Seeds the deterministic half-A (route choice) / half-B (reserved for fit_pigauto()/impute()'s gate calibration and conformal scoring) split of each trait's validation cells under predict_method = "auto" with real validation cells (see “"auto"” above). NULL (default) leaves the split to whatever the ambient RNG state is at call time, matching this package's other optional seed arguments. Ignored when predict_method is "exact" / "per_column", when predict_route is supplied, or when there are no validation cells.

Value

A list with:

mu: Numeric matrix (n_species x p_latent), baseline means in latent scale.
se: Numeric matrix (n_species x p_latent), standard errors.
path: Named character vector, one entry per trait_map entry (name = trait/group name), recording which dispatch branch produced that trait's baseline: "joint_mvn" (Phase 2 continuous-only joint), "threshold_joint" (Phase 3/5 binary + ordinal + continuous-family liability joint), "ovr_categorical" (Phase 6 one-vs-rest K-fits), "per_column_bm" (per-column Brownian-motion conditional MVN, including the lambda-aware path and the per-trait ordinal fallback that beats threshold-joint on held-out val MSE), "multi_proportion_bm" (per-component per_column_bm on the CLR columns of a multi_proportion group – same kernel as "per_column_bm", labelled separately because multi_proportion is architecturally distinct), "label_propagation" (phylogenetic label propagation for binary/categorical/ordinal, and the per-trait ordinal fallback that beats the other two candidates), "zi_gate_lp" (label propagation for a zero-inflated count's gate column), or "zi_mag_constant" (global mean/sd fallback for a zero-inflated count's magnitude column when fewer than 5 non-zero observations exist). For zi_count traits (2 latent columns: gate then magnitude), path reports the GATE column's dispatch – the magnitude column can independently land on "joint_mvn", "threshold_joint", "per_column_bm", or "zi_mag_constant" and is not separately surfaced here.
lambda_per_trait: Named numeric vector (length p_latent, names = colnames(data$X_scaled)), the lambda actually used for each latent column: 1 for every column not eligible for lambda estimation (all discrete columns; continuous-family columns under lambda_mode %in% c("cv", "bayes")), else the estimated or fixed value.
lambda_block: Numeric scalar, the shared lambda used internally by whichever joint fit ran for its Sigma M-step / opt-in exact conditional / opt-in EM refine (see R/joint_mvn_solver.R's $lambda_block); NA when no joint fit ran.
lambda_mode: Character, echoes the resolved lambda_mode argument.
predict_method_used: Character scalar, "exact", "per_column", or "auto". Under a concrete predict_method ("exact" / "per_column"), aggregated across every joint fit that ran: "exact" only if every joint fit achieved exact; "per_column" if any fell back, if none ran (a pure per-column baseline), or if every joint fit used joint_solver = "rphylopars" (no exact/per_column dichotomy there). Under predict_method = "auto" (or a forced predict_route), always "auto"; see $predict_method_by_trait for the resolved per-trait routes.
predict_method_by_trait: Named character vector (names = trait names), one entry per trait_map entry, giving the route ("exact" or "per_column") actually used for that trait's cells. Under a concrete predict_method, every entry equals predict_method_used. Under "auto" (or a forced predict_route), the per-trait validation choice (see predict_method's "auto" entry above).

Details

When splits is supplied the val and test cells are masked to NA before fitting, so the baseline is evaluated under the same conditions as fit_pigauto.

Continuous-family columns use Brownian-motion conditional MVN baselines on the phylogenetic correlation matrix, either independently or through the joint MVN path when the data and optional dependencies support it. Binary, ordinal, categorical, and zero-inflated gate columns use the appropriate label-propagation or threshold/liability baseline candidates, with per-column fallbacks when a joint path is not available.

Covariates and the joint baseline (P1-8): data$covariates is only used by the per-column BM path (bm_impute_col_with_cov()). The joint MVN and threshold-joint (Rphylopars) baselines do not accept a covariate design matrix, so when a joint path is selected (BM-eligible columns >= 2, or binary/ordinal cols present, with Rphylopars available) any supplied covariates are ignored for the BASELINE and a warning is emitted; covariates still reach the GNN correction via fit_pigauto regardless of which baseline path fires.

Per-type lambda dispatch (arc/lambda-per-type; corrected in S4/S6, feat/joint-lambda-default): lambda_mode only ever governs the baseline for CONTINUOUS-FAMILY columns (continuous, count, proportion, zi_count magnitude) – NOT ordinal. Binary, ordinal, categorical, and zero-inflated gate columns always stay at lambda = 1 in EVERY path (threshold-joint, OVR-categorical, label propagation, and the per-trait ordinal path-selection fallback below), regardless of lambda_mode – there is no discrete-trait analogue of Pagel's lambda, and previously forcing these columns onto label propagation any time lambda_mode != "fixed_1" cost 19pp of Trophic.Level accuracy on AVONET (0.789 -> 0.600; see docs/dev-log/2026-08-16-external-comparison-results.md). When lambda_mode != "fixed_1" and the threshold-joint baseline fires for a dataset with binary/ordinal AND continuous-family columns, the joint liability fit's own continuous-column OUTPUT is now USED (each continuous-family column gets its own lambda_k via the joint solver's lambda_cols), rather than being discarded for a separate lambda-aware per-column BM re-fit as in the pre-S4 design. Ordinal liability columns inside that same joint fit, and the BM-via-MVN / K-class-OVR-LP alternatives the "Per-trait ordinal path selection" block below compares against it, all stay at lambda = 1.

Examples

# \donttest{
data(avonet300, tree300, package = "pigauto")
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
traits <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                     c("Mass", "Wing.Length"), drop = FALSE]
rownames(traits) <- tree$tip.label
pd     <- preprocess_traits(traits, tree)
splits <- make_missing_splits(pd$X_scaled, trait_map = pd$trait_map)
bl     <- fit_baseline(pd, tree, splits)
# }