## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
if (!"package:drmTMB" %in% search()) {
  library(drmTMB)
}

## ----one-call-example---------------------------------------------------------
set.seed(20260801)
n <- 160L
habitat_score <- seq(-1.5, 1.5, length.out = n)
shared_tendency <- rnorm(n)

paired_data <- data.frame(
  habitat_score = habitat_score,
  bred = rbinom(
    n,
    size = 1,
    prob = plogis(-0.3 + 0.25 * habitat_score + 0.7 * shared_tendency)
  ),
  offspring = rnbinom(
    n,
    mu = exp(0.7 + 0.2 * habitat_score + 0.6 * shared_tendency),
    size = 2
  )
)

assoc <- biv_associate(
  bf(mu = bred ~ habitat_score),
  bf(mu = offspring ~ habitat_score, sigma = ~ 1),
  family = list(binomial(), nbinom2()),
  data = paired_data,
  association = ~ 1
)

association(assoc)
alpha_se <- sqrt(diag(vcov(assoc)))
alpha_se
suppressWarnings(confint(assoc))
suppressWarnings(confint(assoc, type = "eta"))

## ----staged-diagram, echo = FALSE, fig.width = 8, fig.height = 3.1, fig.alt = "Flow diagram showing one complete paired data set splitting into a binary margin and a count margin in stage 1, then joining in stage 2 where both margins are frozen and eta is estimated.", fig.cap = "The frozen-margin workflow. The two response models are fitted first on the same complete rows. Their fitted probabilities or distributions are held fixed while stage 2 estimates a latent-normal association, eta. For every admitted route, the public Godambe covariance propagates uncertainty from both fitted margins into alpha-scale Wald intervals when fit-specific diagnostics pass."----
old_par <- par(mar = c(0, 0, 1, 0))
plot.new()
plot.window(xlim = c(0, 12), ylim = c(0, 4))
boxes <- list(
  c(0.3, 1.3, 2.7, 2.7), c(3.5, 2.2, 6.2, 3.5),
  c(3.5, 0.5, 6.2, 1.8), c(8.2, 1.3, 11.6, 2.7)
)
cols <- c("#E7F3F5", "#F2E8C9", "#F2E8C9", "#DCEAD7")
for (i in seq_along(boxes)) {
  rect(boxes[[i]][1], boxes[[i]][2], boxes[[i]][3], boxes[[i]][4],
       col = cols[i], border = "#112638", lwd = 1.2)
}
text(1.5, 2.35, "One complete\npaired data set", cex = 0.9)
text(4.85, 3.05, "Stage 1\nBinary margin\nP(breeding = 1)", cex = 0.84)
text(4.85, 1.15, "Stage 1\nCount margin\nE(offspring)", cex = 0.84)
text(9.9, 2.35, "Stage 2\nFreeze both margins\nEstimate eta (or eta_i)", cex = 0.88)
arrows(2.7, 2.35, 3.45, 2.85, length = 0.08, lwd = 1.2)
arrows(2.7, 1.65, 3.45, 1.15, length = 0.08, lwd = 1.2)
arrows(6.2, 2.85, 8.15, 2.45, length = 0.08, lwd = 1.2)
arrows(6.2, 1.15, 8.15, 2.05, length = 0.08, lwd = 1.2)
par(old_par)

## ----association-interval-example, eval = FALSE-------------------------------
# association(assoc)       # bounded latent-normal eta point estimate
# sqrt(diag(vcov(assoc)))  # Godambe SE for alpha
# confint(assoc)           # Wald interval for alpha
# confint(assoc, type = "eta") # transformed interval for constant eta

## ----unresolved-example, eval = FALSE-----------------------------------------
# if (identical(assoc$status, "boundary_unresolved")) {
#   assoc$diagnostics
#   # Report that no association estimate was returned. Do not clip or repair it.
# }

## ----association-slope-example, eval = FALSE----------------------------------
# assoc_by_habitat <- biv_associate(
#   bf(mu = bred ~ habitat_score),
#   bf(mu = offspring ~ habitat_score, sigma = ~ season),
#   family = list(binomial(), nbinom2()),
#   data = paired_data,
#   association = ~ habitat_score
# )
# 
# association(assoc_by_habitat)
# association(assoc_by_habitat, type = "fitted")
# 
# new_habitats <- data.frame(habitat_score = c(-1, 0, 1))
# predict(
#   assoc_by_habitat,
#   newdata = new_habitats,
#   type = "eta",
#   se.fit = TRUE,
#   interval = "confidence"
# )