Reproducible workflow for generating the main visualizations

This appendix provides the code used to generate the main figures in the manuscript. It includes setup, data preparation, model-based prediction workflows, plot construction, and export code for the final multi-panel figures.

Most chunks are included as a reproducible workflow record and are not executed during rendering of the website.

Plotting utilities

This section defines reusable plotting objects to ensure a consistent visual style across figures.

Custom theme

Show custom plotting theme
custom_theme <- function(show_legend = FALSE) {
  # Set legend position string
  legend_pos <- if (show_legend) "bottom" else "none"
  
  # Use theme_classic() as the base
  theme_classic(
    base_size = 10,
    base_family = "Arial"
  ) +
  theme(
    # --- Add minor ticks back in ---
    # This is the line that fixes your problem
    axis.ticks.minor = element_line(color = "black", linewidth = 0.25), 
    
    # --- Legend and margins ---
    legend.position = legend_pos,
    plot.margin = margin(5, 5, 5, 5),
    
    # --- Backgrounds (optional, theme_classic is not transparent) ---
    panel.background = element_rect(fill = 'transparent', colour = NA),
    plot.background = element_rect(fill = 'transparent', colour = NA),
    
    # --- Text and Titles (your original customizations) ---
    axis.title.x = element_text(margin = margin(t = 5)),
    axis.title.y = element_text(margin = margin(r = 5)),
    axis.text.x = element_text(margin = margin(t = 2)),
    axis.text.y = element_text(hjust = 1, margin = margin(r = 2)),
    plot.title = element_text(hjust = 0.5, face = "bold", size = rel(1.2), margin = margin(b = 8)),
    plot.subtitle = element_text(hjust = 0.5, size = rel(1.0), margin = margin(b = 8)),
    
    # --- Axis Lines (already in theme_classic, but good to keep) ---
    axis.line = element_line(color = "black"),
    axis.ticks = element_line(color = "black")
  )
}

Shared x-axis scale

Show shared x-axis scale
# Create a reusable x-axis scale for all plots
my_x_scale <- scale_x_continuous(
  name = "Time (years)", # Set the axis title here
  breaks = seq(2003, 2023, 2),
  minor_breaks = seq(2003, 2023, 1),
  labels = function(breaks) {
    # This function shows the label only if the year is the first one (2003)
    # or is a multiple of 4 years from the start.
    ifelse((breaks - 2003) %% 4 == 0, as.character(breaks), "")
  },
  guide = guide_axis(minor.ticks = TRUE)
)

Load data and fitted models

Data and pre-trained models were loaded from project files prior to figure generation.

Show bloom data and model loading
bloom_dat_raw <- read.csv(here("..","data", "bloom_dates.csv"))
bloom_model <- readRDS(here("..","Rdata", "model_delta5.rds"))
Show female data and model loading
female_dat_raw <- read.csv(here("..","data", "Lay_date_mismatches_250702.csv"))
mismatch_modelf <- readRDS(here("..","Rdata", "female_fit_m11.rds"))
fledging_modelf <- readRDS(here("..","Rdata", "female_fit_Q17.rds"))
Show male data and model loading
male_dat_raw <- read.csv(here("..","data", "Lay_date_mismatches_250702.csv"))
mismatch_modelm <- readRDS(here("..","Rdata", "male_fit_m11.rds"))
fledging_modelm <- readRDS(here("..","Rdata", "male_fit_Q17.rds"))

Data preprocessing

Raw data were transformed to create the variables used for plotting and model-based prediction.

Bloom timing data

Show bloom preprocessing
# --- Process Bloom Data ---
bloom_dat <- bloom_dat_raw %>%
  mutate(
    bloom_start_thr5_as_date = as.Date(bloom_start_thr5, format = "%d/%m/%Y"),
    winter_start_date = as.Date(paste0(season_year - 1, "-12-21")),
    delta5 = as.numeric(difftime(bloom_start_thr5_as_date, winter_start_date, units = "days")),
    season_year = as.numeric(season_year),
    TIME = as.numeric(scale(season_year, center = TRUE, scale = FALSE))
  )
Show female preprocessing
# --- Process Female Mismatch & Fledging Data ---
dat_female_mismatch <- female_dat_raw %>%
  filter(SEX == 1) %>%
  mutate(
    ori_age = as.numeric(AGE),
    ori_time = as.numeric(YEAR),
    TIME = scale(YEAR, center = TRUE, scale = FALSE),
    AGE = scale(as.numeric(AGE), center = TRUE, scale = FALSE),
    AFR = scale(as.numeric(AFR), center = TRUE, scale = FALSE),
    LONGEVITY = scale(as.numeric(LONGEVITY), center = TRUE, scale = FALSE),
    YEAR = as.factor(YEAR),
    RING = as.factor(RING)
  )

dat_female_fledging <- female_dat_raw %>%
  filter(SEX == 1) %>%
  mutate(
    ori_age = as.numeric(AGE),
    ori_time = as.numeric(YEAR),
    TIME = scale(YEAR, center = TRUE, scale = FALSE),
    ZX5MISMATCH = scale(as.numeric(X5MISMATCH), center = TRUE, scale = FALSE),
    AGE = scale(as.numeric(AGE), center = TRUE, scale = FALSE),
    AFR = scale(as.numeric(AFR), center = TRUE, scale = FALSE),
    LONGEVITY = scale(as.numeric(LONGEVITY), center = TRUE, scale = FALSE),
    FLEDS = factor(FLEDS, levels = c(0, 1, 2, 3), ordered = TRUE)
  )
Show male preprocessing
dat_male_mismatch <- male_dat_raw %>%
  filter(SEX == 0) %>%
  mutate(
    ori_age = as.numeric(AGE),
    ori_time = as.numeric(YEAR),
    TIME = scale(YEAR, center = TRUE, scale = FALSE),
    AGE = scale(as.numeric(AGE), center = TRUE, scale = FALSE),
    AFR = scale(as.numeric(AFR), center = TRUE, scale = FALSE),
    LONGEVITY = scale(as.numeric(LONGEVITY), center = TRUE, scale = FALSE),
    YEAR = as.factor(YEAR),
    RING = as.factor(RING)
  )

dat_male_fledging <- male_dat_raw %>%
  filter(SEX == 0) %>%
  mutate(
    ori_age = as.numeric(AGE),
    ori_time = as.numeric(YEAR),
    TIME = scale(YEAR, center = TRUE, scale = FALSE),
    ZX5MISMATCH = scale(as.numeric(X5MISMATCH), center = TRUE, scale = FALSE),
    AGE = scale(as.numeric(AGE), center = TRUE, scale = FALSE),
    AFR = scale(as.numeric(AFR), center = TRUE, scale = FALSE),
    LONGEVITY = scale(as.numeric(LONGEVITY), center = TRUE, scale = FALSE),
    FLEDS = factor(FLEDS, levels = c(0, 1, 2, 3), ordered = TRUE))

Figure 3: The Effect of Female/males Age

This figure summarizes age-related patterns in phenological mismatch and fledging success, together with the association between mismatch and fledging outcomes.

mean_of_original_age_female <- mean(dat_female_mismatch$ori_age)

age_draws_female <- mismatch_modelf %>%
  epred_draws(
    newdata = expand.grid(
      AGE = seq(min(dat_female_mismatch$AGE), max(dat_female_mismatch$AGE), length.out = 100),
      TIME = mean(dat_female_mismatch$TIME),
      AFR = mean(dat_female_mismatch$AFR),
      LONGEVITY = mean(dat_female_mismatch$LONGEVITY)
    ),
    re_formula = NA,
    dpar = TRUE
  ) %>%
  mutate(actual_age = AGE + mean_of_original_age_female)

dat_count_mismatch_age_female <- dat_female_mismatch %>%
  group_by(ori_age, X5MISMATCH) %>%
  summarize(n = n(), .groups = 'drop')

dat_count_sigma_age_female <- dat_with_residuals_female %>%
  group_by(ori_age, raw_sigma) %>%
  summarize(n = n(), .groups = 'drop')
mean_of_original_age_male <- mean(dat_male_mismatch$ori_age)
age_draws_male <- mismatch_modelm %>%
  epred_draws(
    newdata = expand.grid(
      AGE = seq(min(dat_male_mismatch$AGE), max(dat_male_mismatch$AGE), length.out = 100),
      TIME = mean(dat_male_mismatch$TIME),
      AFR = mean(dat_male_mismatch$AFR),
      LONGEVITY = mean(dat_male_mismatch$LONGEVITY)
    ),
    re_formula = NA,
    dpar = TRUE
  ) %>%
  mutate(actual_age = AGE + mean_of_original_age_male)

dat_count_mismatch_age_male <- dat_male_mismatch %>%
  group_by(ori_age, X5MISMATCH) %>%
  summarize(n = n(), .groups = 'drop')


dat_count_sigma_age_male <- dat_with_residuals_male %>%
  group_by(ori_age, raw_sigma) %>%
  summarize(n = n(), .groups = 'drop')

To visualize individual variation, predictions for a random sample of 20 individual birds are generated, showing their unique age-related trajectories.

set.seed(42) # For reproducible random sampling
num_to_sample <- 50
rings_to_plot_female <- sample(unique(dat_female_mismatch$RING), num_to_sample)

newdata_grid_female <- expand.grid(
  AGE = seq(min(dat_female_mismatch$AGE), max(dat_female_mismatch$AGE), length.out = 100),
  RING = rings_to_plot_female,
  TIME = mean(dat_female_mismatch$TIME),
  AFR = mean(dat_female_mismatch$AFR),
  LONGEVITY = mean(dat_female_mismatch$LONGEVITY)
)

individual_predictions_female <- mismatch_modelf %>%
  add_linpred_draws(newdata = newdata_grid_female, re_formula = ~(1 + AGE | RING)) %>%
  group_by(RING, AGE) %>%
  summarise(.value = median(.linpred), .groups = "drop") %>%
  mutate(actual_age = AGE + mean_of_original_age_female)
set.seed(42) 
num_to_sample <- 50
rings_to_plot_male <- sample(unique(dat_male_mismatch$RING), num_to_sample)


newdata_grid_male <- expand.grid(
  AGE = seq(min(dat_male_mismatch$AGE), max(dat_male_mismatch$AGE), length.out = 100),
  RING = rings_to_plot_male,
  TIME = mean(dat_male_mismatch$TIME),
  AFR = mean(dat_male_mismatch$AFR),
  LONGEVITY = mean(dat_male_mismatch$LONGEVITY)
)

individual_predictions_male <- mismatch_modelm %>%
  add_linpred_draws(newdata = newdata_grid_male, re_formula = ~(1 + AGE | RING)) %>%
  group_by(RING, AGE) %>%
  summarise(.value = median(.linpred), .groups = "drop") %>%
  mutate(actual_age = AGE + mean_of_original_age_male)

The next two chunks create Plot G (mismatch vs. age) and Plot H (consistency vs. age).

# --- Panel G: Mismatch vs. Age ---
plot_G <- ggplot() +
  geom_point(data = dat_count_mismatch_age_female, aes(x = ori_age, y = X5MISMATCH, size = n),
             shape = 21, alpha = 0.2, color = "#9096b1", fill = "#d2bfb8") +
  stat_lineribbon(data = age_draws_female, aes(x = actual_age, y = .epred),
                  .width = 0.95, alpha = 0.4, fill = "#a0786e") +
  #geom_line(data = individual_predictions_female, aes(x = actual_age, y = .value, group = RING),
            #color = "#D95F02", linewidth = 0.5, alpha = 0.6)+
  stat_lineribbon(data = age_draws_female, aes(x = actual_age, y = .epred),
                  fill = NA, color = "black", linewidth = 1) +
  scale_x_continuous(breaks = seq(1, 23, 2),limits = c(1,23),guide = guide_axis(minor.ticks = TRUE)) +
  scale_y_continuous(breaks = seq(-90, 210, 30),minor_breaks = seq(-90,210,15), limits = c(-90, 210),guide = guide_axis(minor.ticks = TRUE)) +
  labs(
    subtitle = "\u03B2 = 0.334, 95% CI = [0.291, 0.377]",
    y = "Phenological mismatch (days)",
    x = "Age (years)"
  ) +
  custom_theme() + 
  theme(legend.position = "none")+theme(plot.subtitle = element_text(face="bold"))

# --- Panel H: Consistency vs. Age ---
plot_H <- ggplot() +
  geom_jitter(data = dat_count_sigma_age_female, aes(x = ori_age, y = raw_sigma), width = 0.2, height = 0, shape = 21, alpha = 0.2,stroke=0.1, fill = "#b1ab90", colour = "#b8cbd2",size=1) +
  stat_lineribbon(data = age_draws_female, aes(x = actual_age, y = sigma),
                  .width = 0.95, alpha = 0.4, fill = "#cbba62") +
  stat_lineribbon(data = age_draws_female, aes(x = actual_age, y = sigma),
                  fill = NA, color = "black", linewidth = 1) +
  scale_x_continuous(breaks = seq(1, 23, 2),limits = c(1,23),guide = guide_axis(minor.ticks = TRUE)) +
  scale_y_continuous(breaks = seq(0, 150, 30), limits = c(0, 150),guide = guide_axis(minor.ticks = TRUE)) +
  labs(
    subtitle = "\u03B2 = 0.003, 95% CI = [0.001, 0.004]",
    y = "Inconsistency (\u03C3)",
    x = "Age (years)"
  ) +
  custom_theme() +
  theme(legend.position = "none")+theme(plot.subtitle = element_text(face="bold"))
# --- Panel G: Mismatch vs. Age ---
plot_Gm <- ggplot() +
  geom_point(data = dat_count_mismatch_age_male, aes(x = ori_age, y = X5MISMATCH, size = n),
             shape = 21, alpha = 0.2, color = "#9096b1", fill = "#d2bfb8") +
  stat_lineribbon(data = age_draws_female, aes(x = actual_age, y = .epred),
                  .width = 0.95, alpha = 0.4, fill = "#a0786e") +
 # geom_line(data = individual_predictions_female, aes(x = actual_age, y = .value, group = RING),
  #          color = "#D95F02", linewidth = 0.5, alpha = 0.6)+
  stat_lineribbon(data = age_draws_male, aes(x = actual_age, y = .epred),
                  fill = NA, color = "black", linewidth = 1) +
  scale_x_continuous(breaks = seq(1, 23, 2),limits = c(1,23),guide = guide_axis(minor.ticks = TRUE)) +
  scale_y_continuous(breaks = seq(-90, 210, 30),minor_breaks = seq(-90,210,15), limits = c(-90, 210),guide = guide_axis(minor.ticks = TRUE)) +
  labs(
    subtitle = "\u03B2 = 0.283, 95% CI = [0.235, 0.332]",
    y = "Phenological mismatch (days)",
    x = "Age (years)"
  ) +
  custom_theme() + 
  theme(legend.position = "none")+theme(plot.subtitle = element_text(face="bold"))

# --- Panel H: Consistency vs. Age ---
plot_Hm <- ggplot() +
  geom_jitter(data = dat_count_sigma_age_male, aes(x = ori_age, y = raw_sigma), width = 0.2, height = 0, shape = 21, alpha = 0.2,stroke=0.1, fill = "#b1ab90", colour = "#b8cbd2",size=1) +
  stat_lineribbon(data = age_draws_male, aes(x = actual_age, y = sigma),
                  .width = 0.95, alpha = 0.4, fill = "#cbba62") +
  stat_lineribbon(data = age_draws_male, aes(x = actual_age, y = sigma),
                  fill = NA, color = "black", linewidth = 1) +
  scale_x_continuous(breaks = seq(1, 23, 2),limits = c(1,23),guide = guide_axis(minor.ticks = TRUE)) +
  scale_y_continuous(breaks = seq(0, 150, 30), limits = c(0, 150),guide = guide_axis(minor.ticks = TRUE)) +
  labs(
    subtitle = "\u03B2 = 0.004, 95% CI = [0.003, 0.005]",
    y = "Inconsistency (\u03C3)",
    x = "Age (years)"
  ) +
  custom_theme() +
  theme(legend.position = "none")

Predictions are now generated from the fledgling success model across the range of female ages.

mean_of_original_fledging_agef <- mean(dat_female_fledging$ori_age)

fledging_age_drawsf <- fledging_modelf %>%
  epred_draws(
    newdata = expand.grid(
      AGE = seq(min(dat_female_fledging$AGE), max(dat_female_fledging$AGE), length.out = 100),
      TIME = mean(dat_female_fledging$TIME),
      ZX5MISMATCH = mean(dat_female_fledging$ZX5MISMATCH),
      AFR = mean(dat_female_fledging$AFR),
      LONGEVITY = mean(dat_female_fledging$LONGEVITY)
    ),
    re_formula = NA,
    dpar = TRUE 
  ) %>%
  mutate(actual_age = AGE + mean_of_original_fledging_agef)

fledging_success_draws_agef <- fledging_age_drawsf %>%
  filter(.category != "0") %>%
  group_by(across(-c(.epred, .category))) %>% # Groups by all cols except .epred & .category
  summarize(prob_success = sum(.epred), .groups = 'drop')

fledging_success_props_agef <- dat_female_fledging%>%
  mutate(SUCCESS = ifelse(FLEDS == "0", "Failure (0)", "Success (1+)")) %>%
  count(ori_age, SUCCESS) %>%
  tidyr::complete(ori_age, 
                  SUCCESS = c("Failure (0)", "Success (1+)"), 
                  fill = list(n = 0)) %>%
  group_by(ori_age) %>%
  mutate(
         total_n_year = sum(n),
         proportion = ifelse(total_n_year > 0, n / total_n_year, 0)) %>%
  filter(SUCCESS == "Success (1+)") %>%
  ungroup()
mean_of_original_fledging_agem <- mean(dat_male_fledging$ori_age)

fledging_age_drawsm <- fledging_modelm %>%
  epred_draws(
    newdata = expand.grid(
      AGE = seq(min(dat_male_fledging$AGE), max(dat_male_fledging$AGE), length.out = 100),
      TIME = mean(dat_male_fledging$TIME),
      ZX5MISMATCH = mean(dat_male_fledging$ZX5MISMATCH),
      AFR = mean(dat_male_fledging$AFR),
      LONGEVITY = mean(dat_male_fledging$LONGEVITY)
    ),
    re_formula = NA,
    dpar = TRUE 
  ) %>%
  mutate(actual_age = AGE + mean_of_original_fledging_agem)


fledging_success_draws_agem <- fledging_age_drawsm %>%
  filter(.category != "0") %>%
  group_by(across(-c(.epred, .category))) %>% # Groups by all cols except .epred & .category
  summarize(prob_success = sum(.epred), .groups = 'drop')

fledging_success_props_agem <- dat_male_fledging%>%
  mutate(SUCCESS = ifelse(FLEDS == "0", "Failure (0)", "Success (1+)")) %>%
  count(ori_age, SUCCESS) %>%
  tidyr::complete(ori_age, 
                  SUCCESS = c("Failure (0)", "Success (1+)"), 
                  fill = list(n = 0)) %>%
  group_by(ori_age) %>%
  mutate(
         total_n_year = sum(n),
         proportion = ifelse(total_n_year > 0, n / total_n_year, 0)) %>%
  filter(SUCCESS == "Success (1+)") %>%
  ungroup()

The next chunks build the bottom row of the figure: fledgling probabilities vs. age.

# --- PanelA: Fledging Probability vs. Age ---
plot_fledging_success_age <- ggplot() +

  geom_point(data = fledging_success_props_agef, 
             aes(x = ori_age, y = proportion,size=total_n_year), 
             shape = 21, alpha = 0.4, stroke = 1, fill = "#9dbcbe", colour = "#a68551") +
  stat_lineribbon(data = fledging_success_draws_agef, 
                  aes(x = actual_age, y = prob_success), 
                  .width = 0.95, alpha = 0.4, fill = "#5172a6") + 
  stat_lineribbon(data = fledging_success_draws_agef, 
                  aes(x = actual_age, y = prob_success), 
                  fill = NA, color = "black", linewidth = 1) +
  

  scale_x_continuous(
    breaks = seq(1, 23, 2),
    minor_breaks = seq(1, 23, 1),
    guide = guide_axis(minor.ticks = TRUE)
  ) +
  scale_y_continuous(
    breaks = seq(0, 1, 0.2), 
    minor_breaks = seq(0, 1, 0.1),
    limits = c(0, 1),
    guide = guide_axis(minor.ticks = TRUE)
  ) +

  labs(
    subtitle = "\u03B2 = -0.075, 95% CI = [-0.156, -0.031]",
    y = "Probability of nest success",
    x = "Age (years)"
  ) +

  custom_theme(show_legend = FALSE)+theme(plot.subtitle = element_text(face="bold"))
plot_fledging_success_agem <- ggplot() +

  geom_point(data = fledging_success_props_agem, 
             aes(x = ori_age, y = proportion, size=total_n_year), 
             shape = 21, alpha = 0.4, stroke = 1, fill = "#9dbcbe", colour = "#a68551") +
  stat_lineribbon(data = fledging_success_draws_agem, 
                  aes(x = actual_age, y = prob_success), 
                  .width = 0.95, alpha = 0.4, fill = "#5172a6") + 
  stat_lineribbon(data = fledging_success_draws_agem, 
                  aes(x = actual_age, y = prob_success), 
                  fill = NA, color = "black", linewidth = 1) +
  

  scale_x_continuous(
    breaks = seq(1, 23, 2),
    minor_breaks = seq(1, 23, 1),
    guide = guide_axis(minor.ticks = TRUE)
  ) +
  scale_y_continuous(
    breaks = seq(0, 1, 0.2), 
    minor_breaks = seq(0, 1, 0.1),
    limits = c(0, 1),
    guide = guide_axis(minor.ticks = TRUE)
  ) +

  labs(
    subtitle = "\u03B2 = -0.039, 95% CI = [-0.088, -0.014]",
    y = "Probability of nest success",
    x = "Age (years)"
  ) +
  

  custom_theme(show_legend = FALSE)+theme(plot.subtitle = element_text(face="bold"))

This panel shows the direct relationship between phenological mismatch and the probability of fledgling success.

Model predictions are generated for fledgling success across the full range of observed mismatch values.

mean_of_original_mismatchf <- mean(dat_female_fledging$X5MISMATCH, na.rm = TRUE)
# --- 1. Generate Raw Predictions ---
mismatch_drawsf <- fledging_modelf %>%
  epred_draws(
    newdata = expand.grid(
      ZX5MISMATCH = seq(min(dat_female_fledging$ZX5MISMATCH), max(dat_female_fledging$ZX5MISMATCH), length.out = 100),
      TIME = mean(dat_female_fledging$TIME),
      AGE = mean(dat_female_fledging$AGE),
      AFR = mean(dat_female_fledging$AFR),
      LONGEVITY = mean(dat_female_fledging$LONGEVITY)
    ),
    re_formula = NA,
    dpar = TRUE
  ) %>%
  mutate(actual_mismatch = ZX5MISMATCH + mean_of_original_mismatchf)


# --- 2. Aggregate Predictions into 3 Categories (0, 1-2, 3) ---
fledging_grouped_drawsf <- mismatch_drawsf %>%
  # Create the new grouping variable based on .category
  mutate(fled_group = case_when(
    .category == "0" ~ "0",
    .category %in% c("1", "2") ~ "1-2", # Grouping 1 and 2 to see the quadratic effect
    .category == "3" ~ "3",
    TRUE ~ NA_character_ # Just in case there are >3, or handle accordingly
  )) %>%
  # Filter out any categories we aren't interested in (if any exist)
  filter(!is.na(fled_group)) %>%
  # Sum the probabilities (.epred) within the new groups for every draw
  group_by(actual_mismatch, .draw, fled_group) %>% 
  summarize(prob_val = sum(.epred), .groups = 'drop') 

# --- 3. Calculate Observed Proportions for the 3 Categories ---
fledging_grouped_propsf <- dat_female_fledging %>%
  filter(!is.na(X5MISMATCH)) %>%
  # Create the same grouping variable in the raw data
  mutate(
    fled_group = case_when(
      FLEDS == "0" ~ "0",
      FLEDS %in% c("1", "2") ~ "1-2",
      FLEDS == "3" ~ "3",
      TRUE ~ "Other"
    )
  ) %>%
  filter(fled_group %in% c("0", "1-2", "3")) %>%
  # Count occurrences per mismatch value and group
  count(X5MISMATCH, fled_group) %>%
  # Ensure all groups exist for all mismatch values (fill with 0 if missing)
  tidyr::complete(X5MISMATCH, 
                  fled_group = c("0", "1-2", "3"), 
                  fill = list(n = 0)) %>%
  # Calculate proportions
  group_by(X5MISMATCH) %>%
  mutate(total_n_year = sum(n),
         proportion = ifelse(total_n_year > 0, n / total_n_year, 0)) %>%
  ungroup()

# --- 4. Generate Label Data for Plotting ---
group_mismatch_label_dataf <- fledging_grouped_drawsf %>%
  filter(actual_mismatch == min(actual_mismatch) | actual_mismatch == max(actual_mismatch)) %>%
  group_by(actual_mismatch, fled_group) %>%
  summarise(mean_prob_perc = mean(prob_val) * 100, .groups = 'drop')

group_mismatch_label_dataf
mean_of_original_mismatchm <- mean(dat_male_fledging$X5MISMATCH, na.rm = TRUE)

mismatch_drawsm <- fledging_modelm %>%
  epred_draws(
    newdata = expand.grid(
      ZX5MISMATCH = seq(min(dat_male_fledging$ZX5MISMATCH), max(dat_male_fledging$ZX5MISMATCH), length.out = 100),
      TIME = mean(dat_male_fledging$TIME),
      AGE = mean(dat_male_fledging$AGE),
      AFR = mean(dat_male_fledging$AFR),
      LONGEVITY = mean(dat_male_fledging$LONGEVITY)
    ),
    re_formula = NA,
    dpar = TRUE
  ) %>%
  mutate(actual_mismatch = ZX5MISMATCH + mean_of_original_mismatchm)


# --- 2. Aggregate Predictions into 3 Categories (0, 1-2, 3) ---
fledging_grouped_drawsm <- mismatch_drawsm %>%
  # Create the new grouping variable based on .category
  mutate(fled_group = case_when(
    .category == "0" ~ "0",
    .category %in% c("1", "2") ~ "1-2", # Grouping 1 and 2 to see the quadratic effect
    .category == "3" ~ "3",
    TRUE ~ NA_character_ # Just in case there are >3, or handle accordingly
  )) %>%

  filter(!is.na(fled_group)) %>%
  # Sum the probabilities (.epred) within the new groups for every draw
  group_by(actual_mismatch, .draw, fled_group) %>% 
  summarize(prob_val = sum(.epred), .groups = 'drop') 

# --- 3. Calculate Observed Proportions for the 3 Categories ---
fledging_grouped_propsm <- dat_male_fledging %>%
  filter(!is.na(X5MISMATCH)) %>%
  # Create the same grouping variable in the raw data
  mutate(
    fled_group = case_when(
      FLEDS == "0" ~ "0",
      FLEDS %in% c("1", "2") ~ "1-2",
      FLEDS == "3" ~ "3",
      TRUE ~ "Other"
    )
  ) %>%
  filter(fled_group %in% c("0", "1-2", "3")) %>%
  # Count occurrences per mismatch value and group
  count(X5MISMATCH, fled_group) %>%
  # Ensure all groups exist for all mismatch values (fill with 0 if missing)
  tidyr::complete(X5MISMATCH, 
                  fled_group = c("0", "1-2", "3"), 
                  fill = list(n = 0)) %>%
  # Calculate proportions
  group_by(X5MISMATCH) %>%
  mutate(total_n_year = sum(n),
         proportion = ifelse(total_n_year > 0, n / total_n_year, 0)) %>%
  ungroup()

# --- 4. Generate Label Data for Plotting ---
group_mismatch_label_datam <- fledging_grouped_drawsm %>%
  filter(actual_mismatch == min(actual_mismatch) | actual_mismatch == max(actual_mismatch)) %>%
  group_by(actual_mismatch, fled_group) %>%
  summarise(mean_prob_perc = mean(prob_val) * 100, .groups = 'drop')

group_mismatch_label_datam
plot_mismatch_prob_altf <- ggplot() + 
  # 1. Observed Data Points (binned proportions)
  # Mapped fill to fled_group to see observed 0s vs 1-2s vs 3s
  geom_point(data = fledging_grouped_propsf, 
             aes(x = X5MISMATCH, y = proportion, size = total_n_year, fill = fled_group), 
             shape = 21, alpha = 0.2, stroke = 0.5, color = "black") +

  # 2. Posterior Ribbons (95% CI)
  stat_lineribbon(data = fledging_grouped_drawsf, 
                  aes(x = actual_mismatch, y = prob_val, fill = fled_group), 
                  .width = 0.95, alpha = 0.3, color = NA) +

  # 3. Posterior Mean Lines
  # We use .width = 0 to draw just the line without another ribbon on top
  stat_lineribbon(data = fledging_grouped_drawsf, 
                  aes(x = actual_mismatch, y = prob_val, color = fled_group), 
                  .width = 0, fill = NA, linewidth = 1) +

  # 4. Scales and Colors
  scale_y_continuous(
    breaks = seq(0, 1, 0.2), 
    minor_breaks = seq(0, 1, 0.1),  
    limits = c(0, 1),
    guide = guide_axis(minor.ticks = TRUE) 
  ) +
  scale_x_continuous(
    breaks = seq(-90, 210, by = 30), 
    minor_breaks = seq(-90, 210, by = 15), 
    limits = c(-90, 210),
    guide = guide_axis(minor.ticks = TRUE) 
  ) +
  
  # Manual colors: Adjust these hex codes to your preference
  scale_fill_manual(values = c("0" = "#999999", "1-2" = "#E69F00", "3" = "#56B4E9"), name = "Number of fledglings") +
  scale_color_manual(values = c("0" = "#999999", "1-2" = "#E69F00", "3" = "#56B4E9"), name = "Number of fledglings") +

  # 5. Labels
  labs(
    subtitle = "\u03B2 = -0.0003, 95% CI = [-0.0007, -0.0001]",
    y = "Probability of fledging outcomes",
    x = "Phenological mismatch (days)"
  ) +
  guides(size = "none",
    # title.position = "left" puts the title next to the keys
    # nrow = 1 forces all keys to be in one line
    fill = guide_legend(title.position = "left", nrow = 1),
    color = guide_legend(title.position = "left", nrow = 1))+
  custom_theme() + theme(
    plot.subtitle = element_text(face = "bold", color = "black"),
    axis.text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    
    # Move legend to bottom
    legend.position = "bottom",
    # Ensure items are not too spaced out
    legend.spacing.x = unit(0.5, "cm")
  )
plot_mismatch_prob_altm <- ggplot() + 
  # 1. Observed Data Points (binned proportions)
  # Mapped fill to fled_group to see observed 0s vs 1-2s vs 3s
  geom_point(data = fledging_grouped_propsm, 
             aes(x = X5MISMATCH, y = proportion, size = total_n_year, fill = fled_group), 
             shape = 21, alpha = 0.2, stroke = 0.5, color = "black") +

  # 2. Posterior Ribbons (95% CI)
  stat_lineribbon(data = fledging_grouped_drawsm, 
                  aes(x = actual_mismatch, y = prob_val, fill = fled_group), 
                  .width = 0.95, alpha = 0.3, color = NA) +

  # 3. Posterior Mean Lines
  # We use .width = 0 to draw just the line without another ribbon on top
  stat_lineribbon(data = fledging_grouped_drawsm, 
                  aes(x = actual_mismatch, y = prob_val, color = fled_group), 
                  .width = 0, fill = NA, linewidth = 1) +

  # 4. Scales and Colors
  scale_y_continuous(
    breaks = seq(0, 1, 0.2), 
    minor_breaks = seq(0, 1, 0.1),  
    limits = c(0, 1),
    guide = guide_axis(minor.ticks = TRUE) 
  ) +
  scale_x_continuous(
    breaks = seq(-90, 210, by = 30), 
    minor_breaks = seq(-90, 210, by = 15), 
    limits = c(-90, 210),
    guide = guide_axis(minor.ticks = TRUE) 
  ) +
  
  # Manual colors: Adjust these hex codes to your preference
  scale_fill_manual(values = c("0" = "#999999", "1-2" = "#E69F00", "3" = "#56B4E9"), name = "Number of fledglings") +
  scale_color_manual(values = c("0" = "#999999", "1-2" = "#E69F00", "3" = "#56B4E9"), name = "Number of fledglings") +

  # 5. Labels
  labs(
    subtitle = "\u03B2 = -0.0002, 95% CI = [-0.0005, -0.0000]",
    y = "Probability of fledging outcomes",
    x = "Phenological mismatch (days)"
  ) +
  guides(size = "none",
    # title.position = "left" puts the title next to the keys
    # nrow = 1 forces all keys to be in one line
    fill = guide_legend(title.position = "left", nrow = 1),
    color = guide_legend(title.position = "left", nrow = 1))+
  custom_theme() + theme(
    plot.subtitle = element_text(face = "bold", color = "black"),
    axis.text = element_text(color = "black"),
    axis.title = element_text(color = "black"),
    
    # Move legend to bottom
    legend.position = "bottom",
    # Ensure items are not too spaced out
    legend.spacing.x = unit(0.5, "cm")
)
figure3_female <- (plot_G | plot_H) /
                (plot_fledging_success_age  | plot_mismatch_prob_altf ) +
  plot_annotation(
    tag_levels = 'A'
  )

# Print the final figure
figure3_female
figure3_male <- (plot_Gm | plot_Hm) /
                (plot_fledging_success_agem  | plot_mismatch_prob_altm ) +
  plot_annotation(
    tag_levels = 'A'
  )

# Print the final figure
figure3_male

The assembled figure is then saved as both a PDF and JPG file.

ggsave(filename = here("..","Plots", "panel_fig3_females.pdf"),
       plot = figure3_female, width = 7,  
       height = 9, 
       dpi = 300, device = cairo_pdf)

ggsave(filename = here("..","Plots", "panel_fig3_females.jpg"),
       plot = figure3_female,
       width = 7,  
       height = 9,  
       dpi = 300)




ggsave(filename = here("..","Plots", "panel_fig3_males.pdf"),
       plot = figure3_male,
       dpi = 300, device = cairo_pdf)

ggsave(filename = here("..","Plots", "panel_fig3_males.jpg"),
       plot = figure3_male, 
       width = 7,  
       height = 9,  
       dpi = 300)
R version 4.5.2 (2025-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_Guernsey.utf8  LC_CTYPE=English_Guernsey.utf8   
[3] LC_MONETARY=English_Guernsey.utf8 LC_NUMERIC=C                     
[5] LC_TIME=English_Guernsey.utf8    

time zone: America/Edmonton
tzcode source: internal

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

other attached packages:
 [1] ggtext_0.1.2        extrafont_0.20      cmdstanr_0.9.0.9000
 [4] bayesplot_1.15.0    ggthemes_5.2.0      emmeans_2.0.1      
 [7] rstan_2.32.7        StanHeaders_2.32.10 brms_2.23.0        
[10] Rcpp_1.1.1          tidybayes_3.0.7     lubridate_1.9.4    
[13] forcats_1.0.1       stringr_1.6.0       purrr_1.2.1        
[16] readr_2.1.6         tidyr_1.3.2         tibble_3.3.1       
[19] tidyverse_2.0.0     arm_1.14-4          lme4_1.1-38        
[22] Matrix_1.7-4        patchwork_1.3.2     ggdist_3.3.3       
[25] here_1.0.2          broom.mixed_0.2.9.6 ggplot2_4.0.1      
[28] dplyr_1.1.4         MASS_7.3-65        

loaded via a namespace (and not attached):
 [1] Rdpack_2.6.5          gridExtra_2.3         inline_0.3.21        
 [4] sandwich_3.1-1        rlang_1.1.7           magrittr_2.0.4       
 [7] multcomp_1.4-29       furrr_0.3.1           otel_0.2.0           
[10] matrixStats_1.5.0     compiler_4.5.2        loo_2.9.0            
[13] vctrs_0.7.1           pkgconfig_2.0.3       arrayhelpers_1.1-0   
[16] fastmap_1.2.0         backports_1.5.0       rmarkdown_2.30       
[19] tzdb_0.5.0            ps_1.9.1              nloptr_2.2.1         
[22] xfun_0.56             jsonlite_2.0.0        broom_1.0.12         
[25] parallel_4.5.2        R6_2.6.1              stringi_1.8.7        
[28] RColorBrewer_1.1-3    extrafontdb_1.1       parallelly_1.46.1    
[31] boot_1.3-32           estimability_1.5.1    knitr_1.51           
[34] zoo_1.8-15            pacman_0.5.1          splines_4.5.2        
[37] timechange_0.3.0      tidyselect_1.2.1      rstudioapi_0.18.0    
[40] abind_1.4-8           yaml_2.3.12           codetools_0.2-20     
[43] processx_3.8.6        pkgbuild_1.4.8        listenv_0.10.0       
[46] lattice_0.22-7        withr_3.0.2           bridgesampling_1.2-1 
[49] S7_0.2.1              posterior_1.6.1       coda_0.19-4.1        
[52] evaluate_1.0.5        future_1.69.0         survival_3.8-6       
[55] RcppParallel_5.1.11-1 xml2_1.5.2            pillar_1.11.1        
[58] tensorA_0.36.2.1      checkmate_2.3.3       stats4_4.5.2         
[61] reformulas_0.4.3.1    distributional_0.6.0  generics_0.1.4       
[64] rprojroot_2.1.1       hms_1.1.4             rstantools_2.6.0     
[67] scales_1.4.0          minqa_1.2.8           globals_0.18.0       
[70] xtable_1.8-4          glue_1.8.0            tools_4.5.2          
[73] mvtnorm_1.3-3         grid_4.5.2            Rttf2pt1_1.3.14      
[76] QuickJSR_1.9.0        rbibutils_2.4.1       nlme_3.1-168         
[79] cli_3.6.5             svUnit_1.0.8          Brobdingnag_1.2-9    
[82] gtable_0.3.6          digest_0.6.39         TH.data_1.1-5        
[85] htmlwidgets_1.6.4     farver_2.1.2          htmltools_0.5.9      
[88] lifecycle_1.0.5       gridtext_0.1.5