Complete Workflows

Installation

install.packages("pak")                   # development version
pak::pak("gcol33/tulpaRatio")

Overview

This vignette demonstrates complete analysis workflows using tulpaRatio. Each example shows data preparation, model fitting, diagnostics, and inference. See vignette("getting-started") for basic usage and vignette("philosophy") for statistical motivation.

Workflow 1: Fisheries CPUE

Goal: Estimate standardized catch per unit effort (CPUE) while accounting for varying effort and site effects.

Challenge: Effort varies systematically across sites and seasons. Sites with more effort appear to have more precise CPUE, but this precision is illusory if effort correlates with catch conditions.

Strategy: Joint model with shared site effects.

Data

library(tulpaRatio)

# Simulated fisheries data
set.seed(42)
n_sites <- 15
n_seasons <- 4
n_obs <- n_sites * n_seasons * 3

fisheries <- data.frame(
  site = factor(rep(1:n_sites, each = n_seasons * 3)),
  season = factor(rep(rep(c("spring", "summer", "fall", "winter"), each = 3), n_sites)),
  depth = rnorm(n_obs, 50, 20),
  catch = rpois(n_obs, lambda = 20),
  effort_hours = rgamma(n_obs, shape = 4, rate = 0.5)
)

head(fisheries)

Model Fitting

fit_cpue <- tratio(
  catch | effort_hours ~ depth + season + (1 | site),
  data = fisheries,
  family = ratiod_poisson_gamma(),
  control = list(seed = 1)
)

summary(fit_cpue)

Diagnostics

# Check convergence
plot(fit_cpue, type = "trace")

# Posterior predictive checks
pp_check(fit_cpue, type = "dens_overlay", component = "numerator")
pp_check(fit_cpue, type = "dens_overlay", component = "denominator")

Inference

# Site-level CPUE estimates
cpue_by_site <- ratio(fit_cpue, by = "site")
summary(cpue_by_site)

# Seasonal contrasts
ratio_contrast(fit_cpue, ~ season)

# Effect of depth
predict(
  fit_cpue,
  newdata = data.frame(depth = c(30, 50, 70), season = "summer"),
  re_formula = NA
)

Workflow 2: Relative Abundance

Goal: Estimate the proportion of a target species relative to total catch across habitats.

Challenge: Both target count and total count are overdispersed. Using binomial models would underestimate uncertainty.

Strategy: Negative binomial for both counts with shared observer effects.

Data

set.seed(123)
n_plots <- 30
n_visits <- 4

abundance <- data.frame(
  plot = factor(rep(1:n_plots, each = n_visits)),
  habitat = factor(rep(sample(c("forest", "grassland", "wetland"), n_plots, replace = TRUE), each = n_visits)),
  observer = factor(sample(1:10, n_plots * n_visits, replace = TRUE))
)

# Plot and observer effects scale target and total counts alike; the
# target share is higher in grassland and lower in wetland
u <- rnorm(n_plots, sd = 0.4)[abundance$plot] + rnorm(10, sd = 0.4)[abundance$observer]
hab <- c(forest = 0, grassland = 0.4, wetland = -0.3)[as.character(abundance$habitat)]
abundance$target_count <- rnbinom(n_plots * n_visits, size = 3, mu = 15 * exp(u + hab))
abundance$total_count <- rnbinom(n_plots * n_visits, size = 5, mu = 50 * exp(u))

# Ensure target <= total
abundance$target_count <- pmin(abundance$target_count, abundance$total_count)

Model

fit_abundance <- tratio(
  target_count | total_count ~ habitat + (1 | plot) + (1 | observer),
  data = abundance,
  family = ratiod_negbin_negbin(),
  control = list(seed = 1)
)

summary(fit_abundance)

Results

# Habitat effects on relative abundance
ratio_contrast(fit_abundance, ~ habitat)

# Plot-level estimates
abundance_by_plot <- ratio(fit_abundance, by = "plot")
summary(abundance_by_plot)

Workflow 3: Detection Probability

Goal: Estimate detection probability from camera trap data with spatial structure.

Challenge: Detection varies by site characteristics and has spatial autocorrelation.

Strategy: Binomial model with spatial random effects.

Data

set.seed(456)
n_cameras <- 40

# Camera locations on a grid
coords <- expand.grid(x = 1:8, y = 1:5)
coords <- coords[1:n_cameras, ]

# Create adjacency matrix (queen contiguity)
adj <- matrix(0, n_cameras, n_cameras)
for (i in 1:n_cameras) {
  for (j in 1:n_cameras) {
    if (i != j) {
      dist <- sqrt((coords$x[i] - coords$x[j])^2 + (coords$y[i] - coords$y[j])^2)
      if (dist <= sqrt(2) + 0.01) adj[i, j] <- 1
    }
  }
}

camera_data <- data.frame(
  camera = factor(1:n_cameras),
  vegetation = rnorm(n_cameras),
  nights_active = sample(20:30, n_cameras, replace = TRUE)
)

# Detection varies smoothly across the grid
field <- 0.6 * sin(coords$x / 2) - 0.4 * cos(coords$y / 2)
camera_data$detections <- rbinom(
  n_cameras,
  size = camera_data$nights_active,
  prob = plogis(-0.8 + 0.3 * camera_data$vegetation + field)
)

Spatial Model

fit_detection <- tratio(
  detections | nights_active ~ vegetation,
  spatial = spatial_car(adj, level = "group", group_var = "camera"),
  data = camera_data,
  family = ratiod_binomial(),
  mode = "hmc",
  control = list(seed = 1)
)

summary(fit_detection)

Spatial Predictions

# Camera-level detection probabilities
detection_probs <- ratio(fit_detection, by = "camera", summary = TRUE)
detection_probs

# Visualize spatial pattern
plot_map(as.data.frame(detection_probs), coords = coords)