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,
  discrete_lambda = c("estimate", "fixed_1"),
  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 a joint baseline. "hard" (default) thresholds binary proportions at 0.5 and uses the most frequent categorical class. "soft" preserves species-level proportions and combines truncated-Gaussian liability expectations, so intermediate class frequencies retain fractional evidence. Only relevant to multi-observation binary/categorical joint baselines.
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) 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; their lambda is governed by discrete_lambda (below), and is 1 when discrete_lambda = "fixed_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". When lambda_fixed = NULL, requesting either unsupported mode issues a warning and fits these columns at lambda = 1. An explicit lambda_fixed takes precedence.
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. To replay a previous fit's per-trait lambda without re-estimating it, supply that fit's $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.
discrete_lambda: character, "estimate" (default) or "fixed_1". With "estimate", the binary, ordinal and zero-inflated-count gate liability columns of the threshold-joint baseline, and the one-vs-rest column of each categorical class, get their own Pagel's lambda, so a weak-signal discrete trait shrinks toward its prevalence instead of copying close relatives. "fixed_1" restores the previous behaviour (binary and ordinal at lambda = 1, or at the shared joint lambda on the "exact" route; categorical at 1). In the reported simulated Brownian-motion and OU comparison (n = 100 to 1000, 30 percent missing; see the changelog), estimating discrete lambda alone added about 0.01 to 0.04 discrete accuracy at Pagel's lambda of 0.3 or below; together with safety_floor = FALSE and phylo_signal_gate = FALSE (also the defaults) the gain was 0.04 to 0.12, and binary and categorical accuracy matched or beat BACE's. Costs: with lambda near 1 and n >= 1000, accuracy -0.002 to -0.005 and Brier score +0.010 to +0.016; on AVONET300's categorical traits, Brier +0.014 to +0.020 with accuracy unchanged or better. Ordinal accuracy still trails BACE's at lambda = 0.3. The categorical one-vs-rest lambdas are not reported in $lambda_per_trait. Recorded in the returned list and in a fit's model_config; binary and ordinal lambdas are replayed by lambda_fixed.
em_iterations: integer. Maximum passes through each liability-prior refinement loop for binary, ordinal, and one-vs-rest categorical baselines. Default 0L disables iterative refinement and uses a single-pass baseline. 1L also permits only the initial pass; 2L permits at most one feedback pass. After each pass, estimated liability variances inform the next prior. Refinement can stop early on convergence or a later fit failure. A categorical pass may require K one-vs-rest fits; predict_method = "auto" may fit two route candidates. Continuous-only joint MVN fits are unaffected.
em_tol: numeric. Relative-Frobenius convergence tolerance for liability-prior refinement (default 1e-3). The binary/ordinal threshold-joint loop compares the variance diagonal by default, or the full covariance matrix when em_offdiag = TRUE. Categorical one-vs-rest refinement compares the vector of class variances.
em_offdiag: logical. When TRUE and em_iterations >= 2L, binary/ordinal liability priors use conditional-MVN means and SDs based on the other traits and the full covariance matrix from the preceding pass. With the default FALSE, refinement updates per-trait prior SDs only. This option does not change categorical one-vs-rest refinement.
joint_solver: character. Solver for the joint MVN, threshold-joint, and categorical one-vs-rest baselines. "inhouse" (default) uses pigauto's joint solver. "rphylopars" delegates to Rphylopars::phylopars(), using its lambda model under lambda_mode = "estimate", with automatic fallback to "inhouse" on failure or implausible output. This does not affect routes outside these joint fits. "cv"/"bayes" lambda modes instead use per-column numeric fits. See fit_baseline for dispatch and fallbacks.
predict_method: character. Prediction route for the in-house joint solver. "auto" (default) 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 for this outer exact-versus-per-column choice; half B is reserved for gate calibration and conformal scoring. This separation applies only to that outer choice. For ordinal traits, each candidate fit's inner threshold-joint-versus-BM selector has already used the full validation split, including half B. Below 38 rows each half would keep fewer than the 19 cells needed for a 95 percent split-conformal interval, 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 requests "exact". With splits = NULL, no validation comparison is possible, so every trait requests exact prediction; unsupported or unusable exact fits may still fall back to per-column prediction. 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 so the mean and covariance models use the same phylogenetic correlation. 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. Maximum passes of optional cross-trait refinement of joint-baseline cell imputations using the estimated covariance matrix. 0L disables this refinement. The covariance step must shrink each iteration; otherwise the last accepted iterate is restored and $diverged is set. This has no effect on traits predicted by the "exact" route, which returns the analytic conditional directly. Under predict_method = "auto", it can therefore affect only traits routed to "per_column". Evaluate this opt-in setting on held-out cells from the intended data.
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) with route-specific uncertainty values. Brownian-motion conditional-MVN routes return conditional standard deviations; the zero-inflated-count constant magnitude fallback returns its observed-sample SD. Decoded threshold, one-vs-rest categorical, and label-propagation routes may store zero placeholders, which do not indicate certainty. A zero can also be genuine conditional uncertainty for a numeric observation. These baseline values are distinct from the type-specific uncertainties in pred$se.
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" (continuous-family joint MVN), "threshold_joint" (binary, ordinal, and continuous-family liability joint), "ovr_categorical" (K one-vs-rest categorical fits), "per_column_bm" (per-column Brownian-motion conditional MVN, including lambda-aware fits and an ordinal candidate selected by fit-specific validation MSE against threshold-joint), "multi_proportion_bm" (component-wise BM on CLR columns), "label_propagation" (phylogenetic label propagation for eligible discrete columns), "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 with one entry per latent column, named by colnames(data$X_scaled). Entries report values available for each column; exact joint prediction also uses the shared lambda_block; categorical one-vs-rest fit lambdas are not reported in this vector. For covariate-aware columns where lambda could not be estimated, the effective fallback is 1 and the latent column names are recorded in the vector's lambda_unestimated attribute.
lambda_block: Numeric scalar, the shared lambda used internally by the joint fit for covariance estimation and cross-trait conditional calculations; 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 meet its route requirements. The default in-house joint solver does not require Rphylopars; that package is an optional alternative solver. 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 baseline routes: data$covariates are used by the per-column BM path (a separate Brownian-motion fit for each eligible numeric trait). Joint MVN paths fit eligible numeric traits together; threshold-joint paths add discrete liabilities to a joint fit. Joint MVN and threshold-joint fits both omit the covariate design and emit a warning. One-versus-rest categorical fits omit covariates but do not trigger the separate baseline warning themselves. In a mixed fit, another selected joint route can still trigger it. Each category gets its own class-versus-rest fit. In a mixed fit, per-column numeric BM fallbacks can still use covariates for the columns they handle. The warning applies to the selected joint baseline only. GNN covariate use is separate and occurs only when a GNN correction is enabled in fit_pigauto(). These inputs do not guarantee improved predictions.

Per-type lambda dispatch: lambda_mode controls the continuous-family columns (continuous, count, proportion, and zi_count magnitude). Under "estimate", the selected joint or per-column numeric route estimates each eligible column's lambda; "cv" and "bayes" use the per-column route and are unsupported by the covariate-aware per-column fallback, which uses lambda 1 and warns when no explicit lambda_fixed was supplied. Ordinal, binary, and zi_count gate liabilities are handled by the discrete-liability rules instead of this numeric-column setting. discrete_lambda controls their route-specific behavior: estimated values may be used in threshold-joint fits, while exact joint predictions share the fit's lambda_block. Label-propagation and categorical one-vs-rest routes have additional reporting limits. Covariance estimates alone do not imply that returned predictions used cross-trait conditioning; see predict_method below.

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