Introduction


Progress of the task: 3% complete (3 of 100)



This document provides a conceptual and methodological walkthrough of papers on parental coordination methodology. For each paper, the relevant methods are extracted, summarized, and—when possible—visually reconstructed to clarify how the analysis was performed. Based on these descriptions, function/s are then developed to process data in the same way as in the original studies.

Other methodological details are extracted from each paper, including: basic bibliographic information and methodological details as specified in the two tables below.


Bibliographic information for each paper included in the methodological review of parental coordination.



Methodological information from each paper included in the methodological review of parental coordination



#001 Mariette & Griffith 2012

Mariette MM, Griffith SC. 2012. Nest visit synchrony is high and correlates with reproductive success in the wild Zebra finch Taeniopygia guttata. Journal of Avian Biology. 43(2):131–140. https://doi.org/10.1111/j.1600-048X.2012.05555.x


Methods excerpt

During nestling rearing, when parents arrived at the box together, the male and female usually took turns to enter the box individually rather than entering at the same time, and left the nest area at the same time (unpubl.). Typically, when one bird was inside the nest, its partner sat outside the nest, but within about a meter (Zann 1996, unpubl.) and during this time partners maintain contact vocally (Elie et al. 2010). Males and females each visited the nest on average 12.57 times a day during nestling rearing (see Results; hereafter ‘nest visit rate’; by contrast we refer to the combined nest visit rate by both parents (sum of visits by male and female partners) as ‘pair nest visit rate’), for 4.76 (SD 3.06) min per visit, as measured by the PIT-tag detection system. To account for the transition time between the exit of one parent and the entry of its partner, we considered that the two partners were present together at the nest (i.e. synchronous visit) if they individually triggered the decoder within two minutes of each other. Using a period of less than two minutes would have meant that on some occasions when the birds actually arrived at the box together the second bird would not have had a chance to enter the box after the first bird had been in (see Mariette et al. 2011 for further details). Increasing this threshold by 5 min only increased the proportion of visits where both parents were together by 5% and did not qualitatively change any of the results (including the test of observed versus expected nest visit synchrony), so we only present data for the 2-min transition time. Pairs were active for 13 h (i.e. 780 min) d 1 throughout the season (SeptemberDecember) so that we estimated that each parent was present at the nest for (2 4.76) 12.57/780 100 11% of the day (i.e. for one adult: (2 min transition time mean visit duration) daily visit rate/active day length 100). We conservatively added 2 min to the average visit duration in the above calculation to account for the time the bird may be waiting in the nesting tree before or after it entered the nest-box (i.e. ‘transition time’). If the male and female visited the nest independently of each other, and randomly, we would expect them to be recorded together at the nest 11 11% 1.2% of the time. Based on the above logic, we calculated the expected time together at the nest for each pair separately, using the individual’s specific nest visit rate and nest visit duration, to compare it with the observed proportion of time that partners spent at the nest together per day using a paired t test (where values were paired per pair).

We investigated whether parents actively coordinated their visits to the nest (as opposed to nest synchrony merely resulting from the fact that both parents were independently following similar schedule) by testing whether synchrony was maintained after one parent visited the nest alone. If synchrony were passive, a visit by only one of the parents would indicate that partners are out of phase and synchrony is disrupted. Partners should then be less synchronised after a visit alone than before, unless the parent that skipped a visit readjusted its next visit time on that of its partner. We therefore compared (using a paired t test) the interval between partners’ visits for the visit preceding and that following a visit where only one parent was present and expected the latter to be longer that the former if nest visit synchrony were passive. To do so, we measured the (square-root transformed) interval between the time the first parent last exited the box and the time its partner first entered the box. When parental visits overlapped (i.e. the first partner re-entered the box after the second partner entered), the interval between partners was set to zero for that visit.


Breakdown the logic

Observed synchrony was calculated from nest‑entry records, assuming that two visits were synchronous when they overlapped or, if they did not overlap, when their entry times were separated by no more than 2 minutes.


Then to calculated the expected values, there were two steps:

Step 1 — computation of % of day each parent was at the nest, using this formula

\(p_i = \frac{(\text{visit duration} + \text{transition time}) \cdot \text{visit rate}}{\text{active day length}}\)

when:

visit duration = 4.76 min

transition time = 2 min

visit rate = 12.57 visits/day

active day = 780 min


Then:

\(p = \frac{(4.76 + 2)\cdot 12.57}{780} \approx 0.11 = 11\%\)


Step 2 — calculation of expected synchrony if visits are independent

If male and female visit independently:

\(\text{Expected synchrony} = p_{\text{male}} \cdot p_{\text{female}}\)

With both ≈ 0.11:

0\(0.11 \times 0.11 = 0.0121 = 1.2\%\)

Then observed proportion of time that partners spent at the nest together per day, compared with the expected proportion if they were independent, using a paired t test (where values were paired per pair)

Visualization of the passive synchronization analysis - Interval after solo expected to be longer than before


# passive synchronization

df_passive <- tibble(
  event = c(
    "A_exit_before", 
    "B_entry_before",
    "A_entry_solo", 
    "A_exit_after", 
    "B_entry_after"
  ),
  time = c(1, 2, 4, 6, 10),
  parent = c("A", "B", "A", "A", "B")
)

intervals_passive <- tibble(
  interval = c("before", "after"),
  xstart = c(1, 6),
  xend   = c(2, 10),
  y      = c(1.3, 1.3)
)

rects_passive <- tibble(
  label = c("synchrony", "solo"),
  xmin  = c(0.8, 3.5),
  xmax  = c(2.2, 4.5),
  ymin  = c(0.85, 0.85),
  ymax  = c(1.15, 1.15)
)

ggplot() +
  geom_rect(data = rects_passive,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = label),
            alpha = 0.25) +
  geom_segment(aes(x = 0.5, xend = 11, y = 1, yend = 1),
               linewidth = 0.4, colour = "grey60") +
  geom_point(data = df_passive,
             aes(x = time, y = 1, colour = parent),
             size = 4) +
  geom_text(data = df_passive,
            aes(x = time, y = 1.18, label = event),
            size = 4, angle = 90, vjust = 0) +
  geom_segment(data = intervals_passive %>% filter(interval == "before"),
               aes(x = xstart, xend = xend, y = y, yend = y),
               colour = "forestgreen", linewidth = 2,
               arrow = arrow(length = unit(0.25, "cm"))) +
  geom_segment(data = intervals_passive %>% filter(interval == "after"),
               aes(x = xstart, xend = xend, y = y, yend = y),
               colour = "firebrick", linewidth = 2,
               arrow = arrow(length = unit(0.25, "cm"))) +
  annotate("text", x = 1.5, y = 1.45, label = "Interval BEFORE", colour = "forestgreen") +
  annotate("text", x = 8, y = 1.45, label = "Interval AFTER", colour = "firebrick") +
  scale_colour_manual(values = c("A" = "#1b9e77", "B" = "#d95f02")) +
  scale_fill_manual(values = c("Synchronous visit" = "skyblue", 
                               "Solo visit" = "orange")) +
  theme_minimal(base_size = 14) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        panel.grid = element_blank(),
        legend.title = element_blank(),
        legend.position = "bottom")

Visualization of the active synchronization analysis - Interval after solo expected to be similar to those before:


# active synchroniztion

df_active <- tibble(
  event = c(
    "A_exit_before", 
    "B_entry_before",
    "A_entry_solo", 
    "A_exit_after", 
    "B_entry_after"
  ),
  time = c(1, 2, 4, 6, 7),   # partner returns quickly
  parent = c("A", "B", "A", "A", "B")
)

intervals_active <- tibble(
  interval = c("before", "after"),
  xstart = c(1, 6),
  xend   = c(2, 7),
  y      = c(1.3, 1.3)
)

rects_active <- tibble(
  label = c("synchrony", "solo"),
  xmin  = c(0.8, 3.5),
  xmax  = c(2.2, 4.5),
  ymin  = c(0.85, 0.85),
  ymax  = c(1.15, 1.15)
)

ggplot() +
  geom_rect(data = rects_active,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = label),
            alpha = 0.25) +
  geom_segment(aes(x = 0.5, xend = 11, y = 1, yend = 1),
               linewidth = 0.4, colour = "grey60") +
  geom_point(data = df_active,
             aes(x = time, y = 1, colour = parent),
             size = 4) +
  geom_text(data = df_active,
            aes(x = time, y = 1.18, label = event),
            size = 4, angle = 90, vjust = 0) +
  geom_segment(data = intervals_active %>% filter(interval == "before"),
               aes(x = xstart, xend = xend, y = y, yend = y),
               colour = "forestgreen", linewidth = 2,
               arrow = arrow(length = unit(0.25, "cm"))) +
  geom_segment(data = intervals_active %>% filter(interval == "after"),
               aes(x = xstart, xend = xend, y = y, yend = y),
               colour = "firebrick", linewidth = 2,
               arrow = arrow(length = unit(0.25, "cm"))) +
  annotate("text", x = 1.5, y = 1.45, label = "Interval BEFORE", colour = "forestgreen") +
  annotate("text", x = 6.5, y = 1.45, label = "Interval AFTER", colour = "firebrick") +
  scale_colour_manual(values = c("A" = "#1b9e77", "B" = "#d95f02")) +
  scale_fill_manual(values = c("Synchronous visit" = "skyblue", 
                               "Solo visit" = "orange")) +
  theme_minimal(base_size = 14) +
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        panel.grid = element_blank(),
        legend.title = element_blank(),
        legend.position = "bottom")


Developed functions

synchrony_mariette2()

# Prop of overalping time in respect to total time (Mariette and Griffith 2012)
synchrony_mariette2 <- function(df, transition_sec = 120, active_day_sec = 780 * 60) {
  
  # split male and female
  male   <- df %>% filter(sex == "male")
  female <- df %>% filter(sex == "female")
  
  # count visits
  V_m <- nrow(male)
  V_f <- nrow(female)
  
  # mean visit durations
  mean_dur_m <- if (V_m > 0) mean(as.numeric(male$departure_time - male$arrival_time, units = "secs")) else 0
  mean_dur_f <- if (V_f > 0) mean(as.numeric(female$departure_time - female$arrival_time, units = "secs")) else 0
  
  # presence proportion
  # p = (transition + mean_duration) * visit_rate / active_day_length
  p_m <- if (V_m > 0) (transition_sec + mean_dur_m) * V_m / active_day_sec else 0
  p_f <- if (V_f > 0) (transition_sec + mean_dur_f) * V_f / active_day_sec else 0
  
  # expected synchrony under independence
  expected_sync <- p_m * p_f
  
  # observed synchrony
  total_overlap <- 0
  
  for (i in seq_len(nrow(male))) {
    m_start <- male$arrival_time[i]
    m_end   <- male$departure_time[i]
    
    overlaps <- female %>%
      filter(arrival_time <= m_end & departure_time >= m_start)
    
    if (nrow(overlaps) > 0) {
      for (j in seq_len(nrow(overlaps))) {
        f_start <- overlaps$arrival_time[j]
        f_end   <- overlaps$departure_time[j]
        
        overlap_start <- max(m_start, f_start)
        overlap_end   <- min(m_end, f_end)
        
        if (overlap_end > overlap_start) {
          total_overlap <- total_overlap +
            as.numeric(overlap_end - overlap_start, units = "secs")
        }
      }
    }
  }
  
  observed_sync <- total_overlap / active_day_sec
  
  # Interpretation ----
  message("\n--- Nest Visit Synchrony Interpretation (seconds-based) ---\n")
  message(sprintf("Observed synchrony: %.4f (%.2f%% of the active day)",
                  observed_sync, observed_sync * 100))
  message(sprintf("Expected synchrony under independence: %.4f (%.2f%% of the active day)",
                  expected_sync, expected_sync * 100))
  
  if (observed_sync > expected_sync) {
    message("\nInterpretation: Observed synchrony is HIGHER than expected under independent visits.")
    message("This suggests active coordination between parents.\n")
  } else if (observed_sync < expected_sync) {
    message("\nInterpretation: Observed synchrony is LOWER than expected under independent visits.")
    message("This suggests parents may be out of phase or avoiding overlap.\n")
  } else {
    message("\nInterpretation: Observed synchrony matches expected synchrony.")
    message("This is consistent with independent, uncoordinated visits.\n")
  }
  
  # Return numerical results
  return(list(
    male_visits      = V_m,
    female_visits    = V_f,
    mean_dur_m_sec   = mean_dur_m,
    mean_dur_f_sec   = mean_dur_f,
    p_m              = p_m,
    p_f              = p_f,
    expected_sync    = expected_sync,
    observed_sync    = observed_sync,
    overlap_seconds  = total_overlap
  ))
}

synchrony_disruption_mariette()

# Extra funs from Mariette and Griffith 2012
synchrony_disruption_mariette <- function(df) {
  
  # helper: compute interval between two visits (seconds)
  # interval = time between first parent exiting and partner entering
  # if overlap → interval = 0
  compute_interval <- function(exit_time, entry_time) {
    if (entry_time <= exit_time) return(0)
    as.numeric(entry_time - exit_time, units = "secs")
  }
  
  # identify visits where only one parent was present
  # "solo visit" = no overlap with the other parent
  is_overlapping <- function(v1_start, v1_end, v2_start, v2_end) {
    (v1_start <= v2_end) & (v2_start <= v1_end)
  }
  
  df$solo <- FALSE
  
  for (i in seq_len(nrow(df))) {
    focal <- df[i, ]
    partner <- df %>% filter(sex != focal$sex)
    
    overlaps <- partner %>%
      filter(is_overlapping(arrival_time, departure_time,
                            focal$arrival_time, focal$departure_time))
    
    if (nrow(overlaps) == 0) df$solo[i] <- TRUE
  }
  
  # extract solo visits
  solo_visits <- df %>% filter(solo)
  if (nrow(solo_visits) < 1) {
    message("No solo visits detected — cannot compute synchrony disruption test.")
    return(NULL)
  }
  
  before_intervals <- c()
  after_intervals  <- c()
  
  # for each solo visit, compute BEFORE and AFTER intervals
  for (i in seq_len(nrow(solo_visits))) {
    solo <- solo_visits[i, ]
    idx  <- which(df$arrival_time == solo$arrival_time)
    
    # BEFORE interval: partner entry relative to previous focal exit
    if (idx > 1) {
      prev <- df[idx - 1, ]
      if (prev$sex != solo$sex) {
        interval_before <- compute_interval(prev$departure_time, solo$arrival_time)
        before_intervals <- c(before_intervals, sqrt(interval_before))
      }
    }
    
    # AFTER interval: partner entry relative to focal exit
    if (idx < nrow(df)) {
      next_visit <- df[idx + 1, ]
      if (next_visit$sex != solo$sex) {
        interval_after <- compute_interval(solo$departure_time, next_visit$arrival_time)
        after_intervals <- c(after_intervals, sqrt(interval_after))
      }
    }
  }
  
  # ensure equal lengths for paired t-test
  n <- min(length(before_intervals), length(after_intervals))
  before_intervals <- before_intervals[1:n]
  after_intervals  <- after_intervals[1:n]
  
  # paired t-test
  t_res <- t.test(after_intervals, before_intervals, paired = TRUE)
  
  # Interpretation
  message("\n--- Synchrony Disruption Test ---\n")
  message(sprintf("Number of usable solo events: %d", n))
  message(sprintf("Mean sqrt(interval BEFORE): %.3f", mean(before_intervals)))
  message(sprintf("Mean sqrt(interval AFTER):  %.3f", mean(after_intervals)))
  
  if (t_res$p.value < 0.05) {
    if (mean(after_intervals) > mean(before_intervals)) {
      message("\nInterpretation: AFTER intervals are significantly LONGER than BEFORE intervals.")
      message("This suggests synchrony is PASSIVE and is disrupted by a solo visit.\n")
    } else {
      message("\nInterpretation: AFTER intervals are significantly SHORTER than BEFORE intervals.")
      message("This suggests parents actively re-align their visits after a solo visit.\n")
    }
  } else {
    message("\nInterpretation: No significant difference between BEFORE and AFTER intervals.")
    message("This suggests synchrony is MAINTAINED even after a solo visit.\n")
  }
  
  return(list(
    before_intervals_sqrt = before_intervals,
    after_intervals_sqrt  = after_intervals,
    t_test = t_res
  ))
}


#002 Mariette & Griffith 2015

Mariette MM, Griffith SC. 2012. The adaptive significance of provisioning and foraging coordination between breeding partners. The American Naturalist 185, 270–280. https://doi.org/10.1086/679441


Methods excerpt

Nest and foraging synchrony were measured as the square root arcsine–transformed proportion of visits where partners were together. Equity (i.e., similarity) in partners’ work rate was measured as 1 minus the absolute difference between male and female visit rates divided by the pair visit rate, so partners that made exactly the same number of visits had an equity score of 1.


Breakdown the logic

Synchrony (nest visits based) measured as the proportion of visits where both partners were together is:

\[ \text{prop\_together} = \frac{\text{together\_visits}}{\text{total male visits} + \text{total female visits}} \]

The square root arcsine transformation is then applied:

\[ \text{synchrony} = \arcsin\!\left(\sqrt{\text{prop\_together}}\right) \]

Equity in patners work rate measured using this formula:

\[ \text{equity} = 1 - \frac{\left| \text{male visits} - \text{female visits} \right|}{\text{male visits} + \text{female visits}} \]


Developed functions

synchrony_mariette1()

synchrony_mariette1 <- function(df, threshold_sec = 120) {
  
  # male-female split
  male_visits   <- df %>% filter(sex == "male")
  female_visits <- df %>% filter(sex == "female")
  
  together_count <- 0
  
  # bouts overlap
  for (i in seq_len(nrow(male_visits))) {
    m_start <- male_visits$arrival_time[i]
    m_end   <- male_visits$departure_time[i]
    
    overlaps <- female_visits %>%
      filter(arrival_time <= m_end & departure_time >= m_start)
    
    if (nrow(overlaps) > 0) {
      for (j in seq_len(nrow(overlaps))) {
        f_start <- overlaps$arrival_time[j]
        f_end   <- overlaps$departure_time[j]
        
        overlap_start <- max(m_start, f_start)
        overlap_end   <- min(m_end, f_end)
        
        overlap_duration <- as.numeric(overlap_end - overlap_start, units = "secs")
        
        if (overlap_duration >= threshold_sec) {
          together_count <- together_count + 1
          break
        }
      }
    }
  }
  
  # proportion of overlapping bouts (partners together)
  total_visits <- nrow(male_visits) + nrow(female_visits)
  prop_together <- if (total_visits > 0) together_count / total_visits else 0
  
  # arcsin(sqrt(x)) transformation
  arcsin_sqrt <- asin(sqrt(prop_together))
  
  # --- Interpretation ----
  message("\n--- Nest Visit Synchrony Interpretation (visit-based) ---\n")
  message(sprintf("Together visits: %d out of %d total visits", together_count, total_visits))
  message(sprintf("Proportion together: %.4f (%.2f%% of visits)", prop_together, prop_together * 100))
  message(sprintf("Arcsin-sqrt transformed synchrony: %.4f", arcsin_sqrt))
  
  if (prop_together > 0.5) {
    message("\nInterpretation: Synchrony is relatively HIGH — partners are often together.")
    message("This suggests strong coordination or pair bonding.\n")
  } else if (prop_together > 0) {
    message("\nInterpretation: Synchrony is MODERATE — partners sometimes overlap.")
    message("This suggests partial coordination but also independent activity.\n")
  } else {
    message("\nInterpretation: Synchrony is LOW — partners never overlapped above threshold.")
    message("This suggests independent, uncoordinated visits.\n")
  }
  
  # outputs
  return(list(
    proportion_together = prop_together,
    arcsin_sqrt = arcsin_sqrt,
    total_male_visits = nrow(male_visits),
    total_female_visits = nrow(female_visits),
    together_visits = together_count
  ))
}

equity_mariette()

equity_mariette <- function(df) {
  
  # Split male and female visits
  male_visits   <- df %>% filter(sex == "male")
  female_visits <- df %>% filter(sex == "female")
  
  # Count visits
  V_m <- nrow(male_visits)
  V_f <- nrow(female_visits)
  V_total <- V_m + V_f
  
  # Equity formula: 1 - |male - female| / total
  equity <- if (V_total > 0) {
    1 - abs(V_m - V_f) / V_total
  } else {
    NA
  }
  
  # --- Interpretation ----
  message("\n--- Equity in Work Rate Interpretation ---\n")
  message(sprintf("Male visits: %d", V_m))
  message(sprintf("Female visits: %d", V_f))
  message(sprintf("Total visits: %d", V_total))
  message(sprintf("Equity score: %.4f", equity))
  
  if (is.na(equity)) {
    message("\nInterpretation: No visits recorded, equity cannot be calculated.\n")
  } else if (equity == 1) {
    message("\nInterpretation: PERFECT equity — both partners contributed equally.\n")
  } else if (equity >= 0.75) {
    message("\nInterpretation: HIGH equity — partners contributed fairly similarly.\n")
  } else if (equity >= 0.5) {
    message("\nInterpretation: MODERATE equity — some imbalance in contributions.\n")
  } else {
    message("\nInterpretation: LOW equity — one partner contributed much more than the other.\n")
  }
  
  # Return numerical results
  return(list(
    male_visits   = V_m,
    female_visits = V_f,
    total_visits  = V_total,
    equity_score  = equity
  ))
}


#003 Wojczulanis-Jakubas et al 2023

Wojczulanis-Jakubas K, Płóciennik J, Guinebretiere A, Hałupka L. 2012. Cooperative parental performance at chick provisioning in a small passerine, the Reed Warbler Acrocephalus scirpaceus. Behavioural Ecology and Sociobiology, 77, 123. https://doi.org/10.1007/s00265-023-03397-5


Methods excerpt

To test whether parents tend to temporally synchronise their feeding visits, we first calculated observed proportion of nest visits when partners met. Then, we applied Monte Carlo randomisation to test how the observed proportion is different from the one that could be expected by chance. Before any calculation, to time interval of each nest visit we added a time margin of 30 s at each side of the interval, to account for a possible partners’ encounter outside the camera view. This two-side margin makes in total of 1 min, and follows recommendations of Leniowski and Węgrzyn (2018a). As in the system of blackcaps (Leniowski and Węgrzyn 2018a), we had numerous observations that RW parents often arrive together with food, but then one of them perches in close vicinity of the nest waiting for the other to complete feeding. Thus, the added time buffer accounts for subsequent visits of the partners, performed shortly one after the other. In other words, the visit was considered to be synchronised if two partners met at the nest or paid their visit within 30 s time interval in respect to each other. For randomisation, we first shuffled the observed sequences of nest visits and inter-visit intervals (with the added time margins To test whether parents tend to temporally synchronise their feeding visits, we first calculated observed proportion of nest visits when partners met. Then, we applied Monte Carlo randomisation to test how the observed proportion is different from the one that could be expected by chance. Before any calculation, to time interval of each nest visit we added a time margin of 30 s at each side of the interval, to account for a possible partners’ encounter outside the camera view. This two-side margin makes in total of 1 min, and follows recommendations of Leniowski and Węgrzyn (2018a). As in the system of blackcaps (Leniowski and Węgrzyn 2018a), we had numerous observations that RW parents often arrive together with food, but then one of them perches in close vicinity of the nest waiting for the other to complete feeding. Thus, the added time buffer accounts for subsequent visits of the partners, performed shortly one after the other. In other words, the visit was considered to be synchronised if two partners met at the nest or paid their visit within 30 s time interval in respect to each other. For randomisation, we first shuffled the observed sequences of nest visits and inter-visit intervals (with the added time margins


Breakdown the logic

Step 1. Observed Synchrony
- For each nest, calculate the observed proportion of visits where both partners encountered each other.
- Encounters are defined as:
- Partners arriving together at the nest, or
- Visits occurring within ±30 seconds of each other.
- A time buffer of 30 seconds is added to both sides of each visit interval (total 1 minute) to account for:
- Possible encounters outside camera view.
- Cases where one parent waits nearby while the other feeds.

Step 2. Randomisation (Monte Carlo)
- To test whether synchrony differs from chance:
- Shuffle the observed sequences of visits and inter-visit intervals (with the added buffer).
- Recalculate the synchrony proportion for each randomised dataset.
- Repeat this procedure 1000 times to generate a distribution of randomised synchrony proportions.

Step 3. Statistical Testing
- Compare the observed synchrony proportion to the randomised distribution. - Calculate p-values:
- Right-tail: probability observed ≥ randomised.
- Left-tail: probability observed ≤ randomised.
- Two-tailed: probability observed differs from expected mean by at least as much as observed.

Step 4. Synchronisation Index
- Define the synchronisation index as:

\[ \text{Index} = \frac{\text{obs} - \text{mean(exp)}}{\text{mean(exp)}} \]

  • Interpretation:
    • Positive values → Partners synchronise more than expected by chance.
    • Zero values → Partners behave randomly, consistent with independence.
    • Negative values → Partners avoid synchronisation, rarely meeting at the nest.

Step 5. Overall p-value
- For multiple nests, combine individual p-values using Fisher’s method (via the R package metap).
- This provides an overall test of synchronisation across the dataset.

Step 6. Interpretation
- Observed synchrony higher than expected → evidence of active coordination.
- Observed synchrony lower than expected → evidence of avoidance or out-of-phase behaviour.
- Synchrony index magnitude reflects the intensity of synchronisation.


Developed functions

synchrony_wojczulanis1()

synchrony_wojczulanis1 <- function(df, buffer_sec = 30, n_sim = 1000) {
  # df: data.frame with columns sex, arrival_time, departure_time (POSIXct), single nest
  
  library(ggplot2)
  
  # --- prep ---
  df <- df[order(df$arrival_time), ]
  sexes <- sort(unique(df$sex))
  
  # convert to numeric seconds relative to min time
  t0 <- min(df$arrival_time)
  arr_sec  <- as.numeric(difftime(df$arrival_time, t0, units = "secs"))
  departure_sec <- as.numeric(difftime(df$departure_time, t0, units = "secs"))
  
  # add ±30s buffer to each visit (total 1 min margin)
  df$start_ext <- arr_sec - buffer_sec
  df$end_ext   <- departure_sec + buffer_sec
  
  # --- helper: synchrony proportion from extended intervals ---
  synchrony_prop <- function(df) {
    s1 <- sexes[1]
    s2 <- sexes[2]
    has_partner <- logical(nrow(df))
    for (i in seq_len(nrow(df))) {
      this_sex <- df$sex[i]
      opp_sex  <- if (this_sex == s1) s2 else s1
      opp_df   <- df[df$sex == opp_sex, ]
      st <- df$start_ext[i]
      en <- df$end_ext[i]
      overlap <- (pmax(st, opp_df$start_ext) < pmin(en, opp_df$end_ext))
      has_partner[i] <- any(overlap)
    }
    mean(has_partner)
  }
  
  # --- helper: get visit and gap durations per sex ---
  get_durations <- function(df_sex) {
    df_sex <- df_sex[order(df_sex$arrival_time), ]
    arr <- as.numeric(difftime(df_sex$arrival_time, t0, units = "secs"))
    departure <- as.numeric(difftime(df_sex$departure_time, t0, units = "secs"))
    
    # add buffer to each visit duration
    visit_dur <- (departure - arr) + 2 * buffer_sec
    if (length(arr) > 1) {
      gap_dur <- arr[-1] - departure[-length(departure)]
    } else {
      gap_dur <- numeric(0)
    }
    list(visit_dur = visit_dur, gap_dur = gap_dur, n_visits = length(arr))
  }
  
  # --- helper: simulate sequence for one sex ---
  simulate_sex <- function(durs) {
    k <- durs$n_visits
    if (k == 0) return(data.frame(start_ext = numeric(0), end_ext = numeric(0)))
    
    v <- sample(durs$visit_dur, size = k, replace = FALSE)
    g <- if (k > 1) sample(durs$gap_dur, size = k - 1, replace = FALSE) else numeric(0)
    
    start_ext <- numeric(k)
    end_ext   <- numeric(k)
    start_ext[1] <- 0
    end_ext[1]   <- start_ext[1] + v[1]
    
    if (k > 1) {
      for (i in 2:k) {
        start_ext[i] <- end_ext[i - 1] + g[i - 1]
        end_ext[i]   <- start_ext[i] + v[i]
      }
    }
    data.frame(start_ext = start_ext, end_ext = end_ext)
  }
  
  # --- observed synchrony ---
  obs_sync <- synchrony_prop(df)
  
  # --- durations per sex ---
  d_sex1 <- get_durations(df[df$sex == sexes[1], ])
  d_sex2 <- get_durations(df[df$sex == sexes[2], ])
  
  # --- Monte Carlo randomisation ---
  sim_sync <- numeric(n_sim)
  for (r in seq_len(n_sim)) {
    sim1 <- simulate_sex(d_sex1)
    sim2 <- simulate_sex(d_sex2)
    
    sim_df <- rbind(
      data.frame(sex = sexes[1], start_ext = sim1$start_ext, end_ext = sim1$end_ext),
      data.frame(sex = sexes[2], start_ext = sim2$start_ext, end_ext = sim2$end_ext)
    )
    sim_sync[r] <- synchrony_prop(sim_df)
  }
  
  # --- p-values ---
  p_right <- mean(sim_sync >= obs_sync)
  p_left  <- mean(sim_sync <= obs_sync)
  p_two   <- mean(abs(sim_sync - mean(sim_sync)) >= abs(obs_sync - mean(sim_sync)))
  
  # --- synchrony index ---
  exp_mean <- mean(sim_sync)
  sync_index <- (obs_sync - exp_mean) / exp_mean
  
  # --- interpretation ----
  message("\n--- Synchrony Monte Carlo Interpretation ---\n")
  message(sprintf("Observed synchrony proportion: %.4f", obs_sync))
  message(sprintf("Expected synchrony (mean randomised): %.4f", exp_mean))
  message(sprintf("Synchrony index: %.4f", sync_index))
  message(sprintf("Monte Carlo p-values: right = %.4f, left = %.4f, two-tailed = %.4f",
                  p_right, p_left, p_two))
  
  if (sync_index > 0) {
    message("\nInterpretation: Positive synchrony index — partners show MORE synchronisation than expected by chance.\n")
  } else if (sync_index < 0) {
    message("\nInterpretation: Negative synchrony index — partners show LESS synchronisation than expected, suggesting avoidance.\n")
  } else {
    message("\nInterpretation: Synchrony index near zero — partners behave as if independent (random performance).\n")
  }
  
  # --- plot ---
  sync_plot <- ggplot(data.frame(sim_sync = sim_sync), aes(x = sim_sync)) +
    geom_histogram(color = "black", fill = "skyblue", bins = 30) +
    geom_vline(xintercept = obs_sync, color = "red", linewidth = 1.2) +
    labs(
      title = "Randomised Synchrony Distribution",
      x = "Synchrony proportion",
      y = "Count"
    ) +
    theme_minimal()
  
  print(sync_plot)
  
  # --- return ---
  list(
    observed_synchrony = obs_sync,
    expected_synchrony = exp_mean,
    synchrony_index = sync_index,
    simulated_synchronies = sim_sync,
    p_right = p_right,
    p_left = p_left,
    p_two_tailed = p_two,
    plot = sync_plot
  )
}