Prior-sensitivity sweep for the phylogenetic correlation penalty

Refits a penalized (MAP) phylogenetic model across a range of cor_sd values so you can see whether a weakly identified coupling is data-informed or prior-shaped. There is no universal cor_sd: a coupling that is stable across the sweep is data-informed, while one that tracks cor_sd is prior-shaped. This is the sweep the penalized/MAP workflow asks you to run; it is most informative for coupled location-scale or bivariate phylogenetic models that actually estimate a phylogenetic correlation.

Usage

drm_phylo_penalty_sweep(
  formula,
  data,
  family = gaussian(),
  cor_sd = c(0.25, 0.5, 1),
  sd_u = 1,
  sd_alpha = 0.05,
  control = drm_control(),
  ...
)

Arguments

formula, data, family, control

Passed to drmTMB().

cor_sd

Numeric vector of positive correlation-penalty SDs to sweep.

sd_u, sd_alpha

Penalty SD-prior parameters; see drm_phylo_penalty().

...

Further arguments passed to drmTMB().

Value

A list with $summary – a data frame with one row per cor_sd giving convergence, pdHess, logLik, and any fit error – and $fits – the fitted objects, named by cor_sd, for extracting corpars(), coef(), and other couplings.

Examples

# \donttest{
if (requireNamespace("ape", quietly = TRUE)) {
  set.seed(20260601)
  n_tip <- 10
  tree <- ape::rcoal(n_tip)
  tree$tip.label <- paste0("sp_", seq_len(n_tip))
  A <- ape::vcv(tree, corr = TRUE)
  u <- as.vector(t(chol(A)) %*% rnorm(n_tip)) * 0.6
  species <- factor(rep(tree$tip.label, each = 2), levels = tree$tip.label)
  x <- rnorm(length(species))
  dat <- data.frame(
    y = 0.3 + 0.5 * x + u[rep(seq_len(n_tip), each = 2)] +
      rnorm(length(species), sd = 0.5),
    x = x,
    species = species
  )

  # A coupled location-scale phylo model has two phylogenetic SDs, so the
  # cor_sd sweep is informative (see drm_phylo_penalty()).
  out <- drm_phylo_penalty_sweep(
    bf(
      y ~ x + phylo(1 | species, tree = tree),
      sigma ~ phylo(1 | species, tree = tree)
    ),
    data = dat,
    family = gaussian(),
    cor_sd = c(0.5, 1)
  )
  out$summary
}
#>   cor_sd convergence pdHess    logLik error
#> 1    0.5           0   TRUE -12.67675  <NA>
#> 2    1.0           0   TRUE -12.62564  <NA>
# }