Custom Priors and Mechanistic Constraints • because

Introduction

One of the strengths of Bayesian modeling is the ability to incorporate prior knowledge into your analysis. In the because package, we generally use “uninformative” or weakly informative priors by default to let the data speak for itself. However, there are many cases where you might want to inject specific information:

Mechanistic Knowledge: You know a parameter (like a growth rate) must be positive or lies within a biologically feasible range.
Previous Studies: You have estimates from a previous meta-analysis or experiment that you want to use as a starting point.
Measurement Error: You know the measurement error of an instrument and want to fix or constrain the residual variance.

This vignette demonstrates how to use the priors argument in because() to override default priors and inject your own mechanistic constraints.

Example 1: Comparing Default and Custom Priors

For this example, let’s simulate Absolute Growth Rate in grams per day (g/day) for hypothetical juvenile birds responding to ambient temperature.

library(because)

set.seed(42)
N <- 30
# Ambient temperature in Celsius (Mean 20C)
Temp_Raw <- rnorm(N, mean = 20, sd = 5)

# Center it so '0' represents the mean (20C).
# This ensures the intercept represents growth at average conditions,
# rather than at 0°C (where the bird would be frozen solid!).
Temp_Centered <- Temp_Raw - 20

# True relationship:
# Baseline growth at mean temp is 10 g/day.
# Every degree of warming above mean adds 0.5 g/day.
Growth_g_day <- 0.5 * Temp_Centered + 10 + rnorm(N, sd = 2)

df <- data.frame(Temp_Centered, Growth_g_day)

Default Model

Let’s fit a standard model first. By default, because assigns weakly informative Gaussian priors (dnorm(0, 0.01)) to intercepts (alpha) and slopes (beta). This corresponds to a Standard Deviation of 10. These anchors are broad enough to be uninformative for scientific data while preventing the numerical instability and convergence failure ( $Rhat > 1.1$ ) often seen with extremely “flat” priors in complex hierarchical models.

For Gaussian response variables, the residual standard deviation gets a half-uniform prior: sigma_e ~ dunif(0, 100), from which the precision is derived as tau_e <- 1 / (sigma_e * sigma_e). This is scale-invariant and works safely whether data are standardised or on their original (raw) scale.

fit_default <- because(
    equations = list(Growth_g_day ~ Temp_Centered),
    data = df
)

summary(fit_default)

Custom Prior Model

Now, suppose we have strong prior knowledge from scaling theory or previous literature: 1. Intercept (Alpha): At mean temperatures (20°C), metabolic constraints suggest growth should be around 10 g/day. 2. Slope (Beta): We expect a positive effect of temperature, likely around 0.5 g/day/°C.

The parameter names usually follow this convention:

Intercepts: alpha_{Response} (e.g., alpha_Growth_g_day)

Slopes: beta_{Response}_{Predictor} (e.g., beta_Growth_g_day_Temp_Centered)

Residual SD (Gaussian): sigma_e_{Response} (e.g., sigma_e_Growth_g_day) — this is what you override in priors for Gaussian responses.

Residual precision (derived): tau_e_{Response} — set automatically from sigma_e; do not set directly unless you need a legacy dgamma override.

# Define our custom priors
my_priors <- list(
    # Strong prior on intercept: Mean 10, Precision 100 (SD = 0.1)
    alpha_Growth_g_day = "dnorm(10, 100)",

    # Informative prior on slope: Mean 0.5, Precision 400 (SD = 0.05)
    beta_Growth_g_day_Temp_Centered = "dnorm(0.5, 400)"
)

fit_custom <- because(
    equations = list(Growth_g_day ~ Temp_Centered),
    data = df,
    priors = my_priors
)

summary(fit_custom)

Notice how the credible intervals for the custom model will be tighter and centered closer to our priors, especially if the data were sparse (small N).

Visualizing the Impact

We can plot the posterior estimates together to see the “shrinkage” effect of our informative priors.

# We can use the helper function plot_posterior() to compare models.
# By passing a list of models, they are overlaid on the same plot.
# The 'parameter' argument uses partial matching (regex), so "Growth_g_day"
# will match both the intercept (alphaGrowth_g_day) and the slope (beta_Growth_g_day...).

plot_posterior(
    model = list("Default" = fit_default, "Custom" = fit_custom),
    parameter = "Growth_g_day"
)

Example 2: Mechanistic Constraints

In many biological contexts, a negative effect is not just unlikely, it is physically impossible.

For example, consider the relationship between Body Mass and Metabolic Rate (Kleiber’s Law). A larger animal simply cannot consume less energy than a smaller one (a negative slope), all else equal. The physics of life requires costs to scale positively with mass.

However, if our sample size is small and measurement error is huge, we might accidentally estimate a negative slope. We can prevent this by enforcing a positive prior.

# Simulate data following Kleiber's Law: MR = a * Mass^0.75
# Taking logs: log(MR) = log(a) + 0.75 * log(Mass)
set.seed(42)
Mass <- runif(30, 10, 100)
# True scaling exponent is 0.75
# We add enough noise that the estimated slope might be negative by chance
Log_Mass <- log(Mass)
Log_MR <- 1.5 + 0.75 * Log_Mass + rnorm(30, sd = 2.5)

df_kleiber <- data.frame(Log_Mass, Log_MR)

# Prior: The scaling exponent (slope) must be positive
# Physics dictates that metabolic cost increases with mass.
positive_prior <- list(
    beta_Log_MR_Log_Mass = "dnorm(0, 1) T(0, )"
)

# 1. Fit Default Model (Unconstrained) - Log-Log regression
fit_default_kleiber <- because(
    equations = list(Log_MR ~ Log_Mass),
    data = df_kleiber
)

# 2. Fit Mechanistic Model (Constrained)
fit_mech_kleiber <- because(
    equations = list(Log_MR ~ Log_Mass),
    data = df_kleiber,
    priors = positive_prior
)

# Compare Estimates
summary(fit_default_kleiber)
summary(fit_mech_kleiber)

# Visualize: Unconstrained vs. Truncated
plot_posterior(
    list(Unconstrained = fit_default_kleiber, Mechanistic = fit_mech_kleiber),
    parameter = "beta_Log_MR_Log_Mass",
    density_args = list(Mechanistic = list(from = 0))
)

Example 3: Informing Residual Variance Components

We usually estimate the residual standard deviation (sigma_e) from the data. However, you might have prior knowledge about the expected noise level — for example, a laboratory instrument with known precision, or a measurement protocol with a documented coefficient of variation.

[!NOTE] Distinction from variability argument: The variability argument in because() is designed for Measurement Error in predictors (Error-in-Variables) or when you have repeated measures per individual.

Custom priors on sigma_e_*, shown below, are for providing information about the Residual Variance of the response variable (which includes both Process Error and unrecognised Observation Error).

[!TIP] Why sigma_e and not tau_e? Since version >= 0.5, because priors Gaussian residual error on the standard deviation scale (sigma_e ~ dunif(0, 100)). This is scale-invariant: it works correctly whether your response variable is in grams or tonnes. The precision tau_e is derived deterministically as tau_e <- 1 / (sigma_e * sigma_e) and cannot be given a prior directly — use sigma_e_* as the override key.

For example, if field protocol experience tells you that measurement noise in growth measurements is typically 1–3 g/day, you can constrain the prior on sigma_e accordingly:

# Simulate small, noisy dataset
set.seed(42)
N_small <- 15
Temp_Small <- rnorm(N_small, 20, 5)
# Actual residual SD = 3 (quite noisy)
Growth_Small <- 0.5 * Temp_Small + rnorm(N_small, sd = 3)
df_small <- data.frame(Growth_g_day = Growth_Small, Temp_Centered = Temp_Small - 20)

# 1. Default Model
# dunif(0, 100) on sigma_e: completely open to any residual SD up to 100.
# With only 15 observations the posterior will be wide.
fit_default_var <- because(
    equations = list(Growth_g_day ~ Temp_Centered),
    data = df_small,
    quiet = TRUE
)

# 2. Constrained Model
# Prior knowledge: residual SD should be small, around 1 g/day.
# dunif(0, 2) restricts sigma_e to [0, 2], strongly informative.
variance_prior <- list(
    sigma_e_Growth_g_day = "dunif(0, 2)"
)

fit_constrained_var <- because(
    equations = list(Growth_g_day ~ Temp_Centered),
    data = df_small,
    priors = variance_prior,
    quiet = TRUE
)

# Visualize: The constrained posterior for sigma will be shifted left (towards 1) and sharper
plot_posterior(
    list(Default = fit_default_var, Constrained = fit_constrained_var),
    parameter = "sigma"
)

[!NOTE] Legacy tau_e_* overrides: If you previously used priors = list(tau_e_Response = "dgamma(...)"), this still works — because checks for a tau_e_* entry in priors first and uses it verbatim, bypassing the sigma parameterisation entirely. This ensures backward compatibility.

Default Priors Reference

It is important to know what you are overriding. because aims to use weakly informative defaults that provide minimal regularization while allowing the data to dominate the posterior.

Parameter Type	Parameter Name	Default Prior	Description
Intercepts	`alpha_*`	`dnorm(0, 0.01)`	Weakly Informative SD = 10
Coefficients	`beta_*`	`dnorm(0, 0.01)`	Weakly Informative SD = 10
Residual SD (Gaussian)	`sigma_e_*`	`dunif(0, 100)`	Half-uniform on σ; scale-invariant, robust to raw or standardised data. Override with `sigma_e_*` in `priors`.
Residual Precision (Gaussian)	`tau_e_*`	derived	Set deterministically as `1/(sigma_e * sigma_e)`. Not assigned a prior directly.
Precision (other families)	`tau_*`	`dgamma(1, 1)`	Gamma(1,1) for non-Gaussian precision components (random effects).
Phylo Signal	`lambda_*`	`dunif(0, 1)`	Uniform on [0,1]
Zero-Inflation	`psi_*`	`dbeta(1, 1)`	Beta(1,1) (Uniform on [0,1])
Ordinal Cutpoints	`cutpoint_*`	`dnorm(0, 0.01)`	SD = 10
NegBinomial Size	`r_*`	`dgamma(0.01, 0.01)`	Wide Gamma for dispersion parameter

Precision vs Variance: JAGS uses precision $\tau = 1/\sigma^2$ . A prior of dnorm(0, 0.01) means a normal distribution with mean 0 and precision $0.01$ , which corresponds to a variance of $100$ (SD = 10). This is uninformative for most standardized biological effects but remains numerically stable during adaptation.

Parameter Names Reference

To use custom priors, you need to know the exact internal name of the parameter in the JAGS code.

Common patterns:

alpha_{Response}: Intercept
beta_{Response}_{Predictor}: Regression coefficient
sigma_e_{Response}: Residual standard deviation for Gaussian responses — use this key in priors to override the default dunif(0, 100)
tau_e_{Response}: Residual precision (derived from sigma_e; not directly set via priors — use sigma_e_* instead)
lambda_{Response}: Phylogenetic signal (0–1)
psi_{Response}: Zero-inflation probability

If you are unsure, run a quick model with n.adapt=0, n.iter=0 (just compilation) and check the generated model code:

fit_check <- because(
    equations = list(Growth_g_day ~ Temp_Centered),
    data = df,
    n.iter = 0, quiet = TRUE
)

# Print the JAGS model
fit_check$model

```