A small count example
The example below uses two small trees and a count response. The same marker can also enter a Gaussian location model or an ordinary NB2 mean model.
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)
#> plant pollinator floral_density visits
#> 1 plant_1 poll_1 -0.4671070 3
#> 1.1 plant_1 poll_1 0.2959469 0
#> 1.2 plant_1 poll_1 0.3627907 8
#> 1.3 plant_1 poll_1 -1.0481202 2
#> 1.4 plant_1 poll_1 -1.7112862 4
#> 2 plant_2 poll_1 -0.4815308 2Fit the pair-level phylogenetic interaction with one marker:
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)
#> <drm_check: 14 checks>
#> ok: 13; notes: 1; warnings: 0; errors: 0
#> check status
#> optimizer_convergence ok
#> optimizer_budget ok
#> finite_objective ok
#> logsigma_clamp_active ok
#> fixed_gradient ok
#> sdreport_status ok
#> hessian_positive_definite ok
#> standard_errors_finite ok
#> standard_errors_inflated ok
#> dropped_rows ok
#> positive_scale ok
#> random_effect_sd_boundary ok
#> interval_reliability_scope note
#> fixed_effect_design_size ok
#> value
#> 0
#> iterations=13; function=19; gradient=14
#> 141.3
#> <NA>
#> max=0.000001537; component=beta_mu[2]
#> ok
#> TRUE
#> range=[0.09184,0.1935]
#> n_inflated=0; max_se=0.1935; median_se=0.1427
#> nobs=80; dropped=0
#> min=1.000
#> min=0.4569; boundary=0.0001000; term=mu.phylo_interaction(1 | plant:pollinator)
#> sd_targets=1; assessed_here=0
#> total_mb=0.006866; max_cols=2; largest=mu; largest_class=matrix; largest_density=1.000
#> message
#> nlminb convergence code is 0.
#> Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#> Objective and log-likelihood are finite.
#> The log(sigma) clamp does not apply to this family.
#> Maximum absolute fixed gradient is <= 0.001; largest component is beta_mu[2].
#> TMB::sdreport() completed successfully.
#> sdreport reports a positive-definite Hessian.
#> All fixed-effect standard errors are finite.
#> No fixed-effect standard error is inflated relative to the others.
#> No rows were dropped by model-frame or known-covariance filtering.
#> All fitted scale values are finite and positive.
#> All fitted random-effect standard deviations are finite, positive, and above the requested lower-boundary warning threshold.
#> This fit has 1 random-effect standard-deviation target. `check_drm()` assesses the fit, not interval reliability: a target can pass every check above and still return an interval that `confint()` warns about at a variance boundary. Before reporting an interval, call `confint()` and read `conf.status` and any boundary warning.
#> Dense fixed-effect design matrices are modest for this fit.The fitted standard deviation is the scale of the latent pair
deviations on the linear predictor scale. For this Poisson model, that
is the log-count scale. Use ranef() for conditional pair
deviations and profile_targets() for the pair-SD estimate
and its uncertainty status. These are different outputs: the deviations
describe this fitted network, whereas the SD says how much pair
deviation the model estimates across comparable pairs.
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")
]
#> parm estimate scale
#> 3 sd:mu:phylo_interaction(1 | plant:pollinator) 0.4568795 response
#> profile_ready profile_note
#> 3 TRUE ready
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 conditional_log_rate_deviation
#> 4 plant_4:poll_1 -0.69065604
#> 10 plant_4:poll_2 -0.55844860
#> 28 plant_4:node6 -0.41636821
#> 20 plant_2:poll_4 -0.39245709
#> 3 plant_3:poll_1 -0.38738209
#> 6 node7:poll_1 -0.35934604
#> 19 plant_1:poll_4 -0.33024318
#> 9 plant_3:poll_2 -0.32961024
#> 14 plant_2:poll_3 -0.30169199
#> 12 node7:poll_2 -0.29601961
#> 23 node6:poll_4 -0.24090009
#> 27 plant_3:node6 -0.23899744
#> 32 plant_2:node7 -0.23138303
#> 30 node7:node6 -0.21845522
#> 31 plant_1:node7 -0.12360352
#> 35 node6:node7 -0.11832885
#> 17 node6:poll_3 -0.11408646
#> 13 plant_1:poll_3 -0.04056738
#> 7 plant_1:poll_2 0.01665583
#> 11 node6:poll_2 0.06059367
#> 29 node6:node6 0.12038470
#> 8 plant_2:poll_2 0.16512519
#> 25 plant_1:node6 0.17735915
#> 26 plant_2:node6 0.18379496
#> 36 node7:node7 0.26034263
#> 5 node6:poll_1 0.30056043
#> 33 plant_3:node7 0.31274844
#> 18 node7:poll_3 0.32572115
#> 15 plant_3:poll_3 0.33320190
#> 2 plant_2:poll_1 0.38625968
#> 24 node7:poll_4 0.45530675
#> 34 plant_4:node7 0.46827945
#> 1 plant_1:poll_1 0.51542162
#> 21 plant_3:poll_4 0.60504343
#> 16 plant_4:poll_3 0.64396153
#> 22 plant_4:poll_4 0.76087682The table is a reportable point-estimate summary of
which observed pairs depart above or below the fixed-effect log rate. It
is not an interval for an individual pair. Before reporting the pair SD
as an uncertainty interval, inspect profile_ready and
profile_note above, then check the exact capability
evidence for this family and target; a computable profile is not by
itself a calibration claim.
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")
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.