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

## ----bivariate-coscale-simulate-----------------------------------------------
set.seed(1)
n <- 180
dat <- data.frame(
  food = rnorm(n),
  temperature = rnorm(n),
  disturbance = rnorm(n)
)
mu1 <- 0.2 + 0.5 * dat$food + 0.2 * dat$temperature
mu2 <- -0.1 + 0.4 * dat$food
sigma1 <- exp(-0.2 + 0.2 * dat$food - 0.1 * dat$temperature)
sigma2 <- exp(0.1 - 0.2 * dat$food)
eta_rho12 <- -0.1 + 0.4 * dat$disturbance
rho12 <- tanh(eta_rho12)
e1 <- rnorm(n)
e2 <- rho12 * e1 + sqrt(1 - rho12^2) * rnorm(n)
dat$activity <- mu1 + sigma1 * e1
dat$boldness <- mu2 + sigma2 * e2

## ----bivariate-coscale-fit----------------------------------------------------
fit_biv <- drmTMB(
  drm_formula(
    mu1 = activity ~ food + temperature,
    mu2 = boldness ~ food,
    sigma1 = ~ food + temperature,
    sigma2 = ~ food,
    rho12 = ~ disturbance
  ),
  family = c(gaussian(), gaussian()),
  data = dat
)

## ----bivariate-coscale-output-------------------------------------------------
check_drm(fit_biv)
summary(fit_biv)

## ----bivariate-coscale-rho12--------------------------------------------------
coef(fit_biv, "rho12")
head(rho12(fit_biv))

## ----bivariate-coscale-corpairs-----------------------------------------------
corpairs(fit_biv)[
  c("level", "from_response", "to_response", "class", "parameter",
    "estimate", "min", "max", "modelled")
]

## ----bivariate-coscale-interpretation-table-----------------------------------
newdat <- data.frame(
  food = 0,
  temperature = 0,
  disturbance = c(-1, 0, 1)
)
rho_table <- data.frame(
  disturbance = newdat$disturbance,
  sigma_activity = predict(fit_biv, newdata = newdat, dpar = "sigma1"),
  sigma_boldness = predict(fit_biv, newdata = newdat, dpar = "sigma2"),
  rho12 = rho12(fit_biv, newdata = newdat)
)
rho_table$residual_covariance <- with(
  rho_table,
  rho12 * sigma_activity * sigma_boldness
)
rho_table$residual_variance_activity <- rho_table$sigma_activity^2
rho_table$residual_variance_boldness <- rho_table$sigma_boldness^2
round(rho_table, 3)

## ----bivariate-coscale-report-table-------------------------------------------
report_grid <- data.frame(
  food = 0,
  temperature = 0,
  disturbance = c(-1.5, 0, 1.5)
)
report_table <- data.frame(
  disturbance = report_grid$disturbance,
  eta_rho12 = predict(
    fit_biv,
    newdata = report_grid,
    dpar = "rho12",
    type = "link"
  ),
  rho12 = rho12(fit_biv, newdata = report_grid),
  sigma_activity = predict(fit_biv, newdata = report_grid, dpar = "sigma1"),
  sigma_boldness = predict(fit_biv, newdata = report_grid, dpar = "sigma2")
)
report_table$residual_covariance <- with(
  report_table,
  rho12 * sigma_activity * sigma_boldness
)
report_table$residual_variance_activity <- report_table$sigma_activity^2
report_table$residual_variance_boldness <- report_table$sigma_boldness^2
round(report_table, 3)

## ----bivariate-coscale-raw-vs-residual----------------------------------------
round(data.frame(
  raw_activity_boldness_correlation = stats::cor(dat$activity, dat$boldness),
  mean_fitted_residual_rho12 = mean(rho12(fit_biv))
), 3)

## ----bivariate-coscale-rho12-curve, fig.width = 6, fig.height = 4, fig.cap = "Fitted residual `rho12` over disturbance after modelling response-specific means and residual SDs; the ribbon is a 95% Wald confidence interval computed for each supplied row (coverage not certified) and the dotted line marks zero residual correlation.", fig.alt = "Line plot of fitted residual rho12 correlation over disturbance with a pale 95 percent Wald confidence ribbon and a dotted zero line, increasing from negative to positive values and crossing zero near average disturbance."----
rho_grid <- data.frame(
  food = 0,
  temperature = 0,
  disturbance = seq(-2, 2, length.out = 80)
)
rho_pred <- predict_parameters(
  fit_biv,
  newdata = rho_grid,
  dpar = "rho12",
  conf.int = TRUE
)
rho_grid$rho12 <- rho_pred$estimate
rho_grid$conf.low <- rho_pred$conf.low
rho_grid$conf.high <- rho_pred$conf.high

if (requireNamespace("ggplot2", quietly = TRUE)) {
  ggplot2::ggplot(rho_grid, ggplot2::aes(disturbance, rho12)) +
    ggplot2::geom_hline(yintercept = 0, linetype = "dotted", colour = "grey60") +
    ggplot2::geom_ribbon(
      ggplot2::aes(ymin = conf.low, ymax = conf.high),
      fill = "#006D77",
      alpha = 0.18,
      colour = NA
    ) +
    ggplot2::geom_line(linewidth = 0.9, colour = "#006D77") +
    ggplot2::coord_cartesian(ylim = c(-1, 1)) +
    ggplot2::labs(
      x = "Disturbance",
      y = "Fitted residual correlation (rho12)",
      title = "Residual coupling changes with disturbance",
      subtitle = "Ribbon is a 95% Wald interval; dotted line marks rho12 = 0"
    ) +
    ggplot2::theme_minimal(base_size = 11) +
    ggplot2::theme(
      panel.grid.minor = ggplot2::element_blank(),
      plot.title = ggplot2::element_text(face = "bold"),
      plot.subtitle = ggplot2::element_text(colour = "grey30"),
      plot.background = ggplot2::element_rect(fill = "white", colour = NA),
      panel.background = ggplot2::element_rect(fill = "white", colour = NA)
    )
} else {
  old_par <- graphics::par(bg = "white")
  on.exit(graphics::par(old_par), add = TRUE)
  plot(
    rho12 ~ disturbance,
    data = rho_grid,
    type = "l",
    lwd = 2,
    ylim = c(-1, 1),
    xlab = "Disturbance",
    ylab = "Fitted residual correlation (rho12)"
  )
  abline(h = 0, lty = 3, col = "grey60")
}

## ----bivariate-coscale-group-simulate-----------------------------------------
set.seed(98)
n_ID <- 55
n_each <- 7
ID <- factor(rep(seq_len(n_ID), each = n_each))
id_index <- as.integer(ID)
n <- length(ID)

food <- rnorm(n)
disturbance <- rnorm(n)

sigma_activity <- 0.35
sigma_boldness <- 0.45
sd_activity_ID <- 0.65
sd_boldness_ID <- 0.65
rho_ID <- 0.55
rho_residual <- 0.15

u_activity <- rnorm(n_ID)
u_boldness <- rho_ID * u_activity + sqrt(1 - rho_ID^2) * rnorm(n_ID)

e_activity <- rnorm(n)
e_boldness <- rho_residual * e_activity +
  sqrt(1 - rho_residual^2) * rnorm(n)

dat_group <- data.frame(ID = ID, food = food, disturbance = disturbance)
dat_group$activity <- 0.25 + 0.45 * food + 0.15 * disturbance +
  sd_activity_ID * u_activity[id_index] + sigma_activity * e_activity
dat_group$boldness <- -0.15 + 0.35 * food +
  sd_boldness_ID * u_boldness[id_index] + sigma_boldness * e_boldness

## ----bivariate-coscale-group-fit----------------------------------------------
fit_group <- drmTMB(
  bf(
    mu1 = activity ~ food + disturbance + (1 | p | ID),
    mu2 = boldness ~ food + (1 | p | ID),
    sigma1 = ~1,
    sigma2 = ~1,
    rho12 = ~1,
    corpair(ID, level = "group", block = "p", from = "mu1", to = "mu2") ~ 1
  ),
  family = biv_gaussian(),
  data = dat_group,
  control = list(eval.max = 500, iter.max = 500)
)

## ----bivariate-coscale-group-check--------------------------------------------
group_checks <- check_drm(fit_group)
group_checks[
  group_checks$check %in% c(
    "optimizer_convergence",
    "random_effect_sd_boundary",
    "rho12_boundary",
    "biv_mu_random_effect_covariance"
  ),
  c("check", "status", "value", "message")
]

## ----bivariate-coscale-group-corpairs-----------------------------------------
pair_table <- corpairs(
  fit_group
)
pair_table[
  ,
  c(
    "level", "group", "block", "from_dpar", "to_dpar", "class",
    "parameter", "estimate", "conf.status", "interval_source", "modelled"
  )
]

## ----bivariate-coscale-group-corpairs-plot, fig.width = 6.6, fig.height = 3.4, fig.cap = "Point display separating residual `rho12` from the individual-level mean-mean random-intercept correlation; the dotted line marks zero correlation. The group row is not shown with an interval because coverage-backed validation remains planned.", fig.alt = "Point plot comparing residual rho12 correlation with the individual-level mean-mean random-intercept correlation, with hollow point estimates and a dotted zero line on a correlation scale from minus one to one."----
if (requireNamespace("ggplot2", quietly = TRUE)) {
  pair_plot_table <- pair_table
  pair_plot_table$display_label <- ifelse(
    pair_plot_table$level == "residual",
    "Residual\nrho12",
    "Individual\nmean-mean"
  )
  plot_corpairs(pair_plot_table, label = "display_label", facet = NULL) +
    ggplot2::labs(
      title = "Residual and individual-level correlations are separate",
      subtitle = "Point estimates only; dotted line marks zero correlation",
      colour = NULL
    ) +
    ggplot2::theme_minimal(base_size = 11) +
    ggplot2::theme(
      panel.grid.minor = ggplot2::element_blank(),
      plot.title = ggplot2::element_text(face = "bold"),
      plot.subtitle = ggplot2::element_text(colour = "grey30"),
      legend.position = "none"
    )
}

## ----bivariate-coscale-group-covariance-table---------------------------------
group_covariance <- summary(fit_group)$covariance
report_group_covariance <- group_covariance[
  ,
  c(
    "level", "group", "block", "from_response", "to_response", "class",
    "correlation", "from_sd", "to_sd", "covariance", "from_scale",
    "to_scale", "covariance_conf.status"
  )
]
numeric_columns <- vapply(report_group_covariance, is.numeric, logical(1))
report_group_covariance[numeric_columns] <- lapply(
  report_group_covariance[numeric_columns],
  round,
  3
)
report_group_covariance

## ----bivariate-coscale-group-targets------------------------------------------
group_targets <- profile_targets(fit_group)
group_cor_dpar <- 'corpair(ID, level = "group", block = "p", from = "mu1", to = "mu2")'
group_target_names <- c(
  "sd:mu:mu1:(1 | p | ID)",
  "sd:mu:mu2:(1 | p | ID)",
  paste0("fixef:", group_cor_dpar, ":(Intercept)")
)
group_targets[
  match(group_target_names, group_targets$parm),
  c("parm", "profile_ready", "profile_note")
]