| Title: | Bayesian Hierarchical Models for Ratios, Rates, and Proportions |
|---|---|
| Description: | An opinionated framework for Bayesian hierarchical modelling of ratios, rates, and proportions. Ratios are treated as derived quantities, with all inference performed on the latent processes generating the numerator and denominator. Supports shared random effects, latent factors for unmeasured confounders, spatial structure ('CAR'/'BYM2'/'GP'), spatiotemporal interaction (Knorr-Held Type I-IV), zero-inflation variants, and full uncertainty propagation to ratio posteriors. Provides multiple inference backends including custom 'HMC'/'NUTS' (default), 'Laplace' approximation, and 'Polya-Gamma' Gibbs sampling. |
| Authors: | Gilles Colling [aut, cre, cph] (ORCID: <https://orcid.org/0000-0003-3070-6066>) |
| Maintainer: | Gilles Colling <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.7.3 |
| Built: | 2026-10-13 12:05:29 UTC |
| Source: | https://github.com/gcol33/tulparatio |
Convert a fitted ratio model to a format compatible with the posterior and tidybayes packages. Returns draws in a tidy format suitable for further analysis and visualization.
## S3 method for class 'ratiod_fit' as_draws(x, ...)## S3 method for class 'ratiod_fit' as_draws(x, ...)
x |
A |
... |
Additional arguments passed to |
The returned draws object is compatible with:
posterior package functions (summarise_draws, rvar, etc.)
tidybayes package functions (spread_draws, gather_draws, etc.)
bayesplot package visualization functions
To work in tidybayes, extract the draws with this function and pass them on to its verbs.
A draws_df object from the posterior package.
posterior::as_draws_df(), ratio() for ratio-specific extraction
# Generate synthetic data set.seed(111) n <- 48 site <- rep(letters[1:12], each = 4) u <- rnorm(12, sd = 0.5)[match(site, letters[1:12])] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(3 + 0.3 * depth + u)), effort = rgamma(n, shape = 11, rate = exp(-u)), depth = depth, site = site ) # Fit model fit <- tratio( count | effort ~ depth + (1 | site), data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Convert to draws format draws <- as_draws(fit) # Use with posterior package if (requireNamespace("posterior", quietly = TRUE)) { posterior::summarise_draws(draws) }# Generate synthetic data set.seed(111) n <- 48 site <- rep(letters[1:12], each = 4) u <- rnorm(12, sd = 0.5)[match(site, letters[1:12])] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(3 + 0.3 * depth + u)), effort = rgamma(n, shape = 11, rate = exp(-u)), depth = depth, site = site ) # Fit model fit <- tratio( count | effort ~ depth + (1 | site), data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Convert to draws format draws <- as_draws(fit) # Use with posterior package if (requireNamespace("posterior", quietly = TRUE)) { posterior::summarise_draws(draws) }
Posterior Draws of a Fitted Ratio Model as a Data Frame
## S3 method for class 'ratiod_fit' as.data.frame(x, ...)## S3 method for class 'ratiod_fit' as.data.frame(x, ...)
x |
A ratiod_fit object |
... |
Ignored. |
A draws_df data frame: one row per posterior draw, one column per
parameter, plus the .chain, .iteration and .draw columns.
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) draws <- as.data.frame(fit) dim(draws) names(draws)[1:6]sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) draws <- as.data.frame(fit) dim(draws) names(draws)[1:6]
Checks MCMC diagnostics (Rhat, ESS, divergences) and prints actionable warnings with suggestions for fixing issues.
check_diagnostics(fit, quiet = FALSE)check_diagnostics(fit, quiet = FALSE)
fit |
A ratiod_fit object |
quiet |
Logical; if TRUE, suppress output and just return status. |
Invisibly returns a list with diagnostic status:
Logical; TRUE if all diagnostics pass
Number of divergent transitions
Number of parameters with Rhat > 1.01
Number of parameters with ESS < 400
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) check_diagnostics(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) check_diagnostics(fit)
Provides a comprehensive diagnostic report for a ratio model, combining convergence metrics, divergence information, and actionable recommendations.
diagnostic_summary(fit) ## S3 method for class 'ratiod_diagnostic_summary' print(x, ...)diagnostic_summary(fit) ## S3 method for class 'ratiod_diagnostic_summary' print(x, ...)
fit |
A |
x |
A ratiod_diagnostic_summary object. |
... |
Ignored. |
A list with class ratiod_diagnostic_summary containing:
Overall status: "PASS", "WARN", or "FAIL"
Number of divergent transitions
Percentage of divergent transitions
Data frame of parameters with worst Rhat
Data frame of parameters with worst ESS
E-BFMI value (HMC only)
Character vector of recommendations
check_diagnostics(), diagnostics(), plot_diagnostics()
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) diagnostic_summary(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) diagnostic_summary(fit)
Computes per-parameter improved Rhat and bulk / tail effective sample size
by delegating to the engine's tulpa::diagnostics(), the single source of
truth for convergence diagnostics in the tulpa ecosystem. Works with all
backends.
## S3 method for class 'ratiod_fit' diagnostics(fit, pars = NULL, ...) ## S3 method for class 'ratiod_diagnostics' print(x, ...)## S3 method for class 'ratiod_fit' diagnostics(fit, pars = NULL, ...) ## S3 method for class 'ratiod_diagnostics' print(x, ...)
fit |
A ratiod_fit object |
pars |
Character vector of parameter names (default: all). |
... |
Ignored. |
x |
A ratiod_diagnostics object |
Rhat is the improved value of Vehtari et al. (2021): the maximum of rank-normalized split-Rhat and folded split-Rhat. Split-Rhat is defined for a single chain (the chain is split in half), so single-chain fits also get a value.
A data frame with columns: parameter, rhat, ess_bulk, ess_tail.
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) diagnostics(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) diagnostics(fit)
Summarize the draws of the named parameters as a point estimate (median or mean) and a credible interval (quantile or highest-density), one row per parameter per interval width.
draws_interval( object, ..., .width = 0.95, .point = c("median", "mean"), .interval = c("qi", "hdi") ) ## S3 method for class 'ratiod_fit' draws_interval( object, ..., .width = 0.95, .point = c("median", "mean"), .interval = c("qi", "hdi") )draws_interval( object, ..., .width = 0.95, .point = c("median", "mean"), .interval = c("qi", "hdi") ) ## S3 method for class 'ratiod_fit' draws_interval( object, ..., .width = 0.95, .point = c("median", "mean"), .interval = c("qi", "hdi") )
object |
A |
... |
Parameter names to summarize (unquoted). |
.width |
Width(s) of credible intervals. Default 0.95. |
.point |
Point estimate function: "median" (default) or "mean". |
.interval |
Interval type: "qi" for quantile interval (default), "hdi" for highest density interval. |
A data frame with columns for the variable name, point estimate, and interval bounds.
# Generate synthetic data set.seed(444) n <- 48 site <- rep(letters[1:12], each = 4) u <- rnorm(12, sd = 0.5)[match(site, letters[1:12])] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(3 + 0.3 * depth + u)), effort = rgamma(n, shape = 9, rate = exp(-u)), depth = depth, site = site ) # Fit model fit <- tratio( count | effort ~ depth + (1 | site), data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Get point estimates and 95% intervals draws_interval(fit, beta_num, beta_denom) # Multiple interval widths draws_interval(fit, beta_num, .width = c(0.5, 0.9, 0.95))# Generate synthetic data set.seed(444) n <- 48 site <- rep(letters[1:12], each = 4) u <- rnorm(12, sd = 0.5)[match(site, letters[1:12])] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(3 + 0.3 * depth + u)), effort = rgamma(n, shape = 9, rate = exp(-u)), depth = depth, site = site ) # Fit model fit <- tratio( count | effort ~ depth + (1 | site), data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Get point estimates and 95% intervals draws_interval(fit, beta_num, beta_denom) # Multiple interval widths draws_interval(fit, beta_num, .width = c(0.5, 0.9, 0.95))
The long counterpart of draws_wide(): one row per draw per parameter,
with the parameter in .variable and its value in .value.
draws_long(object, ..., regex = FALSE, ndraws = NULL) ## S3 method for class 'ratiod_fit' draws_long(object, ..., regex = FALSE, ndraws = NULL)draws_long(object, ..., regex = FALSE, ndraws = NULL) ## S3 method for class 'ratiod_fit' draws_long(object, ..., regex = FALSE, ndraws = NULL)
object |
A |
... |
Parameter names to extract (unquoted). |
regex |
Logical; if TRUE, treat parameter names as regex patterns. |
ndraws |
Number of draws to return. If NULL, returns all. |
A data frame with columns .chain, .iteration, .draw,
.variable, and .value.
# Generate synthetic data set.seed(333) n <- 48 site <- rep(letters[1:12], each = 4) u <- rnorm(12, sd = 0.5)[match(site, letters[1:12])] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(3 + 0.3 * depth + u)), effort = rgamma(n, shape = 12, rate = exp(-u)), depth = depth, site = site ) # Fit model fit <- tratio( count | effort ~ depth + (1 | site), data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Gather all fixed effects draws_long(fit, beta_num, beta_denom) # Use for plotting if (requireNamespace("ggplot2", quietly = TRUE)) { library(ggplot2) draws_long(fit, beta_num) |> ggplot(aes(x = .value)) + geom_histogram(bins = 40) + facet_wrap(~ .variable) }# Generate synthetic data set.seed(333) n <- 48 site <- rep(letters[1:12], each = 4) u <- rnorm(12, sd = 0.5)[match(site, letters[1:12])] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(3 + 0.3 * depth + u)), effort = rgamma(n, shape = 12, rate = exp(-u)), depth = depth, site = site ) # Fit model fit <- tratio( count | effort ~ depth + (1 | site), data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Gather all fixed effects draws_long(fit, beta_num, beta_denom) # Use for plotting if (requireNamespace("ggplot2", quietly = TRUE)) { library(ggplot2) draws_long(fit, beta_num) |> ggplot(aes(x = .value)) + geom_histogram(bins = 40) + facet_wrap(~ .variable) }
Extract posterior draws for the named parameters into a wide data frame,
one row per draw and one column per parameter, alongside the .chain,
.iteration and .draw bookkeeping columns.
Parameters are named unquoted, and a name matches every indexed element it
prefixes, so beta_num reaches beta_num[1], beta_num[2] and so on.
draws_wide(object, ..., regex = FALSE, ndraws = NULL) ## S3 method for class 'ratiod_fit' draws_wide(object, ..., regex = FALSE, ndraws = NULL)draws_wide(object, ..., regex = FALSE, ndraws = NULL) ## S3 method for class 'ratiod_fit' draws_wide(object, ..., regex = FALSE, ndraws = NULL)
object |
A |
... |
Parameter names to extract (unquoted). Supports patterns like
|
regex |
Logical; if TRUE, treat parameter names as regex patterns. |
ndraws |
Number of draws to return. If NULL, returns all. |
A data frame with columns .chain, .iteration, .draw, and
one column per requested parameter.
# Generate synthetic data set.seed(222) n <- 48 site <- rep(letters[1:12], each = 4) u <- rnorm(12, sd = 0.5)[match(site, letters[1:12])] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(3 + 0.3 * depth + u)), effort = rgamma(n, shape = 10, rate = exp(-u)), depth = depth, site = site ) # Fit model fit <- tratio( count | effort ~ depth + (1 | site), data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Extract fixed effects draws_wide(fit, beta_num, beta_denom) # Subsample draws draws_wide(fit, sigma_re, ndraws = 100)# Generate synthetic data set.seed(222) n <- 48 site <- rep(letters[1:12], each = 4) u <- rnorm(12, sd = 0.5)[match(site, letters[1:12])] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(3 + 0.3 * depth + u)), effort = rgamma(n, shape = 10, rate = exp(-u)), depth = depth, site = site ) # Fit model fit <- tratio( count | effort ~ depth + (1 | site), data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Extract fixed effects draws_wide(fit, beta_num, beta_denom) # Subsample draws draws_wide(fit, sigma_re, ndraws = 100)
Compute posterior distributions of fitted values (expected values) for both the numerator and denominator processes from a fitted ratio model.
## S3 method for class 'ratiod_fit' fitted( object, component = c("all", "numerator", "denominator", "ratio"), summary = TRUE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_fitted' print(x, n = 10, ...)## S3 method for class 'ratiod_fit' fitted( object, component = c("all", "numerator", "denominator", "ratio"), summary = TRUE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_fitted' print(x, n = 10, ...)
object |
A |
component |
Which component to return: "numerator", "denominator", "ratio", or "all" (default). |
summary |
Logical; if TRUE, return summary statistics instead of full posterior draws (default TRUE). |
probs |
Quantiles to compute if |
... |
Ignored |
x |
A ratiod_fitted object |
n |
Number of rows to print |
If summary = TRUE, a data frame with posterior summaries.
If summary = FALSE, a list with matrices of posterior draws.
# Generate synthetic data set.seed(789) n <- 40 df <- data.frame( count = rpois(n, lambda = 15), effort = rgamma(n, shape = 8, rate = 1), depth = rnorm(n) ) # Fit model fit <- tratio( count | effort ~ depth, data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Get fitted value summaries fitted(fit) # Get just numerator fitted values fitted(fit, component = "numerator") # Get full posterior draws fitted(fit, summary = FALSE)# Generate synthetic data set.seed(789) n <- 40 df <- data.frame( count = rpois(n, lambda = 15), effort = rgamma(n, shape = 8, rate = 1), depth = rnorm(n) ) # Fit model fit <- tratio( count | effort ~ depth, data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Get fitted value summaries fitted(fit) # Get just numerator fitted values fitted(fit, component = "numerator") # Get full posterior draws fitted(fit, summary = FALSE)
Performs Geweke's convergence diagnostic, comparing the mean of the first portion of a chain to the last portion. Useful for single-chain diagnostics.
geweke_test(fit, frac1 = 0.1, frac2 = 0.5, pars = NULL) ## S3 method for class 'ratiod_geweke' print(x, ...)geweke_test(fit, frac1 = 0.1, frac2 = 0.5, pars = NULL) ## S3 method for class 'ratiod_geweke' print(x, ...)
fit |
A |
frac1 |
Fraction of chain for first window (default: 0.1). |
frac2 |
Fraction of chain for second window (default: 0.5). |
pars |
Character vector of parameter names (default: all main params). |
x |
A ratiod_geweke object. |
... |
Ignored. |
The Geweke test computes a z-score comparing the means of early and late portions of a chain. Large z-scores (|z| > 2) indicate the chain has not converged.
A data frame with columns: parameter, z_score, p_value.
diagnostics(), check_diagnostics()
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) geweke_test(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) geweke_test(fit)
Checks whether tulpaRatio was compiled with GPU support (CUDA or OpenCL) and whether compatible GPU hardware is detected.
gpu_available()gpu_available()
GPU support requires:
tulpaRatio compiled with GPU support (-DUSE_CUDA or -DUSE_OPENCL)
Compatible GPU hardware (NVIDIA for CUDA, any GPU for OpenCL)
Appropriate drivers installed
The solver = "gpu" option in spatial_gp() will automatically fall back
to PCG if GPU support is unavailable.
Logical; TRUE if GPU acceleration is available, FALSE otherwise.
spatial_gp() for using GPU-accelerated GP models
# Check if GPU is available gpu_available() # Get detailed GPU info if (gpu_available()) { gpu_info() }# Check if GPU is available gpu_available() # Get detailed GPU info if (gpu_available()) { gpu_info() }
Returns information about available GPU devices including name, memory, compute capability, and backend (CUDA or OpenCL).
gpu_info() ## S3 method for class 'tulpaRatio_gpu_info' print(x, ...)gpu_info() ## S3 method for class 'tulpaRatio_gpu_info' print(x, ...)
x |
GPU info object from |
... |
Ignored |
A list with GPU information, or NULL if no GPU is available.
Logical; whether GPU is available
Character; "cuda", "opencl", or "none"
Integer; number of GPU devices
List of device info (name, memory, compute capability)
info <- gpu_info() if (!is.null(info) && info$available) { cat("GPU backend:", info$backend, "\n") cat("Devices:", info$device_count, "\n") }info <- gpu_info() if (!is.null(info) && info$available) { cat("GPU backend:", info$backend, "\n") cat("Devices:", info$device_count, "\n") }
Functions to check GPU availability and capabilities for accelerated GP computations.
gpu_available() ratiod_threads()gpu_available() ratiod_threads()
Hurdle Binomial Specification
hurdle_binomial(formula = NULL)hurdle_binomial(formula = NULL)
formula |
Optional formula for the hurdle probability (P(Y > 0)). If NULL, uses intercept-only model. |
A ratiod_zi object
# Intercept-only hurdle hurdle_binomial()# Intercept-only hurdle hurdle_binomial()
Hurdle Negative Binomial Specification
hurdle_negbin(formula = NULL)hurdle_negbin(formula = NULL)
formula |
Optional formula for the hurdle probability (P(Y > 0)). If NULL, uses intercept-only model. |
A ratiod_zi object
# Intercept-only hurdle hurdle_negbin() # Hurdle probability depends on sampling effort hurdle_negbin(~ log_effort)# Intercept-only hurdle hurdle_negbin() # Hurdle probability depends on sampling effort hurdle_negbin(~ log_effort)
Hurdle Poisson Specification
hurdle_poisson(formula = NULL)hurdle_poisson(formula = NULL)
formula |
Optional formula for the hurdle probability (P(Y > 0)). If NULL, uses intercept-only model. |
A ratiod_zi object
# Intercept-only hurdle hurdle_poisson() # Hurdle probability depends on effort hurdle_poisson(~ effort)# Intercept-only hurdle hurdle_poisson() # Hurdle probability depends on effort hurdle_poisson(~ effort)
Displays information about the inference modes available in tulpaRatio, including their tiers, guarantees, and appropriate use cases.
inference_mode_info()inference_mode_info()
inference_mode_info()inference_mode_info()
Define latent factors for capturing unmeasured structure shared between numerator and denominator. Each observation carries factor scores, and each process loads on them with its own loading, so the factors model the residual covariance of an observation's two linear predictors.
latent_factor( n_factors = 1L, prior = NULL, shared = TRUE, constraint = c("sum_to_zero", "first_zero") ) ## S3 method for class 'ratiod_latent' print(x, ...) ## S3 method for class 'ratiod_fit' latent_factors(fit, summary = TRUE, probs = c(0.025, 0.5, 0.975), ...)latent_factor( n_factors = 1L, prior = NULL, shared = TRUE, constraint = c("sum_to_zero", "first_zero") ) ## S3 method for class 'ratiod_latent' print(x, ...) ## S3 method for class 'ratiod_fit' latent_factors(fit, summary = TRUE, probs = c(0.025, 0.5, 0.975), ...)
n_factors |
Integer; number of latent factors. Default is 1. With
|
prior |
Prior for the diagonal loadings. Default is a PC prior with P(loading > 1) = 0.01, which shrinks toward simpler models. |
shared |
Logical; if TRUE (default), the factors load on both the numerator and the denominator linear predictor. If FALSE, they load on the numerator only. |
constraint |
Identifiability constraint for the factor scores:
|
x |
A ratiod_latent object |
... |
Ignored |
fit |
A |
summary |
Logical; |
probs |
Quantiles for the summary. |
When modeling ratios, the numerator and denominator often share unmeasured drivers. For example, in relative abundance data, both the focal species count and total count might be affected by:
Observer skill (unmeasured)
Local microhabitat conditions (unmeasured)
Weather on sampling day (unmeasured)
Latent factors capture this shared structure without requiring the drivers to be measured.
For observation i with K latent factors:
where:
are factor scores, with
for each k
is the loading matrix, lower triangular
with a positive diagonal
The factors add the covariance between the two
linear predictors of an observation. A 2 x 2 covariance has three free
entries: one factor (two loadings) fills it up to rank one and two factors
(three loadings) fill it completely, so a third factor would not be
identified. Fixing lower triangular with a positive diagonal
makes it the Cholesky factor of that covariance, which is unique; without
it the loadings and scores rotate freely.
The ratio's linear predictor gains
: a factor that loads
equally on both processes cancels in the ratio, one that loads unequally
moves it.
Each diagonal loading carries prior (a PC prior by default). Each
off-diagonal loading carries a normal prior centred at zero, with the
standard deviation that gives the same 99% quantile as the
diagonal prior.
n_factors = 1: one unmeasured driver acting on both processes, with
process-specific strength
n_factors = 2: an unrestricted residual covariance between the two
processes
Latent factors differ from random effects in several ways:
Random effects are grouped (e.g., site-level), factors are observation-level
Random effects require grouping structure, factors don't
Factors capture residual correlation not explained by observed predictors
You can use both together: random effects for known grouping, factors for residual unmeasured structure.
A ratiod_latent object for use in tratio().
tratio() for model fitting, prior_pc() for prior specification
# Basic latent factor (single factor on both processes) latent_factor() # Two latent factors: an unrestricted residual covariance latent_factor(n_factors = 2) # Custom prior (more regularization) latent_factor(n_factors = 1, prior = prior_pc(U = 0.5, alpha = 0.01)) # Numerator-only factor (not shared) latent_factor(n_factors = 1, shared = FALSE) # Use in model fitting set.seed(1) n <- 60 df <- data.frame( species_count = rnbinom(n, size = 5, mu = 15), total_count = rnbinom(n, size = 8, mu = 100), habitat = factor(rep(c("forest", "grassland"), each = n / 2)), site = factor(rep(1:10, each = n / 10)) ) fit <- tratio( species_count | total_count ~ habitat + (1 | site), data = df, family = ratiod_negbin_negbin(), latent = latent_factor(n_factors = 1), control = list(iter = 400, chains = 1, seed = 1) ) head(latent_factors(fit))# Basic latent factor (single factor on both processes) latent_factor() # Two latent factors: an unrestricted residual covariance latent_factor(n_factors = 2) # Custom prior (more regularization) latent_factor(n_factors = 1, prior = prior_pc(U = 0.5, alpha = 0.01)) # Numerator-only factor (not shared) latent_factor(n_factors = 1, shared = FALSE) # Use in model fitting set.seed(1) n <- 60 df <- data.frame( species_count = rnbinom(n, size = 5, mu = 15), total_count = rnbinom(n, size = 8, mu = 100), habitat = factor(rep(c("forest", "grassland"), each = n / 2)), site = factor(rep(1:10, each = n / 10)) ) fit <- tratio( species_count | total_count ~ habitat + (1 | site), data = df, family = ratiod_negbin_negbin(), latent = latent_factor(n_factors = 1), control = list(iter = 400, chains = 1, seed = 1) ) head(latent_factors(fit))
Compute leave-one-out cross-validation using Pareto-smoothed importance sampling (PSIS-LOO).
## S3 method for class 'ratiod_fit' loo(x, ...)## S3 method for class 'ratiod_fit' loo(x, ...)
x |
A |
... |
Additional arguments passed to loo::loo |
The pointwise log-likelihood is reconstructed from the posterior draws and
the model (see ratiod_compare()); it works on any fit produced by the
hmc, ess, sghmc, or vi backend, for every family including the zero-inflated,
one-inflated and hurdle ones. Each observation's likelihood is evaluated at
the full linear predictor: fixed, random and structured effects.
A loo object
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) loo(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) loo(fit)
Number of Divergent Transitions
n_divergent(fit)n_divergent(fit)
fit |
A ratiod_fit object |
Integer count of divergent transitions
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) n_divergent(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) n_divergent(fit)
One-Inflated Binomial Specification
oi_binomial(formula = NULL)oi_binomial(formula = NULL)
formula |
Optional formula for the one-inflation probability (excess at the upper boundary, y = n). If NULL, uses intercept-only model. |
A ratiod_zi object
# Intercept-only OI oi_binomial()# Intercept-only OI oi_binomial()
Creates autocorrelation function (ACF) plots for selected parameters. High autocorrelation indicates slow mixing and low effective sample size.
plot_acf(fit, pars = NULL, lags = 25, n_pars = 6)plot_acf(fit, pars = NULL, lags = 25, n_pars = 6)
fit |
A |
pars |
Character vector of parameter names. If NULL, selects worst-mixing parameters based on ESS. |
lags |
Maximum number of lags to compute (default: 25). |
n_pars |
Maximum number of parameters to plot (default: 6). |
Ideal ACF plots show rapid decay to zero. Slow decay indicates high autocorrelation, which reduces effective sample size and may indicate poor mixing.
A ggplot object (if ggplot2 available) or base R plot (invisible).
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) plot_acf(fit) plot_acf(fit, pars = c("beta_num[2]", "sigma_re"))sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) plot_acf(fit) plot_acf(fit, pars = c("beta_num[2]", "sigma_re"))
Creates visualizations to investigate divergent transitions. Parallel coordinates and scatter plots highlight where in parameter space divergences occur.
plot_divergences(fit, pars = NULL, type = c("parcoord", "scatter"))plot_divergences(fit, pars = NULL, type = c("parcoord", "scatter"))
fit |
A |
pars |
Character vector of parameter names. If NULL, uses variance parameters. |
type |
Plot type: "parcoord" (parallel coordinates) or "scatter". |
Divergent transitions indicate regions of high posterior curvature that the sampler cannot efficiently explore. Common causes:
Very narrow funnels (hierarchical models)
Strong correlations
Multi-modality
A ggplot object or base R plot (invisible).
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) n_divergent(fit) plot_divergences(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) n_divergent(fit) plot_divergences(fit)
Creates overlaid histograms of marginal energy and energy transition, along with the E-BFMI statistic. Low E-BFMI indicates poor exploration of the posterior.
plot_energy(fit)plot_energy(fit)
fit |
A |
E-BFMI (Energy Bayesian Fraction of Missing Information) compares the
distribution of energy levels to energy transitions. Values below 0.3
indicate the sampler may not be exploring the full posterior. The energy
is the Hamiltonian after each iteration's momentum draw; it is centred and
differenced within each chain, and the subtitle reports the lowest
per-chain E-BFMI (tulpa::energy_bfmi()).
A ggplot object (if ggplot2 available) or base R plot.
diagnostic_summary(), check_diagnostics()
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) plot_energy(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) plot_energy(fit)
Creates a visual display of effective sample size (ESS) for all parameters, expressed as a ratio of ESS to total samples. Low ESS indicates high autocorrelation.
plot_ess(fit, type = c("bulk", "tail"), threshold = 400, pars = NULL)plot_ess(fit, type = c("bulk", "tail"), threshold = 400, pars = NULL)
fit |
A |
type |
Type of ESS: "bulk" (default) or "tail". |
threshold |
Minimum acceptable ESS (default: 400). |
pars |
Character vector of parameter names to include. |
ESS/iter ratio interpretation:
Green: ESS >= threshold
Yellow: threshold/2 <= ESS < threshold
Red: ESS < threshold/2
A ggplot object (if ggplot2 available) or base R plot (invisible).
plot_rhat(), diagnostic_summary(), diagnostics()
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) plot_ess(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) plot_ess(fit)
Create publication-ready maps from tulpaRatio spatial predictions. Supports plotting ratio estimates, uncertainty (credible interval width), and individual process predictions (numerator/denominator).
plot_map( x, coords = NULL, what = c("ratio", "numerator", "denominator", "uncertainty"), summary = c("median", "mean", "q2.5", "q97.5", "sd"), newdata = NULL, title = NULL, palette = "viridis", points = FALSE, point_color = "grey30", point_size = 0.5, na_color = "transparent", crs = NULL, legend_title = NULL, ... )plot_map( x, coords = NULL, what = c("ratio", "numerator", "denominator", "uncertainty"), summary = c("median", "mean", "q2.5", "q97.5", "sd"), newdata = NULL, title = NULL, palette = "viridis", points = FALSE, point_color = "grey30", point_size = 0.5, na_color = "transparent", crs = NULL, legend_title = NULL, ... )
x |
A |
coords |
A data frame or matrix with spatial coordinates (columns
named 'x'/'y', 'X'/'Y', 'lon'/'lat', 'longitude'/'latitude', or
'Easting'/'Northing'). Required if |
what |
What to plot: "ratio" (default), "numerator", "denominator", or "uncertainty". |
summary |
Which summary statistic: "median" (default), "mean", "q2.5", "q97.5", or "sd". |
newdata |
Optional data frame with prediction locations and covariates.
If NULL and |
title |
Plot title. If NULL, auto-generated based on |
palette |
Color palette: "viridis" (default), "magma", "plasma", "inferno", "cividis", "mako", "rocket", or a custom vector of colors. |
points |
Logical; if TRUE, overlay observation points. Default FALSE. |
point_color |
Color for observation points. Default "grey30". |
point_size |
Size for observation points. Default 0.5. |
na_color |
Color for NA values. Default "transparent". |
crs |
Coordinate reference system (proj4 string or EPSG code). If NULL, uses planar coordinates. |
legend_title |
Title for the color legend. If NULL, auto-generated. |
... |
Additional arguments passed to ggplot2 theme functions. |
This function provides a streamlined workflow for visualizing spatial predictions from ratio models. It handles:
Extracting predictions from ratiod_fit objects
Converting to appropriate spatial format (stars/sf)
Creating publication-quality maps with sensible defaults
Uncertainty visualization via credible interval width
For custom maps or more control, extract predictions using ratio() or
predict() and use ggplot2 directly with geom_stars() or geom_sf().
A ggplot2 object that can be further customized.
This function requires ggplot2. For raster-style maps, stars and sf
are also needed. Install with:
install.packages(c("ggplot2", "stars", "sf"))
ratio() for extracting ratio posteriors,
predict.ratiod_fit() for predictions at new locations
# Twenty sites on a chain, each visited three times set.seed(123) n_sites <- 20 sites <- data.frame(site = factor(seq_len(n_sites)), x = sort(runif(n_sites)), y = runif(n_sites)) df <- sites[rep(seq_len(n_sites), 3), ] df$count <- rpois(nrow(df), lambda = 8) df$effort <- rgamma(nrow(df), shape = 4, rate = 1) df$elevation <- rnorm(nrow(df)) adj <- matrix(0, n_sites, n_sites) for (i in 1:(n_sites - 1)) adj[i, i + 1] <- adj[i + 1, i] <- 1 fit <- tratio( count | effort ~ elevation, data = df, family = ratiod_poisson_gamma(), spatial = spatial_car(adj, group_var = "site"), mode = "laplace" ) if (requireNamespace("ggplot2", quietly = TRUE)) { plot_map(fit, coords = df[, c("x", "y")]) }# Twenty sites on a chain, each visited three times set.seed(123) n_sites <- 20 sites <- data.frame(site = factor(seq_len(n_sites)), x = sort(runif(n_sites)), y = runif(n_sites)) df <- sites[rep(seq_len(n_sites), 3), ] df$count <- rpois(nrow(df), lambda = 8) df$effort <- rgamma(nrow(df), shape = 4, rate = 1) df$elevation <- rnorm(nrow(df)) adj <- matrix(0, n_sites, n_sites) for (i in 1:(n_sites - 1)) adj[i, i + 1] <- adj[i + 1, i] <- 1 fit <- tratio( count | effort ~ elevation, data = df, family = ratiod_poisson_gamma(), spatial = spatial_car(adj, group_var = "site"), mode = "laplace" ) if (requireNamespace("ggplot2", quietly = TRUE)) { plot_map(fit, coords = df[, c("x", "y")]) }
Create a multi-panel figure with ratio estimate and uncertainty maps side by side.
plot_map_panel(x, newdata = NULL, ncol = 2, ...)plot_map_panel(x, newdata = NULL, ncol = 2, ...)
x |
A |
newdata |
Optional data frame with prediction locations. |
ncol |
Number of columns in the grid. Default 2. |
... |
Additional arguments passed to |
A patchwork object (if patchwork is installed) or a list of ggplots.
# Twenty sites on a chain, each visited three times set.seed(123) n_sites <- 20 sites <- data.frame(site = factor(seq_len(n_sites)), x = sort(runif(n_sites)), y = runif(n_sites)) df <- sites[rep(seq_len(n_sites), 3), ] df$count <- rpois(nrow(df), lambda = 8) df$effort <- rgamma(nrow(df), shape = 4, rate = 1) df$elevation <- rnorm(nrow(df)) adj <- matrix(0, n_sites, n_sites) for (i in 1:(n_sites - 1)) adj[i, i + 1] <- adj[i + 1, i] <- 1 fit <- tratio( count | effort ~ elevation, data = df, family = ratiod_poisson_gamma(), spatial = spatial_car(adj, group_var = "site"), mode = "laplace" ) if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("patchwork", quietly = TRUE)) { plot_map_panel(fit, coords = df[, c("x", "y")], ncol = 1) }# Twenty sites on a chain, each visited three times set.seed(123) n_sites <- 20 sites <- data.frame(site = factor(seq_len(n_sites)), x = sort(runif(n_sites)), y = runif(n_sites)) df <- sites[rep(seq_len(n_sites), 3), ] df$count <- rpois(nrow(df), lambda = 8) df$effort <- rgamma(nrow(df), shape = 4, rate = 1) df$elevation <- rnorm(nrow(df)) adj <- matrix(0, n_sites, n_sites) for (i in 1:(n_sites - 1)) adj[i, i + 1] <- adj[i + 1, i] <- 1 fit <- tratio( count | effort ~ elevation, data = df, family = ratiod_poisson_gamma(), spatial = spatial_car(adj, group_var = "site"), mode = "laplace" ) if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("patchwork", quietly = TRUE)) { plot_map_panel(fit, coords = df[, c("x", "y")], ncol = 1) }
Creates a pairs plot showing bivariate relationships between parameters. Divergent transitions (if present) are highlighted to help identify problematic posterior regions.
plot_pairs( fit, pars = NULL, highlight_divergent = TRUE, n_pars = 5, alpha = 0.3 )plot_pairs( fit, pars = NULL, highlight_divergent = TRUE, n_pars = 5, alpha = 0.3 )
fit |
A |
pars |
Character vector of parameter names. If NULL, selects main variance parameters. |
highlight_divergent |
Logical; highlight divergent transitions in red (default: TRUE). |
n_pars |
Maximum number of parameters (default: 5). |
alpha |
Point transparency (default: 0.3). |
Pairs plots help identify:
Strong correlations between parameters (potential non-identifiability)
Multimodality
Regions where divergences cluster (indicating problematic geometry)
A ggplot object (if ggplot2/GGally available) or base R plot.
plot_divergences(), diagnostics()
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) plot_pairs(fit, pars = c("beta_num[1]", "beta_num[2]", "sigma_re"))sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) plot_pairs(fit, pars = c("beta_num[1]", "beta_num[2]", "sigma_re"))
Creates a visual display of Rhat values for all parameters, with color-coding to highlight convergence issues. Rhat values > 1.01 indicate potential convergence problems.
plot_rhat(fit, threshold = 1.01, pars = NULL)plot_rhat(fit, threshold = 1.01, pars = NULL)
fit |
A |
threshold |
Rhat threshold for warnings (default: 1.01). |
pars |
Character vector of parameter names to include. If NULL (default), includes all main parameters (excludes high-dimensional RE/spatial). |
Color coding:
Green: Rhat < 1.01 (converged)
Yellow: 1.01 <= Rhat < 1.05 (borderline)
Red: Rhat >= 1.05 (not converged)
A ggplot object (if ggplot2 available) or base R plot (invisible).
plot_ess(), diagnostic_summary(), diagnostics()
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) plot_rhat(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) plot_rhat(fit)
Draw replicated numerator and denominator responses from the posterior
predictive distribution: at each posterior draw the two fitted means are
assembled from that draw's parameters, and a response is drawn from the
family at those means. The replicates are what pp_check() compares the
observed data against.
## S3 method for class 'ratiod_fit' posterior_predict(object, ndraws = NULL, ...)## S3 method for class 'ratiod_fit' posterior_predict(object, ndraws = NULL, ...)
object |
A |
ndraws |
Number of posterior draws to replicate. |
... |
Unused. |
A list with y_num_rep and y_denom_rep, each a draws-by-observation
matrix, plus y_num and y_denom, the observed responses. A family with a
fixed denominator (the binomial arms, whose denominator is the trial count)
reports the trials as y_denom_rep.
set.seed(4) n <- 40 df <- data.frame( count = rpois(n, 12), total = rpois(n, 80), x = rnorm(n) ) fit <- tratio(count | total ~ x, data = df, family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 1, seed = 1)) rep <- posterior_predict(fit, ndraws = 50) dim(rep$y_num_rep)set.seed(4) n <- 40 df <- data.frame( count = rpois(n, 12), total = rpois(n, 80), x = rnorm(n) ) fit <- tratio(count | total ~ x, data = df, family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 1, seed = 1)) rep <- posterior_predict(fit, ndraws = 50) dim(rep$y_num_rep)
Generate posterior predictive checks for a fitted ratio model. Compares observed data to replicated data from the posterior predictive distribution.
## S3 method for class 'ratiod_fit' pp_check( object, type = c("dens_overlay", "scatter", "intervals", "stat"), component = c("numerator", "denominator", "both"), stat = mean, ndraws = 50, ... )## S3 method for class 'ratiod_fit' pp_check( object, type = c("dens_overlay", "scatter", "intervals", "stat"), component = c("numerator", "denominator", "both"), stat = mean, ndraws = 50, ... )
object |
A |
type |
Type of check: "dens_overlay", "scatter", "intervals", "stat" |
component |
Which component: "numerator", "denominator", or "both" |
stat |
Function for "stat" type (default: mean) |
ndraws |
Number of posterior draws to replicate (default: 50) |
... |
Additional arguments passed to bayesplot functions |
A ggplot object
# Simulate data and fit model set.seed(123) n <- 50 df <- data.frame( count = rnbinom(n, size = 5, mu = 15), total = rnbinom(n, size = 8, mu = 100), x = rnorm(n), site = factor(rep(1:10, each = 5)) ) fit <- tratio( count | total ~ x + (1 | site), data = df, family = ratiod_negbin_negbin(), control = list(iter = 400, chains = 1, seed = 1) ) # Density overlay of replicated against observed numerators pp_check(fit, type = "dens_overlay")# Simulate data and fit model set.seed(123) n <- 50 df <- data.frame( count = rnbinom(n, size = 5, mu = 15), total = rnbinom(n, size = 8, mu = 100), x = rnorm(n), site = factor(rep(1:10, each = 5)) ) fit <- tratio( count | total ~ x + (1 | site), data = df, family = ratiod_negbin_negbin(), control = list(iter = 400, chains = 1, seed = 1) ) # Density overlay of replicated against observed numerators pp_check(fit, type = "dens_overlay")
Compute posterior predictions for new data from a fitted ratio model.
## S3 method for class 'ratiod_fit' predict( object, newdata, type = c("response", "link"), component = c("ratio", "numerator", "denominator", "all"), summary = TRUE, probs = c(0.025, 0.5, 0.975), re_formula = NULL, allow_new_levels = FALSE, coords.0 = NULL, include_spatial = TRUE, return_spatial = FALSE, n_samples = NULL, ... ) ## S3 method for class 'ratiod_prediction' print(x, n = 10, ...)## S3 method for class 'ratiod_fit' predict( object, newdata, type = c("response", "link"), component = c("ratio", "numerator", "denominator", "all"), summary = TRUE, probs = c(0.025, 0.5, 0.975), re_formula = NULL, allow_new_levels = FALSE, coords.0 = NULL, include_spatial = TRUE, return_spatial = FALSE, n_samples = NULL, ... ) ## S3 method for class 'ratiod_prediction' print(x, n = 10, ...)
object |
A |
newdata |
Data frame with new observations. Must contain all predictor variables used in the original model. |
type |
Type of prediction: "response" for fitted values on the response scale, "link" for linear predictor scale. |
component |
Which component to predict: "ratio" (default), "numerator", "denominator", or "all". |
summary |
Logical; if TRUE, return summary statistics instead of full posterior draws (default TRUE). |
probs |
Quantiles to compute if |
re_formula |
Formula for random effects in prediction.
Use |
allow_new_levels |
Logical; if TRUE, allows new group levels not in the original data (predictions use population mean). |
coords.0 |
Matrix of spatial coordinates for new prediction locations (required for GP/HSGP spatial models). Should be N x 2 matrix with same coordinate system as training data. |
include_spatial |
Logical; if TRUE (default), include spatial random effects in predictions. Set to FALSE for fixed-effects-only predictions. |
return_spatial |
Logical; if TRUE, return spatial effect samples at
new locations as |
n_samples |
Number of posterior samples to draw for backends that don't store samples directly (e.g., Laplace, VI). Default uses all available samples for HMC/PG. |
... |
Ignored |
x |
A ratiod_prediction object |
n |
Number of rows to print |
If summary = TRUE, a data frame with posterior summaries.
If summary = FALSE, a matrix of posterior draws (rows = draws, cols = obs).
If return_spatial = TRUE, returns a list with predictions and w.0.samples.
# A smooth spatial field w shared by count and effort set.seed(321) n <- 60 x <- runif(n) y <- runif(n) w <- 0.6 * sin(2 * pi * x) + 0.6 * cos(2 * pi * y) depth <- rnorm(n) season <- factor(sample(c("summer", "winter"), n, replace = TRUE)) df <- data.frame( count = rpois(n, exp(2.5 + 0.3 * depth + 0.4 * (season == "winter") + w)), effort = rgamma(n, shape = 9, rate = exp(-w)), depth = depth, season = season, x = x, y = y ) fit <- tratio( count | effort ~ depth + season, data = df, family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) # Predict for new data; levels of `season` follow the training data new_df <- data.frame(depth = c(-1, 0, 1), season = "summer") predict(fit, newdata = new_df) # Full posterior draws of the predicted ratio (draws x rows) dim(predict(fit, newdata = new_df, summary = FALSE)$ratio) # Spatial prediction at new locations (kriging) from a GP fit fit_gp <- tratio( count | effort ~ depth, data = df, spatial = spatial_gp(coords = c("x", "y")), family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) new_coords <- matrix(runif(10), ncol = 2) predict(fit_gp, newdata = data.frame(depth = rnorm(5)), coords.0 = new_coords)# A smooth spatial field w shared by count and effort set.seed(321) n <- 60 x <- runif(n) y <- runif(n) w <- 0.6 * sin(2 * pi * x) + 0.6 * cos(2 * pi * y) depth <- rnorm(n) season <- factor(sample(c("summer", "winter"), n, replace = TRUE)) df <- data.frame( count = rpois(n, exp(2.5 + 0.3 * depth + 0.4 * (season == "winter") + w)), effort = rgamma(n, shape = 9, rate = exp(-w)), depth = depth, season = season, x = x, y = y ) fit <- tratio( count | effort ~ depth + season, data = df, family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) # Predict for new data; levels of `season` follow the training data new_df <- data.frame(depth = c(-1, 0, 1), season = "summer") predict(fit, newdata = new_df) # Full posterior draws of the predicted ratio (draws x rows) dim(predict(fit, newdata = new_df, summary = FALSE)$ratio) # Spatial prediction at new locations (kriging) from a GP fit fit_gp <- tratio( count | effort ~ depth, data = df, spatial = spatial_gp(coords = c("x", "y")), family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) new_coords <- matrix(runif(10), ncol = 2) predict(fit_gp, newdata = data.frame(depth = rnorm(5)), coords.0 = new_coords)
print() gives a short description of a fit: how it was inferred, the
family and the data. summary() returns the full report: posterior mean,
sd and credible interval of every main parameter with its convergence
diagnostics, which its print() method lays out. coef() returns the
posterior means of the fixed effects and vcov() their posterior
covariance; coef(summary(fit)) returns the summary table itself.
nobs() and formula() return the number of observations and the model
formula.
## S3 method for class 'ratiod_sghmc' print(x, digits = 3, ...) ## S3 method for class 'ratiod_sgld' print(x, digits = 3, ...) ## S3 method for class 'ratiod_vi' print(x, digits = 3, ...) ## S3 method for class 'ratiod_vi' summary(object, prob = 0.95, ...) ## S3 method for class 'summary.ratiod_vi' print(x, digits = 3, ...) ## S3 method for class 'ratiod_fit' print(x, ...) ## S3 method for class 'ratiod_fit' summary(object, prob = 0.95, ...) ## S3 method for class 'summary.ratiod_fit' print(x, ...) ## S3 method for class 'summary.ratiod_fit' coef(object, ...) ## S3 method for class 'ratiod_fit' coef(object, ...) ## S3 method for class 'ratiod_fit' vcov(object, ...) ## S3 method for class 'ratiod_fit' nobs(object, ...) ## S3 method for class 'ratiod_fit' formula(x, ...) ## S3 method for class 'ratiod_fit' plot(x, pars = NULL, type = c("trace", "dens", "both"), n_col = NULL, ...)## S3 method for class 'ratiod_sghmc' print(x, digits = 3, ...) ## S3 method for class 'ratiod_sgld' print(x, digits = 3, ...) ## S3 method for class 'ratiod_vi' print(x, digits = 3, ...) ## S3 method for class 'ratiod_vi' summary(object, prob = 0.95, ...) ## S3 method for class 'summary.ratiod_vi' print(x, digits = 3, ...) ## S3 method for class 'ratiod_fit' print(x, ...) ## S3 method for class 'ratiod_fit' summary(object, prob = 0.95, ...) ## S3 method for class 'summary.ratiod_fit' print(x, ...) ## S3 method for class 'summary.ratiod_fit' coef(object, ...) ## S3 method for class 'ratiod_fit' coef(object, ...) ## S3 method for class 'ratiod_fit' vcov(object, ...) ## S3 method for class 'ratiod_fit' nobs(object, ...) ## S3 method for class 'ratiod_fit' formula(x, ...) ## S3 method for class 'ratiod_fit' plot(x, pars = NULL, type = c("trace", "dens", "both"), n_col = NULL, ...)
x |
A fit returned by |
digits |
Number of significant digits to print (default 3). |
... |
Ignored |
object |
A ratiod_fit object |
prob |
Probability mass for credible intervals (default 0.95) |
pars |
Character vector of parameter names to plot. If NULL (default), plots the first 12 main parameters (fixed effects, variance components). |
type |
Type of plot: "trace" (default), "dens" (density), or "both". |
n_col |
Number of columns in the plot grid (default: auto). |
An object of class summary.ratiod_vi.
The input x, invisibly.
summary(): a summary.ratiod_fit object. coef(): a named
vector, numerator coefficients prefixed num: and denominator ones
denom:. vcov(): the matching covariance matrix.
Invisibly returns the plot object (base graphics or ggplot2).
# Create simple dataset set.seed(456) n <- 48 site <- factor(rep(1:6, each = 8)) u <- rnorm(6, sd = 0.5)[site] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(2 + 0.3 * depth + u)), effort = rgamma(n, shape = 2, rate = 0.5 * exp(-u)), depth = depth, site = site ) fit <- tratio( count | effort ~ depth + (1|site), data = df, family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 2, seed = 1) ) plot(fit) plot(fit, pars = "beta_num", type = "both")# Create simple dataset set.seed(456) n <- 48 site <- factor(rep(1:6, each = 8)) u <- rnorm(6, sd = 0.5)[site] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(2 + 0.3 * depth + u)), effort = rgamma(n, shape = 2, rate = 0.5 * exp(-u)), depth = depth, site = site ) fit <- tratio( count | effort ~ depth + (1|site), data = df, family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 2, seed = 1) ) plot(fit) plot(fit, pars = "beta_num", type = "both")
Specify a beta prior for a parameter bounded in (0, 1).
prior_beta(alpha = 1, beta = 1)prior_beta(alpha = 1, beta = 1)
alpha |
First shape parameter. Must be positive. |
beta |
Second shape parameter. Must be positive. |
alpha = beta = 1: Uniform(0, 1)
alpha = beta = 2: Symmetric, peaked at 0.5
alpha > beta: Skewed toward 1
alpha < beta: Skewed toward 0
A ratiod_prior object
prior_beta(1, 1) # Uniform prior_beta(2, 2) # Symmetric, centered at 0.5 prior_beta(5, 2) # Skewed toward 1 (for high autocorrelation)prior_beta(1, 1) # Uniform prior_beta(2, 2) # Symmetric, centered at 0.5 prior_beta(5, 2) # Skewed toward 1 (for high autocorrelation)
Specify an exponential prior for a positive parameter.
prior_exponential(rate = 1)prior_exponential(rate = 1)
rate |
Rate parameter (lambda). Must be positive. |
A ratiod_prior object
prior_exponential(1)prior_exponential(1)
Specify a gamma prior for a positive parameter.
prior_gamma(shape = 2, rate = 0.1)prior_gamma(shape = 2, rate = 0.1)
shape |
Shape parameter (alpha). Must be positive. |
rate |
Rate parameter (beta). Must be positive. |
Mean = shape/rate, Variance = shape/rate^2.
A ratiod_prior object
prior_gamma(2, 0.1) # Mean = 20, weakly informative prior_gamma(1, 1) # Exponential(1)prior_gamma(2, 0.1) # Mean = 20, weakly informative prior_gamma(1, 1) # Exponential(1)
Specify a half-Cauchy prior for a positive parameter. Has heavier tails than half-normal, often used for variance parameters.
prior_half_cauchy(scale = 2.5)prior_half_cauchy(scale = 2.5)
scale |
Scale parameter. Must be positive. |
A ratiod_prior object
prior_half_cauchy(2.5)prior_half_cauchy(2.5)
Specify a half-normal (truncated at 0) prior for a positive parameter.
prior_half_normal(sd = 1)prior_half_normal(sd = 1)
sd |
Scale parameter. Must be positive. |
A ratiod_prior object
prior_half_normal(1)prior_half_normal(1)
Specify a normal prior for a parameter.
prior_normal(mean = 0, sd = 2.5)prior_normal(mean = 0, sd = 2.5)
mean |
Prior mean. Default 0. |
sd |
Prior standard deviation. Must be positive. |
A ratiod_prior object
prior_normal(0, 2.5) prior_normal(0, 1)prior_normal(0, 2.5) prior_normal(0, 1)
Specify a PC prior for a positive parameter (typically a standard deviation). PC priors shrink toward simpler models by penalizing deviation from a base model.
prior_pc(U = 1, alpha = 0.01)prior_pc(U = 1, alpha = 0.01)
U |
Upper bound. P(x > U) = alpha. |
alpha |
Tail probability. Default 0.01. |
The PC prior is specified via: P(sigma > U) = alpha
U is the upper bound you consider "large"
alpha is the probability of exceeding it
This implies an exponential prior with rate = -log(alpha) / U.
A ratiod_prior object
Simpson, D., Rue, H., Riebler, A., Martins, T. G., & Sorbye, S. H. (2017). Penalising model component complexity: A principled, practical approach to constructing priors. Statistical Science, 32(1), 1-28.
prior_pc(U = 1, alpha = 0.01) # P(sigma > 1) = 0.01 prior_pc(U = 0.5, alpha = 0.05) # Tighter, P(sigma > 0.5) = 0.05prior_pc(U = 1, alpha = 0.01) # P(sigma > 1) = 0.01 prior_pc(U = 0.5, alpha = 0.05) # Tighter, P(sigma > 0.5) = 0.05
Simulate numerator and denominator responses from the prior predictive
distribution: parameters are drawn from priors, the two linear predictors
are assembled on the design formula and data imply, and a response is
drawn from family at each draw. Comparing the resulting ratios against
what is plausible for the system under study is the check on whether the
priors are reasonable before any data are fitted.
Random-effect terms contribute at a standard deviation drawn from
priors$sigma, one per intercept and one per slope column. Correlation
between slopes within a term is not drawn, so a correlated term simulates
at its marginal scales.
prior_predict(formula, family, data, priors = NULL, n = 100) ## S3 method for class 'ratiod_priorpred' print(x, ...)prior_predict(formula, family, data, priors = NULL, n = 100) ## S3 method for class 'ratiod_priorpred' print(x, ...)
formula |
A model formula, using the combined |
family |
A |
data |
Data frame supplying the predictors, the grouping factors and, for the binomial families, the trial counts. Response columns need not be present; where one is missing a placeholder is used, since nothing about the prior predictive distribution depends on the observed response. |
priors |
A |
n |
Number of prior predictive draws. |
x |
A |
... |
Ignored |
An object of class ratiod_priorpred: a list with the n x N
matrices y_num, y_denom, ratio, mu_num and mu_denom, the drawn
coefficients beta_num and beta_denom, the drawn random-effect scales
sigma, the drawn dispersions phi_num and phi_denom, and the family
and data the draws were taken on.
x, invisibly
sim_ratiod() for simulating from a known parameterisation,
ratiod_priors() for the prior specification
set.seed(1) df <- data.frame( count = rep(0L, 40), total = rep(0L, 40), x = rnorm(40), site = factor(rep(1:8, each = 5)) ) pp <- prior_predict( count | total ~ x + (1 | site), family = ratiod_negbin_negbin(), data = df, n = 50 ) dim(pp$ratio) summary(as.vector(pp$ratio)) # A tighter coefficient prior narrows the predicted ratios pp_tight <- prior_predict( count | total ~ x + (1 | site), family = ratiod_negbin_negbin(), data = df, priors = ratiod_priors(beta = prior_normal(0, 0.5)), n = 50 )set.seed(1) df <- data.frame( count = rep(0L, 40), total = rep(0L, 40), x = rnorm(40), site = factor(rep(1:8, each = 5)) ) pp <- prior_predict( count | total ~ x + (1 | site), family = ratiod_negbin_negbin(), data = df, n = 50 ) dim(pp$ratio) summary(as.vector(pp$ratio)) # A tighter coefficient prior narrows the predicted ratios pp_tight <- prior_predict( count | total ~ x + (1 | site), family = ratiod_negbin_negbin(), data = df, priors = ratiod_priors(beta = prior_normal(0, 0.5)), n = 50 )
Display the priors a fit uses for every component ratiod_priors() leaves
unspecified, as a starting point for customization.
priors_default(family = NULL, spatial = FALSE, temporal = FALSE)priors_default(family = NULL, spatial = FALSE, temporal = FALSE)
family |
A tulpaRatio family object (e.g., |
spatial |
Logical; if TRUE, include the BYM2 mixing prior. |
temporal |
Logical; if TRUE, include the AR1 correlation prior. |
Fixed effects (beta): Normal(0, 10) on the link scale.
Random effect SD (sigma): Half-Cauchy(2.5); the Laplace backend uses a PC prior with P(sigma > 1) = 0.01.
Dispersion (phi): Gamma(2, 2) for the lognormal log-scale SD, Gamma(2, 0.5) for a Gamma shape (gamma_gamma, poisson_gamma), Gamma(2, 0.1) for a negative binomial size or beta-binomial precision.
Temporal correlation (rho): Beta(2, 2) on (rho + 1) / 2, symmetric
about rho = 0.
Spatial mixing (rho_spatial): Beta(1, 1), uniform on (0, 1).
Invisibly, a ratiod_priors object holding the defaults (for
family = NULL, with the negative binomial dispersion default).
ratiod_priors() for creating custom priors
priors_default(ratiod_negbin_negbin()) priors_default(ratiod_poisson_gamma(), spatial = TRUE, temporal = TRUE) priors_default()priors_default(ratiod_negbin_negbin()) priors_default(ratiod_poisson_gamma(), spatial = TRUE, temporal = TRUE) priors_default()
Compute posterior distributions of ratios from a fitted ratio model. Ratios are derived quantities: E[numerator] / E[denominator], computed from each posterior draw, so they carry the full posterior uncertainty.
ratio( object, newdata = NULL, type = c("response", "log", "logit"), by = NULL, summary = FALSE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_ratio' summary(object, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_ratio' print(x, ...) ## S3 method for class 'ratiod_ratio_summary' print(x, n = 10, ...)ratio( object, newdata = NULL, type = c("response", "log", "logit"), by = NULL, summary = FALSE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_ratio' summary(object, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_ratio' print(x, ...) ## S3 method for class 'ratiod_ratio_summary' print(x, n = 10, ...)
object |
A |
newdata |
Optional data frame for prediction. If NULL, uses original data. |
type |
Scale for ratio: "response" (default), "log", or "logit". |
by |
Optional grouping variable name to aggregate ratios by group. |
summary |
Logical; if TRUE, return summary statistics instead of full posterior draws. |
probs |
Quantiles to compute if |
... |
Ignored |
x |
A |
n |
Number of rows to print |
A ratiod_ratio object containing posterior draws or summaries.
# Create a simple family fam <- ratiod_poisson_gamma() # Generate synthetic data: a site effect shared by count and effort set.seed(123) n <- 60 site <- rep(letters[1:6], each = 10) u <- rnorm(6, sd = 0.5)[match(site, letters[1:6])] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(3 + 0.3 * depth + u)), effort = rgamma(n, shape = 10, rate = exp(-u)), depth = depth, site = site ) # Fit model fit <- tratio( count | effort ~ depth + (1 | site), data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Extract all observation-level ratios r <- ratio(fit) summary(r) # Log-ratio scale r_log <- ratio(fit, type = "log") # Group-level ratios r_site <- ratio(fit, by = "site")# Create a simple family fam <- ratiod_poisson_gamma() # Generate synthetic data: a site effect shared by count and effort set.seed(123) n <- 60 site <- rep(letters[1:6], each = 10) u <- rnorm(6, sd = 0.5)[match(site, letters[1:6])] depth <- rnorm(n) df <- data.frame( count = rpois(n, exp(3 + 0.3 * depth + u)), effort = rgamma(n, shape = 10, rate = exp(-u)), depth = depth, site = site ) # Fit model fit <- tratio( count | effort ~ depth + (1 | site), data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Extract all observation-level ratios r <- ratio(fit) summary(r) # Log-ratio scale r_log <- ratio(fit, type = "log") # Group-level ratios r_site <- ratio(fit, by = "site")
Compute posterior distributions of contrasts (differences or ratios) between conditions for the derived ratio quantity.
ratio_contrast( object, contrast, type = c("difference", "ratio"), ref = NULL, ... ) ## S3 method for class 'ratiod_contrast' print(x, ...)ratio_contrast( object, contrast, type = c("difference", "ratio"), ref = NULL, ... ) ## S3 method for class 'ratiod_contrast' print(x, ...)
object |
A |
contrast |
A formula specifying the contrast, e.g., |
type |
Type of contrast: "difference" (default) or "ratio" |
ref |
Reference level for factor contrasts (default: first level) |
... |
Ignored |
x |
A ratiod_contrast object |
A data frame with posterior summaries of contrasts
# Generate synthetic data with grouping variable set.seed(456) n <- 60 df <- data.frame( count = rpois(n, lambda = 25), effort = rgamma(n, shape = 12, rate = 1), season = sample(c("summer", "winter"), n, replace = TRUE) ) # Fit model fit <- tratio( count | effort ~ season, data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Compare CPUE between seasons ratio_contrast(fit, ~ season) # Ratio of ratios (multiplicative comparison) ratio_contrast(fit, ~ season, type = "ratio")# Generate synthetic data with grouping variable set.seed(456) n <- 60 df <- data.frame( count = rpois(n, lambda = 25), effort = rgamma(n, shape = 12, rate = 1), season = sample(c("summer", "winter"), n, replace = TRUE) ) # Fit model fit <- tratio( count | effort ~ season, data = df, family = ratiod_poisson_gamma(), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Compare CPUE between seasons ratio_contrast(fit, ~ season) # Ratio of ratios (multiplicative comparison) ratio_contrast(fit, ~ season, type = "ratio")
Compute model-averaged predictions using stacking or pseudo-BMA weights. Model weights are derived from LOO-CV or WAIC, and predictions are combined accordingly.
ratiod_average( ..., weights = c("loo", "waic", "pbma", "pbma+"), newdata = NULL, type = c("ratio", "numerator", "denominator"), summary = TRUE ) ## S3 method for class 'ratiod_average' print(x, digits = 3, ...) ## S3 method for class 'ratiod_average' fitted(object, ...) ## S3 method for class 'ratiod_average' weights(object, ...)ratiod_average( ..., weights = c("loo", "waic", "pbma", "pbma+"), newdata = NULL, type = c("ratio", "numerator", "denominator"), summary = TRUE ) ## S3 method for class 'ratiod_average' print(x, digits = 3, ...) ## S3 method for class 'ratiod_average' fitted(object, ...) ## S3 method for class 'ratiod_average' weights(object, ...)
... |
Multiple |
weights |
Method for computing weights: "loo" (stacking, default), "waic", "pbma" (pseudo-BMA), or "pbma+" (pseudo-BMA+ with Bayesian bootstrap). |
newdata |
Optional new data for predictions. If NULL, uses fitted values. |
type |
Type of prediction: "ratio" (default), "numerator", or "denominator". |
summary |
Logical; if TRUE (default), return summary statistics. If FALSE, return full posterior draws. |
x |
A ratiod_average object |
digits |
Number of digits for printing |
object |
A ratiod_average object |
Model averaging accounts for model uncertainty by combining predictions from multiple candidate models. The weighting methods are:
stacking (weights = "loo"): Optimal linear combination minimizing
leave-one-out prediction error. Recommended default.
pseudo-BMA (weights = "pbma"): Weights proportional to exp(-ELPD).
Can be unstable with similar models.
pseudo-BMA+ (weights = "pbma+"): Bayesian bootstrap variant that
accounts for estimation uncertainty.
For ratio predictions, averaging is performed on the log scale to respect the multiplicative nature of ratios.
A ratiod_average object containing:
weights: Model weights
predictions: Model-averaged predictions
models: List of model names
Predictions data frame or matrix
Named numeric vector of weights
ratiod_compare() for model comparison, loo::stacking_weights()
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) fit0 <- tratio(y_num | y_denom ~ 1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) avg <- ratiod_average(fit0, fit) avg weights(avg)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) fit0 <- tratio(y_num | y_denom ~ 1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) avg <- ratiod_average(fit0, fit) avg weights(avg)
Utilities for benchmarking ratio model fitting performance across different dataset sizes and configurations.
Fits a model on simulated data and reports timing and performance metrics.
ratiod_benchmark( N = 500, p = 3, n_groups = 0, family = ratiod_negbin_negbin(), n_iter = 1000, n_warmup = 500, n_chains = 1, n_threads = 1, verbose = FALSE ) ## S3 method for class 'ratiod_benchmark' print(x, ...)ratiod_benchmark( N = 500, p = 3, n_groups = 0, family = ratiod_negbin_negbin(), n_iter = 1000, n_warmup = 500, n_chains = 1, n_threads = 1, verbose = FALSE ) ## S3 method for class 'ratiod_benchmark' print(x, ...)
N |
Number of observations |
p |
Number of predictors |
n_groups |
Number of random effect groups (0 for none) |
family |
Model family (default: ratiod_negbin_negbin()) |
n_iter |
Total iterations |
n_warmup |
Warmup iterations |
n_chains |
Number of chains |
n_threads |
Number of threads per chain |
verbose |
Print progress (default: FALSE) |
x |
A ratiod_benchmark object |
... |
Ignored |
A list with timing and diagnostic information
# Quick benchmark bench <- ratiod_benchmark(N = 100, p = 3, n_iter = 200, n_warmup = 100) print(bench)# Quick benchmark bench <- ratiod_benchmark(N = 100, p = 3, n_iter = 200, n_warmup = 100) print(bench)
Run benchmarks across multiple configurations to compare performance.
ratiod_benchmark_compare(configs, n_iter = 500, n_warmup = 250, n_threads = 1)ratiod_benchmark_compare(configs, n_iter = 500, n_warmup = 250, n_threads = 1)
configs |
A list of configuration lists, each with N, p, n_groups |
n_iter |
Iterations per benchmark |
n_warmup |
Warmup iterations |
n_threads |
Number of threads |
A data frame with benchmark results
# Compare benchmarks configs <- list( list(N = 50, p = 2, n_groups = 0), list(N = 100, p = 2, n_groups = 0) ) comparison <- ratiod_benchmark_compare(configs, n_iter = 200, n_warmup = 100) print(comparison)# Compare benchmarks configs <- list( list(N = 50, p = 2, n_groups = 0), list(N = 100, p = 2, n_groups = 0) ) comparison <- ratiod_benchmark_compare(configs, n_iter = 200, n_warmup = 100) print(comparison)
Two-process model where the numerator follows a beta-binomial distribution, allowing for overdispersion in binomial data. The denominator is the known number of trials.
Use cases:
Overdispersed proportions (e.g., infection rates across sites)
Hierarchical binomial with extra-binomial variation
Proportions with correlation within clusters
ratiod_beta_binomial(link = "logit")ratiod_beta_binomial(link = "logit")
link |
Link function for the success probability: "logit" (default), "probit" or "cloglog". |
The beta-binomial distribution arises when the success probability itself follows a beta distribution. This creates overdispersion (variance > binomial variance) and is useful when:
There is unmodelled heterogeneity in success probability
Observations within clusters are correlated
The binomial assumption of independent trials is violated
The overdispersion parameter (intra-class correlation)
ranges from 0 (binomial) to 1 (maximum overdispersion).
A ratiod_family object
# Create family object fam <- ratiod_beta_binomial() print(fam) # Simulate overdispersed proportion data set.seed(123) n <- 40 df <- data.frame( infected = rbinom(n, size = 20, prob = rbeta(n, 2, 5)), tested = rep(20L, n), treatment = factor(rep(c("control", "treated"), each = n/2)), site = factor(rep(1:8, each = n/8)) ) # Fit model fit <- tratio( infected | tested ~ treatment + (1 | site), data = df, family = ratiod_beta_binomial(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create family object fam <- ratiod_beta_binomial() print(fam) # Simulate overdispersed proportion data set.seed(123) n <- 40 df <- data.frame( infected = rbinom(n, size = 20, prob = rbeta(n, 2, 5)), tested = rep(20L, n), treatment = factor(rep(c("control", "treated"), each = n/2)), site = factor(rep(1:8, each = n/8)) ) # Fit model fit <- tratio( infected | tested ~ treatment + (1 | site), data = df, family = ratiod_beta_binomial(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Trial-based model where the numerator represents successes out of a known number of trials, the denominator.
Use cases:
Detection / availability (occupancy-style)
Successes / trials (hierarchical binomial)
Proportions with known denominators
ratiod_binomial(link = "logit")ratiod_binomial(link = "logit")
link |
Link function for the success probability: "logit" (default), "probit" or "cloglog". |
A ratiod_family object
# Create family object fam <- ratiod_binomial() print(fam) # Simulate detection/availability data set.seed(123) n <- 50 df <- data.frame( detections = rbinom(n, size = 10, prob = 0.4), trials = rep(10L, n), effort = rnorm(n), site = factor(rep(1:10, each = 5)) ) # Fit model fit <- tratio( detections | trials ~ effort + (1 | site), data = df, family = ratiod_binomial(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create family object fam <- ratiod_binomial() print(fam) # Simulate detection/availability data set.seed(123) n <- 50 df <- data.frame( detections = rbinom(n, size = 10, prob = 0.4), trials = rep(10L, n), effort = rnorm(n), site = factor(rep(1:10, each = 5)) ) # Fit model fit <- tratio( detections | trials ~ effort + (1 | site), data = df, family = ratiod_binomial(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Compare multiple ratio models using LOO-CV or WAIC.
ratiod_compare(..., criterion = c("loo", "waic"))ratiod_compare(..., criterion = c("loo", "waic"))
... |
Multiple |
criterion |
Comparison criterion: "loo" (default) or "waic" |
A loo_compare object
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) fit0 <- tratio(y_num | y_denom ~ 1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) ratiod_compare(fit0, fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) fit0 <- tratio(y_num | y_denom ~ 1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) ratiod_compare(fit0, fit)
Distribution families for the numerator and denominator processes. Each family specifies the likelihood structure for joint modelling.
## S3 method for class 'ratiod_family' print(x, ...)## S3 method for class 'ratiod_family' print(x, ...)
x |
A ratiod_family object |
... |
Ignored |
ratiod_poisson_gamma() ratiod_binomial()ratiod_poisson_gamma() ratiod_binomial()
Two-process model where both numerator and denominator follow gamma distributions. This is the natural family for ratios of positive continuous quantities.
Use cases:
Biomass ratios (species biomass / total biomass)
Concentration ratios (analyte / reference)
Duration ratios (time in activity / total time)
ratiod_gamma_gamma(link_num = "log", link_denom = "log")ratiod_gamma_gamma(link_num = "log", link_denom = "log")
link_num |
Link function for the numerator mean: "log" (default), "identity" or "inverse". |
link_denom |
Link function for the denominator mean: "log" (default), "identity" or "inverse". |
The gamma distribution is appropriate for positive continuous data with right skew. The ratio of two gamma random variables does not have a simple closed form, which is why ratio models them jointly.
Shape parameters are estimated for both numerator and denominator, allowing different amounts of variability in each process.
A ratiod_family object
# Create family object fam <- ratiod_gamma_gamma() print(fam) # Simulate biomass ratio data set.seed(123) n <- 50 df <- data.frame( species_biomass = rgamma(n, shape = 3, rate = 0.5), total_biomass = rgamma(n, shape = 10, rate = 0.3), habitat = factor(rep(c("forest", "grassland"), each = n/2)), site = factor(rep(1:10, each = n/10)) ) # Fit model fit <- tratio( species_biomass | total_biomass ~ habitat + (1 | site), data = df, family = ratiod_gamma_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create family object fam <- ratiod_gamma_gamma() print(fam) # Simulate biomass ratio data set.seed(123) n <- 50 df <- data.frame( species_biomass = rgamma(n, shape = 3, rate = 0.5), total_biomass = rgamma(n, shape = 10, rate = 0.3), habitat = factor(rep(c("forest", "grassland"), each = n/2)), site = factor(rep(1:10, each = n/10)) ) # Fit model fit <- tratio( species_biomass | total_biomass ~ habitat + (1 | site), data = df, family = ratiod_gamma_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Hurdle model for binomial data where zero counts are modeled separately from positive counts using a truncated binomial.
ratiod_hurdle_binomial(link_num = "logit", link_hurdle = "logit")ratiod_hurdle_binomial(link_num = "logit", link_hurdle = "logit")
link_num |
Link function for the success probability: "logit" (default), "probit" or "cloglog". |
link_hurdle |
Link function for the hurdle probability P(Y > 0): "logit" (default), "probit" or "cloglog". |
A ratiod_family object with hurdle structure
fam <- ratiod_hurdle_binomial() print(fam)fam <- ratiod_hurdle_binomial() print(fam)
Hurdle model where zero vs positive counts are modelled as separate processes.
where is the probability of a positive count and
is the truncated-at-zero negative binomial.
ratiod_hurdle_negbin( link_num = "log", link_denom = "log", link_hurdle = "logit", denom_family = "negbin" )ratiod_hurdle_negbin( link_num = "log", link_denom = "log", link_hurdle = "logit", denom_family = "negbin" )
link_num |
Link function for the count mean: "log". |
link_denom |
Link function for the denominator mean: "log" (default); a gamma denominator also takes "identity" or "inverse". |
link_hurdle |
Link function for the hurdle probability: "logit" (default), "probit" or "cloglog". |
denom_family |
Denominator distribution: "negbin" (default), "gamma" for continuous effort, or "fixed" for a known exposure. A known exposure is data: it scales the count mean, which is then the rate times the exposure, and the ratio is the rate. |
The hurdle model differs from zero-inflation in interpretation:
Zero-inflated: Zeros come from two sources (structural + sampling)
: probability of structural zero
Some zeros also come from the count process
Hurdle: Zeros are all from one source (the "hurdle" process)
: probability of zero
All positive counts from truncated count process
Choose hurdle when zeros represent a distinct biological/sampling state.
A ratiod_family object with hurdle structure
# Create family object fam <- ratiod_hurdle_negbin() print(fam) # Simulate hurdle data (presence/absence + abundance) set.seed(123) n <- 60 presence_prob <- 0.7 df <- data.frame( count = ifelse(runif(n) > presence_prob, 0, rnbinom(n, size = 3, mu = 10) + 1), total = rnbinom(n, size = 5, mu = 50), habitat = factor(rep(c("forest", "grassland"), each = n/2)), site = factor(rep(1:10, each = n/10)) ) # Fit model fit <- tratio( count | total ~ habitat + (1 | site), data = df, family = ratiod_hurdle_negbin(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create family object fam <- ratiod_hurdle_negbin() print(fam) # Simulate hurdle data (presence/absence + abundance) set.seed(123) n <- 60 presence_prob <- 0.7 df <- data.frame( count = ifelse(runif(n) > presence_prob, 0, rnbinom(n, size = 3, mu = 10) + 1), total = rnbinom(n, size = 5, mu = 50), habitat = factor(rep(c("forest", "grassland"), each = n/2)), site = factor(rep(1:10, each = n/10)) ) # Fit model fit <- tratio( count | total ~ habitat + (1 | site), data = df, family = ratiod_hurdle_negbin(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Hurdle model with truncated Poisson for positive counts.
ratiod_hurdle_pois( link_num = "log", link_denom = "log", link_hurdle = "logit", denom_family = "gamma" )ratiod_hurdle_pois( link_num = "log", link_denom = "log", link_hurdle = "logit", denom_family = "gamma" )
link_num |
Link function for the count mean: "log". |
link_denom |
Link function for the denominator mean: "log" (default); a gamma denominator also takes "identity" or "inverse". |
link_hurdle |
Link function for the hurdle probability: "logit" (default), "probit" or "cloglog". |
denom_family |
Denominator distribution: "gamma" (default), "negbin" for a count denominator, or "fixed" for a known exposure. |
A ratiod_family object with hurdle structure
ratiod_hurdle_pois()ratiod_hurdle_pois()
Specify latent factors to capture shared unmeasured confounders between numerator and denominator processes. Latent factors are particularly useful when you suspect that both processes are driven by common unmeasured variables.
latent_factor(n_factors = 2)latent_factor(n_factors = 2)
Two-process model where both numerator and denominator follow log-normal distributions. Appropriate when the underlying process is multiplicative and residuals are log-normally distributed.
Use cases:
Abundance indices (CPUE with continuous effort)
Economic ratios (costs, prices)
Body condition indices
ratiod_lognormal(link_num = "log", link_denom = "log", denom_fixed = FALSE)ratiod_lognormal(link_num = "log", link_denom = "log", denom_fixed = FALSE)
link_num |
Link function for the numerator median: "log". |
link_denom |
Link function for the denominator median: "log". |
denom_fixed |
If |
The log-normal distribution arises when a variable is the product of
many independent positive factors. If where is
normal, then is log-normal.
For ratios, if both numerator and denominator are log-normal, the ratio is also log-normal (difference of normals on log scale). This can simplify interpretation when working on the log scale.
A ratiod_family object
# Create family object fam <- ratiod_lognormal() print(fam) # Simulate body condition data set.seed(123) n <- 60 df <- data.frame( weight = rlnorm(n, meanlog = 3, sdlog = 0.3), length_cubed = rlnorm(n, meanlog = 2, sdlog = 0.2), age = sample(1:5, n, replace = TRUE), sex = factor(rep(c("M", "F"), each = n/2)), cohort = factor(rep(1:6, each = n/6)) ) # Fit model fit <- tratio( weight | length_cubed ~ age + sex + (1 | cohort), data = df, family = ratiod_lognormal(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create family object fam <- ratiod_lognormal() print(fam) # Simulate body condition data set.seed(123) n <- 60 df <- data.frame( weight = rlnorm(n, meanlog = 3, sdlog = 0.3), length_cubed = rlnorm(n, meanlog = 2, sdlog = 0.2), age = sample(1:5, n, replace = TRUE), sex = factor(rep(c("M", "F"), each = n/2)), cohort = factor(rep(1:6, each = n/6)) ) # Fit model fit <- tratio( weight | length_cubed ~ age + sex + (1 | cohort), data = df, family = ratiod_lognormal(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Two-process model where the numerator follows a negative binomial distribution and the denominator (effort/exposure) follows a Gamma distribution. This fills the gap between Poisson-Gamma (equidispersed counts) and NegBin-NegBin (integer denominator): use it when counts are overdispersed and effort is continuous.
Use cases:
Overdispersed catch per unit effort (CPUE)
Overdispersed observations per continuous exposure
Any ratio of overdispersed counts over continuous positive effort
ratiod_negbin_gamma(link_num = "log", link_denom = "log")ratiod_negbin_gamma(link_num = "log", link_denom = "log")
link_num |
Link function for the count rate: "log". |
link_denom |
Link function for the effort mean: "log" (default), "identity" or "inverse". |
A ratiod_family object
# Create family object fam <- ratiod_negbin_gamma() print(fam) # Simulate overdispersed CPUE data set.seed(123) n <- 60 df <- data.frame( catch = rnbinom(n, size = 5, mu = 8), effort_hours = rgamma(n, shape = 4, rate = 1), depth = rnorm(n), vessel = factor(rep(1:6, each = n/6)) ) # Fit model fit <- tratio( catch | effort_hours ~ depth + (1 | vessel), data = df, family = ratiod_negbin_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create family object fam <- ratiod_negbin_gamma() print(fam) # Simulate overdispersed CPUE data set.seed(123) n <- 60 df <- data.frame( catch = rnbinom(n, size = 5, mu = 8), effort_hours = rgamma(n, shape = 4, rate = 1), depth = rnorm(n), vessel = factor(rep(1:6, each = n/6)) ) # Fit model fit <- tratio( catch | effort_hours ~ depth + (1 | vessel), data = df, family = ratiod_negbin_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Two-process model where both numerator and denominator follow negative binomial distributions. This is the core family for count-count ratios where both processes exhibit overdispersion.
Use cases:
Species counts / total counts (relative abundance)
Events / opportunities (both counts)
Any ratio where numerator and denominator are overdispersed counts
ratiod_negbin_negbin(link_num = "log", link_denom = "log")ratiod_negbin_negbin(link_num = "log", link_denom = "log")
link_num |
Link function for the numerator mean: "log". |
link_denom |
Link function for the denominator mean: "log". |
A ratiod_family object
# Create family object fam <- ratiod_negbin_negbin() print(fam) # Simulate data for relative abundance set.seed(123) n <- 60 df <- data.frame( sp_count = rnbinom(n, size = 5, mu = 15), total_count = rnbinom(n, size = 8, mu = 100), habitat = factor(rep(c("forest", "grassland"), each = n/2)), site = factor(rep(1:10, each = n/10)) ) # Fit model fit <- tratio( sp_count | total_count ~ habitat + (1 | site), data = df, family = ratiod_negbin_negbin(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create family object fam <- ratiod_negbin_negbin() print(fam) # Simulate data for relative abundance set.seed(123) n <- 60 df <- data.frame( sp_count = rnbinom(n, size = 5, mu = 15), total_count = rnbinom(n, size = 8, mu = 100), habitat = factor(rep(c("forest", "grassland"), each = n/2)), site = factor(rep(1:10, each = n/10)) ) # Fit model fit <- tratio( sp_count | total_count ~ habitat + (1 | site), data = df, family = ratiod_negbin_negbin(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
One-inflated binomial for proportions with excess ones (100% success). Useful when perfect detection/success has a structural component.
ratiod_oibinomial(link_num = "logit", link_oi = "logit")ratiod_oibinomial(link_num = "logit", link_oi = "logit")
link_num |
Link function for the success probability: "logit" (default), "probit" or "cloglog". |
link_oi |
Link function for the one-inflation probability: "logit" (default), "probit" or "cloglog". |
Use this when excess 100% success rates occur due to:
Perfect detection in highly suitable habitat
Saturation effects
Structural constraints ensuring full success
A ratiod_family object with one-inflation
# Create family object fam <- ratiod_oibinomial() print(fam)# Create family object fam <- ratiod_oibinomial() print(fam)
Two-process model where the numerator follows a Poisson distribution and the denominator (effort/exposure) follows a Gamma distribution. This is the natural family for CPUE-type data.
Use cases:
Catch per unit effort (CPUE)
Observations per hour
Events per continuous exposure
ratiod_poisson_gamma(link_num = "log", link_denom = "log")ratiod_poisson_gamma(link_num = "log", link_denom = "log")
link_num |
Link function for the count rate: "log". |
link_denom |
Link function for the effort mean: "log" (default), "identity" or "inverse". |
A ratiod_family object
# Create family object fam <- ratiod_poisson_gamma() print(fam) # Simulate CPUE data set.seed(123) n <- 60 season <- factor(rep(c("spring", "summer", "fall"), each = n/3)) vessel <- factor(rep(1:10, each = n/10)) u <- rnorm(10, sd = 0.5)[vessel] depth <- rnorm(n) df <- data.frame( catch = rpois(n, exp(1.5 + 0.3 * depth + 0.4 * (season == "summer") + u)), effort_hours = rgamma(n, shape = 4, rate = exp(-u)), depth = depth, season = season, vessel = vessel ) # Fit model fit <- tratio( catch | effort_hours ~ depth + season + (1 | vessel), data = df, family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create family object fam <- ratiod_poisson_gamma() print(fam) # Simulate CPUE data set.seed(123) n <- 60 season <- factor(rep(c("spring", "summer", "fall"), each = n/3)) vessel <- factor(rep(1:10, each = n/10)) u <- rnorm(10, sd = 0.5)[vessel] depth <- rnorm(n) df <- data.frame( catch = rpois(n, exp(1.5 + 0.3 * depth + 0.4 * (season == "summer") + u)), effort_hours = rgamma(n, shape = 4, rate = exp(-u)), depth = depth, season = season, vessel = vessel ) # Fit model fit <- tratio( catch | effort_hours ~ depth + season + (1 | vessel), data = df, family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Specify priors for model parameters. Supports both PC (penalized complexity) priors for variance components and standard distributions for other parameters.
ratiod_priors( beta = NULL, sigma = NULL, phi = NULL, rho_temporal = NULL, rho_spatial = NULL ) ## S3 method for class 'ratiod_priors' print(x, ...) ## S3 method for class 'ratiod_prior' print(x, ...)ratiod_priors( beta = NULL, sigma = NULL, phi = NULL, rho_temporal = NULL, rho_spatial = NULL ) ## S3 method for class 'ratiod_priors' print(x, ...) ## S3 method for class 'ratiod_prior' print(x, ...)
beta |
Prior for fixed effects. Default: |
sigma |
Prior for random effect SDs. Default: |
phi |
Prior for the dispersion parameter. Default: a Gamma prior set
by family: |
rho_temporal |
Prior for a temporal autocorrelation, placed on
|
rho_spatial |
Prior for spatial proportion (BYM2). Default:
|
x |
A |
... |
Ignored |
The MCMC samplers (mode = "hmc", "ess", "pg", "gibbs", "vi",
"sghmc", "sgld") evaluate a Normal prior with mean 0 on fixed effects,
a Half-Cauchy prior on random-effect SDs and a Gamma prior on the dispersion
parameter; a component specified with another family is refused when the
model is fitted. The Laplace backend takes prior_pc() or
prior_half_normal() on sigma. The constructors prior_normal(),
prior_half_normal(), prior_half_cauchy(), prior_gamma(),
prior_beta(), prior_exponential() and prior_pc() build each
component; prior_predict() draws from any of them.
A ratiod_priors object. A component left NULL takes its
default, the same for every sampler and for prior_predict().
Simpson, D., Rue, H., Riebler, A., Martins, T. G., & Sorbye, S. H. (2017). Penalising model component complexity: A principled, practical approach to constructing priors. Statistical Science, 32(1), 1-28.
# Default priors ratiod_priors() # Custom fixed effect prior ratiod_priors(beta = prior_normal(0, 1)) # Tighter random effect prior ratiod_priors(sigma = prior_half_cauchy(0.5)) # PC prior on the random effect SD (Laplace backend, prior_predict()) ratiod_priors(sigma = prior_pc(U = 0.5, alpha = 0.01)) # Informative prior for temporal correlation ratiod_priors(rho_temporal = prior_beta(5, 2)) # Prior mode at ~0.8# Default priors ratiod_priors() # Custom fixed effect prior ratiod_priors(beta = prior_normal(0, 1)) # Tighter random effect prior ratiod_priors(sigma = prior_half_cauchy(0.5)) # PC prior on the random effect SD (Laplace backend, prior_predict()) ratiod_priors(sigma = prior_pc(U = 0.5, alpha = 0.01)) # Informative prior for temporal correlation ratiod_priors(rho_temporal = prior_beta(5, 2)) # Prior mode at ~0.8
Functions to specify spatial random effects for ratio models. Spatial effects are shared between numerator and denominator by default, which helps prevent spurious ratio effects from spatially-structured unmeasured confounders.
## S3 method for class 'ratiod_spatial' print(x, ...)## S3 method for class 'ratiod_spatial' print(x, ...)
x |
A ratiod_spatial object |
... |
Ignored |
adj <- matrix(0, 4, 4) adj[cbind(1:3, 2:4)] <- 1 adj <- adj + t(adj) spatial_car(adj, group_var = "region") spatial_gp(coords = c("x", "y"))adj <- matrix(0, 4, 4) adj[cbind(1:3, 2:4)] <- 1 adj <- adj + t(adj) spatial_car(adj, group_var = "region") spatial_gp(coords = c("x", "y"))
Functions to specify temporal random effects for ratio models. Temporal effects are shared between numerator and denominator by default, which helps prevent spurious ratio effects from temporally-structured unmeasured confounders.
## S3 method for class 'ratiod_temporal' print(x, ...)## S3 method for class 'ratiod_temporal' print(x, ...)
x |
A ratiod_temporal object |
... |
Ignored |
temporal_rw1("year") temporal_ar1("year")temporal_rw1("year") temporal_ar1("year")
Returns the maximum number of OpenMP threads available.
ratiod_threads()ratiod_threads()
Integer number of threads
ratiod_threads()ratiod_threads()
Visual and numerical checks comparing observed data to posterior predictive distributions. Essential for assessing model fit.
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) if (requireNamespace("ggplot2", quietly = TRUE)) { pp_check(fit, type = "dens_overlay") } loo(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) if (requireNamespace("ggplot2", quietly = TRUE)) { pp_check(fit, type = "dens_overlay") } loo(fit)
Distribution families for handling excess zeros in count data. Zero-inflation occurs when there are more zeros than expected under the standard count distribution.
Two model types are available:
Zero-inflated (ZI): Mixture of point mass at zero and count distribution
Hurdle: Separate processes for zero vs positive counts
## S3 method for class 'ratiod_family_zi' print(x, ...)## S3 method for class 'ratiod_family_zi' print(x, ...)
x |
A ratiod_family_zi object |
... |
Ignored |
ratiod_zipois() ratiod_hurdle_negbin()ratiod_zipois() ratiod_hurdle_negbin()
Zero-inflated binomial for the numerator when modeling proportions with excess zeros (e.g., species detection/non-detection).
ratiod_zibinomial(link_num = "logit", link_zi = "logit")ratiod_zibinomial(link_num = "logit", link_zi = "logit")
link_num |
Link function for the success probability: "logit" (default), "probit" or "cloglog". |
link_zi |
Link function for the zero-inflation probability: "logit" (default), "probit" or "cloglog". |
Use this family when modeling proportions (successes/trials) where excess zeros occur beyond what binomial predicts. Common applications:
Species occupancy with false negatives
Survey data with non-response
Epidemiological data with under-reporting
A ratiod_family object with zero-inflation
ratiod_binomial(), ratiod_zinegbin()
# Create family object fam <- ratiod_zibinomial() print(fam)# Create family object fam <- ratiod_zibinomial() print(fam)
Zero-inflated negative binomial for the numerator process. Models excess zeros as a mixture:
where is the zero-inflation probability and
is the negative binomial PMF.
ratiod_zinegbin( link_num = "log", link_denom = "log", link_zi = "logit", denom_family = "negbin" )ratiod_zinegbin( link_num = "log", link_denom = "log", link_zi = "logit", denom_family = "negbin" )
link_num |
Link function for the count mean: "log". |
link_denom |
Link function for the denominator mean: "log" (default); a gamma denominator also takes "identity" or "inverse". |
link_zi |
Link function for the zero-inflation probability: "logit" (default), "probit" or "cloglog". |
denom_family |
Denominator distribution: "negbin" (default), "gamma" for continuous effort, or "fixed" for a known exposure. A known exposure is data: it scales the count mean, which is then the rate times the exposure, and the ratio is the rate. |
The zero-inflation probability can have its own linear predictor via the
zi argument in tratio(). If not specified, a single intercept is estimated.
For ecological data, zero-inflation often represents:
Structural zeros (species truly absent vs not detected)
Sampling zeros (inadequate effort)
False negatives
A ratiod_family object with zero-inflation
ratiod_hurdle_negbin() for hurdle model alternative
# Create family object fam <- ratiod_zinegbin() print(fam) # Simulate zero-inflated count data set.seed(123) n <- 60 zi_prob <- 0.3 habitat <- factor(rep(c("forest", "grassland"), each = n/2)) site <- factor(rep(1:10, each = n/10)) u <- rnorm(10, sd = 0.5)[site] df <- data.frame( count = ifelse(runif(n) < zi_prob, 0, rnbinom(n, size = 3, mu = exp(2 + 0.5 * (habitat == "grassland") + u))), total = rnbinom(n, size = 5, mu = 50 * exp(u)), habitat = habitat, site = site ) # Fit model fit <- tratio( count | total ~ habitat + (1 | site), data = df, family = ratiod_zinegbin(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create family object fam <- ratiod_zinegbin() print(fam) # Simulate zero-inflated count data set.seed(123) n <- 60 zi_prob <- 0.3 habitat <- factor(rep(c("forest", "grassland"), each = n/2)) site <- factor(rep(1:10, each = n/10)) u <- rnorm(10, sd = 0.5)[site] df <- data.frame( count = ifelse(runif(n) < zi_prob, 0, rnbinom(n, size = 3, mu = exp(2 + 0.5 * (habitat == "grassland") + u))), total = rnbinom(n, size = 5, mu = 50 * exp(u)), habitat = habitat, site = site ) # Fit model fit <- tratio( count | total ~ habitat + (1 | site), data = df, family = ratiod_zinegbin(), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Zero-inflated Poisson for the numerator process. Useful when overdispersion is primarily due to excess zeros.
ratiod_zipois( link_num = "log", link_denom = "log", link_zi = "logit", denom_family = "gamma" )ratiod_zipois( link_num = "log", link_denom = "log", link_zi = "logit", denom_family = "gamma" )
link_num |
Link function for the count mean: "log". |
link_denom |
Link function for the denominator mean: "log" (default); a gamma denominator also takes "identity" or "inverse". |
link_zi |
Link function for the zero-inflation probability: "logit" (default), "probit" or "cloglog". |
denom_family |
Denominator distribution: "gamma" (default), "negbin" for a count denominator, or "fixed" for a known exposure. |
Use ratiod_zipois() when:
Excess zeros are the main source of overdispersion
The non-zero counts follow Poisson (no additional overdispersion)
Use ratiod_zinegbin() when:
Both excess zeros AND overdispersion in positive counts
A ratiod_family object with zero-inflation
# Create family object fam <- ratiod_zipois() print(fam) # Simulate zero-inflated CPUE data set.seed(123) n <- 60 zi_prob <- 0.25 df <- data.frame( catch = ifelse(runif(n) < zi_prob, 0, rpois(n, lambda = 5)), effort = rgamma(n, shape = 4, rate = 1), depth = rnorm(n), vessel = factor(rep(1:6, each = n/6)) ) # Fit model fit <- tratio( catch | effort ~ depth + (1 | vessel), data = df, family = ratiod_zipois(denom_family = "gamma"), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create family object fam <- ratiod_zipois() print(fam) # Simulate zero-inflated CPUE data set.seed(123) n <- 60 zi_prob <- 0.25 df <- data.frame( catch = ifelse(runif(n) < zi_prob, 0, rpois(n, lambda = 5)), effort = rgamma(n, shape = 4, rate = 1), depth = rnorm(n), vessel = factor(rep(1:6, each = n/6)) ) # Fit model fit <- tratio( catch | effort ~ depth + (1 | vessel), data = df, family = ratiod_zipois(denom_family = "gamma"), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Zero-and-one inflated binomial for proportions with excess zeros AND ones. Commonly needed in ecological applications where both absence (structural) and perfect detection (saturation) have distinct processes.
ratiod_zoibinomial(link_num = "logit", link_zi = "logit", link_oi = "logit")ratiod_zoibinomial(link_num = "logit", link_zi = "logit", link_oi = "logit")
link_num |
Link function for the success probability: "logit" (default), "probit" or "cloglog". |
link_zi |
Link function for the zero-inflation probability: "logit" (default), "probit" or "cloglog". |
link_oi |
Link function for the one-inflation probability: "logit" (default), "probit" or "cloglog". |
The ZOIB (Zero-and-One Inflated Binomial) model has three components:
Zero process: P(structural zero) = pi_0
One process: P(structural one | not zero) = pi_1
Binomial process: P(Y=y | non-structural) follows binomial
Use when both boundaries (0 and n) have excess observations.
A ratiod_family object with zero-and-one inflation
# Create family object fam <- ratiod_zoibinomial() print(fam)# Create family object fam <- ratiod_zoibinomial() print(fam)
Generate simulated datasets from prior or posterior predictive distributions. Useful for prior predictive checks, simulation-based calibration, and understanding model behavior.
sim_ratiod( n = 100, family = ratiod_negbin_negbin(), formula = NULL, beta_num = NULL, beta_denom = NULL, sigma_re = 0.5, phi_num = 5, phi_denom = 5, n_groups = 10, n_per_group = NULL, spatial = NULL, seed = NULL ) ## S3 method for class 'ratiod_simdata' print(x, ...)sim_ratiod( n = 100, family = ratiod_negbin_negbin(), formula = NULL, beta_num = NULL, beta_denom = NULL, sigma_re = 0.5, phi_num = 5, phi_denom = 5, n_groups = 10, n_per_group = NULL, spatial = NULL, seed = NULL ) ## S3 method for class 'ratiod_simdata' print(x, ...)
n |
Number of observations to simulate. |
family |
A tulpaRatio family object specifying the distributions. |
formula |
Optional formula for generating covariates and structure. If NULL, generates intercept-only data. |
beta_num |
Numeric vector of true numerator coefficients. If NULL, drawn from prior. |
beta_denom |
Numeric vector of true denominator coefficients. If NULL, drawn from prior. |
sigma_re |
Standard deviation of random effects. Default 0.5. |
phi_num |
Overdispersion for numerator (negbin/gamma families). Default 5. |
phi_denom |
Overdispersion for denominator (negbin/gamma families). Default 5. |
n_groups |
Number of random effect groups. Default 10. |
n_per_group |
Observations per group. If NULL, groups are unbalanced. |
spatial |
Optional spatial structure specification. |
seed |
Random seed for reproducibility. |
x |
A ratiod_simdata object |
... |
Ignored |
The simulation generates data according to the ratio model structure:
For negbin_negbin:
y_num ~ NegBin(mu_num, phi_num)
y_denom ~ NegBin(mu_denom, phi_denom)
log(mu_num) = X %*% beta_num + re
log(mu_denom) = X %*% beta_denom + re (shared RE)
For binomial:
y_num ~ Binomial(y_denom, p)
logit(p) = X %*% beta + re
For poisson_gamma:
y_num ~ Poisson(mu_num)
y_denom ~ Gamma(shape, rate) where rate = shape / mu_denom
A list with class ratiod_simdata containing:
Data frame with simulated observations
List of true parameter values used
Family specification
Number of observations
tratio() for model fitting, prior_predict() for prior checks
# Simulate CPUE-type data sim <- sim_ratiod( n = 200, family = ratiod_poisson_gamma(), beta_num = c(2, 0.5), # intercept + depth effect beta_denom = c(1, 0.2), # effort model sigma_re = 0.3, n_groups = 20, seed = 123 ) # View simulated data head(sim$data) # Check true ratio distribution true_ratio <- sim$data$y_num / sim$data$y_denom hist(true_ratio, main = "Simulated ratios") # Simulate binomial proportion data sim_binom <- sim_ratiod( n = 100, family = ratiod_binomial(), beta_num = c(-0.5, 1), # logit-scale coefficients n_groups = 10, seed = 456 ) # Fit model to simulated data fit <- tratio( y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) # Compare estimated vs. true parameters # summary(fit) # sim$true_params# Simulate CPUE-type data sim <- sim_ratiod( n = 200, family = ratiod_poisson_gamma(), beta_num = c(2, 0.5), # intercept + depth effect beta_denom = c(1, 0.2), # effort model sigma_re = 0.3, n_groups = 20, seed = 123 ) # View simulated data head(sim$data) # Check true ratio distribution true_ratio <- sim$data$y_num / sim$data$y_denom hist(true_ratio, main = "Simulated ratios") # Simulate binomial proportion data sim_binom <- sim_ratiod( n = 100, family = ratiod_binomial(), beta_num = c(-0.5, 1), # logit-scale coefficients n_groups = 10, seed = 456 ) # Fit model to simulated data fit <- tratio( y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 400, chains = 1, seed = 1) ) # Compare estimated vs. true parameters # summary(fit) # sim$true_params
Generate multiple simulated datasets for simulation-based calibration (SBC). Each dataset is generated with parameters drawn from the prior.
sim_ratiod_sbc( n_sims = 100, n = 100, family = ratiod_negbin_negbin(), priors = NULL, n_groups = 10, seed = NULL ) ## S3 method for class 'ratiod_sbc_sims' print(x, ...)sim_ratiod_sbc( n_sims = 100, n = 100, family = ratiod_negbin_negbin(), priors = NULL, n_groups = 10, seed = NULL ) ## S3 method for class 'ratiod_sbc_sims' print(x, ...)
n_sims |
Number of simulations to run. |
n |
Number of observations per simulation. |
family |
A tulpaRatio family object. |
priors |
Prior specification for parameter generation. |
n_groups |
Number of random effect groups. |
seed |
Random seed. |
x |
A ratiod_sbc_sims object |
... |
Ignored |
A list of ratiod_simdata objects.
sim_ratiod() for single dataset simulation
# Generate datasets for SBC (returns immediately) sims <- sim_ratiod_sbc( n_sims = 3, n = 30, family = ratiod_negbin_negbin(), seed = 42 ) print(sims) # Fit model to each simulation # ranks <- lapply(sims, function(sim) { # fit <- tratio(y_num | y_denom ~ x1 + (1 | group), # data = sim$data, # family = ratiod_negbin_negbin(), # control = list(iter = 400, chains = 1, seed = 1)) # sum(as.matrix(fit$draws)[, "beta_num[1]"] < sim$true_params$beta_num[1]) # })# Generate datasets for SBC (returns immediately) sims <- sim_ratiod_sbc( n_sims = 3, n = 30, family = ratiod_negbin_negbin(), seed = 42 ) print(sims) # Fit model to each simulation # ranks <- lapply(sims, function(sim) { # fit <- tratio(y_num | y_denom ~ x1 + (1 | group), # data = sim$data, # family = ratiod_negbin_negbin(), # control = list(iter = 400, chains = 1, seed = 1)) # sum(as.matrix(fit$draws)[, "beta_num[1]"] < sim$true_params$beta_num[1]) # })
Specify a Besag-York-Mollie 2 (BYM2) spatial random effect. BYM2 decomposes the spatial effect into a structured (ICAR) component and an unstructured (IID) component, with a mixing parameter controlling the proportion of variance attributable to spatial structure.
BYM2 is preferred over plain CAR when you want to:
Distinguish structured vs unstructured spatial variation
Have an interpretable spatial fraction parameter
Use the scaling from Riebler et al. (2016)
spatial_bym2( adjacency, level = c("group", "obs"), group_var = NULL, shared = TRUE, scale_factor = NULL, parameterization = c("standard", "collapsed") )spatial_bym2( adjacency, level = c("group", "obs"), group_var = NULL, shared = TRUE, scale_factor = NULL, parameterization = c("standard", "collapsed") )
adjacency |
Adjacency matrix (sparse or dense). A symmetric matrix where entry (i,j) is 1 if areas i and j are neighbors, 0 otherwise. |
level |
Level at which spatial structure applies:
|
group_var |
Name of the grouping variable in data (required if
|
shared |
Logical; if TRUE (default), spatial effect enters both numerator and denominator linear predictors identically. |
scale_factor |
Scaling factor for the ICAR component. If NULL (default), computed from the adjacency matrix following Riebler et al. |
parameterization |
Parameterization of the spatial random effect.
Either |
A ratiod_spatial object
Riebler, A., Sorbye, S. H., Simpson, D., & Rue, H. (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research, 25(4), 1145-1165.
# Create adjacency matrix for 10 regions (chain structure) adj <- matrix(0, 10, 10) for (i in 1:9) { adj[i, i+1] <- adj[i+1, i] <- 1 } # Create BYM2 spatial structure bym2 <- spatial_bym2(adj, level = "group", group_var = "region") print(bym2) # Generate synthetic epidemiological data set.seed(456) n_regions <- 10 epi_data <- data.frame( region = factor(1:n_regions), age = rnorm(n_regions, 50, 10), cases = rbinom(n_regions, size = 100, prob = 0.15), population = rep(100L, n_regions) ) # Disease mapping with BYM2 spatial smoothing fit <- tratio( cases | population ~ age, spatial = spatial_bym2(adj, level = "group", group_var = "region"), data = epi_data, family = ratiod_binomial(), mode = "laplace" ) summary(fit)# Create adjacency matrix for 10 regions (chain structure) adj <- matrix(0, 10, 10) for (i in 1:9) { adj[i, i+1] <- adj[i+1, i] <- 1 } # Create BYM2 spatial structure bym2 <- spatial_bym2(adj, level = "group", group_var = "region") print(bym2) # Generate synthetic epidemiological data set.seed(456) n_regions <- 10 epi_data <- data.frame( region = factor(1:n_regions), age = rnorm(n_regions, 50, 10), cases = rbinom(n_regions, size = 100, prob = 0.15), population = rep(100L, n_regions) ) # Disease mapping with BYM2 spatial smoothing fit <- tratio( cases | population ~ age, spatial = spatial_bym2(adj, level = "group", group_var = "region"), data = epi_data, family = ratiod_binomial(), mode = "laplace" ) summary(fit)
Specify a conditional autoregressive (CAR) spatial random effect. Supports both intrinsic CAR (ICAR) and proper CAR variants.
spatial_car( adjacency, level = c("group", "obs"), group_var = NULL, proper = FALSE, shared = TRUE, parameterization = c("standard", "collapsed"), center = TRUE )spatial_car( adjacency, level = c("group", "obs"), group_var = NULL, proper = FALSE, shared = TRUE, parameterization = c("standard", "collapsed"), center = TRUE )
adjacency |
Adjacency matrix (sparse or dense). A symmetric matrix where entry (i,j) is 1 if areas i and j are neighbors, 0 otherwise. |
level |
Level at which spatial structure applies:
|
group_var |
Name of the grouping variable in data (required if
|
proper |
Logical; if FALSE (default), uses Intrinsic CAR (ICAR) with rho = 1 fixed. If TRUE, uses proper CAR with rho estimated. See Details. |
shared |
Logical; if TRUE (default), spatial effect enters both numerator and denominator linear predictors identically. |
parameterization |
Parameterization of the spatial random effect.
Either |
center |
Logical; if TRUE (default), the field's mean is removed on its
way into the linear predictor, so the intercept carries the level and the
field carries deviations from it. The prior still reads the uncentred
field. |
The CAR model specifies that:
where is the number of neighbors of area i and
denotes that j is a neighbor of i.
This leads to a precision matrix:
where D is the diagonal matrix of neighbor counts and W is the adjacency matrix.
ICAR (proper = FALSE, default)
Sets rho = 1 (fixed)
Improper prior (rank-deficient Q)
Requires sum-to-zero constraint for identifiability
Simpler: one fewer parameter to estimate
Standard choice for disease mapping
Equivalent to RW1 on a graph
Proper CAR (proper = TRUE)
Estimates rho in (0, 1) from data
Proper prior (integrates to 1)
rho measures spatial autocorrelation strength
rho -> 0: approaches independence (IID)
rho -> 1: approaches ICAR
More flexible but additional parameter to estimate
Prior: rho ~ Beta(1, 1) (uniform on (0, 1))
A ratiod_spatial object
spatial_bym2() for decomposed spatial + IID effects,
spatial_gp() for continuous spatial effects
# Create adjacency matrix for 10 regions (chain structure) adj <- matrix(0, 10, 10) for (i in 1:9) { adj[i, i+1] <- adj[i+1, i] <- 1 } # ICAR (default, rho = 1 fixed) icar <- spatial_car(adj, level = "group", group_var = "site") print(icar) # Proper CAR (rho estimated) proper_car <- spatial_car(adj, level = "group", group_var = "site", proper = TRUE) print(proper_car) # Generate synthetic data with spatial structure set.seed(123) n_sites <- 10 n_per_site <- 5 df <- data.frame( site = rep(1:n_sites, each = n_per_site), x = rnorm(n_sites * n_per_site), count = rpois(n_sites * n_per_site, 20), effort = rgamma(n_sites * n_per_site, shape = 5, rate = 1) ) # ICAR model (standard disease mapping approach) fit_icar <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), spatial = spatial_car(adj, level = "group", group_var = "site"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Proper CAR model (estimate spatial autocorrelation) fit_car <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), spatial = spatial_car(adj, level = "group", group_var = "site", proper = TRUE), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Extract rho from proper CAR summary(fit_car) # Shows rho_spatial parameter# Create adjacency matrix for 10 regions (chain structure) adj <- matrix(0, 10, 10) for (i in 1:9) { adj[i, i+1] <- adj[i+1, i] <- 1 } # ICAR (default, rho = 1 fixed) icar <- spatial_car(adj, level = "group", group_var = "site") print(icar) # Proper CAR (rho estimated) proper_car <- spatial_car(adj, level = "group", group_var = "site", proper = TRUE) print(proper_car) # Generate synthetic data with spatial structure set.seed(123) n_sites <- 10 n_per_site <- 5 df <- data.frame( site = rep(1:n_sites, each = n_per_site), x = rnorm(n_sites * n_per_site), count = rpois(n_sites * n_per_site, 20), effort = rgamma(n_sites * n_per_site, shape = 5, rate = 1) ) # ICAR model (standard disease mapping approach) fit_icar <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), spatial = spatial_car(adj, level = "group", group_var = "site"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Proper CAR model (estimate spatial autocorrelation) fit_car <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), spatial = spatial_car(adj, level = "group", group_var = "site", proper = TRUE), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Extract rho from proper CAR summary(fit_car) # Shows rho_spatial parameter
Specify a Gaussian Process (GP) spatial random effect using coordinate-based distances. Uses Nearest Neighbor Gaussian Process (NNGP) approximation for computational efficiency with large datasets (scales to millions of observations).
This provides a continuous spatial effect that captures smooth spatial variation, unlike CAR/BYM2 which require discrete areal units.
spatial_gp( coords, cov = c("exponential", "matern", "gaussian", "spherical"), nu = 1.5, nn = 15, shared = TRUE, scale_coords = TRUE, parameterization = c("centered", "noncentered", "collapsed") ) ## S3 method for class 'ratiod_gp' print(x, ...)spatial_gp( coords, cov = c("exponential", "matern", "gaussian", "spherical"), nu = 1.5, nn = 15, shared = TRUE, scale_coords = TRUE, parameterization = c("centered", "noncentered", "collapsed") ) ## S3 method for class 'ratiod_gp' print(x, ...)
coords |
A one-sided formula specifying coordinate columns (e.g.,
|
cov |
Covariance function: |
nu |
Smoothness parameter for Matern covariance. Common values:
|
nn |
Number of nearest neighbors for NNGP approximation. Default 15. Larger values give better approximation but slower computation. |
shared |
Logical; if TRUE (default), spatial effect enters both numerator and denominator. Set to FALSE for process-specific spatial effects (triggers warning about potential confounding). |
scale_coords |
Logical; if TRUE (default), coordinates are scaled to unit variance before computing distances. |
parameterization |
Parameterization of the latent field. One of
|
x |
A ratiod_gp object |
... |
Ignored |
The GP spatial model adds a spatially-correlated random effect to the linear predictor:
where follows a Gaussian process:
The correlation function depends on distance and
range parameter :
Exponential:
Matern:
Gaussian:
Spherical: for
NNGP approximation: For computational tractability, we use the Nearest Neighbor Gaussian Process (Datta et al., 2016), which conditions each location on its k nearest neighbors. This reduces complexity from O(n^3) to O(n*k^2), enabling models with millions of observations.
A ratiod_gp object
Datta, A., Banerjee, S., Finley, A. O., & Gelfand, A. E. (2016). Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, 111(514), 800-812.
spatial_car(), spatial_bym2() for areal spatial effects,
spatial_svc() for spatially-varying coefficients,
spatial_multiscale() for multi-scale spatial effects
# Create GP spatial structure gp <- spatial_gp(~ lon + lat) print(gp) # Generate synthetic spatial data with a smooth field w shared by # count and effort set.seed(789) n <- 40 lon <- runif(n, 0, 10) lat <- runif(n, 0, 10) w <- 0.6 * sin(lon / 2) + 0.6 * cos(lat / 2) depth <- rnorm(n) temp <- rnorm(n) df <- data.frame( lon = lon, lat = lat, depth = depth, temp = temp, count = rpois(n, exp(3 + 0.3 * depth - 0.2 * temp + w)), effort = rgamma(n, shape = 4, rate = exp(-w)) ) # Continuous spatial effect with exponential covariance fit <- tratio( count | effort ~ depth + temp, data = df, family = ratiod_poisson_gamma(), spatial = spatial_gp(~ lon + lat), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1, adapt_delta = 0.95) ) summary(fit) # Smoother spatial field with Matern covariance fit2 <- tratio( count | effort ~ depth + temp, data = df, family = ratiod_poisson_gamma(), spatial = spatial_gp(~ lon + lat, cov = "matern", nu = 1.5), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1, adapt_delta = 0.95) )# Create GP spatial structure gp <- spatial_gp(~ lon + lat) print(gp) # Generate synthetic spatial data with a smooth field w shared by # count and effort set.seed(789) n <- 40 lon <- runif(n, 0, 10) lat <- runif(n, 0, 10) w <- 0.6 * sin(lon / 2) + 0.6 * cos(lat / 2) depth <- rnorm(n) temp <- rnorm(n) df <- data.frame( lon = lon, lat = lat, depth = depth, temp = temp, count = rpois(n, exp(3 + 0.3 * depth - 0.2 * temp + w)), effort = rgamma(n, shape = 4, rate = exp(-w)) ) # Continuous spatial effect with exponential covariance fit <- tratio( count | effort ~ depth + temp, data = df, family = ratiod_poisson_gamma(), spatial = spatial_gp(~ lon + lat), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1, adapt_delta = 0.95) ) summary(fit) # Smoother spatial field with Matern covariance fit2 <- tratio( count | effort ~ depth + temp, data = df, family = ratiod_poisson_gamma(), spatial = spatial_gp(~ lon + lat, cov = "matern", nu = 1.5), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1, adapt_delta = 0.95) )
Specify a Hilbert Space Gaussian Process approximation for spatial effects. HSGP approximates a GP using Laplacian eigenfunctions, providing O(NM^2) complexity with analytical gradients instead of O(Nk^2) with numerical gradients for NNGP.
This gives approximately 50x speedup over standard NNGP while maintaining high accuracy for smooth spatial fields.
spatial_hsgp(coords, m = 6, c = 1.5, shared = TRUE, scale_coords = TRUE) ## S3 method for class 'ratiod_hsgp' print(x, ...)spatial_hsgp(coords, m = 6, c = 1.5, shared = TRUE, scale_coords = TRUE) ## S3 method for class 'ratiod_hsgp' print(x, ...)
coords |
A one-sided formula specifying coordinate columns (e.g.,
|
m |
Number of basis functions per dimension. Total basis functions will be m^2. Default 6. Higher values give better approximation but slower computation. Recommended range: 5-15. |
c |
Boundary factor controlling domain extension beyond data range.
Default 1.5. The domain is extended to |
shared |
Logical; if TRUE (default), spatial effect enters both numerator and denominator. Set to FALSE for process-specific spatial effects (triggers warning about potential confounding). |
scale_coords |
Logical; if TRUE (default), coordinates are scaled to unit variance before computing basis functions. |
x |
A ratiod_hsgp object |
... |
Ignored |
HSGP approximates a GP as:
where:
are Laplacian eigenfunctions (products of sines)
is the spectral density of the squared exponential kernel
are basis coefficients
The hyperparameters are:
: marginal variance (PC prior: P(sigma > 1) = 0.01)
: lengthscale (LogNormal(0, 1) prior)
Advantages over NNGP:
Analytical gradients enable ~50x speedup
Simple parameter interpretation
Works well for smooth spatial fields
Limitations:
Assumes squared exponential (smooth) covariance
Less accurate for rough fields (Matern with low nu)
Approximation quality depends on m and c choices
A ratiod_hsgp object
Riutort-Mayol, G., Bürkner, P. C., Andersen, M. R., Solin, A., & Vehtari, A. (2023). Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming. Statistics and Computing, 33(1), 17.
spatial_gp() for NNGP-based GP, spatial_car() for areal effects
# Create HSGP spatial structure hsgp <- spatial_hsgp(~ lon + lat) print(hsgp) # Generate synthetic spatial data set.seed(42) n <- 100 lon <- runif(n, 0, 10) lat <- runif(n, 0, 10) w <- 0.6 * sin(lon / 2) + 0.6 * cos(lat / 2) x <- rnorm(n) df <- data.frame( lon = lon, lat = lat, x = x, count = rpois(n, exp(3 + 0.3 * x + w)), effort = rgamma(n, shape = 4, rate = exp(-w)) ) # Fast spatial effect with HSGP (much faster than spatial_gp) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), spatial = spatial_hsgp(~ lon + lat, m = 6), mode = "hmc", control = list(iter = 500, warmup = 250, chains = 2, seed = 1) ) summary(fit)# Create HSGP spatial structure hsgp <- spatial_hsgp(~ lon + lat) print(hsgp) # Generate synthetic spatial data set.seed(42) n <- 100 lon <- runif(n, 0, 10) lat <- runif(n, 0, 10) w <- 0.6 * sin(lon / 2) + 0.6 * cos(lat / 2) x <- rnorm(n) df <- data.frame( lon = lon, lat = lat, x = x, count = rpois(n, exp(3 + 0.3 * x + w)), effort = rgamma(n, shape = 4, rate = exp(-w)) ) # Fast spatial effect with HSGP (much faster than spatial_gp) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), spatial = spatial_hsgp(~ lon + lat, m = 6), mode = "hmc", control = list(iter = 500, warmup = 250, chains = 2, seed = 1) ) summary(fit)
Specify a multi-scale spatial random effect that decomposes spatial variation into local (fine-scale) and regional (broad-scale) components. Each scale has its own range and variance parameters, and both are expanded on one Hilbert-space basis (Solin and Sarkka 2020).
spatial_multiscale( coords, scales = c("local", "regional"), range_local = c(0.01, 1), range_regional = c(1, 10), m = 10, c_boundary = 1.5, shared = TRUE, scale_coords = TRUE, sampler = c("auto", "lbfgs") ) ## S3 method for class 'ratiod_multiscale' print(x, ...)spatial_multiscale( coords, scales = c("local", "regional"), range_local = c(0.01, 1), range_regional = c(1, 10), m = 10, c_boundary = 1.5, shared = TRUE, scale_coords = TRUE, sampler = c("auto", "lbfgs") ) ## S3 method for class 'ratiod_multiscale' print(x, ...)
coords |
A one-sided formula specifying coordinate columns (e.g.,
|
scales |
Character vector specifying scale names. Default: |
range_local |
Prior range for local scale as |
range_regional |
Prior range for regional scale as |
m |
Basis functions per dimension, giving |
c_boundary |
Boundary factor: the domain is extended by this multiple before the basis is built. Default 1.5. |
shared |
Logical; if TRUE (default), spatial effects enter both numerator and denominator. |
scale_coords |
Logical; if TRUE (default), coordinates are scaled to unit variance before computing distances. |
sampler |
Warmup mass-matrix adaptation for the two GP blocks:
|
x |
A ratiod_multiscale object |
... |
Ignored |
The multi-scale model decomposes spatial variation additively:
where each component follows an independent Gaussian process with a squared-exponential covariance:
The range priors centre each scale's lengthscale prior.
Each field is the expansion of m^2 basis functions of the domain with
standard normal coefficients scaled by the covariance's spectral density,
so its cost does not grow with the number of locations.
Identifiability: the data inform the sum of the two fields, and the priors on the variances and ranges decide how it is split between them. Non-overlapping range priors give each scale its own band of the spectrum. A field with one value per location for each scale (a nearest-neighbour GP per scale) leaves that split along a ridge in as many dimensions as there are locations, which NUTS does not traverse, so the basis is the only representation offered.
A ratiod_multiscale object
Solin, A. and Sarkka, S. (2020). Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Computing, 30(2), 419-446.
spatial_gp() for single-scale GP, temporal_multiscale() for
multi-scale temporal effects
# Create multi-scale spatial structure ms <- spatial_multiscale( ~ lon + lat, range_local = c(0.1, 0.5), range_regional = c(1, 5) ) print(ms) # A regional trend plus local variation, shared by count and effort set.seed(101) n <- 60 lon <- runif(n, 0, 10) lat <- runif(n, 0, 10) w <- 0.5 * sin(lon / 3) + 0.5 * cos(lat / 3) + 0.3 * sin(2 * lon) * cos(2 * lat) depth <- rnorm(n) temp <- rnorm(n) df <- data.frame( lon = lon, lat = lat, depth = depth, temp = temp, count = rpois(n, exp(3 + 0.3 * depth - 0.2 * temp + w)), effort = rgamma(n, shape = 5, rate = exp(-w)) ) # Multi-scale spatial with local and regional components fit <- tratio( count | effort ~ depth + temp, data = df, family = ratiod_poisson_gamma(), spatial = spatial_multiscale( ~ lon + lat, range_local = c(0.1, 0.5), range_regional = c(1, 5), m = 6 ), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create multi-scale spatial structure ms <- spatial_multiscale( ~ lon + lat, range_local = c(0.1, 0.5), range_regional = c(1, 5) ) print(ms) # A regional trend plus local variation, shared by count and effort set.seed(101) n <- 60 lon <- runif(n, 0, 10) lat <- runif(n, 0, 10) w <- 0.5 * sin(lon / 3) + 0.5 * cos(lat / 3) + 0.3 * sin(2 * lon) * cos(2 * lat) depth <- rnorm(n) temp <- rnorm(n) df <- data.frame( lon = lon, lat = lat, depth = depth, temp = temp, count = rpois(n, exp(3 + 0.3 * depth - 0.2 * temp + w)), effort = rgamma(n, shape = 5, rate = exp(-w)) ) # Multi-scale spatial with local and regional components fit <- tratio( count | effort ~ depth + temp, data = df, family = ratiod_poisson_gamma(), spatial = spatial_multiscale( ~ lon + lat, range_local = c(0.1, 0.5), range_regional = c(1, 5), m = 6 ), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Apply Restricted Spatial Regression to mitigate spatial confounding. RSR orthogonalizes the spatial effect to the covariate space, preventing the spatial random effect from absorbing covariate information.
This is important when covariates are spatially smooth (e.g., climate variables, elevation) because the spatial random effect can "steal" variance from these covariates, leading to biased coefficient estimates.
spatial_rsr(spatial, restrict_to) ## S3 method for class 'ratiod_rsr' print(x, ...)spatial_rsr(spatial, restrict_to) ## S3 method for class 'ratiod_rsr' print(x, ...)
spatial |
A spatial specification ( |
restrict_to |
Formula specifying which covariates to orthogonalize
against (e.g., |
x |
A ratiod_rsr object |
... |
Passed to underlying print method |
The RSR approach (Reich et al., 2006; Hodges & Reich, 2010) modifies the spatial random effect to be orthogonal to the fixed effect design matrix:
where is the projection matrix onto the column
space of X.
When to use RSR:
Covariates are spatially smooth (environmental gradients)
Interested in causal interpretation of covariate effects
Coefficients appear attenuated toward zero
When NOT to use RSR:
Covariates are spatially uncorrelated
Spatial effect is the primary quantity of interest
Prediction is the main goal (not causal inference)
A modified spatial specification with RSR enabled
Reich, B. J., Hodges, J. S., & Zadnik, V. (2006). Effects of residual smoothing on the posterior of the fixed effects in disease-mapping models. Biometrics, 62(4), 1197-1206.
Hodges, J. S., & Reich, B. J. (2010). Adding spatially-correlated errors can mess up the fixed effect you love. The American Statistician, 64(4), 325-334.
# Create RSR spatial structure rsr <- spatial_rsr( spatial_gp(~ lon + lat), restrict_to = ~ depth + temp ) print(rsr) # Generate synthetic spatial data with confounding set.seed(404) n <- 40 lon <- runif(n, 0, 10) lat <- runif(n, 0, 10) # Smooth field w shared by count and effort w <- 0.6 * sin(lon / 2) + 0.6 * cos(lat / 2) # Make depth and temp spatially correlated depth <- lon/5 + rnorm(n, 0, 0.5) # Correlated with lon temp <- lat/5 + rnorm(n, 0, 0.5) # Correlated with lat df <- data.frame( lon = lon, lat = lat, depth = depth, temp = temp, count = rpois(n, exp(2.5 + 0.3 * depth - 0.2 * temp + w)), effort = rgamma(n, shape = 4, rate = exp(-w)) ) # Standard GP (may have spatial confounding) fit1 <- tratio( count | effort ~ depth + temp, data = df, family = ratiod_poisson_gamma(), spatial = spatial_gp(~ lon + lat), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # RSR to protect depth and temp coefficients fit2 <- tratio( count | effort ~ depth + temp, data = df, family = ratiod_poisson_gamma(), spatial = spatial_rsr( spatial_gp(~ lon + lat), restrict_to = ~ depth + temp ), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Compare coefficient estimates summary(fit1) # May be attenuated summary(fit2) # Protected from spatial confounding# Create RSR spatial structure rsr <- spatial_rsr( spatial_gp(~ lon + lat), restrict_to = ~ depth + temp ) print(rsr) # Generate synthetic spatial data with confounding set.seed(404) n <- 40 lon <- runif(n, 0, 10) lat <- runif(n, 0, 10) # Smooth field w shared by count and effort w <- 0.6 * sin(lon / 2) + 0.6 * cos(lat / 2) # Make depth and temp spatially correlated depth <- lon/5 + rnorm(n, 0, 0.5) # Correlated with lon temp <- lat/5 + rnorm(n, 0, 0.5) # Correlated with lat df <- data.frame( lon = lon, lat = lat, depth = depth, temp = temp, count = rpois(n, exp(2.5 + 0.3 * depth - 0.2 * temp + w)), effort = rgamma(n, shape = 4, rate = exp(-w)) ) # Standard GP (may have spatial confounding) fit1 <- tratio( count | effort ~ depth + temp, data = df, family = ratiod_poisson_gamma(), spatial = spatial_gp(~ lon + lat), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # RSR to protect depth and temp coefficients fit2 <- tratio( count | effort ~ depth + temp, data = df, family = ratiod_poisson_gamma(), spatial = spatial_rsr( spatial_gp(~ lon + lat), restrict_to = ~ depth + temp ), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Compare coefficient estimates summary(fit1) # May be attenuated summary(fit2) # Protected from spatial confounding
Specify spatially-varying coefficients for ratio models. SVCs allow regression coefficients to vary smoothly across space using a Gaussian process prior. This captures local effects that may differ from global relationships.
Uses Nearest Neighbor Gaussian Process (NNGP) approximation for computational efficiency with large datasets.
spatial_svc( coords, terms = 1, cov = c("exponential", "matern", "gaussian", "spherical"), nn = 15, shared = TRUE, scale_coords = TRUE, approx = c("nngp", "hsgp"), m = 6, c_boundary = 1.5, range = svc_default_range(), parameterization = c("noncentered", "centered") ) ## S3 method for class 'ratiod_svc' print(x, ...)spatial_svc( coords, terms = 1, cov = c("exponential", "matern", "gaussian", "spherical"), nn = 15, shared = TRUE, scale_coords = TRUE, approx = c("nngp", "hsgp"), m = 6, c_boundary = 1.5, range = svc_default_range(), parameterization = c("noncentered", "centered") ) ## S3 method for class 'ratiod_svc' print(x, ...)
coords |
A one-sided formula specifying coordinate columns (e.g.,
|
terms |
Which coefficients should vary spatially. Options:
|
cov |
Covariance function: |
nn |
Number of nearest neighbors for NNGP approximation. Default 15. Larger values give better approximation but slower computation. |
shared |
Logical; if TRUE (default), SVC effects enter both numerator and denominator. Set to FALSE for process-specific SVCs (triggers warning about potential confounding). |
scale_coords |
Logical; if TRUE (default), coordinates are scaled to unit variance before computing distances. |
approx |
Spatial GP approximation. One of |
m |
Number of basis functions per coordinate for HSGP (used only
when |
c_boundary |
Boundary scaling factor for HSGP. The HSGP domain is
extended to |
range |
Prior range for the SVC range parameter |
parameterization |
Coordinate each NNGP term is sampled in.
|
x |
A ratiod_svc object |
... |
Ignored |
The SVC model extends the linear predictor:
where:
are global (non-spatial) coefficients
is the subset of covariates with SVCs
are spatially-varying adjustments at location s
Each SVC follows an independent Gaussian process:
where is the correlation function with range parameter
.
NNGP approximation: For computational tractability, we use the Nearest Neighbor Gaussian Process (Datta et al., 2016), which conditions each location on its k nearest neighbors. This reduces complexity from O(n^3) to O(n*k^2), enabling models with thousands of locations.
Interpretation: A positive SVC for depth at location s means the
depth effect is stronger at s than the global average. The spatial
variance quantifies how much the effect varies across
space.
A ratiod_svc object
Datta, A., Banerjee, S., Finley, A. O., & Gelfand, A. E. ( 2016). Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, 111(514), 800-812.
spatial_car(), spatial_bym2() for areal spatial effects,
svc() for extracting SVC posteriors
# Create SVC specification svc_spec <- spatial_svc(~ lon + lat, terms = 1) print(svc_spec) # Generate synthetic spatial data set.seed(202) n <- 30 df <- data.frame( lon = runif(n, 0, 10), lat = runif(n, 0, 10), depth = rnorm(n), temp = rnorm(n), count = rpois(n, 20), effort = rgamma(n, shape = 4, rate = 1) ) # Spatially-varying intercept (random spatial field) fit <- tratio( count | effort ~ depth, data = df, family = ratiod_poisson_gamma(), spatial = spatial_svc(~ lon + lat, terms = 1), control = list(iter = 200, chains = 1, seed = 1) ) summary(fit) # Extract and plot the spatially-varying coefficients svc_effects <- svc(fit) plot(svc_effects)# Create SVC specification svc_spec <- spatial_svc(~ lon + lat, terms = 1) print(svc_spec) # Generate synthetic spatial data set.seed(202) n <- 30 df <- data.frame( lon = runif(n, 0, 10), lat = runif(n, 0, 10), depth = rnorm(n), temp = rnorm(n), count = rpois(n, 20), effort = rgamma(n, shape = 4, rate = 1) ) # Spatially-varying intercept (random spatial field) fit <- tratio( count | effort ~ depth, data = df, family = ratiod_poisson_gamma(), spatial = spatial_svc(~ lon + lat, terms = 1), control = list(iter = 200, chains = 1, seed = 1) ) summary(fit) # Extract and plot the spatially-varying coefficients svc_effects <- svc(fit) plot(svc_effects)
Functions to specify spatiotemporal interaction effects for ratio models. These capture dependencies that arise when spatial patterns vary over time, or when temporal trends differ across space.
Specify a spatiotemporal interaction effect for ratio models. The interaction captures structured or unstructured deviation from the additive spatial + temporal model.
spatiotemporal( spatial, temporal, type = c("I", "II", "III", "IV", "iid", "separable", "nonsep_gp"), shared = TRUE, parameterization = c("centered", "noncentered", "collapsed"), nonsep_type = c("product", "gneiting"), coords = NULL, cov_space = NULL, cov_time = NULL, nn = NULL ) ## S3 method for class 'ratiod_spatiotemporal' print(x, ...)spatiotemporal( spatial, temporal, type = c("I", "II", "III", "IV", "iid", "separable", "nonsep_gp"), shared = TRUE, parameterization = c("centered", "noncentered", "collapsed"), nonsep_type = c("product", "gneiting"), coords = NULL, cov_space = NULL, cov_time = NULL, nn = NULL ) ## S3 method for class 'ratiod_spatiotemporal' print(x, ...)
spatial |
A spatial specification from |
temporal |
A temporal specification from |
type |
Interaction type:
|
shared |
Logical; if TRUE (default), spatiotemporal effect enters both numerator and denominator. Set to FALSE for process-specific effects (triggers warning about potential confounding). |
parameterization |
With the HMC backend: |
nonsep_type |
For |
coords |
Required for |
cov_space, cov_time
|
For the GP types: the spatial and temporal
kernels, one of |
nn |
For the GP types: NNGP neighbours per grid cell (default 15,
capped at |
x |
A ratiod_spatiotemporal object |
... |
Ignored |
Spatiotemporal interactions extend the basic additive model:
to include interactions:
where captures space-time interactions.
Interaction Types (following Knorr-Held, 2000):
Type I: Unstructured interaction - IID
Type II: Structured time, unstructured space - temporal structure at each location
Type III: Structured space, unstructured time - spatial structure at each time
Type IV: Structured space AND time - full Kronecker interaction
Separable Models:
Separable: Covariance is Kronecker product
Non-separable: GP with joint space-time metric
Type I (IID)
Independent random effect for each space-time combination:
This is the simplest form, requiring S*T parameters but capturing no structured interaction.
Type II (Temporal structure per location)
Each location has its own temporal random effect:
This captures location-specific temporal trends but assumes independence across locations.
Type III (Spatial structure per time point)
Each time point has its own spatial random effect:
This captures time-specific spatial patterns but assumes independence across time points.
Type IV (Full structure)
Kronecker product of spatial and temporal precision matrices:
This is the most constrained model, assuming the interaction has the same structure as the marginal effects.
Separable
For GP-based effects, assumes separable covariance:
A ratiod_spatiotemporal object
Knorr-Held, L. (2000). Bayesian modelling of inseparable space-time variation in disease risk. Statistics in Medicine, 19(17-18), 2555-2567.
Margossian, C. C., Vehtari, A., Simpson, D. and Agrawal, R. (2020). Hamiltonian Monte Carlo using an adjoint-differentiated Laplace approximation: Bayesian inference for latent Gaussian models and beyond. Advances in Neural Information Processing Systems 33.
Rue, H., Martino, S. and Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society B, 71(2), 319-392.
spatial_car(), spatial_gp(), temporal_rw1(), temporal_ar1()
# Create adjacency matrix for 10 regions adj <- matrix(0, 10, 10) for (i in 1:9) adj[i, i+1] <- adj[i+1, i] <- 1 # Type I: Unstructured interaction st1 <- spatiotemporal( spatial = spatial_car(adj, level = "group", group_var = "region"), temporal = temporal_rw1("year"), type = "I" ) print(st1) # Type IV: Fully structured interaction st4 <- spatiotemporal( spatial = spatial_car(adj, level = "group", group_var = "region"), temporal = temporal_rw1("year"), type = "IV" ) print(st4) # Generate synthetic spatiotemporal data set.seed(123) n_regions <- 10 n_years <- 8 df <- expand.grid( region = 1:n_regions, year = 2015:(2015 + n_years - 1) ) df$x <- rnorm(nrow(df)) df$count <- rpois(nrow(df), lambda = 20) df$effort <- rgamma(nrow(df), shape = 4, rate = 1) # Fit model with spatiotemporal interaction fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), spatiotemporal = spatiotemporal( spatial = spatial_car(adj, level = "group", group_var = "region"), temporal = temporal_rw1("year"), type = "IV" ), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)# Create adjacency matrix for 10 regions adj <- matrix(0, 10, 10) for (i in 1:9) adj[i, i+1] <- adj[i+1, i] <- 1 # Type I: Unstructured interaction st1 <- spatiotemporal( spatial = spatial_car(adj, level = "group", group_var = "region"), temporal = temporal_rw1("year"), type = "I" ) print(st1) # Type IV: Fully structured interaction st4 <- spatiotemporal( spatial = spatial_car(adj, level = "group", group_var = "region"), temporal = temporal_rw1("year"), type = "IV" ) print(st4) # Generate synthetic spatiotemporal data set.seed(123) n_regions <- 10 n_years <- 8 df <- expand.grid( region = 1:n_regions, year = 2015:(2015 + n_years - 1) ) df$x <- rnorm(nrow(df)) df$count <- rpois(nrow(df), lambda = 20) df$effort <- rgamma(nrow(df), shape = 4, rate = 1) # Fit model with spatiotemporal interaction fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), spatiotemporal = spatiotemporal( spatial = spatial_car(adj, level = "group", group_var = "region"), temporal = temporal_rw1("year"), type = "IV" ), control = list(iter = 400, chains = 1, seed = 1) ) summary(fit)
Extract posterior distributions of spatially-varying coefficients (SVCs) from a fitted ratio model with SVC specification.
svc(object, terms = NULL, summary = FALSE, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_fit' svc(object, terms = NULL, summary = FALSE, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_svc_posterior' summary(object, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_svc_posterior' print(x, ...) ## S3 method for class 'ratiod_svc_posterior' plot(x, term = 1, type = "mean", ...)svc(object, terms = NULL, summary = FALSE, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_fit' svc(object, terms = NULL, summary = FALSE, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_svc_posterior' summary(object, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_svc_posterior' print(x, ...) ## S3 method for class 'ratiod_svc_posterior' plot(x, term = 1, type = "mean", ...)
object |
A |
terms |
Which SVC terms to extract. If NULL (default), extracts all. |
summary |
Logical; if TRUE, return summary statistics instead of full posterior draws. |
probs |
Quantiles to compute if |
... |
Ignored |
x |
A ratiod_svc_posterior object |
term |
Which term to plot (name or index). Default: first term. |
type |
Plot type: "mean" (default), "sd", or quantile (e.g., "q50") |
A ratiod_svc_posterior object containing:
draws: Array of posterior draws (draws x locations x terms)
coords: Coordinate matrix
term_names: Names of SVC terms
spatial_svc(), plot.ratiod_svc_posterior()
# Generate synthetic spatial data set.seed(303) n <- 30 lon <- runif(n, 0, 10) lat <- runif(n, 0, 10) depth <- rnorm(n) # Spatially varying intercept w0 and depth slope w1, shared by count and # effort w0 <- 0.5 * sin(lon / 2) + 0.5 * cos(lat / 2) w1 <- 0.4 * cos(lon / 3) df <- data.frame( lon = lon, lat = lat, depth = depth, count = rpois(n, exp(3 + 0.3 * depth + w0 + w1 * depth)), effort = rgamma(n, shape = 4, rate = exp(-(w0 + w1 * depth))) ) # Fit model with SVC fit <- tratio( count | effort ~ depth, data = df, family = ratiod_poisson_gamma(), spatial = spatial_svc(~ lon + lat, terms = c(1, 2)), control = list(iter = 400, chains = 1, seed = 1) ) # Extract SVC posteriors svc_post <- svc(fit) summary(svc_post)# Generate synthetic spatial data set.seed(303) n <- 30 lon <- runif(n, 0, 10) lat <- runif(n, 0, 10) depth <- rnorm(n) # Spatially varying intercept w0 and depth slope w1, shared by count and # effort w0 <- 0.5 * sin(lon / 2) + 0.5 * cos(lat / 2) w1 <- 0.4 * cos(lon / 3) df <- data.frame( lon = lon, lat = lat, depth = depth, count = rpois(n, exp(3 + 0.3 * depth + w0 + w1 * depth)), effort = rgamma(n, shape = 4, rate = exp(-(w0 + w1 * depth))) ) # Fit model with SVC fit <- tratio( count | effort ~ depth, data = df, family = ratiod_poisson_gamma(), spatial = spatial_svc(~ lon + lat, terms = c(1, 2)), control = list(iter = 400, chains = 1, seed = 1) ) # Extract SVC posteriors svc_post <- svc(fit) summary(svc_post)
Extract posterior distributions of temporal effects from a fitted tulpaRatio model with temporal specification.
temporal( object, component = "all", summary = FALSE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_fit' temporal( object, component = "all", summary = FALSE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_temporal_posterior' print(x, ...) ## S3 method for class 'ratiod_temporal_posterior' summary(object, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_temporal_posterior' plot(x, component = NULL, type = "ribbon", ...)temporal( object, component = "all", summary = FALSE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_fit' temporal( object, component = "all", summary = FALSE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_temporal_posterior' print(x, ...) ## S3 method for class 'ratiod_temporal_posterior' summary(object, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_temporal_posterior' plot(x, component = NULL, type = "ribbon", ...)
object |
A |
component |
Which component to extract for multi-scale models:
|
summary |
Logical; if TRUE, return summary statistics instead of full posterior draws. |
probs |
Quantiles to compute if |
... |
Ignored |
x |
A ratiod_temporal_posterior object |
type |
Plot type: "ribbon" (default) or "line" |
A ratiod_temporal_posterior object
temporal_multiscale(), temporal_rw1()
# Fit model with multi-scale temporal effects set.seed(131) year <- 1:40 # Trend plus a 12-year cycle, shared by count and effort eff <- 0.4 * sin(year / 8) + 0.3 * sin(2 * pi * year / 12) x <- rnorm(40) df <- data.frame( year = year, x = x, count = rpois(40, exp(3 + 0.3 * x + eff)), effort = rgamma(40, shape = 4.5, rate = exp(-eff)) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_multiscale("year", trend = "rw2", seasonal = 12), control = list(iter = 400, chains = 1, seed = 1) ) # Extract all temporal effects temp_post <- temporal(fit) summary(temp_post)# Fit model with multi-scale temporal effects set.seed(131) year <- 1:40 # Trend plus a 12-year cycle, shared by count and effort eff <- 0.4 * sin(year / 8) + 0.3 * sin(2 * pi * year / 12) x <- rnorm(40) df <- data.frame( year = year, x = x, count = rpois(40, exp(3 + 0.3 * x + eff)), effort = rgamma(40, shape = 4.5, rate = exp(-eff)) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_multiscale("year", trend = "rw2", seasonal = 12), control = list(iter = 400, chains = 1, seed = 1) ) # Extract all temporal effects temp_post <- temporal(fit) summary(temp_post)
Specify a first-order autoregressive temporal random effect.
AR1 models temporal correlation where each time point depends on
the previous one: phi[t] = rho * phi[t-1] + epsilon[t].
Unlike RW1/RW2, AR1 is a proper (stationary) model with an estimated autocorrelation parameter rho.
temporal_ar1( time_var, group_var = NULL, shared = TRUE, rho_prior = NULL, parameterization = c("centered", "noncentered") )temporal_ar1( time_var, group_var = NULL, shared = TRUE, rho_prior = NULL, parameterization = c("centered", "noncentered") )
time_var |
Name of the time variable in data. Can be a formula
(e.g., |
group_var |
Optional name of grouping variable for panel data. If provided, a separate AR1 process is estimated for each group. |
shared |
Logical; if TRUE (default), temporal effect enters both numerator and denominator linear predictors identically. |
rho_prior |
Prior on the autocorrelation, as a |
parameterization |
Coordinate the field is sampled in. |
The AR1 process has marginal variance sigma^2 / (1 - rho^2) and correlation between time points t and s of rho^|t-s|.
The precision matrix is tridiagonal and full rank, so no constraints are needed.
rho is estimated on (-1, 1), the range every AR1 correlation in the
package is sampled on, so an alternating series is reachable. The sampled
coordinate is logit((rho + 1) / 2), which is what rho_prior places its
Beta on.
A ratiod_temporal object
# Create temporal AR1 specification temporal_ar1("year") temporal_ar1("year", group = "site") # AR1 temporal correlation set.seed(127) df <- data.frame( year = rep(2005:2024, each = 3), x = rnorm(60), count = rpois(60, lambda = 18), effort = rgamma(60, shape = 3.5, rate = 1) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_ar1("year"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Panel data with group-specific AR1 set.seed(128) df_panel <- data.frame( year = rep(2012:2027, each = 3), site = rep(paste0("site", 1:3), times = 16), x = rnorm(48), count = rpois(48, lambda = 12), effort = rgamma(48, shape = 2.5, rate = 1) ) fit2 <- tratio( count | effort ~ x, data = df_panel, family = ratiod_poisson_gamma(), temporal = temporal_ar1("year", group_var = "site"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) )# Create temporal AR1 specification temporal_ar1("year") temporal_ar1("year", group = "site") # AR1 temporal correlation set.seed(127) df <- data.frame( year = rep(2005:2024, each = 3), x = rnorm(60), count = rpois(60, lambda = 18), effort = rgamma(60, shape = 3.5, rate = 1) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_ar1("year"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Panel data with group-specific AR1 set.seed(128) df_panel <- data.frame( year = rep(2012:2027, each = 3), site = rep(paste0("site", 1:3), times = 16), x = rnorm(48), count = rpois(48, lambda = 12), effort = rgamma(48, shape = 2.5, rate = 1) ) fit2 <- tratio( count | effort ~ x, data = df_panel, family = ratiod_poisson_gamma(), temporal = temporal_ar1("year", group_var = "site"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) )
Specify a Gaussian Process (GP) temporal random effect for irregularly-spaced or continuous time points. Unlike RW1/RW2/AR1 which assume equally-spaced observations, GP temporal effects model correlation as a function of time distance.
This is particularly useful for:
Irregularly-spaced time series
Continuous time (e.g., exact timestamps)
Smooth temporal trends with uncertainty
temporal_gp( time_var, cov = c("exponential", "matern", "gaussian", "periodic"), nu = 1.5, period = NULL, group_var = NULL, shared = TRUE, scale_coords = TRUE, parameterization = c("noncentered", "centered") ) ## S3 method for class 'ratiod_temporal_gp' print(x, ...)temporal_gp( time_var, cov = c("exponential", "matern", "gaussian", "periodic"), nu = 1.5, period = NULL, group_var = NULL, shared = TRUE, scale_coords = TRUE, parameterization = c("noncentered", "centered") ) ## S3 method for class 'ratiod_temporal_gp' print(x, ...)
time_var |
Name of the time variable in data. Can be a formula
(e.g., |
cov |
Covariance function: |
nu |
Smoothness parameter for Matern covariance. Common values:
|
period |
Period for periodic covariance (e.g., 12 for monthly, 365 for
daily data with annual cycle). Only used when |
group_var |
Optional name of grouping variable for panel data. If provided, separate GPs are estimated for each group. |
shared |
Logical; if TRUE (default), temporal effect enters both numerator and denominator. |
scale_coords |
Logical; if TRUE (default), time values are scaled to unit variance before computing distances. |
parameterization |
Parameterization for GP effects:
|
x |
A ratiod_temporal_gp object |
... |
Ignored |
The GP temporal model adds a time-correlated random effect:
where follows a Gaussian process:
The correlation function depends on time distance :
Exponential: - continuous but not differentiable
Matern: Smooth with tunable roughness via
Gaussian: - infinitely differentiable
Periodic: - for seasonal data
Implementation: Uses a state-space representation for O(n) computational complexity when possible (exponential, Matern with half-integer nu).
A ratiod_temporal_gp object
temporal_rw1(), temporal_ar1() for equally-spaced temporal effects,
spatial_gp() for spatial GP effects
# Create GP temporal specification temporal_gp("timestamp") temporal_gp("day", cov = "matern", nu = 1.5) temporal_gp("month", cov = "periodic", period = 12) # Irregularly-spaced time series set.seed(140) times <- sort(runif(30, 0, 100)) df <- data.frame( time = times, x = rnorm(30), count = rpois(30, lambda = 20), effort = rgamma(30, shape = 4, rate = 1) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_gp("time"), control = list(iter = 400, chains = 1, seed = 1) )# Create GP temporal specification temporal_gp("timestamp") temporal_gp("day", cov = "matern", nu = 1.5) temporal_gp("month", cov = "periodic", period = 12) # Irregularly-spaced time series set.seed(140) times <- sort(runif(30, 0, 100)) df <- data.frame( time = times, x = rnorm(30), count = rpois(30, lambda = 20), effort = rgamma(30, shape = 4, rate = 1) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_gp("time"), control = list(iter = 400, chains = 1, seed = 1) )
Specify an independent temporal random effect: each time point gets its own
coefficient with phi[t] ~ N(0, sigma^2), no dependence between them.
Use it for a year or season effect that varies without a trend, and as the
unstructured comparison against which temporal_rw1(), temporal_rw2()
and temporal_ar1() add smoothness.
temporal_iid(time_var, group_var = NULL, shared = TRUE)temporal_iid(time_var, group_var = NULL, shared = TRUE)
time_var |
Name of the time variable in data. Can be a formula
(e.g., |
group_var |
Optional name of grouping variable for panel data.
If provided, separate random walks are estimated for each group. Each
group's walk sums to zero over time, so it carries the group's trend and
not its level; add a group intercept such as |
shared |
Logical; if TRUE (default), temporal effect enters both numerator and denominator linear predictors identically. |
The precision matrix is tau * I: diagonal and full rank, so no constraint
is needed and no ordering of the time points is used. A model whose time
points carry no meaningful order fits the same way here as one whose do.
A ratiod_temporal object
# Create temporal IID specification temporal_iid("year") # Unstructured year effects set.seed(129) df <- data.frame( year = rep(2010:2029, each = 2), x = rnorm(40), count = rpois(40, lambda = 22), effort = rgamma(40, shape = 4, rate = 1) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_iid("year"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) )# Create temporal IID specification temporal_iid("year") # Unstructured year effects set.seed(129) df <- data.frame( year = rep(2010:2029, each = 2), x = rnorm(40), count = rpois(40, lambda = 22), effort = rgamma(40, shape = 4, rate = 1) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_iid("year"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) )
Specify a multi-scale temporal random effect that decomposes temporal variation into trend, seasonal, and short-term components. Each component has its own variance parameter and structure.
This is particularly useful for long time series where patterns exist at multiple temporal scales (e.g., annual cycles, decadal trends).
temporal_multiscale( time_var, trend = c("rw2", "rw1", "none"), seasonal = NULL, short_term = c("ar1", "iid", "none"), group_var = NULL, shared = TRUE, parameterization = c("centered", "noncentered") ) ## S3 method for class 'ratiod_temporal_multiscale' print(x, ...)temporal_multiscale( time_var, trend = c("rw2", "rw1", "none"), seasonal = NULL, short_term = c("ar1", "iid", "none"), group_var = NULL, shared = TRUE, parameterization = c("centered", "noncentered") ) ## S3 method for class 'ratiod_temporal_multiscale' print(x, ...)
time_var |
Name of the time variable in data. Can be a formula
(e.g., |
trend |
Type of trend component: |
seasonal |
Period for seasonal component (integer). Set to |
short_term |
Type of short-term residual component: |
group_var |
Optional name of grouping variable for panel data. |
shared |
Logical; if TRUE (default), temporal effects enter both numerator and denominator. |
parameterization |
Coordinate the intrinsic trend and seasonal arms are
sampled in. |
x |
A ratiod_temporal_multiscale object |
... |
Ignored |
The multi-scale temporal model decomposes variation additively:
where:
trend(t): Long-term smooth change via RW1 or RW2
seasonal(t): Repeating pattern with period P via cyclic RW1
short(t): Residual temporal correlation via AR(1) or IID
Each component has a separate variance parameter with PC priors to favor simpler models (fewer active components).
Sum-to-zero constraints are applied to trend and seasonal components for identifiability with the intercept.
A ratiod_temporal_multiscale object
temporal_rw1(), temporal_rw2(), temporal_ar1() for
single-component temporal effects, spatial_multiscale() for multi-scale
spatial effects
# Create multi-scale temporal specification temporal_multiscale("month_id", trend = "rw2", seasonal = 12, short_term = "ar1") temporal_multiscale("year", trend = "rw2", seasonal = NULL, short_term = "none") # Decompose into trend + seasonal + short-term set.seed(129) month_id <- 1:48 # Trend, 12-month cycle and AR(1) noise, shared by count and effort eff <- 0.4 * sin(month_id / 15) + 0.4 * sin(2 * pi * month_id / 12) + as.numeric(arima.sim(list(ar = 0.6), 48, sd = 0.15)) x <- rnorm(48) df <- data.frame( month_id = month_id, x = x, count = rpois(48, exp(3 + 0.3 * x + eff)), effort = rgamma(48, shape = 4, rate = exp(-eff)) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_multiscale( time_var = "month_id", trend = "rw2", seasonal = 12, short_term = "ar1" ), control = list(iter = 1000, chains = 1, seed = 1) ) # Extract and plot components temporal_effects <- temporal(fit) plot(temporal_effects, component = "trend")# Create multi-scale temporal specification temporal_multiscale("month_id", trend = "rw2", seasonal = 12, short_term = "ar1") temporal_multiscale("year", trend = "rw2", seasonal = NULL, short_term = "none") # Decompose into trend + seasonal + short-term set.seed(129) month_id <- 1:48 # Trend, 12-month cycle and AR(1) noise, shared by count and effort eff <- 0.4 * sin(month_id / 15) + 0.4 * sin(2 * pi * month_id / 12) + as.numeric(arima.sim(list(ar = 0.6), 48, sd = 0.15)) x <- rnorm(48) df <- data.frame( month_id = month_id, x = x, count = rpois(48, exp(3 + 0.3 * x + eff)), effort = rgamma(48, shape = 4, rate = exp(-eff)) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_multiscale( time_var = "month_id", trend = "rw2", seasonal = 12, short_term = "ar1" ), control = list(iter = 1000, chains = 1, seed = 1) ) # Extract and plot components temporal_effects <- temporal(fit) plot(temporal_effects, component = "trend")
Apply Restricted Temporal Regression to mitigate temporal confounding. RTR orthogonalizes the temporal effect to the covariate space, preventing the temporal random effect from absorbing covariate information.
This is important when covariates have temporal trends (e.g., increasing over time) because the temporal random effect can "steal" variance from these covariates, leading to biased coefficient estimates.
temporal_rtr(temporal, restrict_to) ## S3 method for class 'ratiod_rtr' print(x, ...)temporal_rtr(temporal, restrict_to) ## S3 method for class 'ratiod_rtr' print(x, ...)
temporal |
A temporal specification ( |
restrict_to |
Formula specifying which covariates to orthogonalize
against (e.g., |
x |
A ratiod_rtr object |
... |
Passed to underlying print method |
The RTR approach (analogous to Restricted Spatial Regression) modifies the temporal random effect to be orthogonal to the fixed effect design matrix:
where is the projection matrix onto the column
space of X.
When to use RTR:
Covariates have temporal trends (e.g., increasing/decreasing over time)
Interested in causal interpretation of covariate effects
Coefficients appear attenuated toward zero
Time and covariates are correlated
When NOT to use RTR:
Covariates are temporally uncorrelated (random sampling over time)
Temporal effect is the primary quantity of interest
Prediction is the main goal (not causal inference)
A modified temporal specification with RTR enabled
Hanks, E. M., Schliep, E. M., Hooten, M. B., & Hoeting, J. A. (2015). Restricted spatial regression in practice: geostatistical models, confounding, and robustness under model misspecification. Environmetrics, 26(4), 243-254.
temporal_rw1(), temporal_rw2(), temporal_gp(),
spatial_rsr() for the spatial analog
# Create RTR temporal structure rtr <- temporal_rtr( temporal_rw2("year"), restrict_to = ~ temperature + precipitation ) print(rtr) # Generate data with temporal confounding set.seed(170) years <- 2000:2023 n_per_year <- 4 df <- data.frame( year = rep(years, each = n_per_year), # Temperature increasing over time (confounded with year) temperature = rep(seq(15, 18, length.out = length(years)), each = n_per_year) + rnorm(length(years) * n_per_year, 0, 0.5), count = rpois(length(years) * n_per_year, lambda = 20), effort = rgamma(length(years) * n_per_year, shape = 4, rate = 1) ) # Standard RW2 (may have temporal confounding) fit1 <- tratio( count | effort ~ temperature, data = df, family = ratiod_poisson_gamma(), temporal = temporal_rw2("year"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # RTR to protect temperature coefficient fit2 <- tratio( count | effort ~ temperature, data = df, family = ratiod_poisson_gamma(), temporal = temporal_rtr( temporal_rw2("year"), restrict_to = ~ temperature ), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Compare coefficient estimates summary(fit1) # May be attenuated summary(fit2) # Protected from temporal confounding# Create RTR temporal structure rtr <- temporal_rtr( temporal_rw2("year"), restrict_to = ~ temperature + precipitation ) print(rtr) # Generate data with temporal confounding set.seed(170) years <- 2000:2023 n_per_year <- 4 df <- data.frame( year = rep(years, each = n_per_year), # Temperature increasing over time (confounded with year) temperature = rep(seq(15, 18, length.out = length(years)), each = n_per_year) + rnorm(length(years) * n_per_year, 0, 0.5), count = rpois(length(years) * n_per_year, lambda = 20), effort = rgamma(length(years) * n_per_year, shape = 4, rate = 1) ) # Standard RW2 (may have temporal confounding) fit1 <- tratio( count | effort ~ temperature, data = df, family = ratiod_poisson_gamma(), temporal = temporal_rw2("year"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # RTR to protect temperature coefficient fit2 <- tratio( count | effort ~ temperature, data = df, family = ratiod_poisson_gamma(), temporal = temporal_rtr( temporal_rw2("year"), restrict_to = ~ temperature ), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Compare coefficient estimates summary(fit1) # May be attenuated summary(fit2) # Protected from temporal confounding
Specify a first-order random walk temporal random effect.
RW1 assumes that the difference between consecutive time points is
normally distributed: phi[t] - phi[t-1] ~ N(0, sigma^2).
This is an intrinsic model (improper prior) that captures smooth temporal trends. It is equivalent to an ICAR model on a 1D chain.
temporal_rw1( time_var, group_var = NULL, cyclic = FALSE, shared = TRUE, parameterization = c("centered", "noncentered") )temporal_rw1( time_var, group_var = NULL, cyclic = FALSE, shared = TRUE, parameterization = c("centered", "noncentered") )
time_var |
Name of the time variable in data. Can be a formula
(e.g., |
group_var |
Optional name of grouping variable for panel data.
If provided, separate random walks are estimated for each group. Each
group's walk sums to zero over time, so it carries the group's trend and
not its level; add a group intercept such as |
cyclic |
Logical; if TRUE, the random walk wraps around (first and last time points are neighbors). Useful for seasonal data. |
shared |
Logical; if TRUE (default), temporal effect enters both numerator and denominator linear predictors identically. |
parameterization |
Coordinate the walk is sampled in. |
The RW1 precision matrix Q has the form (for T time points):
This is a rank-deficient matrix (rank T-1), so a sum-to-zero constraint is applied for identifiability.
A ratiod_temporal object
# Create temporal RW1 specification temporal_rw1("year") temporal_rw1("month", cyclic = TRUE) # Simple temporal trend: a random walk shared by count and effort set.seed(123) year <- rep(2010:2021, each = 4) trend <- cumsum(rnorm(12, sd = 0.25))[year - 2009] x <- rnorm(48) df <- data.frame( year = year, x = x, count = rpois(48, exp(3 + 0.3 * x + trend)), effort = rgamma(48, shape = 5, rate = exp(-trend)) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_rw1("year"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Panel data: separate trends per site around each site's level set.seed(124) site <- rep(paste0("site", 1:5), times = 12) s <- match(site, paste0("site", 1:5)) t <- rep(1:12, each = 5) walks <- apply(matrix(rnorm(12 * 5, sd = 0.2), 12, 5), 2, cumsum) eff <- rnorm(5, sd = 0.5)[s] + walks[cbind(t, s)] x <- rnorm(60) df_panel <- data.frame( year = 2014 + t, site = site, x = x, count = rpois(60, exp(2.5 + 0.3 * x + eff)), effort = rgamma(60, shape = 3, rate = exp(-eff)) ) # The per-site walks carry the trends, (1 | site) the site levels fit2 <- tratio( count | effort ~ x + (1 | site), data = df_panel, family = ratiod_poisson_gamma(), temporal = temporal_rw1("year", group_var = "site"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Cyclic for monthly seasonal data set.seed(125) month <- rep(1:12, times = 4) season <- 0.5 * sin(2 * pi * month / 12) x <- rnorm(48) df_monthly <- data.frame( month = month, x = x, count = rpois(48, exp(2.5 + 0.3 * x + season)), effort = rgamma(48, shape = 2, rate = exp(-season)) ) fit3 <- tratio( count | effort ~ x, data = df_monthly, family = ratiod_poisson_gamma(), temporal = temporal_rw1("month", cyclic = TRUE), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) )# Create temporal RW1 specification temporal_rw1("year") temporal_rw1("month", cyclic = TRUE) # Simple temporal trend: a random walk shared by count and effort set.seed(123) year <- rep(2010:2021, each = 4) trend <- cumsum(rnorm(12, sd = 0.25))[year - 2009] x <- rnorm(48) df <- data.frame( year = year, x = x, count = rpois(48, exp(3 + 0.3 * x + trend)), effort = rgamma(48, shape = 5, rate = exp(-trend)) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_rw1("year"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Panel data: separate trends per site around each site's level set.seed(124) site <- rep(paste0("site", 1:5), times = 12) s <- match(site, paste0("site", 1:5)) t <- rep(1:12, each = 5) walks <- apply(matrix(rnorm(12 * 5, sd = 0.2), 12, 5), 2, cumsum) eff <- rnorm(5, sd = 0.5)[s] + walks[cbind(t, s)] x <- rnorm(60) df_panel <- data.frame( year = 2014 + t, site = site, x = x, count = rpois(60, exp(2.5 + 0.3 * x + eff)), effort = rgamma(60, shape = 3, rate = exp(-eff)) ) # The per-site walks carry the trends, (1 | site) the site levels fit2 <- tratio( count | effort ~ x + (1 | site), data = df_panel, family = ratiod_poisson_gamma(), temporal = temporal_rw1("year", group_var = "site"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Cyclic for monthly seasonal data set.seed(125) month <- rep(1:12, times = 4) season <- 0.5 * sin(2 * pi * month / 12) x <- rnorm(48) df_monthly <- data.frame( month = month, x = x, count = rpois(48, exp(2.5 + 0.3 * x + season)), effort = rgamma(48, shape = 2, rate = exp(-season)) ) fit3 <- tratio( count | effort ~ x, data = df_monthly, family = ratiod_poisson_gamma(), temporal = temporal_rw1("month", cyclic = TRUE), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) )
Specify a second-order random walk temporal random effect.
RW2 penalizes deviations from linearity, resulting in smoother trends
than RW1: phi[t] - 2*phi[t-1] + phi[t-2] ~ N(0, sigma^2).
temporal_rw2( time_var, group_var = NULL, cyclic = FALSE, shared = TRUE, parameterization = c("centered", "noncentered") )temporal_rw2( time_var, group_var = NULL, cyclic = FALSE, shared = TRUE, parameterization = c("centered", "noncentered") )
time_var |
Name of the time variable in data. Can be a formula
(e.g., |
group_var |
Optional name of grouping variable for panel data.
If provided, separate random walks are estimated for each group. Each
group's walk sums to zero over time, so it carries the group's trend and
not its level; add a group intercept such as |
cyclic |
Logical; if TRUE, the random walk wraps around (first and last time points are neighbors). Useful for seasonal data. |
shared |
Logical; if TRUE (default), temporal effect enters both numerator and denominator linear predictors identically. |
parameterization |
Coordinate the walk is sampled in. |
RW2 produces smoother curves than RW1 because it penalizes the second derivative (curvature) rather than the first derivative (slope). It requires at least 3 time points.
The precision matrix is rank T-2 (two constraints needed).
A ratiod_temporal object
# Create temporal RW2 specification temporal_rw2("year") # Smooth temporal trend set.seed(126) df <- data.frame( year = rep(2008:2027, each = 2), x = rnorm(40), count = rpois(40, lambda = 25), effort = rgamma(40, shape = 4, rate = 1) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_rw2("year"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) )# Create temporal RW2 specification temporal_rw2("year") # Smooth temporal trend set.seed(126) df <- data.frame( year = rep(2008:2027, each = 2), x = rnorm(40), count = rpois(40, lambda = 25), effort = rgamma(40, shape = 4, rate = 1) ) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_rw2("year"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) )
Specify temporally-varying coefficients for ratio models. TVCs allow regression coefficients to change smoothly over time using a Gaussian process or random walk prior. This captures how effects evolve temporally.
temporal_tvc( time_var, terms = 1, structure = c("rw1", "rw2", "ar1", "iid", "gp"), group_var = NULL, shared = TRUE ) ## S3 method for class 'ratiod_tvc' print(x, ...)temporal_tvc( time_var, terms = 1, structure = c("rw1", "rw2", "ar1", "iid", "gp"), group_var = NULL, shared = TRUE ) ## S3 method for class 'ratiod_tvc' print(x, ...)
time_var |
Name of the time variable in data. Can be a formula
(e.g., |
terms |
Which coefficients should vary over time. Options:
|
structure |
Temporal structure for the varying coefficients:
|
group_var |
Optional name of grouping variable for panel data. |
shared |
Logical; if TRUE (default), TVC effects enter both numerator and denominator. Set to FALSE for process-specific TVCs (triggers warning about potential confounding). |
x |
A ratiod_tvc object |
... |
Ignored |
The TVC model extends the linear predictor:
where:
are global (non-time-varying) coefficients
is the subset of covariates with TVCs
are time-varying adjustments
Each TVC follows a temporal process:
Interpretation: A positive TVC for depth at time t means the depth
effect is stronger at t than the global average. The temporal variance
quantifies how much the effect varies over time.
Identifiability: TVCs are centered (sum-to-zero constraint) to separate the time-varying component from the global effect.
A ratiod_tvc object
spatial_svc() for spatially-varying coefficients,
temporal_rw1() for temporal random effects,
tvc() for extracting TVC posteriors
# Create TVC specification tvc <- temporal_tvc("year", terms = 1) print(tvc) # Generate synthetic data with time-varying effect set.seed(150) df <- data.frame( year = rep(2010:2025, each = 5), x = rnorm(80), count = rpois(80, lambda = 20), effort = rgamma(80, shape = 4, rate = 1) ) # Time-varying intercept (random temporal field) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_tvc("year", terms = 1), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) summary(fit) # Time-varying effect of x fit2 <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_tvc("year", terms = c("(Intercept)", "x"), structure = "rw2"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Extract and plot the time-varying coefficients tvc_effects <- tvc(fit2) plot(tvc_effects, "x")# Create TVC specification tvc <- temporal_tvc("year", terms = 1) print(tvc) # Generate synthetic data with time-varying effect set.seed(150) df <- data.frame( year = rep(2010:2025, each = 5), x = rnorm(80), count = rpois(80, lambda = 20), effort = rgamma(80, shape = 4, rate = 1) ) # Time-varying intercept (random temporal field) fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_tvc("year", terms = 1), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) summary(fit) # Time-varying effect of x fit2 <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_tvc("year", terms = c("(Intercept)", "x"), structure = "rw2"), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Extract and plot the time-varying coefficients tvc_effects <- tvc(fit2) plot(tvc_effects, "x")
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.
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() )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() )
formula |
Model formula specifying the response and predictors.
Two syntax options: Combined (recommended) |
data |
Data frame containing all variables. |
family |
A tulpaRatio family object specifying the distributions:
|
formula_num |
Optional one-sided formula for additional numerator
predictors: |
formula_denom |
Optional formula for denominator. Required if main
formula has single response. Can be |
shared |
Formula specifying shared random effects structure:
|
spatial |
Optional spatial structure specification.
See |
temporal |
Optional temporal structure specification.
See |
spatiotemporal |
Optional spatiotemporal interaction specification.
See |
zi |
Optional zero-inflation specification.
See |
latent |
Optional latent factor specification for unmeasured confounders.
See |
priors |
Prior specification. See |
mode |
Inference mode or backend. Accepts either: Tier names (epistemic guarantees):
Backend names (specific implementation):
Structures each backend fits.
See |
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):
NUTS / HMC (
Variational inference (
Stochastic gradient backends (
|
A ratiod_fit object containing:
Posterior draws matrix
Parsed formula specification
Model family
Original data
Inference backend used
MCMC diagnostics (if applicable)
ratio() for extracting ratio posteriors
ratio_contrast() for comparing ratios between groups
pp_check() for posterior predictive checks
loo::loo() and loo::waic() for model comparison
# 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 ) # 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) # Compare ratios between groups ratio_contrast(fit, ~ season)# 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 ) # 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) # Compare ratios between groups ratio_contrast(fit, ~ season)
Extract posterior distributions of temporally-varying coefficients (TVCs) from a fitted ratio model with TVC specification.
tvc( object, terms = NULL, group = 1L, summary = FALSE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_fit' tvc( object, terms = NULL, group = 1L, summary = FALSE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_tvc_posterior' print(x, ...) ## S3 method for class 'ratiod_tvc_posterior' summary(object, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_tvc_posterior' plot(x, term = 1, type = "ribbon", ...)tvc( object, terms = NULL, group = 1L, summary = FALSE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_fit' tvc( object, terms = NULL, group = 1L, summary = FALSE, probs = c(0.025, 0.5, 0.975), ... ) ## S3 method for class 'ratiod_tvc_posterior' print(x, ...) ## S3 method for class 'ratiod_tvc_posterior' summary(object, probs = c(0.025, 0.5, 0.975), ...) ## S3 method for class 'ratiod_tvc_posterior' plot(x, term = 1, type = "ribbon", ...)
object |
A |
terms |
Which TVC terms to extract. If NULL (default), extracts all. |
group |
Which level of the temporal grouping variable to report, as an index into its levels; a model without temporal groups has the one level 1. |
summary |
Logical; if TRUE, return summary statistics instead of full posterior draws. |
probs |
Quantiles to compute if |
... |
Ignored |
x |
A ratiod_tvc_posterior object |
term |
Which term to plot (name or index). Default: first term. |
type |
Plot type: "ribbon" (default) or "line" |
A ratiod_tvc_posterior object containing:
draws: Array of posterior draws (draws x times x terms)
time_levels: Time point labels
term_names: Names of TVC terms
temporal_tvc(), plot.ratiod_tvc_posterior()
# Generate synthetic data set.seed(160) df <- data.frame( year = rep(2015:2024, each = 5), x = rnorm(50), count = rpois(50, lambda = 18), effort = rgamma(50, shape = 3.5, rate = 1) ) # Fit model with TVC fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_tvc("year", terms = c(1, 2)), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Extract TVC posteriors tvc_post <- tvc(fit) summary(tvc_post) # Plot temporal evolution plot(tvc_post, "x")# Generate synthetic data set.seed(160) df <- data.frame( year = rep(2015:2024, each = 5), x = rnorm(50), count = rpois(50, lambda = 18), effort = rgamma(50, shape = 3.5, rate = 1) ) # Fit model with TVC fit <- tratio( count | effort ~ x, data = df, family = ratiod_poisson_gamma(), temporal = temporal_tvc("year", terms = c(1, 2)), mode = "hmc", control = list(iter = 400, chains = 1, seed = 1) ) # Extract TVC posteriors tvc_post <- tvc(fit) summary(tvc_post) # Plot temporal evolution plot(tvc_post, "x")
Ensures the fit object was created with the expected inference mode. Useful for enforcing mode requirements in downstream analysis.
validate_mode(fit, expected_mode, error = TRUE)validate_mode(fit, expected_mode, error = TRUE)
fit |
A ratiod_fit object |
expected_mode |
Mode that should have been used: a tier ( |
error |
If TRUE (default), error on mismatch. If FALSE, return logical. |
If error = FALSE, returns TRUE/FALSE. Otherwise errors on mismatch.
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) validate_mode(fit, "hmc") validate_mode(fit, "laplace", error = FALSE)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) validate_mode(fit, "hmc") validate_mode(fit, "laplace", error = FALSE)
Compute Widely Applicable Information Criterion (WAIC).
## S3 method for class 'ratiod_fit' waic(x, ...)## S3 method for class 'ratiod_fit' waic(x, ...)
x |
A |
... |
Additional arguments passed to loo::waic |
Uses the same reconstructed pointwise log-likelihood as loo.ratiod_fit();
supported on the hmc, ess, sghmc, and vi backends and conditional on the
fixed-effect linear predictor.
A waic object
sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) # waic() warns when an observation's p_waic exceeds 0.4; loo() is then # the estimate to report waic(fit) loo(fit)sim <- sim_ratiod(n = 100, family = ratiod_poisson_gamma(), seed = 1) fit <- tratio(y_num | y_denom ~ x1 + (1 | group), data = sim$data, family = ratiod_poisson_gamma(), control = list(iter = 1000, chains = 2, seed = 1)) # waic() warns when an observation's p_waic exceeds 0.4; loo() is then # the estimate to report waic(fit) loo(fit)
Zero-Inflated Binomial Specification
zi_binomial(formula = NULL)zi_binomial(formula = NULL)
formula |
Optional formula for the zero-inflation probability. If NULL, uses intercept-only model for P(structural zero). |
A ratiod_zi object
# Intercept-only ZI zi_binomial() # ZI probability depends on a covariate zi_binomial(~ habitat)# Intercept-only ZI zi_binomial() # ZI probability depends on a covariate zi_binomial(~ habitat)
Zero-Inflated Negative Binomial Specification
zi_negbin(formula = NULL)zi_negbin(formula = NULL)
formula |
Optional formula for the zero-inflation probability. If NULL, uses intercept-only model for P(structural zero). |
A ratiod_zi object
# Intercept-only ZI zi_negbin() # ZI probability depends on multiple covariates zi_negbin(~ habitat + season)# Intercept-only ZI zi_negbin() # ZI probability depends on multiple covariates zi_negbin(~ habitat + season)
Zero-Inflated Poisson Specification
zi_poisson(formula = NULL)zi_poisson(formula = NULL)
formula |
Optional formula for the zero-inflation probability. If NULL, uses intercept-only model for P(structural zero). |
A ratiod_zi object
# Intercept-only ZI zi_poisson() # ZI probability depends on a covariate zi_poisson(~ habitat)# Intercept-only ZI zi_poisson() # ZI probability depends on a covariate zi_poisson(~ habitat)
Functions to specify zero-inflation or hurdle model components for ratio models. These handle excess zeros in count data where the standard count distribution (Poisson, negative binomial) doesn't adequately capture the zero proportion.
Note on naming convention: These functions (zi_poisson(), zi_negbin(),
hurdle_poisson(), hurdle_negbin()) are ZI specification functions for
the zi argument in tratio(). They are distinct from the family functions
(ratiod_zipois(), ratiod_zinegbin(), ratiod_hurdle_pois(), ratiod_hurdle_negbin())
which specify the complete model family. Use ZI specifications when you want to
add zero-inflation to an existing family; use ZI family functions when you want
a pre-configured zero-inflated family.
## S3 method for class 'ratiod_zi' print(x, ...)## S3 method for class 'ratiod_zi' print(x, ...)
x |
A ratiod_zi object |
... |
Additional arguments (ignored) |
Zero-Inflated Models: Zero-inflated models assume two data-generating processes:
A process that generates "structural zeros" with probability
A count process that generates counts (including zeros) with probability
The probability of observing y is:
Hurdle Models: Hurdle models treat the process as:
A binary process determining zero vs non-zero
A truncated count process (conditional on y > 0)
When to use:
zi_poisson(): Excess zeros with equidispersed positive counts
zi_negbin(): Excess zeros with overdispersed positive counts
hurdle_poisson(): Zero hurdle with equidispersed positive counts
hurdle_negbin(): Zero hurdle with overdispersed positive counts
Invisibly returns the input
zi_poisson() hurdle_negbin()zi_poisson() hurdle_negbin()
Zero-and-One-Inflated Binomial (ZOIB) Specification
zoib(zi_formula = NULL, oi_formula = NULL)zoib(zi_formula = NULL, oi_formula = NULL)
zi_formula |
Optional formula for the zero-inflation probability. If NULL, uses intercept-only model. |
oi_formula |
Optional formula for the one-inflation probability. If NULL, uses intercept-only model. |
A ratiod_zi object
# Intercept-only for both zero and one inflation zoib()# Intercept-only for both zero and one inflation zoib()