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(), spatiotemporal_gp().

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

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 during fitting (default TRUE).

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.
riemannian: enable per-trajectory SoftAbs metric retry on divergence, computing a local Hessian-based metric and retrying. NULL (default) enables for BYM2 and ICAR with a dense mass matrix. Costs about (p+1) gradient evaluations per divergent trajectory only.
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)
set.seed(123)
df <- data.frame(
  catch = rpois(50, 10),
  effort = rgamma(50, 2, 0.5),
  depth = rnorm(50),
  season = factor(sample(c("spring", "summer"), 50, replace = TRUE)),
  site = factor(sample(1:5, 50, replace = TRUE))
)

# \donttest{
# CPUE example with combined formula
fit <- tratio(
  catch | effort ~ depth + season + (1 | site),
  data = df,
  family = ratiod_poisson_gamma(),
  control = list(iter = 200, warmup = 100, chains = 1)
)
#> Inference: Exact (Tier 1)
#>   Backend: hmc
#>   Reason: default (full MCMC)
#> Fitting ratio model...
#>   Family: poisson_gamma
#>   Observations: 50
#>   Iterations: 200 (warmup: 100)
#> Running NUTS sampler...
#>   Parameters: 13
#>   Iterations: 200 (warmup: 100)
#>   Chains: 1 (cores: 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 = 200, warmup = 100, chains = 1)
)
#> Inference: Exact (Tier 1)
#>   Backend: hmc
#>   Reason: default (full MCMC)
#> Fitting ratio model...
#>   Family: poisson_gamma
#>   Observations: 50
#>   Iterations: 200 (warmup: 100)
#> Running NUTS sampler...
#>   Parameters: 12
#>   Iterations: 200 (warmup: 100)
#>   Chains: 1 (cores: 1)

# Extract ratio posteriors
cpue <- ratio(fit)
summary(cpue)
#> tulpaRatio ratio summary (response scale)
#> ===========================================
#> 
#>  obs     mean        sd     q2.5      q50    q97.5
#>    1 2.328175 0.4053337 1.578950 2.329428 3.002349
#>    2 2.329070 0.4374686 1.542388 2.339878 3.071103
#>    3 2.964943 0.6549895 2.021628 2.879466 4.172565
#>    4 2.335104 0.3441534 1.719054 2.320751 3.100055
#>    5 2.329905 0.4541541 1.525576 2.333010 3.142023
#>    6 2.328280 0.3856676 1.605752 2.350520 2.951722
#>    7 2.329694 0.3617594 1.648517 2.352288 2.967819
#>    8 3.103446 0.8767279 1.803543 3.011397 4.743132
#>    9 2.966296 0.6018358 2.008515 2.856380 4.106194
#>   10 2.337615 0.3444195 1.741895 2.307394 3.140574
#> 
#> ... and 40 more rows

# Compare ratios between groups
ratio_contrast(fit, ~ season)
#> tulpaRatio ratio contrasts
#> =====================
#> 
#> Type: difference 
#> Variable: season 
#> Reference: 1 
#> 
#>          contrast      mean        sd       q2.5       q50     q97.5
#>  summer vs spring 0.2398869 0.2151743 -0.1304036 0.2816063 0.6367514
#>  prob_positive
#>            0.8
# }

Fit a ratio, rate, or proportion model

Usage

Arguments

Value

See also

Examples