Generate experimental stochastic completion datasets

Run pigauto's full imputation pipeline and return M stochastic completions of the trait matrix instead of a single point estimate. The conformal-width and Brownian/MC-dropout draws returned here are experimental prediction-diagnostic draws. Do not use these datasets for downstream inference or Rubin pooling. For continuous traits, draws_method = "posterior" instead returns proper Bayesian posterior imputations for the analyses described in "Posterior draws" below: fit the analysis to them with with_imputations() and pool the resulting fits with pool_mi(). The separate analysis-aware backend, multi_impute_analysis(), covers its own documented narrow regime.

Usage

multi_impute(
  traits,
  tree,
  m = 100L,
  draws_method = c("auto", "conformal", "mc_dropout", "posterior"),
  species_col = NULL,
  trait_types = NULL,
  multi_proportion_groups = NULL,
  log_transform = TRUE,
  missing_frac = 0.25,
  covariates = NULL,
  epochs = 2000L,
  verbose = TRUE,
  seed = NULL,
  gnn = FALSE,
  lambda_mode = c("estimate", "fixed_1", "cv", "bayes"),
  posterior_control = list(),
  ...
)

Arguments

traits

data.frame with species as rownames and trait columns. Same input format as impute(). Supported column types are numeric, integer, factor, ordered factor, and logical.

tree

object of class phylo aligned with traits.

m

integer. Number of stochastic completion datasets to generate (default 100). Observed cells are identical across all M datasets; only originally-missing cells vary.

draws_method

character. How stochastic draws are generated for missing cells. One of:

"auto": (default; previously "conformal") Use "posterior" when the data fit its requirements (every trait continuous, one row per species, no covariates, no multi_proportion_groups); otherwise use "conformal" and print a message saying why and that the draws cannot be pooled. The chosen method is stored in result$draws_method. To get the previous default, set draws_method = "conformal".
"conformal": Run the model once, then sample each originally-missing cell from a Normal distribution centred on the point estimate with SD = conformal_score / 1.96. Converting a split-conformal residual quantile to a Normal scale is a heuristic; a nominal held-out conformal diagnostic does not establish that these draws are proper multiple imputations. Falls back to BM-SE-based Normal sampling when conformal scores are unavailable, and to Bernoulli / Categorical draws for discrete traits.
"mc_dropout": Run M stochastic GNN forward passes in training mode (dropout active) on top of stochastic Brownian-motion baseline draws. Brownian draws still contribute between-draw variation when a calibrated GNN gate is zero.
"posterior": Continuous traits only, one row per species, no covariates. Fit a Bayesian multivariate phylogenetic mixed model by MCMC and return m completions drawn from the posterior predictive distribution of the missing cells. No GNN is fitted. See "Posterior draws" below. In simulation (continuous traits, 30% missing completely at random, 200 datasets per cell, n = 100 to 1000, lambda 0.3 to 1, trait correlation 0 or 0.5), a Rubin-pooled downstream PGLS slope and phylogenetic correlation had 95% interval coverage of 0.925 to 0.975 in every cell. At lambda = 1 with trait correlation 0.5, slope coverage was 0.955, 0.940 and 0.965 at n = 100, 300 and 1000 (proper frequentist multiple imputation: 0.955, 0.920, 0.935). At n = 100 with correlated traits the pooled slope was biased by about -0.02 towards zero, with coverage still 0.96. The previous residual prior (posterior_control = list(residual_prior = "iw")) attenuated the relationship between correlated traits at lambda near 1: coverage 0.96, 0.88 and 0.75 at the same three n.

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, two traits, 120 replicates per regime), 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 "posterior".

species_col

character or NULL. If set, marks the column in traits containing species identifiers and enables multiple observations per species. See impute() for details.

trait_types

named character vector overriding auto-detected trait types for specific columns. Required for "proportion" and "zi_count". See impute() and preprocess_traits(). Default NULL (auto-detect).

multi_proportion_groups

named list declaring compositional trait groups (rows summing to 1), e.g. list(diet = c("plant", "invert", "vert")). Forwarded to impute() / preprocess_traits(). Default NULL.

log_transform

logical. Auto-log positive continuous columns (default TRUE).

missing_frac

numeric. Fraction of observed cells held out for validation/test during training (default 0.25). Passed through to impute().

covariates

data.frame or matrix of environmental covariates (fully observed, numeric). Passed through to impute(). Default NULL (no covariates).

epochs

integer. Maximum GNN training epochs (default 2000).

verbose

logical. Print progress (default TRUE).

seed

optional integer. When supplied, makes fitting and imputation draws reproducible; the default NULL uses the current RNG stream.

gnn

logical. Passed through to impute() / fit_pigauto(). When TRUE, the usual GNN correction is trained. When FALSE (default), no GNN is used (baseline-only fit; see fit_pigauto()). With draws_method = "mc_dropout" this degrades gracefully to BM-posterior draws (there is no dropout to run) and a one-time message is printed; draws_method = "conformal" is unaffected. User covariates are ignored under gnn = FALSE (with a warning), as in impute(). The default changed from TRUE to FALSE in 0.11.0.9001; see impute().

lambda_mode

character. Pagel-lambda mode for the BM baseline, forwarded to impute() / fit_pigauto(). "estimate" (default) fits a per-trait Pagel's lambda on each continuous-family (BM-eligible) latent column; discrete traits stay at lambda = 1 unless routed to predict_method = "exact", where they share lambda_block. "fixed_1" preserves the pre-lambda Brownian correlation matrix everywhere; "cv" and "bayes" are alternative per-column estimators. See fit_pigauto() for the full contract, including the predict_method = "exact" / joint_refine_iter > 0 interaction with lambda_block.

posterior_control

list of settings for draws_method = "posterior" (ignored otherwise). Elements, with defaults:

n_chains: 4L. Number of MCMC chains. Chain 1 starts at the per-trait REML Pagel's lambda clamped to [0.02, 0.98] (up to 2,000 tips; above that, or if the REML fit fails, at lambda = 0.9), and the other chains start from values dispersed around it. mi$posterior$start holds the REML lambda and chain 1's lambda before clamping.
n_iter: 5000L. Sweeps per chain after burn-in.
burnin: 1000L. Burn-in sweeps per chain (also used to tune the Metropolis step sizes).
thin: Default: n_iter integer-divided by ceiling(keep_draws / n_chains), so that the chains together keep keep_draws sweeps.
keep_draws: 1000L. Target number of kept posterior predictive draws across chains, used for the per-cell intervals. Must be at least m. The default thin keeps at least this many whenever n_iter is at least ceiling(keep_draws / n_chains). If a user-set n_iter or thin keeps fewer (but at least m), a warning is raised and the intervals use the sweeps that were kept; keeping fewer than m is an error.
param_uncertainty: "full" (default), "none" or "both". "none" is an improper plug-in mode for validation only: the covariance matrices are fixed at their posterior means from a full run, and only the missing cells (and the trait means) are drawn. Its result carries mi_workflow = "pigauto_posterior_plugin_diagnostic", and with_imputations() and pool_mi() refuse it. "both" is also for validation only: it runs the sampler once and returns the proper results exactly as "full" does, plus the plug-in draws from the same run in mi$posterior_improper.
seed: Integer or NULL. Defaults to seed.
auto_extend: TRUE. If the chains fail the convergence rule (any split R-hat at least 1.05 or any bulk ESS at most 400), every chain continues from where it stopped for another n_iter sweeps and the diagnostics are recomputed on all sweeps after burn-in; this repeats until the rule is met or max_extend extensions have been made. 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. After k extensions the kept draws are every thin * (k + 1)-th sweep, spread over the whole run, so their number does not change. A fit that meets the rule first time is not affected. FALSE never extends.
max_extend: 3L. Maximum number of extensions, a whole number from 0 to 10. The default allows at most 4 times n_iter sweeps per chain after burn-in.
residual_prior: "sep" (default) or "iw". The prior on the residual covariance. "sep" (the separation strategy) gives each residual standard deviation a half-Cauchy prior scaled by the trait's observed standard deviation and the residual correlation a uniform (LKJ(1)) prior, so the residual covariance can shrink towards 0 when traits are close to Brownian motion. "iw" is the inverse-Wishart prior used before pigauto 0.11.0.9002; it biases a downstream relationship between correlated traits towards zero when lambda is near 1 (see draws_method). Kept for reproducing earlier results.

...

additional arguments forwarded to fit_pigauto() via impute(). See fit_pigauto() for the full list; the "Safety floor" section below describes the relevant new v0.9.1.9002 argument.

Value

An object of class "pigauto_mi" with components:

datasets: A list of length m. Each element is a data.frame with the same shape and column types as the input traits; observed cells are preserved and missing cells are filled with the corresponding stochastic draw. These datasets are for prediction diagnostics, not downstream inference, except with draws_method = "posterior" (unless posterior_control$param_uncertainty = "none"), whose datasets are posterior imputations: fit the analysis to them with with_imputations() and pool the fits with pool_mi() (see posterior below and "Posterior draws").
m: Number of stochastic completion datasets.
pooled_point: A single data.frame whose missing cells are replaced by the MC-averaged point estimate. Convenient for reporting but does not provide a valid downstream MI analysis.
se: Matrix of per-cell uncertainty summaries combining the baseline SE and the between-draw standard deviation. A descriptive spread for reporting and ranking only — it is not a Rubin's-rules pooled standard error (it contains no within-imputation variance component and no small-sample df correction) and must not be used for downstream pooling.
mi_workflow: "pigauto_diagnostic_mi", recording that these are prediction-diagnostic completions and cannot be passed to with_imputations() or pool_mi().
imputed_mask: Logical matrix; TRUE where a cell was originally missing.
fit: The underlying pigauto_fit object, retained for diagnostics and for calls to predict() on new data.
data: The pigauto_data object.
tree: The input phylogeny.
species_col: Passed-through species-column name or NULL.
posterior: draws_method = "posterior" only. A list with cell_interval (data.frame: row (row of traits), trait, lower, upper, median; 95% posterior predictive interval on the original scale), diagnostics (data.frame: parameter, rhat, ess_bulk, with attributes "converged", "n_extensions" (the number of automatic extensions made) and "sweeps_per_chain" (the sweeps after burn-in per chain that the diagnostics use)), params (list: Sigma_P and Sigma_E as K x K x draws arrays, lambda and mu as draws x K matrices, on the latent scale), converged, control, hyper (priors), start (REML lambda and chain-1 start lambda, before clamping to [0.02, 0.98]), draw_index (kept sweeps used for the m datasets), wall_s and sweeps (all sweeps, burn-in included, summed over chains). For this method fit is NULL, se holds the posterior predictive SD of each imputed cell, the class is c("pigauto_posterior_mi", "pigauto_mi", "list") and mi_workflow is "pigauto_posterior_mi_v1" ("pigauto_posterior_plugin_diagnostic" with posterior_control$param_uncertainty = "none", which with_imputations() and pool_mi() refuse).
posterior_improper: Only with posterior_control$param_uncertainty = "both" (validation only). A list with datasets (m completed data.frames drawn with the covariance matrices fixed at their posterior means, from the same chain run as datasets), cell_interval (same columns as posterior$cell_interval), and the fixed Sigma_P and Sigma_E. Not for downstream inference: it ignores parameter uncertainty.

Details

The conformal and MC-dropout draws do not condition on a declared substantive analysis model. Consequently, stochastic variation alone does not make them proper or congenial multiple imputations. The analysis-aware backend requires the analysis model before generating draws and dispatches only across its documented supported model classes.

draws_method = "conformal" (default): Run the model once; missing cells are sampled from $x_{ij}^{(k)} \sim \mathrm{N}(\hat\mu_{ij},\; q_{j}/1.96)$ where $q_j$ is the trait-level split-conformal residual quantile. Dividing this quantile by 1.96 is a pragmatic Normal-scale construction, not an inference consequence of a nominal held-out conformal diagnostic. For discrete traits (binary, categorical) it uses Bernoulli / categorical draws from the estimated probability vector. For zi_count the gate is Bernoulli from P(nonzero) and the conditional magnitude is drawn on the log1p-z latent using the conformal score (or the magnitude latent SE when scores are unavailable) — not the reported SE of expected count. For multi_proportion groups it draws the K CLR latent columns with their BM latent SEs, projects back to sum-zero CLR space, and decodes to the simplex.

draws_method = "mc_dropout": Run M GNN forward passes in training mode (dropout active) on top of stochastic BM baseline draws. When r_cal = 0, the GNN-dropout term disappears but the BM draw still contributes between-draw variance.

draws_method = "posterior": see the "Posterior draws" section.

Nakagawa & Freckleton (2008, 2011) review the consequences of ignoring missing data in ecological and comparative analyses and argue for multiple imputation as the default.

When to use this

With draws_method = "conformal" or "mc_dropout", this function is useful for comparing stochastic prediction behavior from one tree; those draws are not for downstream inference. With draws_method = "posterior" (continuous traits), it returns proper posterior imputations for the analyses described in "Posterior draws": with_imputations() fits the analysis to them and pool_mi() pools those fits. multi_impute_analysis() is the separate analysis-aware backend for its own documented regime.

multi_impute_trees() provides an experimental posterior-tree sensitivity path, but tree uncertainty is not supported by multi_impute_analysis().

Posterior draws

With draws_method = "posterior", the latent (z-scored, optionally log-transformed) traits follow the multivariate phylogenetic mixed model $$\mathrm{vec}(Y) \sim \mathrm{N}(1\mu^\top,\; \Sigma_P \otimes R + \Sigma_E \otimes I_n),$$ with $R$ the phylogenetic correlation matrix of the tree, full $K \times K$ phylogenetic ($\Sigma_P$) and residual ($\Sigma_E$) covariance matrices, and a flat prior on $\mu$. The implied Pagel's lambda of trait $k$ is $\Sigma_P[k,k] / (\Sigma_P[k,k] + \Sigma_E[k,k])$.

Priors: $\Sigma_P = \mathrm{diag}(\alpha)\,\Sigma_W\, \mathrm{diag}(\alpha)$ with $\Sigma_W \sim \mathrm{IW}(K+1, I_K)$ and $\alpha \sim \mathrm{N}(0, 1000\, I_K)$ (parameter expansion as in MCMCglmm; Hadfield 2010; Gelman 2006); $\Sigma_E \sim \mathrm{IW}(K+1,\, 0.01\,\mathrm{diag}(s^2))$, where $s^2$ are the observed variances of the latent traits; and a flat prior on $\mu$. Here $\mathrm{IW}(\nu, S)$ has density proportional to $|\Sigma|^{-(\nu+K+1)/2}\exp\{-\mathrm{tr}(S\Sigma^{-1})/2\}$. MCMCglmm parameterises the inverse-Wishart by $(V, \nu)$ with scale matrix $\nu V$, so its equivalent of the $\Sigma_W$ prior is $V = I_K/(K+1)$, $\nu = K+1$, not MCMCglmm's default prior. The prior settings used are returned in mi$posterior$hyper.

The sampler draws the phylogenetic effects, the trait means and the missing cells as one block with a sparse Cholesky factor of the Hadfield and Nakagawa (2010) node precision, and updates the covariance matrices by parameter-expanded Gibbs steps plus Metropolis moves with the phylogenetic effects integrated out. Convergence is checked with rank-normalised split R-hat and bulk effective sample size (Vehtari et al. 2021); mi$posterior$diagnostics reports them. When any R-hat is at least 1.05 or any ESS is at most 400, the chains are extended automatically (posterior_control$auto_extend and max_extend), and a warning is raised if the rule still fails after the last extension.

The m datasets are posterior predictive draws spaced evenly across the kept sweeps of all chains. Per-cell 95% intervals in mi$posterior$cell_interval come from all kept sweeps (keep_draws of them by default), never from the m datasets.

These 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 imputations for analyses that are linear in the imputed traits on that scale and whose variables are all among the imputed traits, for example a phylogenetic regression of one imputed trait on others (on the log scale for log-transformed traits). Not covered: covariates from outside the imputed traits, because the imputation model does not contain them; nonlinear terms (such as squares) or interactions among imputed traits; and analysing a log-transformed trait on its raw scale. The GNN is not used, and gnn, epochs, missing_frac, lambda_mode and other fitting arguments are ignored (a message lists any that were supplied). Non-continuous traits, multiple observations per species (or any species_col), covariates, and input with no missing cells in the rows of traits (tree tips absent from traits do not count) are errors.

Safety floor (v0.9.1.9002+)

When fit_pigauto() was called with safety_floor = TRUE (the default since v0.9.1.9002), the 3-way blend r_BM * BM + r_GNN * GNN + r_MEAN * MEAN propagates through every imputation draw automatically via the updated predict.pigauto_fit(). For draws_method = "mc_dropout" the mean term contributes no between-draw variance (it is a deterministic scalar per column); between-draw variance comes from the BM-draw and GNN-dropout terms only. For draws_method = "conformal" the blend centre is the 3-way prediction and conformal scores remain calibrated on the blended residuals.

References

Rubin DB (1987). Multiple Imputation for Nonresponse in Surveys. Wiley.

Nakagawa S, Freckleton RP (2008). "Missing inaction: the dangers of ignoring missing data." Trends in Ecology & Evolution 23(11): 592-596.

Gelman A (2006). "Prior distributions for variance parameters in hierarchical models." Bayesian Analysis 1(3): 515-534.

Hadfield JD (2010). "MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package." Journal of Statistical Software 33(2): 1-22.

Hadfield JD, Nakagawa S (2010). "General quantitative genetic methods for comparative biology: phylogenies, taxonomies and multi-trait models for continuous and categorical characters." Journal of Evolutionary Biology 23(3): 494-508.

Vehtari A, Gelman A, Simpson D, Carpenter B, Buerkner P-C (2021). "Rank-normalization, folding, and localization: an improved R-hat for assessing convergence of MCMC." Bayesian Analysis 16(2): 667-718.

Nakagawa S, Freckleton RP (2011). "Model averaging, missing data and multiple imputation: a case study for behavioural ecology." Behavioral Ecology and Sociobiology 65(1): 103-116.

Examples

# \donttest{
library(pigauto)
data(avonet300, tree300)
tree <- ape::keep.tip(tree300, tree300$tip.label[seq_len(30L)])
df <- avonet300[match(tree$tip.label, avonet300$Species_Key),
                c("Mass", "Wing.Length"), drop = FALSE]
rownames(df) <- tree$tip.label
df$Mass[seq_len(3L)] <- NA_real_
mi <- multi_impute(df, tree, m = 2L, epochs = 5L, verbose = FALSE)
#> draws_method = "posterior" uses the phylogenetic mixed model only (no GNN is fitted); ignoring: epochs.
print(mi)
#> pigauto posterior multiple imputation (phylogenetic mixed model)
#>   M        : 2 completed datasets (posterior predictive)
#>   Species  : 30
#>   Traits   : 2 -- Mass, Wing.Length
#>   Cells    : 3 imputed / 60 total (5.0%)
#>   MCMC     : 4 chains x (1000 burn-in + 5000 sweeps, thin 20)
#>   lambda   : Mass = 0.99, Wing.Length = 0.99 (posterior means)
#>   Converged: yes (max R-hat 1.007, needs < 1.05; min bulk ESS 1106, needs > 400)
#> 
#>   Access datasets:          mi$datasets[[i]]
#>   Per-cell 95% intervals:   mi$posterior$cell_interval
#>   Downstream inference:     with_imputations(mi, f) then pool_mi()
lapply(mi$datasets, head)
#> [[1]]
#>                              Mass Wing.Length
#> Nothoprocta_pentlandii   117.2626       148.2
#> Eudromia_formosa         213.7835       223.9
#> Rhea_americana          2776.2849       604.5
#> Ptilopachus_petrosus     193.0000       117.2
#> Bambusicola_fytchii      313.7000       142.8
#> Francolinus_psilolaemus  478.3000       167.3
#> 
#> [[2]]
#>                              Mass Wing.Length
#> Nothoprocta_pentlandii   131.7740       148.2
#> Eudromia_formosa         666.2496       223.9
#> Rhea_americana          9505.8317       604.5
#> Ptilopachus_petrosus     193.0000       117.2
#> Bambusicola_fytchii      313.7000       142.8
#> Francolinus_psilolaemus  478.3000       167.3
#> 
# }