---
title: "Complete Workflows"
author: "Gilles Colling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Complete Workflows}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  eval = identical(Sys.getenv("NOT_CRAN"), "true") ||
    identical(Sys.getenv("IN_PKGDOWN"), "true")
)
```

## Installation

```{r, eval = FALSE}
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

```{r}
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

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

summary(fit_cpue)
```

### Diagnostics

```{r}
# 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

```{r}
# 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

```{r}
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

```{r}
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

```{r}
# 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

```{r}
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

```{r}
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

```{r, message = FALSE}
# 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)
```


# Workflow 4: Temporal Trends

**Goal**: Estimate annual trends in population indices while accounting for temporal autocorrelation.

**Strategy**: Random walk prior on year effects.

### Data

```{r}
set.seed(789)
n_years <- 15
n_sites <- 10

trend_data <- expand.grid(
  year = 1:n_years,
  site = factor(1:n_sites)
)
trend_data$effort <- rgamma(nrow(trend_data), shape = 5, rate = 1)
trend_data$count <- rpois(
  nrow(trend_data),
  lambda = 4 * trend_data$effort * exp(0.3 * sin(trend_data$year / 3))
)
```

### Model with Temporal Structure

A temporal effect shared by numerator and denominator cancels in their ratio, so the walk on the index uses `shared = FALSE`; the warning this prints notes that such a walk can also absorb unmeasured year-to-year drivers of the ratio.

```{r}
fit_trend <- tratio(
  count | effort ~ (1 | site),
  temporal = temporal_rw1(time_var = "year", shared = FALSE),
  data = trend_data,
  family = ratiod_poisson_gamma(),
  control = list(seed = 1)
)

summary(fit_trend)
```

### Extract Trend

```{r}
# Year-level index
trend_by_year <- ratio(fit_trend, by = "year")
summary(trend_by_year)

# Plot trend with uncertainty
plot(temporal(fit_trend))
```


## Summary

| Workflow | Family | Key Feature |
|----------|--------|-------------|
| Fisheries CPUE | `poisson_gamma` | Shared site effects |
| Relative Abundance | `negbin_negbin` | Overdispersed counts |
| Detection Probability | `binomial` | Spatial CAR |
| Temporal Trends | `poisson_gamma` | Random walk |

## See Also

- `vignette("getting-started")` - Basic usage
- `vignette("philosophy")` - Statistical motivation
- `vignette("spatial-temporal")` - Spatial and temporal details