Skip to contents

Fit the phylogenetic baseline

Source: R/fit_baseline.R
fit_baseline.Rd

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("fixed_1", "estimate", "cv", "bayes"),
  em_iterations = 0L,
  em_tol = 0.001,
  em_offdiag = FALSE,
  joint_solver = c("inhouse", "rphylopars"),
  joint_refine_iter = 0L
)

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, ordinal, proportion, zi_count magnitude columns only – see “Per-type lambda dispatch” in Details). "fixed_1" preserves the default Brownian correlation matrix; "estimate", "cv", and "bayes" delegate lambda handling to the per-column BM path. Binary/categorical/zi_gate columns are unaffected by lambda_mode and keep the threshold-joint / OVR-categorical baseline. Covariate caveat: when data$covariates is supplied, the per-column path switches to bm_impute_col_with_cov(), which has no lambda argument and always fits at lambda = 1; lambda_mode != "fixed_1" is then silently ignored for BM-eligible columns and a warning is emitted.

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) and is byte-identical to prior releases. "rphylopars" delegates to Rphylopars::phylopars()'s converged REML fit, which measured 0.14-1.27 lower z-RMSE than the in-house solver on AVONET300 (see docs/dev-log/2026-08-16-continuous-gap-diagnosis.md); on failure or non-finite output it falls back to "inhouse" with a warning. Only affects the joint MVN / threshold-joint / OVR categorical paths above; ignored when those paths don't fire. Note that lambda_mode != "fixed_1" disables the continuous-only joint MVN path (it has no lambda argument) but no longer disables the threshold-joint / OVR-categorical paths – see “Per-type lambda dispatch” in Details.

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.

Where it helps, measured (docs/dev-log/2026-08-17-refinement-results.md): at Pagel's lambda around 0.2, joint_refine_iter = 3L lowered simulated RMSE by ~13\ near 1 it changes essentially nothing: on AVONET300, whose four continuous traits sit at lambda 0.993-0.998, every effect was within one Monte Carlo standard error. Note the reason is that refinement has nothing to add at high signal – not that the guard intervenes; on that data the guard does not fire. Use it when phylogenetic signal is moderate or low; expect no benefit on strongly conserved traits.

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.

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): lambda_mode only ever governs the baseline for CONTINUOUS-FAMILY columns (continuous, count, ordinal, proportion, zi_count magnitude). Binary, ordinal, categorical, and zero-inflated gate columns keep the threshold-joint / OVR-categorical joint baseline (always fit at lambda = 1) 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 still uses the continuous-family columns internally to inform the joint Sigma (and hence the binary/ordinal posteriors); only its continuous-column baseline OUTPUT is discarded in favour of the lambda-aware per-column BM fit.

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