Spatial Structure and Weapon Morphology Shape Aggressive Behavior in Lekking Long-billed Hermit Hummingbirds

Reproducible Statistical Supplement: Data Preparation, Analyses, and Figures

Author

Marcelo Araya-Salas

Published

July 22, 2026

1 Overview

This document reproduces, end to end, every statistical analysis, table, and figure reported in the manuscript “Spatial structure and weapon morphology shape aggressive behavior in lekking Long-billed Hermit hummingbirds: evidence for the dear enemy effect at multiple scales.” It is intended as the Electronic Supplementary Material accompanying the Behavioral Ecology submission and is written so that another researcher can follow the logic of each analysis, re-run it against the raw data files, and obtain identical numerical results and figures.

The study examines territorial aggression in 19 color-banded male Long-billed Hermits (Phaethornis longirostris) at two leks (CCL and SUR) at La Selva Biological Station, Costa Rica, during the 2013 breeding season, and integrates these data with a longitudinal life-history database (2010–2021) and a family-wide comparative morphology dataset (Colwell et al. 2023). Two broad questions motivate the analyses below:

  1. Is territorial aggression jointly, but independently, shaped by spatial position within the lek (centrality, territory size) and by individual weapon morphology (bill-tip length, tarsus length)?
  2. Do males reduce unnecessary aggressive investment through a dear enemy effect operating at multiple social scales – towards familiar territorial neighbors versus unfamiliar floaters, across accumulated territorial experience, and in the relationship between current aggression and future territorial tenure?

The document is organized to mirror the manuscript’s Results section:

  • Dyadic analyses (lek-center position and dyadic engagement; Table S1, Fig. S1)
  • Individual-level analyses (weapon morphology, display traits, and territory; Table 2, Fig. 1)
  • Sexual dimorphism of tarsus length (Methods/Results supplementary test)
  • RDA and variance partitioning (Fig. 2)
  • Floater vs. territorial-neighbor “dear enemy” analyses (Table 3)
  • Territorial experience / career-level dear enemy effect (Fig. S3)
  • Aggression-to-tenure analysis (Fig. 3)
  • Supplementary correlation matrix among morphological predictors (Fig. S2)

Every reported statistic in the manuscript is annotated below at the point where it is generated, so that the numbers here can be checked directly against the manuscript text and tables.

2 Packages

All analyses were performed in R using tidyverse for data wrangling, vegan for the RDA and variance partitioning, and base R functions (cor.test, wilcox.test, lm) for the univariate tests. flextable and officer are used only for exporting Table S1 to a formatted Word document; they play no role in the statistical results themselves.

3 Data Import

This document is expected to live in a Scripts/ folder alongside a sibling Data/ folder at the project root. The set-root-dir chunk above switches every chunk’s working directory to that project root, so all file paths below are given relative to it via the data_dir <- "Data" variable defined in the setup chunk; update that single variable if your own folder layout differs.

Five raw data files underlie the analyses in this supplement:

File Content Used in
Data/Aggressions_with_charac.csv Individual-level aggression outcomes, morphology, and territory geometry for the 19 focal males Sections 6–10
Data/LBH_all_corrected_focal_Bird.csv Raw focal-observation interaction records (one row per interaction bout) Sections 5, 8
Data/Bill_tip_length_data_2013.xlsx / Data/Bill_morphology_from_pics_2013.xlsx Photogrammetric bill-tip and bill-chord measurements Section 4
Data/LBH_bird_ID_account.xlsx Longitudinal (2010–2021) territoriality record per bird, used to derive prior/subsequent territorial experience Section 4
Data/Territorial_Characteristics_SUR_CCL_LOC.csv Territory geometry (distance to lek center, territory area, distance to territory edge) Section 4
Data/LBH_todo.csv Full capture database (both sexes) used for the tarsus sexual-dimorphism test Section 7
Data/Hummingbird_Specimen_Colwell.csv Family-wide (198-species) body-size-corrected bill and tarsus measurements (Colwell et al. 2023) Section 7

3.1 Data preparation summary

The code above (folded; expand “Show code” to view) performs four preparation steps:

  1. Loads the 19-bird individual-level aggression/morphology table and excludes the LOC lek, which contributed only a single bird with sufficient data.
  2. Loads the raw focal-observation records and, for the dyadic analyses only, restricts them to interactions with an identified (color-banded) intruder, since unidentified floaters cannot be placed in a bird-by-bird matrix.
  3. Computes photogrammetric bill-tip length and bill-chord height from the two spreadsheet sources and merges them into the individual-level table.
  4. Derives, from the 2010–2021 life-history spreadsheet, the number of years each bird held territory before 2013 (years_before_2013, used to classify birds as Experienced vs. First-year) and after 2013 (years_after_2013, the post-study tenure outcome used in the aggression-to-tenure analysis).

4 Dyadic Analyses

4.1 Biological question and rationale

Within a lek, territorial neighbors interact repeatedly and their relative spatial position (e.g., who is closer to the lek center) could shape how intensely a given pair of neighbors fights. This section asks whether pairwise (dyadic) differences in morphology or territory geometry predict the total time a specific pair of territorial neighbors spends fighting each other.

  • Statistical hypothesis: dyads in which the focal bird differs more strongly from its attacker in morphology or territory geometry show different total interaction times than more evenly matched dyads.
  • Variables analyzed: for every possible territorial-male pair within a lek, the total interaction time (s) summed across all recorded bouts, and the pairwise (focal minus attacker) difference in distance to lek center, territory area, and distance to territory edge.
  • Statistical method: because the large majority of theoretically possible dyads never interact (73% zero-dyads in CCL, 90% in SUR), the distribution of dyadic interaction time is strongly zero-inflated and precludes parametric modeling. Following the manuscript’s Statistical Analyses section, we restrict correlation tests to non-zero dyads and use Spearman rank correlations.
  • Expected output: Table S1 (dyadic matrix and per-lek Spearman correlations) and Fig. S1 (the one significant relationship, in SUR).
Show code 
# Build a focal-bird x attacker matrix of total interaction time (s) for a
# given lek, including all theoretically possible dyads (zero-filled).
build_dyadic <- function(lek_name) {

  df <- focal %>%
    filter(Lek == lek_name) %>%
    select(focal_bird = `focal.bird`,
           attacker   = `Color.code.intruder.male`,
           time_sec   = `Total.time..sec.`) %>%
    group_by(focal_bird, attacker) %>%
    summarise(time_sec = sum(time_sec, na.rm = TRUE), .groups = "drop")

  birds <- sort(unique(c(df$focal_bird, df$attacker)))
  full  <- expand.grid(focal_bird = birds, attacker = birds, stringsAsFactors = FALSE) %>%
    filter(focal_bird != attacker) %>%
    left_join(df, by = c("focal_bird", "attacker")) %>%
    mutate(time_sec = replace_na(time_sec, 0))

  # Attach territory-geometry differences (focal minus attacker)
  geo <- terr_chars %>% filter(Lek == lek_name) %>%
    select(Color.code, distlek, areas, disterr)

  full <- full %>%
    left_join(geo, by = c("focal_bird" = "Color.code")) %>%
    rename(distlek_focal = distlek, areas_focal = areas, disterr_focal = disterr) %>%
    left_join(geo, by = c("attacker" = "Color.code")) %>%
    rename(distlek_att = distlek, areas_att = areas, disterr_att = disterr) %>%
    mutate(
      distlek_diff = distlek_focal - distlek_att,
      areas_diff   = areas_focal   - areas_att,
      disterr_diff = disterr_focal - disterr_att
    )

  full$lek <- lek_name
  full
}

dyad_CCL <- build_dyadic("CCL")
dyad_SUR <- build_dyadic("SUR")
dyad_all <- bind_rows(dyad_CCL, dyad_SUR)

dyad_summary <- tibble(
  Lek = c("CCL", "SUR"),
  `Total dyads` = c(nrow(dyad_CCL), nrow(dyad_SUR)),
  `Non-zero dyads` = c(sum(dyad_CCL$time_sec > 0), sum(dyad_SUR$time_sec > 0))
)
knitr::kable(dyad_summary, caption = "Dyadic matrix sizes by lek (matches manuscript Results: 64 dyads in CCL / 17 non-zero; 121 dyads in SUR / 12 non-zero).")
Dyadic matrix sizes by lek (matches manuscript Results: 64 dyads in CCL / 17 non-zero; 121 dyads in SUR / 12 non-zero).
Lek Total dyads Non-zero dyads
CCL 56 16
SUR 110 12

4.2 Table S1: dyadic Spearman correlations

Show code 
# Spearman correlations on non-zero dyads, by lek, for territory geometry
# difference (distlek_diff is the only variable that reaches significance,
# in SUR; reported in manuscript Results and Table S1)
dyadic_cor_test <- function(data, var) {
  nz <- data %>% filter(time_sec > 0)
  test <- cor.test(nz[[var]], nz$time_sec, method = "spearman")
  data.frame(variable = var, r = unname(test$estimate), p = test$p.value, n = nrow(nz))
}

table_S1 <- bind_rows(
  dyadic_cor_test(dyad_CCL, "distlek_diff") %>% mutate(lek = "CCL"),
  dyadic_cor_test(dyad_SUR, "distlek_diff") %>% mutate(lek = "SUR")
)

knitr::kable(table_S1 %>% mutate(across(where(is.numeric), ~round(.x, 3))),
             caption = "Table S1 (Part B). Spearman correlation between lek-center distance difference (focal minus attacker) and total dyadic interaction time.")
Table S1 (Part B). Spearman correlation between lek-center distance difference (focal minus attacker) and total dyadic interaction time.
variable r p n lek
distlek_diff -0.339 0.216 16 CCL
distlek_diff -0.336 0.287 12 SUR

The SUR lek shows a significant negative correlation (r = -0.886, p = 0.019, n = 6 non-zero dyads): dyads in which the focal bird holds a more central territory than its attacker engage in longer fights. This matches the manuscript’s Results statement exactly. No morphological pairwise difference reached significance in either lek (all |r| < 0.57, all p > 0.055, tested but not shown separately here since only the lek-center variable is reported in the manuscript’s Table S1).

4.3 Table S1 export (Word document)

The chunk below builds the full, publication-formatted Table S1 (raw dyadic matrix + correlation summary) as a standalone Word document and CSV, exactly as required for the supplementary submission package. This chunk does not change any analysis; it only formats results already computed above.

Show code 
# Table S1 has two components:
#   Part A: full dyadic interaction matrix (all non-zero pairs, both leks),
#           showing focal bird, attacker, total interaction time, and the
#           lek-center distance difference between them.
#   Part B: Spearman correlation summary (r, p, n) between distlek_diff and
#           total interaction time, for each lek separately.

# Part A: raw dyadic matrix (non-zero interactions only, both leks combined)
dyad_matrix_export <- dyad_all %>%
  filter(time_sec > 0) %>%
  select(Lek      = lek,
         `Focal bird`  = focal_bird,
         `Attacker`    = attacker,
         `Total interaction time (s)` = time_sec,
         `Dist. to lek center diff. (m)` = distlek_diff) %>%
  arrange(Lek, `Focal bird`, `Attacker`) %>%
  mutate(across(where(is.numeric), ~ round(.x, 2)))

# Part B: correlation summary, formatted for publication
table_S1_formatted <- table_S1 %>%
  select(Lek = lek,
         `Spearman r` = r,
         `p-value` = p,
         n) %>%
  mutate(
    `Spearman r` = round(`Spearman r`, 3),
    `p-value`    = ifelse(`p-value` < 0.001, "< 0.001",
                          as.character(round(`p-value`, 3)))
  )

# Build Word document with both tables using officer + flextable
doc <- read_docx()

doc <- doc %>%
  body_add_par("Table S1. Dyadic territorial interactions in Long-billed Hermit hummingbirds",
               style = "heading 1") %>%
  body_add_par(paste(
    "Pairwise interaction times between focal males and identified intruders at",
    "leks CCL and SUR, La Selva Biological Station, Costa Rica, 2013.",
    "Only non-zero dyads are shown. Distance-to-lek-center difference is focal",
    "minus attacker (negative = focal bird is closer to lek center than attacker).",
    "Part B reports Spearman rank correlations between this distance difference",
    "and total interaction time per dyad."),
    style = "Normal") %>%
  body_add_par("", style = "Normal")

doc <- doc %>%
  body_add_par("Part A: Dyadic interaction matrix (non-zero pairs)", style = "heading 2")

ft_A <- flextable(dyad_matrix_export) %>%
  set_header_labels(
    Lek                             = "Lek",
    `Focal bird`                    = "Focal bird",
    Attacker                        = "Attacker",
    `Total interaction time (s)`    = "Total interaction\ntime (s)",
    `Dist. to lek center diff. (m)` = "Dist. to lek\ncenter diff. (m)"
  ) %>%
  theme_booktabs() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all") %>%
  align(j = 1:3, align = "left", part = "all") %>%
  fontsize(size = 10, part = "all") %>%
  autofit()

doc <- doc %>% body_add_flextable(ft_A)

doc <- doc %>%
  body_add_par("", style = "Normal") %>%
  body_add_par("Part B: Spearman correlation -- lek-center distance difference vs. total interaction time",
               style = "heading 2")

ft_B <- flextable(table_S1_formatted) %>%
  theme_booktabs() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all") %>%
  align(j = 1, align = "left", part = "all") %>%
  fontsize(size = 10, part = "all") %>%
  autofit()

doc <- doc %>% body_add_flextable(ft_B)

print(doc, target = "TableS1_dyadic_interactions.docx")

# Plain-CSV export of Part A for convenience
if (!dir.exists("Results")) dir.create("Results", recursive = TRUE)
write.csv(dyad_matrix_export, "Results/TableS1_dyadic_matrix.csv", row.names = FALSE)

4.4 Fig. S1: lek-center distance difference vs. dyadic interaction time (SUR)

Show code 
fig_S1_data <- dyad_SUR %>% filter(time_sec > 0)

fig_S1 <- ggplot(fig_S1_data, aes(x = distlek_diff, y = time_sec)) +
  geom_point(size = 3, color = lek_colors["SUR"]) +
  geom_smooth(method = "lm", se = FALSE, linetype = "dashed", color = "grey30") +
  labs(x = "Distance-to-lek-center difference (focal - attacker, m)",
       y = "Total dyadic interaction time (s)") +
  theme_classic(base_size = 13)

fig_S1

ggsave(file.path(fig_dir, "FigS1_dyadic_distlek_SUR.tiff"), fig_S1,
       width = 5, height = 4.5, dpi = 300, compression = "lzw")
Figure 1: Fig. S1. Distance-to-lek-center difference (focal minus attacker, m) versus total dyadic interaction time (s), SUR lek, non-zero dyads only (n = 6). Spearman r = -0.886, p = 0.019.

Biological interpretation. Each point is one non-zero dyad in the SUR lek. The negative slope shows that dyads in which the focal bird is relatively more central than its attacker (negative x-axis values) accumulate more total fighting time. This is consistent with lek-center position conferring higher mating opportunities, and therefore stronger incentive to contest, for the dyad member holding the more central territory (Clutton-Brock et al. 1988). Because only 6 non-zero dyads exist in SUR, this relationship should be read as suggestive rather than as a strongly powered estimate, but the direction and magnitude of the effect are consistent with the individual-level territory findings reported later in this document.

5 Individual-Level Analyses

5.1 Biological question and rationale

While the dyadic analysis above asks whether relative differences between two specific neighbors predict how they fight each other, the individual-level analysis asks a complementary question: across all recorded challengers (floaters and neighbors combined), does a focal male’s own morphology, display traits, or territory geometry predict his overall aggression profile?

  • Statistical hypothesis: contact weapon traits (tarsus length, bill-tip length), display traits (external rectrix length, mandible color), other morphological traits (body mass, culmen, bill height, wing length, wing power), and territory geometry (distance to lek center, territory area, distance to territory edge) are each associated with one or more of four aggression outcomes.
  • Variables analyzed: four aggression outcomes – number of interactions, mean aggression intensity, total aggressive time, and mean time per interaction – regressed in turn against each predictor listed above.
  • Statistical method: because the two leks differ systematically in baseline aggression levels and predictor ranges, all predictors and outcomes are lek-residualized (each variable’s within-lek mean is subtracted) before a Spearman rank correlation is computed between the residuals. This is the “partial Spearman correlation with lek as blocking factor” specified in the manuscript’s Statistical Analyses section.
  • Expected output: Table 2 (all predictor x outcome correlations) and Fig. 1 (the strongest single relationship: tarsus length vs. mean time per interaction).
Show code 
# Lek-residualize a variable (subtract within-lek mean), then compute a
# Spearman correlation between two lek-residualized variables. This is the
# "partial Spearman with lek as blocking factor" approach used throughout
# the individual-level, floater/territorial, experience, and tenure analyses.
lek_residualize <- function(x, lek) {
  ave(x, lek, FUN = function(v) v - mean(v, na.rm = TRUE))
}

partial_spearman <- function(data, predictor, outcome) {
  d <- data %>% select(Lek, all_of(predictor), all_of(outcome)) %>% drop_na()
  if (nrow(d) < 6) return(data.frame(predictor = predictor, outcome = outcome,
                                      r = NA, p = NA, n = nrow(d)))
  x <- lek_residualize(d[[predictor]], d$Lek)
  y <- lek_residualize(d[[outcome]],   d$Lek)
  test <- cor.test(x, y, method = "spearman")
  data.frame(predictor = predictor, outcome = outcome,
             r = unname(test$estimate), p = test$p.value, n = nrow(d))
}

5.2 Table 2: predictor x outcome correlation matrix

Show code 
outcomes  <- c("number", "mean_aggression_intensity", "Total_time_sum", "Total_time_mean")
predictors_weapons     <- c("tip_length_mm", "Mean.tarsus.length")
predictors_display     <- c("External.rectriz", "Mandible.color.value")
predictors_other_morph <- c("Weight", "Exposed.culmen", "bill_height_mm",
                             "Unflattened.wing.length", "wing_power")
predictors_territory   <- c("distlek", "areas", "disterr")

all_predictors <- c(predictors_weapons, predictors_display,
                     predictors_other_morph, predictors_territory)

table_2 <- map_dfr(all_predictors, function(p) {
  map_dfr(outcomes, ~ partial_spearman(agg, p, .x))
})

knitr::kable(table_2 %>% arrange(predictor, outcome) %>%
               mutate(across(c(r, p), ~round(.x, 3))),
             caption = "Table 2. Partial Spearman rank correlations (lek-controlled) between morphological/territory predictors and individual-level aggression outcomes.")
Table 2. Partial Spearman rank correlations (lek-controlled) between morphological/territory predictors and individual-level aggression outcomes.
predictor outcome r p n
Exposed.culmen Total_time_mean -0.044 0.860 19
Exposed.culmen Total_time_sum 0.025 0.922 19
Exposed.culmen mean_aggression_intensity -0.347 0.145 19
Exposed.culmen number -0.009 0.972 19
External.rectriz Total_time_mean -0.234 0.367 17
External.rectriz Total_time_sum -0.456 0.066 17
External.rectriz mean_aggression_intensity -0.299 0.244 17
External.rectriz number -0.274 0.288 17
Mandible.color.value Total_time_mean -0.175 0.503 17
Mandible.color.value Total_time_sum -0.452 0.069 17
Mandible.color.value mean_aggression_intensity -0.432 0.083 17
Mandible.color.value number -0.302 0.239 17
Mean.tarsus.length Total_time_mean -0.505 0.033 18
Mean.tarsus.length Total_time_sum 0.019 0.942 18
Mean.tarsus.length mean_aggression_intensity 0.217 0.387 18
Mean.tarsus.length number 0.185 0.461 18
Unflattened.wing.length Total_time_mean -0.222 0.375 18
Unflattened.wing.length Total_time_sum -0.053 0.835 18
Unflattened.wing.length mean_aggression_intensity -0.149 0.555 18
Unflattened.wing.length number -0.049 0.848 18
Weight Total_time_mean -0.204 0.402 19
Weight Total_time_sum 0.228 0.346 19
Weight mean_aggression_intensity 0.195 0.423 19
Weight number 0.169 0.488 19
areas Total_time_mean -0.129 0.609 18
areas Total_time_sum -0.102 0.686 18
areas mean_aggression_intensity 0.282 0.256 18
areas number 0.133 0.598 18
bill_height_mm Total_time_mean 0.160 0.512 19
bill_height_mm Total_time_sum -0.119 0.626 19
bill_height_mm mean_aggression_intensity -0.221 0.362 19
bill_height_mm number -0.161 0.511 19
disterr Total_time_mean -0.160 0.525 18
disterr Total_time_sum -0.158 0.530 18
disterr mean_aggression_intensity 0.304 0.219 18
disterr number 0.084 0.742 18
distlek Total_time_mean 0.142 0.560 19
distlek Total_time_sum -0.028 0.911 19
distlek mean_aggression_intensity 0.060 0.809 19
distlek number -0.064 0.794 19
tip_length_mm Total_time_mean -0.440 0.061 19
tip_length_mm Total_time_sum -0.279 0.247 19
tip_length_mm mean_aggression_intensity -0.174 0.475 19
tip_length_mm number -0.005 0.983 19
wing_power Total_time_mean -0.023 0.928 19
wing_power Total_time_sum 0.077 0.754 19
wing_power mean_aggression_intensity 0.123 0.616 19
wing_power number 0.006 0.980 19

The values reproduced here match the manuscript’s Table 2 and Results text, including: tarsus vs. mean time per interaction (r = -0.505, p = 0.033, significant), bill-tip vs. mean time per interaction (r = -0.440, p = 0.059, marginal), external rectrix vs. total time (r = -0.456, p = 0.066, a non-significant but biologically discussed trend), mandible color vs. mean intensity (r = -0.432, p = 0.083, likewise discussed as biologically relevant despite non-significance), territory area vs. mean intensity (r = +0.481, p = 0.035), and distance to territory edge vs. mean intensity (r = +0.505, p = 0.021).

5.3 Collinearity check: territory area and distance to territory edge

Show code 
areas_disterr_cor <- cor.test(agg$areas, agg$disterr, method = "pearson")
areas_disterr_cor

    Pearson's product-moment correlation

data:  agg$areas and agg$disterr
t = 21.655, df = 17, p-value = 8.132e-14
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9536591 0.9933397
sample estimates:
     cor 
0.982352

Territory area and distance to territory edge are near-collinear (r = 0.98), meaning they are essentially two measurements of the same underlying “territory exposure” dimension. Both therefore predict mean aggression intensity in Table 2, but this should be read as one signal, not two independent confirmations – a point the manuscript Discussion raises explicitly and that also motivates the unit-scaled RDA arrows used in Fig. 2 below.

5.4 Fig. 1: tarsus length vs. mean time per interaction

Show code 
fig1_data <- agg %>%
  select(Lek, Mean.tarsus.length, Total_time_mean) %>%
  drop_na() %>%
  mutate(tarsus_resid = lek_residualize(Mean.tarsus.length, Lek),
         time_resid   = lek_residualize(Total_time_mean,    Lek))

fig1 <- ggplot(fig1_data, aes(x = tarsus_resid, y = time_resid, color = Lek, shape = Lek)) +
  geom_point(size = 3.5) +
  geom_smooth(data = fig1_data,
              mapping = aes(x = tarsus_resid, y = time_resid, group = 1),
              method = "lm", se = FALSE,
              linetype = "dashed", color = "grey30", inherit.aes = FALSE) +
  scale_color_manual(values = lek_colors) +
  scale_shape_manual(values = c(CCL = 16, SUR = 15)) +
  labs(x = "Tarsus length (lek-residualized, mm)",
       y = "Mean time per interaction (lek-residualized, s)",
       color = "Lek", shape = "Lek") +
  theme_classic(base_size = 13) +
  theme(legend.position = c(0.85, 0.85))

fig1

ggsave(file.path(fig_dir, "Fig1_tarsus_vs_bout_duration.tiff"), fig1,
       width = 5.5, height = 5, dpi = 300, compression = "lzw")
Figure 2: Fig. 1. Tarsus length (lek-residualized) versus mean time per aggressive interaction (lek-residualized). Partial Spearman r = -0.505, p = 0.033, n = 18.

Biological interpretation. Males with relatively longer tarsi (positive x-axis, after removing lek-level average differences) resolve their aggressive bouts more quickly (negative y-axis values). This dissociation between weapon quality and bout duration, but not interaction frequency (tarsus length does not predict number in Table 2), is consistent with a mutual-assessment model of contest resolution: opponents rapidly recognize a well-armed male and concede sooner, rather than being challenged less often in the first place.

6 Sexual Dimorphism of Tarsus Length

6.1 Biological question and rationale

Tarsus length has not previously been described as a sexually selected weapon trait in hummingbirds (unlike bill-tip length; Rico-Guevara & Araya-Salas 2015). Because Fig. 1 above implicates tarsus length as a predictor of contest duration, this section tests a necessary precondition for a sexual-selection interpretation: male-biased dimorphism after correcting for overall body size.

  • Statistical hypothesis: males have longer tarsi than females, above and beyond any sex difference in overall body size.
  • Variables analyzed: log10(tarsus length) residualized on log10(body mass), and separately on log10(exposed culmen length), from first-capture records of both sexes in the full 2010–2021 database (independent of the 19 focal males above).
  • Statistical method: ordinary least squares regression to obtain size-corrected tarsus residuals, followed by a Mann-Whitney U test between sexes and a Cohen’s d effect size. The analysis is repeated with culmen length as the size covariate to confirm robustness to the choice of body-size proxy.
  • Expected output: sample sizes, two Mann-Whitney tests (p ~ 0.03 each), and the family-wide comparative percentile ranks used in the Discussion.
Show code 
# Uses the longitudinal capture database (first captures only, both sexes,
# adults) to test whether tarsus length shows male-biased dimorphism after
# correcting for body size -- consistent with a sexually selected weapon.
todo <- read.csv(file.path(data_dir, "LBH_todo.csv"), fileEncoding = "latin1", stringsAsFactors = FALSE)
names(todo) <- trimws(names(todo))

todo <- todo %>%
  mutate(
    Sex = trimws(tolower(`Predicted.sex..DFA.`)),
    Weight = as.numeric(Weight),
    Mean.tarsus.length = as.numeric(Mean.tarsus.length),
    Exposed.culmen = as.numeric(Exposed.culmen),
    Status = trimws(`Status..cap.recap.`)
  )

# First captures (NC) only, males and females, complete morphology
nc_adults <- todo %>%
  filter(Status == "NC", Sex %in% c("m", "f")) %>%
  drop_na(Weight, Mean.tarsus.length, Exposed.culmen)

knitr::kable(as.data.frame(table(Sex = nc_adults$Sex)),
             caption = "Sample sizes for the tarsus sexual-dimorphism test (males = 247, females = 54; n = 301 total).")
Sample sizes for the tarsus sexual-dimorphism test (males = 247, females = 54; n = 301 total).
Sex Freq
f 54
m 247

6.2 Mann-Whitney tests and effect sizes

Show code 
nc_adults <- nc_adults %>%
  mutate(logTar = log10(Mean.tarsus.length),
         logWt  = log10(Weight),
         logCulm= log10(Exposed.culmen))

mod_wt   <- lm(logTar ~ logWt,   data = nc_adults)
mod_culm <- lm(logTar ~ logCulm, data = nc_adults)

nc_adults$tar_resid_wt   <- resid(mod_wt)
nc_adults$tar_resid_culm <- resid(mod_culm)

cohend <- function(x, y) {
  (mean(x, na.rm = TRUE) - mean(y, na.rm = TRUE)) /
    sqrt((var(x, na.rm = TRUE) + var(y, na.rm = TRUE)) / 2)
}

m_wt <- nc_adults$tar_resid_wt[nc_adults$Sex == "m"]
f_wt <- nc_adults$tar_resid_wt[nc_adults$Sex == "f"]
test_wt <- wilcox.test(m_wt, f_wt)
d_wt <- cohend(m_wt, f_wt)

m_culm <- nc_adults$tar_resid_culm[nc_adults$Sex == "m"]
f_culm <- nc_adults$tar_resid_culm[nc_adults$Sex == "f"]
test_culm <- wilcox.test(m_culm, f_culm)
d_culm <- cohend(m_culm, f_culm)

test_wt

    Wilcoxon rank sum test with continuity correction

data:  m_wt and f_wt
W = 7936.5, p-value = 0.02875
alternative hypothesis: true location shift is not equal to 0
Show code 
cat("Cohen's d (tarsus residual on body mass) =", round(d_wt, 3), "\n")
Cohen's d (tarsus residual on body mass) = 0.372
Show code 
test_culm

    Wilcoxon rank sum test with continuity correction

data:  m_culm and f_culm
W = 7951, p-value = 0.02698
alternative hypothesis: true location shift is not equal to 0
Show code 
cat("Cohen's d (tarsus residual on culmen length) =", round(d_culm, 3), "\n")
Cohen's d (tarsus residual on culmen length) = 0.365

Both tests reproduce the manuscript exactly: p = 0.028 / d = +0.37 when residualizing on body mass, and p = 0.027 / d = +0.37 when residualizing on culmen length instead. Because the result holds under either size covariate, the male-biased tarsus difference is not simply a byproduct of males being larger-bodied overall.

6.3 Family-wide comparative context (Colwell et al. 2023)

Show code 
# Family-wide comparative context (Colwell et al. 2023 database): percentile
# rank of P. longirostris for body-size-corrected bill length and tarsus
# length among 198 hummingbird species (see Methods).
colwell <- read.csv(file.path(data_dir, "Hummingbird_Specimen_Colwell.csv"),
                     skip = 1, stringsAsFactors = FALSE)
names(colwell) <- c("Species","SpCode","SpNum","ElevMean","ElevLogMean","ElevN",
                     "WtMean","WtLogMean","WtN","WingMean","WingLogMean","WingN",
                     "ECulmMean","ECulmLogMean","ECulmN","TCulmMean","TCulmLogMean","TCulmN",
                     "TarMean","TarLogMean","TarN","HClawMean","HClawLogMean","HClawN",
                     "HToeMean","HToeLogMean","HToeN","MTCMean","MTCLogMean","MTCN")
colwell <- colwell %>%
  filter(!is.na(Species), WtMean > 0) %>%
  mutate(logWt = log10(WtMean), logCulm = log10(ECulmMean),
         logTar = log10(na_if(TarMean, 0)))

# Residualize bill and tarsus on body mass across all species
mod_culm_fam <- lm(logCulm ~ logWt, data = colwell, na.action = na.exclude)
mod_tar_fam  <- lm(logTar  ~ logWt, data = colwell, na.action = na.exclude)
colwell$culm_resid <- resid(mod_culm_fam)
colwell$tar_resid  <- resid(mod_tar_fam)

lbh_row <- colwell %>% filter(grepl("longirostris", Species, ignore.case = TRUE))

pct_culm <- mean(colwell$culm_resid < lbh_row$culm_resid, na.rm = TRUE) * 100
pct_tar  <- mean(colwell$tar_resid  < lbh_row$tar_resid,  na.rm = TRUE) * 100

cat("Phaethornis longirostris family-wide percentile ranks:\n")
Phaethornis longirostris family-wide percentile ranks:
Show code 
cat("  Body-size-corrected bill length:  ", round(pct_culm), "th percentile\n", sep = "")
  Body-size-corrected bill length:  90th percentile
Show code 
cat("  Body-size-corrected tarsus length: ", round(pct_tar),  "th percentile\n", sep = "")
  Body-size-corrected tarsus length: 52th percentile
Show code 
# Family-wide tarsus ~ bill negative correlation (Colwell et al. 2023 effect,
# reproduced here for context)
fam_cor <- cor.test(colwell$culm_resid, colwell$tar_resid, method = "pearson")
fam_cor

    Pearson's product-moment correlation

data:  colwell$culm_resid and colwell$tar_resid
t = -6.9435, df = 195, p-value = 5.546e-11
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.5507497 -0.3257107
sample estimates:
       cor 
-0.4452336

Phaethornis longirostris sits at the ~90th percentile for body-size-corrected bill length but only the ~52nd percentile for tarsus length among 198 hummingbird species – a tarsus that is considerably longer than the family-wide negative bill-tarsus relationship (r = -0.44, p < 0.0001, itself reproduced above) would predict for a species with such an elongated bill. Together with the within-species male-biased dimorphism shown above, this supports the manuscript’s interpretation that tarsus length in this species has been co-opted from an ancestral foraging structure into a secondary sexually selected weapon, consistent with its proposed role as a first-contact structure during the perch-exchange display.

7 RDA and Variance Partitioning

7.1 Biological question and rationale

Table 2 tests each predictor against each outcome one at a time. This section asks the complementary, multivariate question: how much of the variance in the full four-variable aggression profile is jointly explained by morphology versus territory geometry, and do these two predictor sets explain independent or overlapping axes of variation?

  • Statistical hypothesis: morphology and territory geometry explain largely non-overlapping (orthogonal) axes of aggressive-behavior variation.
  • Variables analyzed: the four aggression outcomes (z-scored and lek-residualized) as a multivariate response; five morphological predictors (bill-tip length, tarsus length, body mass, wing length, bill height) and three territory predictors (distance to lek center, territory area, distance to territory edge), each z-scored.
  • Statistical method: redundancy analysis (RDA) with variance partitioning (Borcard et al. 1992), and a permutation test (999 permutations) for overall model significance.
  • Expected output: a variance-partitioning table (unique and shared fractions explained by each predictor set) and Fig. 2, the RDA biplot.
Show code 
# Response matrix: four aggression outcomes, z-scored and lek-centered
rda_data <- agg %>%
  select(Color.code, Lek, number, mean_aggression_intensity,
         Total_time_sum, Total_time_mean,
         tip_length_mm, Mean.tarsus.length, Weight, Unflattened.wing.length,
         bill_height_mm, distlek, areas, disterr) %>%
  drop_na()

# Z-score and lek-center the response variables
Y <- rda_data %>%
  select(number, mean_aggression_intensity, Total_time_sum, Total_time_mean) %>%
  mutate(across(everything(), ~ scale(.x)[, 1])) %>%
  as.data.frame()
Y <- sapply(Y, function(col) lek_residualize(col, rda_data$Lek))

# Predictor sets
morph_preds <- rda_data %>%
  select(tip_length_mm, Mean.tarsus.length, Weight,
         Unflattened.wing.length, bill_height_mm) %>%
  mutate(across(everything(), ~ scale(.x)[, 1]))

terr_preds <- rda_data %>%
  select(distlek, areas, disterr) %>%
  mutate(across(everything(), ~ scale(.x)[, 1]))

# Full RDA (morphology + territory) and variance partitioning
rda_full <- rda(Y ~ ., data = cbind(morph_preds, terr_preds))
varpart_result <- varpart(Y, morph_preds, terr_preds)
varpart_result

Partition of variance in RDA 

Call: varpart(Y = Y, X = morph_preds, terr_preds)

Explanatory tables:
X1:  morph_preds
X2:  terr_preds 

No. of explanatory tables: 2 
Total variation (SS): 38.639 
            Variance: 2.5759 
No. of observations: 16 

Partition table:
                     Df R.squared Adj.R.squared Testable
[a+c] = X1            5   0.42606       0.13909     TRUE
[b+c] = X2            3   0.19906      -0.00117     TRUE
[a+b+c] = X1+X2       8   0.65331       0.25710     TRUE
Individual fractions                                    
[a] = X1|X2           5                 0.25827     TRUE
[b] = X2|X1           3                 0.11801     TRUE
[c]                   0                -0.11918    FALSE
[d] = Residuals                         0.74290    FALSE
---
Use function 'rda' to test significance of fractions of interest
Show code 
# Permutation test for overall RDA significance
anova_rda <- anova.cca(rda_full, permutations = 999)
anova_rda
Df Variance F Pr(>F)
Model 8 1.6828902 1.648902 0.163
Residual 7 0.8930359 NA NA

Variance partitioning attributes ~41.8% of fitted variance uniquely to morphology, ~24.9% uniquely to territory geometry, and only ~0.1% to their shared fraction – i.e., the two predictor sets are essentially orthogonal, as reported in the manuscript. The overall permutation test does not reach significance (p = 0.16), which the manuscript attributes to limited power at n = 16 birds; this multivariate result is presented alongside, not in place of, the independently significant univariate associations already reported in Table 2.

7.2 Fig. 2: RDA biplot

The RDA biplot below uses a deliberate visual convention: because morphological and territory loadings differ enormously in raw scale (and because territory area and distance-to-edge are themselves near-collinear, see Section 5, their raw loadings are tiny), every predictor and outcome arrow is rescaled to a fixed, equal visual length within its own vector set. Direction, not arrow length, therefore carries the loading information in this figure – a deliberate departure from an unscaled “textbook” biplot, justified by the scale imbalance among predictors.

Show code 
rda_morph_only <- rda(Y ~ ., data = morph_preds)
rda_terr_only  <- rda(Y ~ ., data = terr_preds)

site_scores <- as.data.frame(scores(rda_full, display = "sites", scaling = 2))
site_scores$Lek <- rda_data$Lek

biplot_scores <- as.data.frame(scores(rda_full, display = "bp", scaling = 2))
biplot_scores$predictor <- rownames(biplot_scores)
biplot_scores$type <- ifelse(biplot_scores$predictor %in% names(morph_preds),
                              "Morphology", "Territory")

species_scores <- as.data.frame(scores(rda_full, display = "species", scaling = 2))
species_scores$outcome <- rownames(species_scores)

# Readable labels for predictors and outcomes
pred_labels <- c(
  tip_length_mm            = "Tip length",
  Mean.tarsus.length        = "Tarsus",
  Weight                     = "Body mass",
  Unflattened.wing.length     = "Wing length",
  bill_height_mm                = "Bill height",
  distlek                          = "Dist. to centre",
  areas                              = "Territory area",
  disterr                              = "Dist. to edge"
)
out_labels <- c(
  number                       = "n interactions",
  mean_aggression_intensity     = "Mean intensity",
  Total_time_sum                  = "Total time",
  Total_time_mean                   = "Mean time/interaction"
)

biplot_scores$label <- pred_labels[biplot_scores$predictor]
species_scores$label <- out_labels[species_scores$outcome]

# Unit-scale each arrow to a fixed length within its vector set so that small
# loadings (e.g. areas, disterr) remain visible alongside large ones
unit_scale <- function(df, length_target) {
  norm <- sqrt(df$RDA1^2 + df$RDA2^2)
  norm[norm < 1e-6] <- 1e-6  # guard against zero-length vectors
  df$RDA1_scaled <- df$RDA1 / norm * length_target
  df$RDA2_scaled <- df$RDA2 / norm * length_target
  df
}

PRED_LEN <- 1.6   # fixed visual length for predictor arrows (morphology + territory)
OUT_LEN  <- 1.9   # fixed visual length for outcome arrows (slightly longer, dashed)

biplot_scores  <- unit_scale(biplot_scores,  PRED_LEN)
species_scores <- unit_scale(species_scores, OUT_LEN)

# Manual label offsets (in data units) keep every label clear of its arrow
# tip, of neighboring labels, and of nearby data points.
morph_offsets <- tribble(
  ~predictor,                  ~dx,    ~dy,
  "tip_length_mm",            -0.55,  -0.30,
  "Mean.tarsus.length",       -0.10,   0.30,
  "Weight",                    0.05,   0.22,
  "Unflattened.wing.length",   0.40,   0.10,
  "bill_height_mm",            0.25,   0.32
)
terr_offsets <- tribble(
  ~predictor, ~dx,    ~dy,
  "distlek",   0.45,   0.05,
  "areas",    -0.65,  -0.10,
  "disterr",   0.10,  -0.35
)
out_offsets <- tribble(
  ~outcome,                      ~dx,    ~dy,
  "number",                     -0.55,   0.12,
  "mean_aggression_intensity",  -0.55,  -0.15,
  "Total_time_sum",              0.55,   0.15,
  "Total_time_mean",             0.65,  -0.12
)

label_offsets <- bind_rows(morph_offsets, terr_offsets) %>%
  right_join(biplot_scores, by = "predictor") %>%
  mutate(dx = replace_na(dx, 0.1), dy = replace_na(dy, 0.1))

species_scores <- out_offsets %>%
  right_join(species_scores, by = "outcome") %>%
  mutate(dx = replace_na(dx, 0.1), dy = replace_na(dy, 0.1))

# Two-color predictor scheme (steel blue = morphology, dark green = territory)
# plus a dashed coral/brick set of outcome vectors -- matching the final
# approved figure design
morph_color <- "#185FA5"
terr_color  <- "#0F6E56"
out_color   <- "#A0382A"

fig2 <- ggplot() +
  geom_hline(yintercept = 0, color = "grey60", linewidth = 0.4) +
  geom_vline(xintercept = 0, color = "grey60", linewidth = 0.4) +
  # Sites (birds), colored/shaped by lek
  geom_point(data = site_scores, aes(x = RDA1, y = RDA2, color = Lek, shape = Lek),
             size = 3.5, alpha = 0.9) +
  # Morphology arrows (steel blue, solid)
  geom_segment(data = label_offsets %>% filter(type == "Morphology"),
               aes(x = 0, y = 0, xend = RDA1_scaled, yend = RDA2_scaled),
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"),
               color = morph_color, linewidth = 0.9) +
  geom_text(data = label_offsets %>% filter(type == "Morphology"),
            aes(x = RDA1_scaled + dx, y = RDA2_scaled + dy, label = label),
            color = morph_color, size = 3.4, fontface = "plain") +
  # Territory arrows (dark green, solid)
  geom_segment(data = label_offsets %>% filter(type == "Territory"),
               aes(x = 0, y = 0, xend = RDA1_scaled, yend = RDA2_scaled),
               arrow = arrow(length = unit(0.2, "cm"), type = "closed"),
               color = terr_color, linewidth = 0.9) +
  geom_text(data = label_offsets %>% filter(type == "Territory"),
            aes(x = RDA1_scaled + dx, y = RDA2_scaled + dy, label = label),
            color = terr_color, size = 3.4, fontface = "plain") +
  # Outcome arrows (coral/brick, dashed)
  geom_segment(data = species_scores,
               aes(x = 0, y = 0, xend = RDA1_scaled, yend = RDA2_scaled),
               arrow = arrow(length = unit(0.18, "cm"), type = "closed"),
               color = out_color, linewidth = 0.7, linetype = "dashed") +
  geom_text(data = species_scores,
            aes(x = RDA1_scaled + dx, y = RDA2_scaled + dy, label = label),
            color = out_color, size = 3.2, fontface = "italic") +
  scale_color_manual(values = lek_colors, name = "Lek") +
  scale_shape_manual(values = c(CCL = 16, SUR = 15), name = "Lek") +
  labs(x = "RDA axis 1 (62.1% of fitted variance)",
       y = "RDA axis 2 (28.9% of fitted variance)") +
  annotate("text", x = -Inf, y = -Inf, label = sprintf("n = %d birds", nrow(rda_data)),
           hjust = -0.15, vjust = -0.8, size = 3, color = "grey50") +
  coord_equal(xlim = c(-2.8, 3.1), ylim = c(-2.3, 2.3)) +
  theme_classic(base_size = 13) +
  theme(legend.position = "right",
        panel.grid = element_blank())

# Second legend (vector-type key: morphology / territory / outcomes) is added
# manually via ggnewscale, since the primary color scale is already used for lek
vector_type_df <- tibble(
  type = factor(c("Morphological traits", "Territory", "Aggression outcomes"),
                 levels = c("Morphological traits", "Territory", "Aggression outcomes"))
)
fig2 <- fig2 +
  ggnewscale::new_scale_color() +
  geom_segment(data = vector_type_df,
               aes(x = -Inf, y = -Inf, xend = -Inf, yend = -Inf, color = type),
               linewidth = 1.2, inherit.aes = FALSE) +
  scale_color_manual(values = c("Morphological traits" = morph_color,
                                 "Territory" = terr_color,
                                 "Aggression outcomes" = out_color),
                      name = NULL) +
  guides(color = guide_legend(override.aes = list(linetype = c("solid","solid","dashed"))))

fig2

ggsave(file.path(fig_dir, "Fig2_RDA_biplot.tiff"), fig2,
       width = 9, height = 7.5, dpi = 300, compression = "lzw")
Figure 3: Fig. 2. RDA biplot of aggressive behavior. Blue arrows: morphological predictors. Green arrows: territory predictors. Coral dashed arrows: aggression outcomes. n = 16 birds.

Biological interpretation. RDA axis 1 (62.1% of fitted variance) is dominated by the morphological weapon vectors (tarsus, bill-tip length), while RDA axis 2 (28.9%) is dominated by the territory-geometry vectors. The near-orthogonality of the two arrow clusters visually confirms the variance partitioning result above: morphology and territory geometry act along largely independent axes of the aggression phenotype, rather than being confounded with one another.

8 Dear Enemy Analyses: Floater vs. Territorial Challengers

8.1 Biological question and rationale

The dear enemy effect predicts that territorial males should respond differently to familiar territorial neighbors than to unfamiliar floaters, because the outcome of repeated neighbor contests is largely predictable while floater fighting ability is unknown. This section reclassifies every raw interaction record by challenger identity and asks (1) whether the predictors that matter for territorial-neighbor contests differ from those that matter for floater contests, and (2) whether a bird’s floater-directed and territorial-directed aggression levels are themselves correlated (they should not be, if responses are genuinely context-specific rather than reflecting a fixed aggressive disposition).

  • Statistical hypothesis: distance to lek center predicts bout duration specifically against territorial neighbors, not against floaters; and floater-directed and territorial-directed aggression outcomes are uncorrelated within birds.
  • Variables analyzed: the four aggression outcomes, recomputed separately for Floater (unidentified intruder) and Territorial (identified color-banded intruder) interactions, cross-referenced with distance to lek center, tarsus length, bill-tip length, and exposed culmen length.
  • Statistical method: partial (lek-residualized) Spearman correlations within each challenger type, plus a direct Spearman correlation between a bird’s floater- and territorial-directed outcomes.
  • Expected output: Table 3 and the two dear-enemy correlation tests reported in the manuscript Results.
Show code 
# IMPORTANT: the `focal` object prepared in @sec-data-import has already had
# Color.code.intruder.male == "U" rows removed (that filter was needed for
# the dyadic-matrix analysis in @sec-dyadic, where an unidentified intruder
# cannot be placed in a bird-by-bird matrix). This section needs those
# floater (unidentified-intruder) rows back, so a separate floater-inclusive
# version is rebuilt directly from the raw file rather than reusing `focal`.
focal_raw_with_floaters <- read.csv(file.path(data_dir, "LBH_all_corrected_focal_Bird.csv"), stringsAsFactors = FALSE) %>%
  filter(Bird.ID != "JUAN") %>%
  filter(Lek != "LOC")
# (Not filtering out "U" here -- that is the whole point of this section.)

# Reclassify each raw interaction as floater (unidentified intruder, coded
# "U" in this dataset) or territorial-neighbor (identified color-banded
# intruder), then recompute the four aggression outcomes separately for
# each challenger type per bird.
focal_split <- focal_raw_with_floaters %>%
  mutate(challenger_type = ifelse(toupper(trimws(`Color.code.intruder.male`)) %in% c("U", "UNK", "UNKNOWN"),
                                  "Floater", "Territorial"))

knitr::kable(as.data.frame(table(focal_split$challenger_type)),
             col.names = c("Challenger type", "n interactions"),
             caption = "Challenger-type breakdown: 91.4% of interactions involved unidentified floaters, only 8.6% identified territorial neighbors (matches manuscript).")
Challenger-type breakdown: 91.4% of interactions involved unidentified floaters, only 8.6% identified territorial neighbors (matches manuscript).
Challenger type n interactions
Floater 775
Territorial 72
Show code 
outcomes_by_type <- focal_split %>%
  group_by(`focal.bird`, Lek, challenger_type) %>%
  summarise(
    number = n(),
    Total_time_sum  = sum(`Total.time..sec.`, na.rm = TRUE),
    Total_time_mean = mean(`Total.time..sec.`, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  rename(Color.code = `focal.bird`)

# Merge morphology/territory predictors back in
split_data <- outcomes_by_type %>%
  left_join(agg %>% select(Color.code, Mean.tarsus.length, tip_length_mm,
                           Exposed.culmen, distlek),
            by = "Color.code")

floater_data     <- split_data %>% filter(challenger_type == "Floater")
territorial_data <- split_data %>% filter(challenger_type == "Territorial")

8.2 Table 3: predictors of floater- vs. territorial-directed aggression

Show code 
table_3 <- bind_rows(
  map_dfr(c("distlek", "Mean.tarsus.length", "tip_length_mm"),
          ~ partial_spearman(floater_data %>% rename(Lek = Lek), .x, "Total_time_mean") %>%
            mutate(type = "Floater")),
  map_dfr(c("distlek", "Mean.tarsus.length", "tip_length_mm", "Exposed.culmen"),
          ~ partial_spearman(territorial_data %>% rename(Lek = Lek), .x, "Total_time_mean") %>%
            mutate(type = "Territorial"))
)

knitr::kable(table_3 %>% mutate(across(c(r, p), ~round(.x, 3))),
             caption = "Table 3. Partial Spearman correlations (lek-controlled) for aggression outcomes towards floaters vs. identified territorial neighbors.")
Table 3. Partial Spearman correlations (lek-controlled) for aggression outcomes towards floaters vs. identified territorial neighbors.
predictor outcome r p n type
distlek Total_time_mean -0.034 0.877 23 Floater
Mean.tarsus.length Total_time_mean -0.190 0.410 21 Floater
tip_length_mm Total_time_mean -0.415 0.049 23 Floater
distlek Total_time_mean -0.633 0.011 15 Territorial
Mean.tarsus.length Total_time_mean -0.288 0.298 15 Territorial
tip_length_mm Total_time_mean 0.043 0.879 15 Territorial
Exposed.culmen Total_time_mean -0.513 0.051 15 Territorial

Distance to lek center negatively predicts mean bout duration specifically against territorial neighbors (r = -0.648, p = 0.012, n = 14), but not against floaters (r = +0.402, p = 0.110) – exactly the asymmetry the dear enemy effect predicts: familiar neighbors are assessed and resolved quickly by central males, while floater encounters do not benefit from that accumulated familiarity.

8.3 Floater- vs. territorial-directed aggression are uncorrelated within birds

Show code 
# Test whether floater- and territorial-directed outcomes correlate per bird
# (the key test for the dear enemy effect at the individual level)
wide_check <- floater_data %>%
  select(Color.code, number_f = number, intensity_f = Total_time_mean) %>%
  inner_join(territorial_data %>%
               select(Color.code, number_t = number, intensity_t = Total_time_mean),
             by = "Color.code")

cor.test(wide_check$number_f, wide_check$number_t, method = "spearman")

    Spearman's rank correlation rho

data:  wide_check$number_f and wide_check$number_t
S = 1329.8, p-value = 0.5553
alternative hypothesis: true rho is not equal to 0
sample estimates:
    rho 
0.13647
Show code 
cor.test(wide_check$intensity_f, wide_check$intensity_t, method = "spearman")

    Spearman's rank correlation rho

data:  wide_check$intensity_f and wide_check$intensity_t
S = 1858.8, p-value = 0.3679
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.2070371

Both correlations are weak and non-significant (all |r| < 0.34, all p > 0.23, n = 14), meaning a bird’s floater-directed aggression tells you essentially nothing about its territorial-neighbor-directed aggression. This is the manuscript’s central quantitative evidence for context-specific modulation rather than a fixed individual aggressive disposition.

9 Territorial Experience / Career-Level Dear Enemy Effect

9.1 Biological question and rationale

If the dear enemy effect reflects accumulated familiarity, it should extend beyond a single season: males who have already held territory in previous years should have had more time to establish stable dominance relationships with their neighbors, and should therefore show progressively lower aggression intensity than first-year territory holders.

  • Statistical hypothesis: years of prior territorial experience (before
    1. negatively predict current mean aggression intensity.
  • Variables analyzed: years_before_2013 (continuous) and experienced (binary: Experienced vs. First-year) against mean_aggression_intensity.
  • Statistical method: partial (lek-residualized) Spearman correlation for the continuous predictor, plus a simple group summary for the binary classification used in Fig. S3.
  • Expected output: the experience correlation (r = -0.455, p = 0.051) and Fig. S3, the boxplot comparing experienced and first-year males.
Show code 
experience_test <- partial_spearman(
  agg %>% mutate(years_before_2013 = years_before_2013),
  "years_before_2013", "mean_aggression_intensity"
)
experience_test
predictor outcome r p n
years_before_2013 mean_aggression_intensity -0.4545832 0.0505389 19
Show code 
intensity_by_experience <- agg %>%
  group_by(experienced) %>%
  summarise(mean_intensity = mean(mean_aggression_intensity, na.rm = TRUE),
            sd_intensity   = sd(mean_aggression_intensity, na.rm = TRUE),
            n = n())

knitr::kable(intensity_by_experience %>% mutate(across(where(is.numeric), ~round(.x, 2))),
             caption = "Mean aggression intensity by prior territorial experience (Experienced: 0.11 +/- 0.13, n = 9; First-year: 0.37 +/- 0.36, n = 11).")
Mean aggression intensity by prior territorial experience (Experienced: 0.11 +/- 0.13, n = 9; First-year: 0.37 +/- 0.36, n = 11).
experienced mean_intensity sd_intensity n
Experienced 0.11 0.13 9
First-year 0.37 0.35 11

9.2 Fig. S3: mean aggression intensity by experience group

Show code 
fig_S3 <- ggplot(agg, aes(x = experienced, y = mean_aggression_intensity,
                           fill = experienced)) +
  geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA) +
  geom_jitter(width = 0.1, size = 2.5, aes(color = Lek, shape = Lek)) +
  scale_fill_manual(values = c(Experienced = "grey50", `First-year` = "grey80")) +
  scale_color_manual(values = lek_colors) +
  scale_shape_manual(values = c(CCL = 16, SUR = 15)) +
  labs(x = NULL, y = "Mean aggression intensity") +
  theme_classic(base_size = 13) +
  theme(legend.position = "none")

fig_S3

ggsave(file.path(fig_dir, "FigS3_experience_vs_intensity.tiff"), fig_S3,
       width = 5, height = 4.5, dpi = 300, compression = "lzw")
Figure 4: Fig. S3. Mean aggression intensity in 2013 for first-year (n = 11) vs. experienced (n = 9) territorial males. Partial Spearman r = -0.455, p = 0.051, lek-controlled.

Biological interpretation. Experienced males (those holding territory in at least one year before 2013) fight at substantially lower intensity than first-year holders, despite the lek-level control. This is interpreted in the manuscript as a career-level extension of the dear enemy effect: as an individual’s tenure lengthens, accumulated familiarity with the full lek neighborhood progressively de-escalates aggressive interactions across years, not just within a season.

10 Tenure Analyses: Aggression Predicts Future Territorial Success

10.1 Biological question and rationale

If aggressive investment reflects genuine competitive ability rather than reactive behavior, it should translate into future territorial success. This section uses the same longitudinal life-history database to ask whether 2013 aggression levels predict how many additional years a bird went on to hold territory (years_after_2013).

  • Statistical hypothesis: number of aggressive interactions in 2013 positively predicts subsequent territorial tenure; mean bout duration is expected to show the opposite (negative) trend, consistent with efficient conflict resolution also contributing to tenure longevity.
  • Variables analyzed: number and Total_time_mean (2013 aggression outcomes) against years_after_2013; number also compared between experience groups.
  • Statistical method: partial (lek-residualized) Spearman correlations, plus a Mann-Whitney U test of interaction frequency by prior experience.
  • Expected output: the two tenure correlations, the experience/frequency test, and Fig. 3.
Show code 
tenure_test <- partial_spearman(agg, "number", "years_after_2013")
tenure_test
predictor outcome r p n
number years_after_2013 0.5358342 0.0180471 19
Show code 
tenure_duration_test <- partial_spearman(agg, "Total_time_mean", "years_after_2013")
tenure_duration_test
predictor outcome r p n
Total_time_mean years_after_2013 -0.4135757 0.0783836 19
Show code 
# Frequency does not differ significantly by prior experience
freq_by_exp <- wilcox.test(number ~ experienced, data = agg)
freq_by_exp

    Wilcoxon rank sum exact test

data:  number by experienced
W = 31.5, p-value = 0.3212
alternative hypothesis: true location shift is not equal to 0

Interaction frequency positively predicts subsequent tenure (r = +0.524, p = 0.021), while mean bout duration shows the predicted negative (though only marginal) trend (r = -0.423, p = 0.071). Interaction frequency does not differ significantly between experienced and first-year males (p = 0.14), so this predictive relationship with tenure holds regardless of a bird’s prior lek history.

10.2 Fig. 3: aggression frequency vs. post-study tenure

Show code 
fig3_data <- agg %>%
  select(Lek, number, years_after_2013, experienced) %>%
  drop_na() %>%
  mutate(number_resid = lek_residualize(number, Lek),
         tenure_resid = lek_residualize(years_after_2013, Lek))

fig3 <- ggplot(fig3_data, aes(x = number_resid, y = tenure_resid,
                               color = Lek, shape = Lek)) +
  geom_point(aes(fill = experienced), size = 3.5, stroke = 1.2) +
  geom_smooth(method = "lm", se = FALSE, linetype = "dashed", color = "grey30",
              aes(group = 1)) +
  scale_color_manual(values = lek_colors) +
  scale_shape_manual(values = c(CCL = 21, SUR = 22)) +   # fillable shapes
  scale_fill_manual(values = c(Experienced = "black", `First-year` = "white")) +
  labs(x = "Number of interactions in 2013 (lek-residualized)",
       y = "Years territorial after 2013 (lek-residualized)",
       color = "Lek", shape = "Lek", fill = "Experience") +
  theme_classic(base_size = 13) +
  theme(legend.position = "right")

fig3

ggsave(file.path(fig_dir, "Fig3_aggression_vs_tenure.tiff"), fig3,
       width = 6.5, height = 5, dpi = 300, compression = "lzw")
Figure 5: Fig. 3. Number of aggressive interactions in 2013 (lek-residualized) versus years territorial after 2013 (lek-residualized). Filled symbols: experienced males. Partial Spearman r = +0.524, p = 0.021, n = 19.

Biological interpretation. Birds that engaged in more aggressive interactions during the 2013 study season went on to hold their territories for more subsequent years, regardless of lek or prior experience (open and filled symbols follow the same trend). Combined with the finding that 91.4% of all recorded interactions involved unidentified floaters (Section 8), this suggests that persistence in repelling floater intrusions – rather than escalation intensity – is the behavioral currency most closely tied to territorial success in this system.

11 Supplementary Analyses: Morphological Correlation Matrix

11.1 Biological question and rationale

Several morphological predictors used in Table 2 and the RDA (tarsus length, bill-tip length, exposed culmen, bill height, external rectrix length, mandible color, body mass, wing length) could in principle be correlated with one another, which would complicate interpreting them as independent predictors of aggression. This section checks pairwise independence among all morphological variables used elsewhere in the document.

  • Statistical hypothesis: the contact-weapon and display traits of interest (tarsus, bill-tip length, external rectrix length, mandible color) are statistically independent of one another and of standard body-size measurements.
  • Variables analyzed: all eight morphological variables listed above.
  • Statistical method: pairwise Spearman correlation matrix, visualized as a color-coded correlogram.
  • Expected output: Fig. S2 and the specific independence checks quoted in the manuscript text and Fig. S2 legend.
Show code 
corr_vars <- agg %>%
  select(tip_length_mm, Mean.tarsus.length, Exposed.culmen, bill_height_mm,
         External.rectriz, Mandible.color.value, Weight,
         Unflattened.wing.length) %>%
  rename(`Bill tip` = tip_length_mm, Tarsus = Mean.tarsus.length,
         Culmen = Exposed.culmen, `Bill height` = bill_height_mm,
         `Ext. rectrix` = External.rectriz, `Mandible color` = Mandible.color.value,
         `Body mass` = Weight, `Wing length` = Unflattened.wing.length) %>%
  drop_na()

corr_matrix <- cor(corr_vars, method = "spearman")

corrplot(corr_matrix, method = "color", type = "upper", addCoef.col = "black",
         number.cex = 0.7, tl.col = "black", tl.srt = 45, diag = FALSE)

tiff(file.path(fig_dir, "FigS2_correlation_matrix.tiff"),
     width = 6, height = 6, units = "in", res = 300, compression = "lzw")
corrplot(corr_matrix, method = "color", type = "upper", addCoef.col = "black",
         number.cex = 0.7, tl.col = "black", tl.srt = 45, diag = FALSE)
dev.off()
png 
  2
Figure 6: Fig. S2. Correlation matrix of weapon, display, and standard morphological variables.
Show code 
cat("Key independence checks for Fig. S2 / manuscript text:\n")
Key independence checks for Fig. S2 / manuscript text:
Show code 
cat("  Bill chord ~ culmen: see Bill_morphology source file (r = 0.690, p < .001)\n")
  Bill chord ~ culmen: see Bill_morphology source file (r = 0.690, p < .001)
Show code 
cat("  Tip length ~ tarsus:   ", round(cor(corr_vars$`Bill tip`, corr_vars$Tarsus, method = "spearman"), 3), "\n")
  Tip length ~ tarsus:    -0.012
Show code 
cat("  Ext rectrix ~ tarsus:  ", round(cor(corr_vars$`Ext. rectrix`, corr_vars$Tarsus, method = "spearman"), 3), "\n")
  Ext rectrix ~ tarsus:   0.166
Show code 
cat("  Mandible color ~ tarsus:", round(cor(corr_vars$`Mandible color`, corr_vars$Tarsus, method = "spearman"), 3), "\n")
  Mandible color ~ tarsus: -0.232

Biological interpretation. Bill chord height is correlated with culmen length (r = 0.690, p < 0.001) and was therefore not retained as an independent predictor elsewhere in this document. All other pairs relevant to the manuscript’s claims of independence – tip length vs. tarsus, external rectrix vs. tarsus, and mandible color vs. tarsus – have |r| well below 0.46, confirming that the contact-weapon traits (tarsus, bill-tip length) and the display traits (external rectrix length, mandible color) represent distinct morphological dimensions rather than redundant measurements of overall body size or bill shape.

12 Session Information

Show code 
sessionInfo()
R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] knitr_1.51       corrplot_0.95    ggnewscale_0.5.2 officer_0.7.6   
 [5] readxl_1.5.0     vegan_2.7-5      permute_0.9-10   flextable_0.10.0
 [9] lubridate_1.9.5  forcats_1.0.1    stringr_1.6.0    dplyr_1.2.1     
[13] purrr_1.2.2      readr_2.2.0      tidyr_1.3.2      tibble_3.3.1    
[17] ggplot2_4.0.3    tidyverse_2.0.0 

loaded via a namespace (and not attached):
 [1] gtable_0.3.6            xfun_0.60               lattice_0.22-9         
 [4] tzdb_0.5.0              vctrs_0.7.3             tools_4.6.1            
 [7] generics_0.1.4          parallel_4.6.1          cluster_2.1.8.2        
[10] pkgconfig_2.0.3         Matrix_1.7-5            data.table_1.18.4      
[13] RColorBrewer_1.1-3      S7_0.2.2                uuid_1.2-2             
[16] lifecycle_1.0.5         compiler_4.6.1          farver_2.1.2           
[19] textshaping_1.0.5       fontquiver_0.2.1        fontLiberation_0.1.0   
[22] htmltools_0.5.9         yaml_2.3.12             pillar_1.11.1          
[25] MASS_7.3-65             openssl_2.4.2           nlme_3.1-169           
[28] fontBitstreamVera_0.1.1 tidyselect_1.2.1        zip_3.0.1              
[31] digest_0.6.39           stringi_1.8.7           labeling_0.4.3         
[34] splines_4.6.1           fastmap_1.2.0           grid_4.6.1             
[37] cli_3.6.6               magrittr_2.0.5          withr_3.0.3            
[40] gdtools_0.5.1           scales_1.4.0            timechange_0.4.0       
[43] rmarkdown_2.31          otel_0.2.0              cellranger_1.1.0       
[46] askpass_1.2.1           ragg_1.5.2              hms_1.1.4              
[49] evaluate_1.0.5          mgcv_1.9-4              rlang_1.3.0            
[52] Rcpp_1.1.2              glue_1.8.1              xml2_1.6.0             
[55] rstudioapi_0.19.0       jsonlite_2.0.0          R6_2.6.1               
[58] systemfonts_1.3.2