Fit a Ratio, Rate, or Proportion Model

Fit a Bayesian hierarchical model for ratios, rates, or proportions. The model jointly estimates the numerator and denominator processes with optional shared latent structure.

tratio() is tulpaRatio's front door, the ratio-model member of the tulpa* family of fitting verbs (tulpa::tulpa() for the engine, tulpaObs::tobs() for observation processes).

Core principle: Ratios are derived quantities, computed post hoc from the posterior of the joint model. The ratio is never modelled directly.

Usage

tratio(
  formula,
  data,
  family = ratiod_negbin_negbin(),
  formula_num = NULL,
  formula_denom = NULL,
  shared = NULL,
  spatial = NULL,
  temporal = NULL,
  spatiotemporal = NULL,
  zi = NULL,
  latent = NULL,
  priors = NULL,
  mode = c("auto", "exact", "structured", "optimized", "hmc", "ess", "pg", "gibbs",
    "sghmc", "sgld", "laplace", "vi"),
  control = list()
)

Arguments

formula

Model formula specifying the response and predictors. Two syntax options: Combined (recommended) num | denom ~ predictors + (1|group) where both processes share the same predictors and random effects; Separate num ~ predictors with the formula_denom argument.

data

Data frame containing all variables.

family

A tulpaRatio family object specifying the distributions: ratiod_negbin_negbin() for both count processes (default), ratiod_binomial() for successes/trials, ratiod_poisson_gamma() for count/continuous effort (CPUE).

formula_num

Optional one-sided formula for additional numerator predictors: ~ extra_terms. Added to the main formula.

formula_denom

Optional formula for denominator. Required if main formula has single response. Can be denom ~ predictors.

shared

Formula specifying shared random effects structure: NULL (default) infers from matching random effects in both processes; ~ (1 | group) for explicit shared structure; ~ 0 for independence assumption (triggers warning).

spatial

Optional spatial structure specification. See spatial_car(), spatial_bym2().

temporal

Optional temporal structure specification. See temporal_rw1(), temporal_rw2(), temporal_ar1().

spatiotemporal

Optional spatiotemporal interaction specification. See spatiotemporal().

zi

Optional zero-inflation specification. See zi_poisson(), zi_negbin(), hurdle_poisson(), hurdle_negbin().

latent

Optional latent factor specification for unmeasured confounders. See latent_factor().

priors

Prior specification. See ratiod_priors().

mode

Inference mode or backend. Accepts either:

Tier names (epistemic guarantees):

"auto" (default): Selects between Exact and Structured based on model. Never silently chooses Optimized.
"exact" (Tier 1): Asymptotically correct posterior inference. Credible intervals are interpretable as posterior uncertainty.
"structured" (Tier 2): Accurate inference conditional on structural assumptions (latent Gaussian, conditional independence).
"optimized" (Tier 3): No general correctness guarantee. Point estimates usually good, uncertainty often underestimated. Requires explicit opt-in.

Backend names (specific implementation):

"hmc" (Tier 1): Native HMC/NUTS sampler
"ess" (Tier 1): Elliptical Slice Sampling
"pg" (Tier 1): Pólya-Gamma Gibbs (binomial only)
"sghmc" (Tier 1): Stochastic Gradient HMC (large datasets)
"sgld" (Tier 1): Stochastic Gradient Langevin Dynamics (large datasets)
"laplace" (Tier 2): Laplace approximation
"vi" (Tier 3): Variational Inference

Structures each backend fits. "hmc" fits every structure. The rest carry a subset, and a model naming a structure its backend cannot fit is refused with that structure named:

"ess", "sghmc", "sgld", "vi": spatial, temporal, zi
"gibbs": spatial, and temporal as temporal_tvc()
"pg", "laplace": spatial, or temporal as temporal_multiscale(), one of the two per model

"auto" selects a backend that fits the whole model, so it reaches "laplace" or "gibbs" only for models those carry.

See inference_mode_info() for details on the tier system.

control

Named list of perf / numerical / tuning knobs. Statistical arguments stay in the signature above; everything that tunes how the fit runs lives here. Unknown names are an error, not a silent no-op.

Sampling budget (all backends):

chains: number of MCMC chains (default 4).
iter: total iterations per chain (default 2000).
warmup: warmup iterations per chain (default floor(iter / 2)).
thin: thinning interval (default 1).
cores: cores for parallel chains (default getOption("mc.cores", chains)).
seed: random seed for reproducibility.
verbose: print progress messages and sampler adaptation traces during fitting (default FALSE).

NUTS / HMC (mode = "hmc"):

adapt_delta: target average acceptance probability for step-size adaptation. Higher values give smaller steps, reducing divergences but slowing sampling. NULL (default) selects on model complexity: 0.80 simple, 0.85 ICAR, 0.90 BYM2 and correlated random slopes. Range 0.80–0.99.
max_treedepth: maximum NUTS tree depth. NULL (default) uses 10. Raise if you see many max-treedepth warnings; costs computation per iteration.
metric: mass matrix. "auto" (default) picks on model complexity; "dense" full matrix with Ledoit-Wolf shrinkage (correlated posteriors: random slopes, BYM2, HSGP, TVC); "diag" diagonal (faster per step, may need deeper trees); "block_diag" per parameter group. Auto falls back to diagonal when p > 2000.
walnuts: run each NUTS trajectory as a WALNUTS transition (Bou-Rabee et al. 2025, tulpa's <tulpa/walnuts.h>), which integrates each macro step at the coarsest leapfrog subdivision that keeps the joint log density within tolerance and keeps it only if the reversed step picks the same one. The kernel is exact, and its step follows the local curvature, which is what a funnel-shaped posterior (a hierarchical scale near zero) needs. FALSE (default) runs NUTS.
gradient_mode: "auto" (default) selects the fastest available; "H" hand-coded analytical (fastest); "A_r" arena reverse-mode autodiff (fast, O(N)); "A" forward-mode (O(p*N), thread-safe); "A_t" tape-based reverse-mode (slow, legacy); "N" numerical finite differences (slowest, always works).
re_param: random-effect parameterization. "noncentered" (default) stores z ~ N(0, 1) and computes re = sigma * z (or diag(sigma) * L * z when correlated), suiting weakly-informed effects or small groups; "centered" stores re directly under re ~ N(0, sigma^2), suiting strongly-informed effects or large groups.
L: fixed leapfrog steps. Default 0 leaves NUTS to choose.

Variational inference (mode = "vi"):

vi_variant: "auto" (default) selects on parameter count; "meanfield" diagonal covariance (fastest); "lowrank" low-rank plus diagonal (balanced); "fullrank" full Cholesky (best for small models).

Stochastic gradient backends (mode = "sghmc" / "sgld"):

batch_size: minibatch size (default ceiling(sqrt(n))).
epsilon: step size.
alpha: friction coefficient ("sghmc" only).
schedule_a, schedule_b, schedule_gamma, use_schedule: step-size decay schedule ("sgld" only).

Value

A ratiod_fit object containing:

draws: Posterior draws matrix
formula: Parsed formula specification
family: Model family
data: Original data
backend: Inference backend used
diagnostics: MCMC diagnostics (if applicable)

Examples

# Quick example: create model object (no fitting)
# A site effect u shared by catch and effort
set.seed(123)
n <- 60
site <- factor(rep(1:10, each = 6))
u <- rnorm(10, sd = 0.5)[site]
depth <- rnorm(n)
season <- factor(sample(c("spring", "summer"), n, replace = TRUE))
df <- data.frame(
  catch = rpois(n, exp(1.5 + 0.3 * depth + 0.4 * (season == "summer") + u)),
  effort = rgamma(n, shape = 2, rate = 0.5 * exp(-u)),
  depth = depth,
  season = season,
  site = site
)

# \donttest{
# CPUE example with combined formula
fit <- tratio(
  catch | effort ~ depth + season + (1 | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = list(iter = 400, chains = 1, seed = 1)
)

# Different predictors for each process
fit2 <- tratio(
  catch | effort ~ (1 | site),
  formula_num = ~ depth + season,
  formula_denom = ~ season,
  data = df,
  family = ratiod_poisson_gamma(),
  control = list(iter = 400, chains = 1, seed = 1)
)

# Extract ratio posteriors
cpue <- ratio(fit)
summary(cpue)
#> tulpaRatio ratio summary (response scale)
#> ===========================================
#> 
#>  obs      mean         sd      q2.5       q50     q97.5
#>    1 3.7436438 0.60610916 2.6959559 3.7189720 5.1930981
#>    2 2.2063724 0.28036492 1.7564337 2.1717247 2.8373256
#>    3 2.2619399 0.28851618 1.7900067 2.2267626 2.9056182
#>    4 1.2532252 0.20099361 0.9561902 1.2181848 1.6838281
#>    5 1.2712857 0.19049725 0.9594541 1.2827839 1.6490979
#>    6 5.3053128 1.08025255 3.4702414 5.3036686 7.7003296
#>    7 2.3994758 0.31007898 1.8858403 2.3754053 3.0740724
#>    8 0.3620751 0.08935861 0.2281328 0.3525959 0.5697255
#>    9 2.7164337 0.36686958 2.0611329 2.6919919 3.5191970
#>   10 1.3359713 0.19437160 1.0178089 1.3416604 1.7155200
#> 
#> ... and 50 more rows

# Compare ratios between groups
ratio_contrast(fit, ~ season)
#> tulpaRatio ratio contrasts
#> =====================
#> 
#> Type: difference 
#> Variable: season 
#> Reference: spring 
#> 
#>          contrast     mean        sd      q2.5      q50     q97.5 prob_positive
#>  summer vs spring 0.565998 0.2071027 0.2071584 0.577519 0.9712434             1
# }

Fit a Ratio, Rate, or Proportion Model

Usage

Arguments

Value

See also

Examples