Bayesian hierarchical preservation model — morie_taphonomy

A conjugate Gaussian-linear Bayesian model for a continuous preservation outcome, with informative Normal priors on the coefficients – so domain knowledge (e.g. quicklime's desiccant effect) enters as a prior and the data updates it to a posterior. The coefficient posterior is closed-form (no MCMC): with prior $\beta \sim N(m_0, \mathrm{diag}(s_0^2))$ and noise variance $\sigma^2$ estimated from the OLS residuals (empirical Bayes), $$\Sigma = (X^\top X / \sigma^2 + \Lambda_0)^{-1}, \quad \mu = \Sigma (X^\top y / \sigma^2 + \Lambda_0 m_0).$$

Usage

morie_taphonomy_bhm(
  data,
  outcome = "preservation_score",
  covariates = NULL,
  group = NULL,
  priors = NULL,
  prior_sd_default = 10,
  backend = c("conjugate", "cmdstanr", "brms", "rstanarm"),
  chains = 4L,
  iter = 1000L,
  seed = 42L
)

Arguments

data: A non-empty data.frame.
outcome: Continuous outcome column (default "preservation_score").
covariates: Character vector of predictors. Defaults to the schema covariates + measurements present in data.
group: Optional column giving a grouping factor (e.g. burial context) for partial-pooled random intercepts.
priors: Optional named list mapping a coefficient name to list(mean=, sd=) – its informative Normal prior. Unlisted coefficients get a diffuse prior (sd = prior_sd_default). Use e.g. list(lime_treatment = list(mean = 0.3, sd = 0.1)).
prior_sd_default: Diffuse prior sd for unlisted coefficients (default 10).
backend: One of "conjugate" (default; closed-form posterior, no deps), "cmdstanr" (full-Bayes HMC/NUTS via cmdstanr + a built CmdStan), "brms" (brms formula sampler), or "rstanarm" (rstanarm stan_glmer) – the same informative-prior hierarchical model, sampled either way. Any HMC backend additionally returns the fitted stanfit object.
chains, iter, seed: HMC settings for backend = "cmdstanr" (chains, warmup = sampling iterations per chain, RNG seed).

Value

A named list (all estimates double): coefficients (data.frame: term, post_mean, post_sd, ci_lower, ci_upper, prob_positive), sigma, group_effects (NULL unless group), fitted (posterior-predictive mean per row), n, and a plain-language interpretation.

Details

When group is supplied, a second level is added: group intercepts are partially pooled toward the grand mean by empirical-Bayes (normal-normal) shrinkage $\lambda_j = \tau^2 / (\tau^2 + \sigma^2 / n_j)$, giving a genuine two-level hierarchical model.

For full hierarchical inference by HMC/NUTS, fit rstanarm (stan_glmer) or brms instead; this function is the dependency-free conjugate core.

References

Gelman A, et al. (2013). Bayesian Data Analysis (3rd ed.), Ch. 5 (hierarchical models) & Ch. 14 (conjugate regression). CRC.

Examples

# \donttest{
df <- data.frame(preservation_score = rnorm(20),
                 lime_treatment = rbinom(20, 1, 0.5))
morie_taphonomy_bhm(df, covariates = "lime_treatment",
  priors = list(lime_treatment = list(mean = 0.3, sd = 0.1)))$coefficients
#>             term post_mean    post_sd    ci_lower  ci_upper prob_positive
#> 1    (Intercept) 0.1395269 0.21790972 -0.28756831 0.5666221     0.7390102
#> 2 lime_treatment 0.2910739 0.09763312  0.09971652 0.4824313     0.9985649
# }