## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(drmTMB)

## ----simulate-network---------------------------------------------------------
balanced_partner_tree <- function(n_tip = 4L, prefix = "sp") {
  stopifnot(n_tip >= 2L, log2(n_tip) == floor(log2(n_tip)))
  edges <- matrix(integer(), ncol = 2L)
  edge_lengths <- numeric()
  next_node <- n_tip + 1L

  build <- function(tips) {
    if (length(tips) == 1L) {
      return(tips)
    }
    node <- next_node
    next_node <<- next_node + 1L
    mid <- length(tips) / 2L
    left <- build(tips[seq_len(mid)])
    right <- build(tips[seq.int(mid + 1L, length(tips))])
    edges <<- rbind(edges, c(node, left), c(node, right))
    edge_lengths <<- c(edge_lengths, 1, 1)
    node
  }

  build(seq_len(n_tip))
  structure(
    list(
      edge = edges,
      edge.length = edge_lengths,
      tip.label = paste0(prefix, "_", seq_len(n_tip)),
      Nnode = n_tip - 1L
    ),
    class = "phylo"
  )
}

set.seed(20260531)
plant_tree <- balanced_partner_tree(4, "plant")
pollinator_tree <- balanced_partner_tree(4, "poll")

pair_grid <- expand.grid(
  plant = plant_tree$tip.label,
  pollinator = pollinator_tree$tip.label,
  KEEP.OUT.ATTRS = FALSE
)

plant_clade <- ifelse(pair_grid$plant %in% plant_tree$tip.label[1:2], -1, 1)
pollinator_clade <- ifelse(
  pair_grid$pollinator %in% pollinator_tree$tip.label[1:2],
  -1,
  1
)
pair_signal <- 0.45 * plant_clade * pollinator_clade
names(pair_signal) <- paste(pair_grid$plant, pair_grid$pollinator, sep = ":")

network_dat <- pair_grid[rep(seq_len(nrow(pair_grid)), each = 5), ]
network_dat$floral_density <- rnorm(nrow(network_dat))
eta <- 0.6 +
  0.25 * network_dat$floral_density +
  pair_signal[paste(network_dat$plant, network_dat$pollinator, sep = ":")]
network_dat$visits <- rpois(nrow(network_dat), lambda = exp(eta))
head(network_dat)

## ----fit-pair-count-----------------------------------------------------------
fit_pair <- drmTMB(
  bf(
    visits ~ floral_density +
      phylo_interaction(
        1 | plant:pollinator,
        tree1 = plant_tree,
        tree2 = pollinator_tree
      )
  ),
  family = poisson(link = "log"),
  data = network_dat
)
check_drm(fit_pair)

## ----inspect-pair-count-------------------------------------------------------
pair_term <- "phylo_interaction(1 | plant:pollinator)"
pair_targets <- profile_targets(fit_pair)
pair_targets[
  pair_targets$parm == paste0("sd:mu:", pair_term),
  c("parm", "estimate", "scale", "profile_ready", "profile_note")
]

pair_dev <- ranef(fit_pair, "phylo_interaction_mu")$terms[[pair_term]]
pair_table <- data.frame(
  pair = names(pair_dev),
  conditional_log_rate_deviation = unname(pair_dev),
  row.names = NULL
)
pair_table[order(pair_table$conditional_log_rate_deviation), ]

## ----pair-deviation-figure, fig.width = 7.2, fig.height = 4.2, fig.cap = "Conditional pair deviations for the fitted two-tree Poisson model. Positive values indicate pairs with a higher fitted log visit rate than the fixed effects predict; negative values indicate lower fitted rates. The display is point-only and is not an interval plot.", fig.alt = "A horizontal dot plot of 16 plant-pollinator pairs ordered by conditional log-rate deviation. Points lie on both sides of a vertical zero line, with the caption stating that they are point estimates rather than intervals."----
pair_table <- pair_table[order(pair_table$conditional_log_rate_deviation), ]
graphics::dotchart(
  pair_table$conditional_log_rate_deviation,
  labels = pair_table$pair,
  xlab = "Conditional log-rate deviation",
  main = "Fitted pair deviations"
)
graphics::abline(v = 0, lty = 2, col = "grey40")

## ----independent-pairs, eval = FALSE------------------------------------------
# network_dat$pair_id <- interaction(
#   network_dat$plant,
#   network_dat$pollinator,
#   drop = TRUE
# )
# 
# fit_independent_pair <- drmTMB(
#   visits ~ floral_density + (1 | pair_id),
#   family = poisson(link = "log"),
#   data = network_dat
# )

## ----relmat-pair-route, eval = FALSE------------------------------------------
# fit_pair_relmat <- drmTMB(
#   visits ~ floral_density +
#     relmat(1 | plant_pollinator, Q = Q_pair),
#   family = poisson(link = "log"),
#   data = network_dat
# )

## ----additive-future, eval = FALSE--------------------------------------------
# visits ~ floral_density +
#   phylo(1 | plant, tree = plant_tree) +
#   phylo(1 | pollinator, tree = pollinator_tree) +
#   phylo_interaction(
#     1 | plant:pollinator,
#     tree1 = plant_tree,
#     tree2 = pollinator_tree
#   )