Skip to contents

Fit a global single-regime multivariate GLS model to an empirical dataset and convert the fitted object into a reusable calibration template for null and shift-recovery simulation studies.

Usage

createSimulationTemplate(
  baseline_tree,
  trait_data,
  formula = "trait_data ~ 1",
  response_columns = NULL,
  predictor_columns = NULL,
  ...
)

Arguments

baseline_tree

A rooted phylogenetic tree of class phylo (or coercible via ape::as.phylo()) whose unique, non-empty tip labels match the row names of trait_data.

trait_data

A numeric matrix or data frame containing the empirical data used to parameterize the simulations. Row names must be unique, non-empty, and match baseline_tree$tip.label. Column names may be omitted for matrix inputs; when present, they must be unique and non-empty.

formula

Formula specification passed to mvMORPH::mvgls(). May be a single character string or a formula object. Indexed formulas that reference trait_data directly are supported (for example, "trait_data ~ 1" or "trait_data[, 1:12] ~ trait_data[, 13]"), as are named-column formulas such as cbind(y1, y2) ~ size + grp.

response_columns

Optional column specification identifying the multivariate response. May be integer positions, character column names, or a mix resolvable against trait_data. Because R coerces mixed atomic vectors to character, integer-like strings are treated as positions when they are not exact column names. Numeric positions must be whole numbers. For intercept-only workflows this defaults to all columns. For any workflow where the response is a subset of trait_data, including formulas such as "trait_data[, 1:12] ~ 1", this should be supplied explicitly.

predictor_columns

Optional column specification identifying predictor columns used by the global calibration model. For named-column formulas these are usually inferred from the formula itself. For non-intercept workflows where predictor_columns is omitted, the complement of response_columns is used unless the formula identifies raw predictor columns directly.

...

Additional arguments passed to mvMORPH::mvgls() for the global empirical fit (for example method = "LL" or error = TRUE).

Value

A list of class bifrost_simulation_template containing the aligned inputs, fitted global model, empirical mean structure, response/predictor column metadata, stored fit settings, the empirical residual covariance and degrees of freedom, and covariance summaries used for simulation.

Details

The returned template stores:

  • the aligned empirical tree and calibration data,

  • the fitted global mvgls model,

  • evaluated copies of key fit settings (fit_method, fit_error) that can be reused safely by downstream simulation-study wrappers,

  • the fitted response mean structure,

  • the full empirical residual covariance and residual degrees of freedom used by the empirical Wishart generator, and

  • summaries of the diagonal and off-diagonal elements (variance_mean, variance_sd, covariance_mean, covariance_sd) retained by the manuscript-reproduction generator.

The template is a calibration object, not a full simulation-study design. Its formula defines the one global mean model used to estimate fitted values and residual covariance structure. Downstream simulation studies then regenerate the response block around those fitted means and evaluate intercept-only shift-search behavior on the simulated responses. In other words, a richer calibration model can still feed the response-only trait_data ~ 1 simulation workflow. Predictor columns are retained in the template for calibration provenance, but simulated replicates returned by the study helpers contain the regenerated response block only.

Internally, all supported formulas are normalized into a single simulation formula specification. Indexed trait_data formulas are rewritten onto named columns, formula objects are accepted alongside character strings, predictor schemas record numeric/logical/factor/ordered predictors, and unsupported forms such as transformed responses, . shorthand, and character predictors are rejected early.

Examples

if (FALSE) { # \dontrun{
set.seed(1)
tr <- ape::rtree(20)
X <- matrix(rnorm(20 * 3), ncol = 3)
rownames(X) <- tr$tip.label

tmpl <- createSimulationTemplate(
  baseline_tree = tr,
  trait_data = X,
  formula = "trait_data ~ 1",
  method = "LL"
)

tmpl

size <- rnorm(20)
X_with_size <- cbind(X, size = size)

covariate_tmpl <- createSimulationTemplate(
  baseline_tree = tr,
  trait_data = X_with_size,
  formula = "trait_data[, 1:3] ~ trait_data[, 4]",
  response_columns = 1:3,
  predictor_columns = 4,
  method = "LL"
)

covariate_tmpl$search_formula
} # }