## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
has_penguins <- requireNamespace("palmerpenguins", quietly = TRUE)

## ----penguin-data, eval = has_penguins----------------------------------------
penguins <- subset(
  palmerpenguins::penguins,
  complete.cases(species, sex, year, flipper_length_mm, body_mass_g) &
    flipper_length_mm > 0 & body_mass_g > 0
)
penguins$year_c <- penguins$year - mean(penguins$year)
nrow(penguins)

## ----penguin-raw-scale, eval = has_penguins, fig.width = 6, fig.height = 4, fig.cap = "Raw log-scale flipper length and body mass for analysed complete rows. This descriptive relationship is not the fitted residual association."----
plot(
  log(body_mass_g) ~ log(flipper_length_mm), data = penguins,
  pch = 16, col = grDevices::adjustcolor("#087f8c", alpha.f = 0.45),
  xlab = "log flipper length (mm)", ylab = "log body mass (g)"
)
abline(stats::lm(log(body_mass_g) ~ log(flipper_length_mm), data = penguins),
  col = "#14344a", lwd = 2
)

## ----penguin-direct-lognormal, eval = has_penguins, warning = FALSE, message = FALSE----
fit_log <- drmTMB::drmTMB(
  drmTMB::bf(mu1 = flipper_length_mm ~ species + sex + year_c,
             mu2 = body_mass_g ~ species + sex + year_c,
             sigma1 = ~ 1, sigma2 = ~ 1, rho12 = ~ 1),
  family = drmTMB::biv_lognormal(), data = penguins
)

drmTMB::rho12(fit_log)[1]
drmTMB::check_drm(fit_log)

# All three target the log-residual rho12, not a raw-scale correlation.
ci_wald <- confint(fit_log, parm = "rho12", method = "wald")
ci_profile <- confint(
  fit_log, parm = "rho12", method = "profile", profile_engine = "endpoint"
)
ci_bootstrap <- confint(
  fit_log, parm = "rho12", method = "bootstrap", R = 99, seed = 20260724
)
interval_row <- function(x, label) {
  data.frame(
    label = label,
    estimate = unname(drmTMB::rho12(fit_log)[1]),
    lower = x$lower,
    upper = x$upper,
    status = x$conf.status,
    note = x$profile.message,
    stringsAsFactors = FALSE
  )
}
intervals <- rbind(
  interval_row(ci_wald, "Wald"),
  interval_row(ci_profile, "Profile"),
  interval_row(
    ci_bootstrap,
    sprintf(
      "Bootstrap (%d/%d retained)",
      ci_bootstrap$bootstrap.n,
      ci_bootstrap$bootstrap.n + ci_bootstrap$bootstrap.failed
    )
  )
)
bootstrap_total <- ci_bootstrap$bootstrap.n + ci_bootstrap$bootstrap.failed
cat(sprintf(
  "Bootstrap diagnostic: %d/%d full refits retained. Failed refits remain a diagnostic, not hidden precision.\n",
  ci_bootstrap$bootstrap.n, bootstrap_total
))

## ----penguin-intervals, eval = has_penguins, echo = FALSE, fig.width = 7, fig.height = 3.8, fig.cap = "Reported 95% confidence intervals for the fitted direct log-residual association in the penguin model. Each eye spans the reported interval and is centred on its estimate; the taper is a visual interval cue, not a likelihood, sampling density, or posterior distribution. The retained bootstrap count is a diagnostic, not hidden precision.", fig.alt = "Three 95 percent confidence intervals for the direct log-residual correlation rho12. Wald, profile, and bootstrap estimates are all about 0.35. The bootstrap interval is annotated with its retained full-refit count."----
at <- rev(seq_len(nrow(intervals)))
old_par <- par(no.readonly = TRUE)
par(mar = c(4.2, 7.4, 0.8, 2.6), bg = "white")
x_limits <- range(c(0, intervals$lower, intervals$upper)) + c(-0.08, 0.18)
plot(
  NA, xlim = x_limits, ylim = c(0.5, 3.5), yaxt = "n", bty = "n",
  xlab = expression(rho[12]), ylab = ""
)
abline(v = 0, lty = 2, lwd = 1, col = "grey65")

draw_interval_eye <- function(lower, estimate, upper, y) {
  x <- seq(lower, upper, length.out = 101)
  left_width <- max(estimate - lower, .Machine$double.eps)
  right_width <- max(upper - estimate, .Machine$double.eps)
  taper <- ifelse(
    x <= estimate,
    (x - lower) / left_width,
    (upper - x) / right_width
  )
  half_height <- 0.18 * sqrt(pmax(taper, 0))
  polygon(
    c(x, rev(x)), c(y + half_height, rev(y - half_height)),
    col = grDevices::adjustcolor("#087f8c", alpha.f = 0.28), border = NA
  )
  points(
    estimate, y,
    pch = 21, cex = 1.50, lwd = 1.80,
    col = "#087f8c", bg = "white"
  )
}

for (i in seq_len(nrow(intervals))) {
  draw_interval_eye(intervals$lower[i], intervals$estimate[i], intervals$upper[i], at[i])
}
axis(2, at = at, labels = c("Wald", "Profile", "Bootstrap"), las = 1, tick = FALSE)
text(
  intervals$upper[3] + 0.025, at[3],
  labels = sprintf("%d/%d retained", ci_bootstrap$bootstrap.n, bootstrap_total),
  adj = c(0, 0.5), cex = 0.82, col = "#4b5563"
)
par(old_par)

## ----student-example, eval = FALSE--------------------------------------------
# fit_t <- drmTMB(
#   bf(mu1 = activity ~ habitat, mu2 = boldness ~ habitat,
#      sigma1 = ~ 1, sigma2 = ~ 1, nu = ~ 1, rho12 = ~ 1),
#   family = biv_student(), data = behaviour
# )
# 
# rho12(fit_t)

## ----staged-example, eval = FALSE---------------------------------------------
# paired_data <- na.omit(data.frame(
#   bred, offspring, habitat, habitat_score, season
# ))
# 
# assoc <- biv_associate(
#   bf(mu = bred ~ habitat),
#   bf(mu = offspring ~ habitat, sigma = ~ season),
#   family = list(binomial(), nbinom2()),
#   data = paired_data,
#   association = ~ 1
# )
# 
# association(assoc)
# vcov(assoc)
# confint(assoc)
# confint(assoc, type = "eta")