Quick Start

Installation

pak::pak("gcol33/tulpaRatio")

library(tulpaRatio)

What tulpaRatio Does

tulpaRatio fits Bayesian hierarchical models for ratios, rates, and proportions. Instead of modelling ratios directly (which is statistically problematic), tulpaRatio jointly models the numerator and denominator processes, then derives ratios with full uncertainty propagation.

See vignette("philosophy") for the statistical motivation.

The Three Families

tulpaRatio provides three model families for common ratio types:

Family	Numerator	Denominator	Use Case
`ratiod_negbin_negbin()`	Counts	Counts	Relative abundance, species ratios
`ratiod_binomial()`	Successes	Trials	Detection rates, survival
`ratiod_poisson_gamma()`	Counts	Continuous	CPUE, catch per effort

Basic Example: CPUE Analysis

Catch per unit effort with site-level variation:

# Simulated trawl data
set.seed(123)
n_sites <- 20
n_obs <- 200

trawl_data <- data.frame(
  site = factor(rep(1:n_sites, each = n_obs/n_sites)),
  depth = rnorm(n_obs),  # standardised depth
  season = factor(sample(c("spring", "summer", "fall"), n_obs, replace = TRUE)),
  weather = factor(sample(c("calm", "rough"), n_obs, replace = TRUE))
)
# A site effect raises both effort and catch; catch per unit effort
# falls with depth
site_effect <- rnorm(n_sites, 0, 0.4)[trawl_data$site]
trawl_data$effort <- rgamma(n_obs, shape = 5, rate = 5 / exp(1 + site_effect))
trawl_data$catch <- rpois(n_obs, trawl_data$effort *
                            exp(1 + site_effect - 0.3 * trawl_data$depth))

# Fit model
fit <- tratio(
  catch | effort ~ depth + season + (1 | site),
  data = trawl_data,
  family = ratiod_poisson_gamma(),
  control = list(seed = 1)
)

# View results
summary(fit)

Extract CPUE Posteriors

# Observation-level CPUE
cpue <- ratio(fit)
summary(cpue)

# Compare CPUE between seasons
ratio_contrast(fit, ~ season)

Formula Syntax

Combined Formula (Recommended)

Use num | denom ~ predictors when both processes share predictors:

fit <- tratio(
  catch | effort ~ depth + season + (1 | site),
  data = trawl_data,
  family = ratiod_poisson_gamma(),
  control = list(seed = 1)
)

Process-Specific Predictors

Add terms that affect only numerator or denominator:

fit_split <- tratio(
  catch | effort ~ (1 | site),           # Shared structure
  formula_num = ~ depth + season,         # Catch-specific

  formula_denom = ~ weather,              # Effort-specific
  data = trawl_data,
  family = ratiod_poisson_gamma(),
  control = list(seed = 1)
)
coef(fit_split)

Random Effects

The examples in this section use simulated counts y out of totals n, with a covariate x, plots nested in sites, and survey years crossed with sites:

set.seed(42)
counts <- expand.grid(plot = 1:3, site = 1:12, year = 1:10)
site_eff <- rnorm(12, 0, 0.4)
year_eff <- rnorm(10, 0, 0.3)
plot_eff <- rnorm(36, 0, 0.2)
slope_eff <- rnorm(12, 0, 0.2)
eta <- site_eff[counts$site] + year_eff[counts$year] +
  plot_eff[(counts$site - 1) * 3 + counts$plot]
counts$x <- rnorm(nrow(counts))
counts$n <- rnbinom(nrow(counts), size = 10, mu = 60 * exp(eta))
counts$y <- rnbinom(nrow(counts), size = 10,
                    mu = 12 * exp(eta + (0.3 + slope_eff[counts$site]) * counts$x))
counts[c("plot", "site", "year")] <- lapply(counts[c("plot", "site", "year")], factor)
ctrl <- list(seed = 1)

Random Intercepts

# Single grouping factor
fit_re <- tratio(y | n ~ x + (1 | site), data = counts, control = ctrl)

# Multiple (crossed) grouping factors
fit_re <- tratio(y | n ~ x + (1 | site) + (1 | year), data = counts, control = ctrl)

# Nested grouping factors
fit_re <- tratio(y | n ~ x + (1 | site/plot), data = counts, control = ctrl)

Random Slopes

# Correlated intercept and slope
fit_re <- tratio(y | n ~ x + (1 + x | site), data = counts, control = ctrl)

# Uncorrelated (independent) intercept and slope
fit_re <- tratio(y | n ~ x + (1 + x || site), data = counts, control = ctrl)

Model Checking

# Posterior predictive checks
pp_check(fit, type = "dens_overlay")

# Model comparison
loo(fit)
waic(fit)

# Compare models
fit1 <- tratio(catch | effort ~ depth + (1 | site), data = trawl_data,
               family = ratiod_poisson_gamma(), control = ctrl)
fit2 <- tratio(catch | effort ~ depth + season + (1 | site), data = trawl_data,
               family = ratiod_poisson_gamma(), control = ctrl)
ratiod_compare(fit1, fit2)

The |elpd_diff| < 4 flag marks a difference too small to separate the two models: adding season does not improve the expected predictive accuracy.

Both criteria print warnings here. Each site effect is estimated from ten observations, so leaving one observation out moves it noticeably: waic() then reports observations whose p_waic exceeds 0.4 and recommends loo(), and loo() flags the observations whose Pareto \(k\) is too high for importance sampling to follow. Their contributions to elpd_loo are less reliable than the rest, which the standard errors do not show.

Quick Reference

Function	Purpose
`tratio()`	Fit model
`ratio()`	Extract ratio posteriors
`ratio_contrast()`	Compare ratios between groups
`pp_check()`	Posterior predictive checks
`loo()`, `waic()`	Model comparison
`ratiod_compare()`	Compare multiple models