Persistent Individual Differences and Age Shape Post-Breeding Molt Dynamics in a Partially Migratory Shorebird

Authors

Affiliations

Luke Eberhart-Hertel

Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany

Bashar Jarayseh

Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany

Ailsa Howard

South Bay Banded Dotterel Project, Kaikoura, New Zealand

Ted Howard

South Bay Banded Dotterel Project, Kaikoura, New Zealand

Tony Habraken

Port Waikato Banded Dotterel Project, Port Waikato, New Zealand

Emma Williams

Department of Conservation, Christchurch, New Zealand

Colin O`Donnell

Department of Conservation, Christchurch, New Zealand

Bart Kempenaers

Department of Ornithology, Max Planck Institute for Biological Intelligence, Seewiesen, Germany

Data import and sample summary

We used repeated photographs from an individually marked breeding population at Kaikoura to ask how the expression of the rufous breast band and the timing of its post-breeding molt varies with sex, age, migratory strategy, and breeding phenology. Our full dataset spanned three breeding seasons and was linked to an individual’s breeding records and (if known) their age and migratory status. This section loads the full dataset used throughout the markdown and gives a brief summary of seasonal coverage and sample size before proceeding to the hypotheses and predictions, analysis, and results.

dat_breeding <- 
  readRDS(
    here("data/moult_breeding_KK_data.rds")
  )

We first summarise the number of photographs and the date range covered in each season.

# A tibble: 3 × 4
  season n_obs min_date   max_date  
  <fct>  <int> <date>     <date>    
1 2021     799 2021-08-02 2022-04-22
2 2022     722 2022-07-30 2023-04-14
3 2023     764 2023-08-01 2024-04-23

We then count the number of individuals represented in the full photographic dataset.

[1] 100

Finally, we check how many individuals were sampled in more than one season, because repeated sampling across years is what allows us to estimate persistent individual differences.

# A tibble: 3 × 2
  n_seasons `n()`
      <int> <int>
1         1    34
2         2    28
3         3    38

total, mean, and SD of photos for each individual per season

# A tibble: 1 × 3
  total_photo mean_photos_per_ind_per_season sd_photos_per_ind_per_season
        <int>                          <dbl>                        <dbl>
1        1646                           8.07                         7.48

max and min number of photos per individual per season

# A tibble: 1 × 2
  max_photos_ind_season min_photos_ind_season
                  <int>                 <int>
1                    33                     1

Table of scored photos

Note, in the analysed dataset, for a given individual, each date had only one score, hence the larger sample shown in this table.

Variable definitions

• sourcefile: the directory path to the original photo scored from Ailsa’s catalogue
• rings_comb: the 4-band color combination identifying an individual
• score: the score of the rufous breast band according to the 1-7 scheme shown below
• date: the date of the photo (extracted from the photo’s metadata)
• date_J: the Julian date of the photo, shifted to have an July 1 origin for the autral summer
• sex: sex based on plumage (M or F)
• season: breeding season of the given photo
• migratory_status: Resident (R) or migrant (M) classified from satellite tracking and/or winter resightings
• age: known age in years at the time of the photo (only available for those that were banded locally as chicks and recruited)

Summary table of usable photos for each individual across the three seasons

Scoring scheme

The “score” value is based on a scale from 1 to 7 describing the immaculate-ness and full-ness of the rufous breast band seen in the photos

Sampling distribution of Ailsa’s photos across the seasons

2021-2022 season data

2022-2023 season data

2023-2024 season data

Hypotheses and predictions

We tested the following five hypotheses. The ‘Sex’ and ‘Breeding phenology’ models used the full dataset, while the ‘Age class’, ‘Continuous age’, and ‘Migratory status’ models used subsets of the full dataset for individuals of known age or migratory status (see Table 1 for sample sizes).

Sex

• H1: sexual dichromatism in the rufous band score and sex-specific molt timing
• P1a: males tend to have higher rufous band scores than females, following the general patter of moderate sexual dichromatism in other Charadriidae plovers
• P1b: no difference in post-breeding molt timing due to the socially monogamous and bi-parental care mating system

Breeding phenology

• H2: post-breeding molt is dependent upon breeding phenology
• P2a: individuals that breed earlier molt earlier due a generally more advanced reproductive schedule
• P2b: individuals that are still engaged in parental care late in the season molt later
• P2c: individuals that initiate late nests molt later

Age class

• H3: age-class structures the molt schedule and rufous band expression
• P3a: first-year birds express lower rufous band scores in comparison to after-first-year birds due to more complete pigmentation during feather growth in older birds
• P3b: first-year birds molt earlier than after-first-year birds because the flight feathers of first-years are older (and more worn) than after-first-years, which triggers an earlier complete post-breeding molt program

Continuous age

• H4: beyond the first-year transition continuous age structures the molt schedule and rufous band expression
• P4a: no age-dependent variation in rufous band expression (either within or between individuals)
• P4b: no age-dependent variation in post-breeding molt timing (either within or between individuals)

Migratory strategy

• H5: migratory strategy structures the molt schedule and rufous band expression
• P5a: no strategy-dependent variation in rufous band expression
• P5b: migrants molt earlier than residents due to their more constrained fall migration schedule

Non-linear mixed effects model

We scored the rufous breast band on a seven-point scale, so the response is bounded and declines through the season as birds transition from their breeding (“alternate”) plumage into their non-breeding plumage (“basic”) via the complete post-breeding molt.

We fit the same nonlinear trajectory for all five models and let the predictors explain variation in the starting score and the timing of the decline.

We modeled the decline in rufous band score with a nonlinear mixed-effects curve. Because we centered and scaled Julian date within season before fitting, we write the seasonal date term as x. For observation j on bird i,

y_{ij} \sim \mathrm{Normal}(\mu_{ij}, \sigma).

\mu_{ij} = 1 + (\beta_{ij} - 1) \left( 1 - \frac{1}{1 + \exp[-k(x_{ij} - t_{0,ij})]} \right).

The parameter \beta is the rufous band score prior to the post-breeding molt, t_0 is the midpoint of the decline to the non-breeding plumage, and k controls how steep/fast the molt progresses. Early in the season (x \ll t_0), expected scores are close to \beta; late in the season (x \gg t_0), expected scores approach 1 (i.e. full non-breeding plumage). This function attempts to match the biology of post-breeding moult, in which birds move from full breeding plumage toward full non-breeding plumage over time.

In the sex-specific model, we allowed both \beta and t_0 to depend on sex and to vary among birds and bird-by-season combinations, while keeping k common across birds (i.e., our photographic sampling was too coarse to robustly assess sources of variation in molt speed):

\begin{aligned} \beta_{ij} &= \alpha_{\beta} + \mathbf{X}_{ij}\boldsymbol{\gamma}_{\beta} + b_{\beta,i} + b_{\beta,is}, \\ t_{0,ij} &= \alpha_{t_0} + \mathbf{X}_{ij}\boldsymbol{\gamma}_{t_0} + b_{t_0,i} + b_{t_0,is}, \\ k &= \alpha_k . \end{aligned}

Here, \mathbf{X} denotes the relevant fixed effects for a given model. In the sex-only model, \mathbf{X} contains sex; in the other models, we added breeding phenology, age class, continuous age terms (i.e., within- vs. between-individual), or migratory status. We kept the same nonlinear form throughout so that differences among models reflect changes in predictors of \beta and t_0, not changes in the underlying molt curve.

We used an informative prior on \beta (i.e., normal(5.5, 1); centering the prior near the upper half of the 1-7 scoring scale and keeping it positive), weakly informative priors on t_0 and k, and exponential priors on the group-level standard deviations. This combination stabilized estimation while still allowing substantial variation among birds and among bird-by-season combinations when the data supported it.

The Figure 1 schematic below shows how \beta, t_0, and k shape the curve that we fit in every model. We use representative values here so that the roles of the three parameters are easy to interpret.

Figure 1 - schematic of the non-linear model

# Draw a schematic version of the nonlinear curve used in all models.
# The values of beta, t0, and k are chosen for illustration only.

# Set the schematic parameters.
baseline <- 6
k <- 3.31
t0 <- 1.45

# Define the expected trajectory from the nonlinear model.
moult_mu <- function(x, baseline, k, t0) {
  1 + (baseline - 1) * (1 - 1 / (1 + exp(-k * (x - t0))))
}

# Build a clean x-range around the inflection point.
x <- seq(t0 - 2.5, t0 + 2.5, length.out = 600)

df <- tibble(
  x = x,
  mu = moult_mu(x, baseline = baseline, k = k, t0 = t0)
)

# At x = t0 the curve is halfway between the upper and lower asymptotes.
y_t0 <- moult_mu(t0, baseline = baseline, k = k, t0 = t0)

# The slope at t0 determines the tangent line used to illustrate k.
slope_t0 <- -(baseline - 1) * k / 4
angle_k <- atan(slope_t0) * 180 / pi

y_k <- 7
x_k <- t0 + (y_k - y_t0) / slope_t0

fig1 <- 
  ggplot(df, aes(x = x, y = mu)) +
  geom_hline(yintercept = baseline, linetype = "dashed", linewidth = 1, color = "grey70") +
  ggtext::geom_richtext(
    data = data.frame(x = min(df$x), y = baseline),
    aes(x = x, y = y - 0.1),
    label = "baseline (<span style='font-family:NotoSans; font-style:italic;'>β</span>)",
    family = "Lato", size = 3,
    hjust = 0, vjust = -0.6,
    fill = NA, label.color = NA
  ) +
  geom_vline(xintercept = t0, linetype = "dashed", linewidth = 1, color = "grey70") +
  geom_point(aes(x = t0, y = y_t0), size = 2) +
  annotate(
    "text", size = 3,
    family = "Lato",
    x = t0, y = min(df$mu),
    label = as.expression(bquote("molt midpoint" ~ "(" * italic(t)[0] * ")" == .(round(t0, 2)))),
    parse = FALSE,
    hjust = 1, vjust = 1.5, angle = 270
  ) +
  {
    slope_t0 <- -(baseline - 1) * k / 4
    geom_abline(
      intercept = y_t0 - slope_t0 * t0,
      slope = slope_t0,
      linetype = "dashed", linewidth = 1, color = "grey70"
    )
  } +
  annotate(
    "text", size = 3,
    family = "Lato",
    x = x_k, y = y_k,
    label = as.expression(bquote("molt rate (" * italic(k) * ") = " * .(round(k, 2)))),
    hjust = 0.5,
    vjust = -0.5,
    angle = angle_k + 5
  ) +
  geom_line(linewidth = 2, color = "grey30") +
  geom_point(aes(x = t0, y = y_t0), size = 5, shape = 21, fill = "white", color = "black") +
  geom_point(aes(x = -0.5, y = 6), size = 5, shape = 23, fill = "white", color = "black") +
  labs(
    x = "Scaled date (mean-centered Julian date)",
    y = "Rufous band score"
  ) +
  coord_cartesian(ylim = c(1, 8)) +
  theme_minimal() +
  scale_y_continuous(limits = c(1, 8), breaks = c(1:7)) +
  theme(
    text = element_text(family = "Lato"),
    panel.grid.minor = element_blank(),
    axis.text = element_text(color = "grey30")
  )

fig1

Shared priors

We use the same prior set in all five non-linear moult models so that differences among models reflect changes in predictors rather than changes in prior assumptions. We define that shared prior set once here and reuse it in each brm() call below.

# Shared priors for all non-linear moult models.
moult_priors <- c(
  # The baseline score sits on the 1-7 scoring scale, so we center its
  # prior near the upper half of that range while keeping it positive.
  prior(normal(5.5, 1), nlpar = "baseline", lb = 0),

  # The inflection point is estimated on the season-standardized date scale.
  prior(normal(0, 0.5), nlpar = "t0"),

  # Positive values of k enforce the expected one-way seasonal decline.
  prior(normal(1, 0.5), nlpar = "k", lb = 0),

  # Group-level SDs are regularized, but can still be large if the data
  # support substantial heterogeneity among birds or bird-seasons.
  prior(exponential(2), nlpar = "baseline", class = "sd", group = "rings_comb"),
  prior(exponential(2), nlpar = "baseline", class = "sd", group = "rings_comb:season"),
  prior(exponential(2), nlpar = "t0", class = "sd", group = "rings_comb"),
  prior(exponential(2), nlpar = "t0", class = "sd", group = "rings_comb:season")
)

Analysis: Sex-specific model

We first fit a model to the full dataset and asked whether males and females differed in \beta and t_0. We used this as the reference model because sex was known for all birds and it gives the cleanest estimate of the post-breeding molt trajectory and baseline expression of the rufous band.

dataset

# wrangle data
dat_sex_mod <- dat_breeding %>%
  # subset to columns of interest
  select(score, date_J, rings_comb, sex, season) %>%
  # reduce to distinct data
  distinct() %>% 
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>% 
  ungroup() %>%
  # drop NAs
  drop_na(score, date_J, sex, rings_comb, season) %>% 
  # convert to dataframe
  as.data.frame()

model formula specification

mod_sex <- bf(
  score ~ 1 + (baseline - 1) * (1 - 1 / (1 + exp(-k * (date_J_scaled - t0)))),
  baseline ~ sex + (1 | i | rings_comb) + (1 | rings_comb:season),
  t0 ~ sex + (1 | i | rings_comb) + (1 | rings_comb:season),
  k ~ 1,
  nl = TRUE
)

fit model

fit_mod_sex <- brm(
  formula = mod_sex,
  data = dat_sex_mod,
  family = gaussian(),
  prior = moult_priors,
  control = list(adapt_delta = 0.999, max_treedepth = 15),
  chains = 4,
  cores = 4,
  iter = 10000,
  warmup = 5000,
  backend = "cmdstanr",
  seed = 123
)

Analysis: Breeding Phenology model

We then retained sex in the model and added first breeding date, last breeding date, and last nest initiation date as predictors of molt timing (t_0). Any carry-over effect of breeding should appear first in the timing of the decline, so here we focus on whether later breeding activity shifts the molt midpoint.

datasets

Here we wrangle the full breeding phenology dataset (dat_breeding_mod), which is the version used in the main analysis and reported in the paper. Two more conservative alternatives (dat_breeding_mod2 and dat_breeding_mod3) are treated as auxiliary sensitivity analyses and are presented in the Auxiliary analyses section below.

# wrangle data
dat_breeding_mod <- dat_breeding %>%
  # subset to columns of interest
  select(score, date_J, rings_comb, sex, season, 
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum) %>% 
  # reduce to distinct data
  distinct() %>% 
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>% 
  ungroup() %>%
  # drop NAs
  drop_na(score, date_J, sex, rings_comb, season, 
          last_breeding_date_J_snum, 
          first_breeding_date_J_snum, 
          last_nest_first_obs_J_snum) %>% 
  # convert to dataframe
  as.data.frame()

model formula specification

mod_sex_all_3_breed <- bf(
   score ~ 1 + (baseline - 1) * (1 - 1 / (1 + exp(-k * (date_J_scaled - t0)))),
   baseline ~ sex + (1 | i | rings_comb) + (1 | rings_comb:season),
   t0 ~ sex + first_breeding_date_J_snum + last_nest_first_obs_J_snum + last_breeding_date_J_snum + (1 | i | rings_comb) + (1 | rings_comb:season),
   k ~ 1,
   nl = TRUE
)

fit model

We fit the main breeding phenology model to the full dataset here. The two more conservative sensitivity versions are presented in the Auxiliary analyses section below.

fit_mod_sex_all_3_breed <- brm(
  formula = mod_sex_all_3_breed,
  data = dat_breeding_mod,
  family = gaussian(),
  prior = moult_priors,
  control = list(adapt_delta = 0.999, max_treedepth = 15),
  chains = 4,
  cores = 4,
  iter = 10000,
  warmup = 5000,
  backend = "cmdstanr",
  seed = 123
)

Analysis: Age Class model

We next incorporated age class (first-year versus after-first-year) to the sex-specific model to test on a subser of known aged birds whether individuals in the first molt cycle differ from after-first-year birds in baseline plumage and molt timing.

dataset

# wrangle data
dat_age_class_mod <- dat_breeding %>%
  # subset to columns of interest
  select(score, date_J, rings_comb, sex, season, age_class) %>%
  # reduce to distinct data
  distinct() %>%
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>% 
  ungroup() %>%
  # drop NAs
  drop_na(score, date_J, sex, rings_comb, season, age_class) %>% 
  # convert to dataframe
  as.data.frame()

model formula specification

mod_sex_age_class <- bf(
  score ~ 1 + (baseline - 1) * (1 - 1 / (1 + exp(-k * (date_J_scaled - t0)))),
  baseline ~ sex + age_class + (1 | i | rings_comb) + (1 | rings_comb:season),
  t0 ~ sex + age_class + (1 | i | rings_comb) + (1 | rings_comb:season),
  k ~ 1,
  nl = TRUE
)

fit model

fit_mod_sex_age_class <- brm(
  formula = mod_sex_age_class,
  data = dat_age_class_mod,
  family = gaussian(),
  prior = moult_priors,
  control = list(adapt_delta = 0.999, max_treedepth = 15),
  chains = 4,
  cores = 4,
  iter = 10000,
  warmup = 5000,
  backend = "cmdstanr",
  save_pars = save_pars(all = TRUE),
  seed = 123
)

Analysis: Continuous Age model

Following Van de Pol & Verhulst (2006), we decomposed age (x_{ist}) into each bird’s mean age across the three-year study period (\bar{x}_i; the between-individual effect, implemented as avg_age) and the deviation from that mean for a given observation (x_{ist} - \bar{x}_i; the within-individual effect, implemented as age_dev). This allowed us to test whether any age effect extended beyond the first-year transition from the complex molt cycle rather than simply reflecting a first-year contrast.

datasets

Here we use the full continuous age dataset (dat_age_mod), including individuals with age 1. A second version that excludes age-1 observations (dat_age2_mod) is treated as an auxiliary sensitivity analysis and is presented in the Auxiliary analyses section below.

# calculate average age in the dataset of individuals of known age
avg_ages <- 
  dat_breeding %>%
  distinct(rings_comb, age) %>% 
  filter(!is.na(age)) %>% 
  group_by(rings_comb) %>% 
  summarise(avg_age = mean(age, na.rm = TRUE))

# wrangle data
dat_age_mod <- dat_breeding %>%
  # subset to columns of interest
  select(score, date_J, rings_comb, sex, season, age) %>%
  # reduce to distinct data
  distinct() %>% 
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>%
  # center age
  filter(!is.na(age)) %>% 
  left_join(., avg_ages, by = "rings_comb") %>% 
  # age deviance
  ungroup() %>%
  mutate(age_dev = age - avg_age) %>% 
  # drop NAs
  drop_na(score, date_J, sex, rings_comb, 
          season, avg_age, age_dev) %>% 
  # convert to dataframe
  as.data.frame()

model formula specification

mod_sex_age_avg <- bf(
  score ~ 1 + (baseline - 1) * (1 - 1 / (1 + exp(-k * (date_J_scaled - t0)))),
  baseline ~ sex + age_dev + avg_age + (1 | i | rings_comb) + (1 | rings_comb:season),
  t0 ~ sex + age_dev + avg_age + (1 | i | rings_comb) + (1 | rings_comb:season),
  k ~ 1,
  nl = TRUE
)

fit model

We fit the main continuous-age model to the full known-age dataset here. The sensitivity version excluding age-1 observations is presented in the Auxiliary analyses section below.

fit_mod_sex_age_avg <- brm(
  formula = mod_sex_age_avg,
  data = dat_age_mod,
  family = gaussian(),
  prior = moult_priors,
  control = list(adapt_delta = 0.999, max_treedepth = 15),
  chains = 4,
  cores = 4,
  iter = 10000,
  warmup = 5000,
  backend = "cmdstanr",
  save_pars = save_pars(all = TRUE),
  seed = 123
)

Analysis: Migratory Strategy model

We finally tested whether migratory strategy explained variation in \beta and t_0 after accounting for sex. If migration imposes tighter late-season constraints, we would expect migrants to move through the post-breeding transition earlier than residents.

datasets

# wrangle data
dat_migratory_status_mod <- 
  dat_breeding %>%
  # subset to columns of interest
  select(score, date_J, rings_comb, sex, season, migratory_status) %>%
  # reduce to distinct data
  distinct() %>% 
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>% 
  ungroup() %>%
  # drop NAs
  drop_na(score, date_J, sex, rings_comb, season, migratory_status) %>% 
  # convert to dataframe
  as.data.frame()

model formula specification

mod_sex_migratory_status <- bf(
  score ~ 1 + (baseline - 1) * (1 - 1 / (1 + exp(-k * (date_J_scaled - t0)))),
  baseline ~ sex + migratory_status + (1 | i | rings_comb) + (1 | rings_comb:season),
  t0 ~ sex + migratory_status + (1 | i | rings_comb) + (1 | rings_comb:season),
  k ~ 1,
  nl = TRUE
)

fit model

fit_mod_sex_migratory_status <- brm(
  formula = mod_sex_migratory_status,
  data = dat_migratory_status_mod,
  family = gaussian(),
  prior = moult_priors,
  control = list(adapt_delta = 0.999, max_treedepth = 15),
  chains = 4,
  cores = 4,
  iter = 10000,
  warmup = 5000,
  backend = "cmdstanr",
  seed = 123
)

Results

Using 1,647 scored photographs from 100 birds across three seasons, we found that 1) baseline rufous band score and molt timing both showed persistent differences among birds, 2) age explained variation both in terms of the first-year vs. after-first year contrast, and in the continuous age effect, and 3) there was little evidence that breeding phenology or migratory strategy shifted the molt schedule. The results below start from the saved model objects in out/, so this section can be rerun without refitting the time-consuming Bayesian models above. Helper chunks are hidden in the HTML to keep the shared document readable, while the visible chunks show only the short code needed to reproduce the reported tables.

Diagnostics

Report compact convergence summaries for each fitted model rather than printing full brms output

Model	Max R-hat	Parameters with R-hat > 1.01	Min bulk ESS	Min tail ESS
Sex-only moult model	1.001	0	1015	2409
Sex + breeding phenology model	1.002	0	1808	3466
Sex + age-class model	1.001	0	2783	2319
Sex + continuous-age model	1.003	0	2931	4567
Sex + migratory-status model	1.002	0	2577	1649

Compare the two age parameterisations using the saved 10-fold cross-validation objects.

age_model_comparison <- loo_compare(kfc_age_class, kfc_age_avg) %>%
  as.data.frame() %>%
  rownames_to_column("model") %>%
  mutate(
    model = recode(
      model,
      kfc_age_class = "Age-class model",
      kfc_age_avg = "Continuous-age model"
    ),
    across(where(is.numeric), \(x) round(x, 2))
  )

age_model_comparison %>%
  gt() %>%
  cols_label(
    model = "Model",
    elpd_diff = "elpd diff",
    se_diff = "SE diff"
  )

Model	elpd diff	SE diff	elpd_kfold	se_elpd_kfold	p_kfold	se_p_kfold
sex_age_class_mod$fitted_model	0.0	0.00	-176.58	17.10	95.40	10.96
sex_age_avg_mod$fitted_model	-1.9	3.35	-178.48	16.57	95.74	10.36

Results: Sex-specific model

We use the sex-only model as the main reference because it uses the full dataset. It shows clear sex differences in baseline ornament expression, little evidence for a sex difference in molt midpoint, and substantial individual consistency in both expression and timing.

Parameter	Posterior median [95% CrI]
Baseline alternate plumage score (β), females (reference)	4.69 [4.51, 4.87]
Sex effect on baseline alternate plumage score (β): males − females	1.15 [0.89, 1.42]
Pre-basic moult midpoint (t0), females (reference)	Feb-04 [Jan-29, Feb-10]
Sex effect on pre-basic moult midpoint (t0): males − females	-1.68 [-9.8, 6.25]
Pre-basic moult rate at inflection (days·score⁻¹), females (reference)	18.45 [16.99, 20.04]
Sex effect on pre-basic moult rate at inflection (days·score⁻¹): males − females	-4.37 [-5.52, -3.33]
Between-year repeatability of alternate plumage	0.76 [0.7, 0.81]
Within-year repeatability of alternate plumage	0.51 [0.42, 0.61]
Between-year repeatability of pre-basic moult mid-point date	0.5 [0.41, 0.58]
Within-year repeatability of pre-basic moult mid-point date	0.47 [0.36, 0.56]
Marginal R²	0.68 [0.65, 0.71]
Conditional R²	0.92 [0.91, 0.92]
Between-individual variance in baseline plumage score	0.29 [0.18, 0.45]
Between-individual variance in pre-basic moult midpoint date	0.01 [0, 0.05]
Correlation between baseline alternate plumage score and pre-basic moult midpoint date	-0.47 [-0.97, 0.7]
Residual variance	0.14 [0.13, 0.15]

Results: Breeding Phenology model

The effect of first breeding date, last breeding date, and last nest initiation date were all small on the day scale. Specifically, later breeding activity did not translate into a clear delay in the post-breeding molt midpoint.

Parameter	Posterior median [95% CrI]
Δ pre-basic moult midpoint (t0; days) per +1 SD (within-season) delay in first breeding date	-0.03 [-0.11, 0.05]
Δ pre-basic moult midpoint (t0; days) per +1 SD (within-season) delay in last breeding date	0.03 [-0.08, 0.14]
Δ pre-basic moult midpoint (t0; days) per +1 SD (within-season) delay in last nest date	0.09 [-0.03, 0.21]

Results: Age Class model

The age-class model showed a clear contrast between first-year and after-first-year birds. Older birds tended to have slightly higher baseline scores and later molt midpoints.

Parameter	Posterior median [95% CrI]
Age-class effect on baseline alternate plumage score (β): AFY − FY	0.19 [0.01, 0.66]
Age-class effect on pre-basic moult mid-point date (t0): AFY − FY	28.01 [7.33, 48.76]

Results: Continuous Age model

Treating age as a continuous trait gave the same broad result as the age class model, but demonstrated that the age effect extended beyond the transition from the first-year cycle. Older birds tended to show higher baseline scores and later molt midpoints, with evidence for both between-individual and within-individual age effects.

Parameter	Posterior median [95% CrI]
Between-individual age effect on baseline alternate plumage score (β)	0.04 [0, 0.16]
Within-individual age effect on baseline alternate plumage score (β)	0.18 [0.04, 0.35]
Between-individual age effect on pre-basic moult midpoint (t0)	6.73 [1.77, 12.03]
Within-individual age effect on pre-basic moult midpoint (t0)	7.64 [-2.77, 17.94]

Results: Migratory Strategy model

Residents and migrants showed similar molt timing. Residents tended to have slightly higher baseline scores, but the overlap was broad and the timing contrast was weak.

Parameter	Posterior median [95% CrI]
Migratory-strategy effect on baseline alternate plumage score (β): Resident − Migrant	0.3 [0.03, 0.7]
Migratory-strategy effect on pre-basic moult midpoint (t0): Resident − Migrant	6.93 [-6.47, 21.11]

Table 1

Consolidated table of results

	Posterior median	[95% CrI]	Nind	Model
Fixed effects on baseline breeding plumage score (β)
Intercept, females (reference; score)	4.69	[4.51, 4.87]	100	Sex-specific
Sex effect (males – females; Δ score)	1.15	[0.89, 1.42]	100	Sex-specific
Age-class effect (AFY – FY; Δ score)	0.19	[0.01, 0.66]	30	Age-class
Between-individual age effect (Δ score)	0.04	[0.0, 0.16]	30	Within-individual centering
Within-individual age effect (Δ score)	0.18	[0.04, 0.35]	30	Within-individual centering
Migratory-strategy effect (Resident – Migrant; Δ score)	0.3	[0.03, 0.7]	41	Migratory status
Fixed effects on post-breeding molt midpoint (t0): sex, age, and migration strategy
Intercept, females (reference; date)	4-Feb	[29 Jan, 10 Feb]	100	Sex-specific
Sex effect (males – females; Δ days)	-1.68	[-9.8, 6.25]	100	Sex-specific
Age-class effect (AFY – FY; Δ days)	28.01	[7.33, 48.76]	30	Age-class
Between-individual age effect (Δ days)	6.73	[1.77, 12.03]	30	Within-individual centering
Within-individual age effect (Δ days)	7.64	[-2.77, 17.94]	30	Within-individual centering
Migratory-strategy effect (Resident – Migrant; Δ days)	6.93	[-6.47, 21.11]	41	Migratory status
Fixed effects on post-breeding molt midpoint (t0): effect of +1 SD (within-season) delay in breeding phenology trait
First date observed breeding (Δ days)	-0.03	[-0.11, 0.05]	81	Breeding phenology
Last date observed breeding (Δ days)	0.03	[-0.08, 0.14]	81	Breeding phenology
Last nest initiation date (Δ days)	0.09	[-0.03, 0.21]	81	Breeding phenology
Derived fixed effects: post-breeding molt rate at t0
Intercept, females (reference; days·score⁻¹)	18.45	[16.99, 20.04]	100	Sex-specific
Sex effect (males – females; days·score⁻¹)	-4.37	[-5.52, -3.33]	100	Sex-specific
Repeatability of molt traits (R)
Between-year repeatability of β	0.76	[0.7, 0.81]	100	Sex-specific
Within-year repeatability of β	0.51	[0.42, 0.61]	100	Sex-specific
Between-year repeatability of t0	0.5	[0.41, 0.58]	100	Sex-specific
Within-year repeatability of t0	0.47	[0.36, 0.56]	100	Sex-specific
Model fit
Marginal R²	0.68	[0.65, 0.71]	100	Sex-specific
Conditional R²	0.92	[0.91, 0.92]	100	Sex-specific
Variance components and residual structure
Between-individual variance in β	0.29	[0.18, 0.45]	100	Sex-specific
Between-individual variance in t0	0.01	[0.0, 0.05]	100	Sex-specific
Correlation between β and t0	-0.47	[-0.97, 0.7]	100	Sex-specific
Residual variance	0.14	[0.13, 0.15]	100	Sex-specific
Note: Δ denotes an additive difference relative to the reference category (for categorical predictors) or the expected change in the response per unit change in the predictor (for continuous covariates).

Figures

Figure 2 - raw data overview

breeding_data_all <- readRDS(here("data/breeding_KK_data.rds"))

#### breeding histograms
# 1) Get per-season centering/scaling constants from the ORIGINAL (unscaled) shifted-Julian variables
stats <- breeding_data_all %>%
  group_by(season) %>%
  summarise(
    mu_last  = mean(last_breeding_date_J, na.rm = TRUE),
    sd_last  = sd(last_breeding_date_J,  na.rm = TRUE),
    mu_first = mean(first_breeding_date_J, na.rm = TRUE),
    sd_first = sd(first_breeding_date_J,  na.rm = TRUE),
    mu_nest  = mean(last_nest_first_obs_J, na.rm = TRUE),
    sd_nest  = sd(last_nest_first_obs_J,  na.rm = TRUE),
    .groups = "drop"
  )

# 2) Join constants, invert transform, then convert shifted-Julian (1 = Jul 1) to a common calendar date
ref_start <- as.Date("2001-07-01")  # any non-leap year works well

out_dates <- dat_breeding %>%
  distinct(rings_comb, season,
           last_breeding_date_J_snum,
           first_breeding_date_J_snum,
           last_nest_first_obs_J_snum) %>%
  drop_na() %>%
  left_join(stats, by = "season") %>%
  mutate(
    # If _snum is CENTER-ONLY (scale = 0): x = centered + mean
    last_J  = last_breeding_date_J_snum  + mu_last,
    first_J = first_breeding_date_J_snum + mu_first,
    nest_J  = last_nest_first_obs_J_snum + mu_nest,

    # Convert shifted-Julian to a common calendar date (same year for all seasons)
    last_date_common  = ref_start + round(last_J  - 1),
    first_date_common = ref_start + round(first_J - 1),
    nest_date_common  = ref_start + round(nest_J  - 1),

    # Optional: month-day labels (often what you want for plotting)
    last_md  = format(last_date_common,  "%b-%d"),
    first_md = format(first_date_common, "%b-%d"),
    nest_md  = format(nest_date_common,  "%b-%d")
  ) %>%
  select(rings_comb, season,
         last_J, first_J, nest_J,
         last_date_common, first_date_common, nest_date_common,
         last_md, first_md, nest_md)

out_dates_long <- out_dates %>%
  select(rings_comb, season, last_date_common, first_date_common, nest_date_common) %>%
  pivot_longer(
    cols = c(last_date_common, first_date_common, nest_date_common),
    names_to = "phenology_trait",
    values_to = "date_common"
  ) %>%
  mutate(
    phenology_trait = sub("_date_common$", "_date", phenology_trait),
    phenology_trait = sub("_common$", "", phenology_trait)
  )

# over lapping histograms
nesting_hist <- 
  ggplot(
      out_dates_long %>%
        filter(phenology_trait %in% c("first_date", "last_date")),
      aes(x = date_common, fill = phenology_trait)
    ) +
  geom_histogram(binwidth = 7, position = "identity", alpha = 0.55, colour = NA) +
  theme_minimal() +
  theme(
    text = element_text(family = "Lato"),
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_blank(),
    panel.grid = element_blank(),
    plot.background = element_rect(fill = "white", color = NA),
    legend.position = "none",
    plot.margin = margin(0.2, 0.2, -15, 0.2)
  ) +
  scale_y_continuous(position = "left", breaks = c(0, 10, 20), limits = c(0, 25)) +
  scale_x_date(
    date_labels = "%B", expand = c(0.01, 0.01),
    date_breaks = "2 month",
    limits = c(as.Date("2001-07-01"), as.Date("2002-05-01"))
  ) +
  scale_fill_manual(
    values = c(first_date = "grey70", last_date = "grey30")
  ) +
  annotate(
    "text",
    x = as.Date("2001-07-25"),
    y = 15,
    label = "First\nbreeding\ndates",
    fontface = "italic",
    family = "Lato",
    size = 3,
    hjust = 0.5,
    vjust = 0,
    colour = "grey30"
  ) +
  annotate(
    "text",
    x = as.Date("2002-01-25"),
    y = 15,
    label = "Last\nbreeding\ndates",
    fontface = "italic",
    family = "Lato",
    size = 3,
    hjust = 0.5,
    vjust = 0,
    colour = "grey30"
  ) +
  annotate(
    "text",
    x = as.Date("2002-03-01"),
    y = 20,
    label = "20",
    fontface = "italic",
    family = "Lato",
    size = 3,
    hjust = 0.5,
    vjust = 0,
    colour = "grey30"
  ) +
  annotate(
    "text",
    x = as.Date("2002-03-01"),
    y = 10,
    label = "10",
    fontface = "italic",
    family = "Lato",
    size = 3,
    hjust = 0.5,
    vjust = 0,
    colour = "grey30"
  ) +
  annotate(
    "text",
    x = as.Date("2002-03-01"),
    y = 0,
    label = "0",
    fontface = "italic",
    family = "Lato",
    size = 3,
    hjust = 0.5,
    vjust = 0,
    colour = "grey30"
  )

#### raw data and example individual trajectories
nests_data <- readRDS(here("data/nests_KK_data.rds"))

dates_for_plot <-
  data.frame(date_ = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             date_J = c(1:365))

# 1) Build a clean prediction grid for the two individuals across dates and seasons
newdata_ind <- tidyr::expand_grid(
  date_J_scaled = seq(
    from = min(sex_mod$data$date_J_scaled),
    to   = max(sex_mod$data$date_J_scaled),
    length.out = max(sex_mod$data$date_J) - min(sex_mod$data$date_J) + 1
  ),
  rings_comb = c("RRGO", "RGGL"),
  season     = levels(sex_mod$data$season)   # e.g., c("2021","2022","2023")
) %>%
  # assign sex per individual from your observed data (safer than expanding sex)
  left_join(
    sex_mod$data %>%
      distinct(rings_comb, sex) %>%
      filter(rings_comb %in% c("RRGO", "RGGL")),
    by = "rings_comb"
  ) %>%
  # if you want Julian date (unscaled) + calendar date for plotting, join your helper tables
  left_join(
    data.frame(
      date_J = seq(min(sex_mod$data$date_J), max(sex_mod$data$date_J), length.out =
                     max(sex_mod$data$date_J) - min(sex_mod$data$date_J) + 1),
      date_J_scaled = seq(min(sex_mod$data$date_J_scaled), max(sex_mod$data$date_J_scaled), length.out =
                            max(sex_mod$data$date_J) - min(sex_mod$data$date_J) + 1)
    ),
    by = "date_J_scaled"
  ) %>%
  left_join(dates_for_plot, by = "date_J")

# 2) Extract expected-value draws (includes group-level effects by default)
epred_draws_ind <- add_epred_draws(
  sex_mod$fitted_model,
  newdata = newdata_ind,
  ndraws = 50
)

photo_data <-
  dat_breeding %>%
  filter(filename %in% c("IMG_0053.JPG", "IMG_0117.JPG", "IMG_0146.JPG", "IMG_0722.JPG", "IMG_2065.JPG", "IMG_0411.JPG")) %>%
  select(rings_comb, sex, migratory_status, age_at_banding, age, score, date_J, date) %>%
  left_join(., dates_for_plot, by = "date_J") %>%
  mutate(score = as.numeric(score))

found_dates_all <- 
  nests_data %>% 
  filter(!is.na(found_date)) %>%
  mutate(found_date_J = as.integer(
    found_date - as.Date(paste0(year(add_with_rollback(found_date, -months(6))), "-07-01"))
  ) + 1)

found_dates_for_plot <-
  data.frame(date_ = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))),
             found_date_J = c(1:365))

found_dates_all_ <-
  found_dates_all %>%
  left_join(., found_dates_for_plot, by = "found_date_J") %>%
  # select(-rings_comb) %>%
  distinct()

sex_score_plot <-
  ggplot() +
  # Add the faint box first
  geom_rect(aes(xmin = min(found_dates_all_$date_), xmax = max(found_dates_all_$date_),
                ymin = -Inf, ymax = Inf),
            fill = "grey70", alpha = 0.2) +  # adjust fill and alpha as needed +
  geom_jitter(data = dat_breeding %>% 
                left_join(., dates_for_plot, by = "date_J") %>% 
                mutate(score = as.numeric(score)) %>% 
                filter(!is.na(sex)),
              aes(x = date_, y = score), fill = "black",
              shape = 16, size = 3, alpha = 0.2,
              height = 0.15, width = 0) +
  geom_line(data = epred_draws_ind %>% 
                filter((rings_comb == "RGGL" & season == 2023) | 
                         (rings_comb == "RRGO" & season == 2022)), 
            aes(x = date_, y = .epred, color = sex,
                group = interaction(rings_comb, .draw, season)), alpha = 0.1) +
  geom_point(data = dat_breeding %>% 
                left_join(., dates_for_plot, by = "date_J") %>% 
                mutate(score = as.numeric(score)) %>% 
                filter(!is.na(sex)) %>% 
                filter((rings_comb == "RGGL" & season == 2023) | 
                         (rings_comb == "RRGO" & season == 2022)), 
             aes(x = date_, y = score, fill = sex),
             shape = 21, color = "white", size = 2, stroke = 0.5) +
  geom_point(data = photo_data %>% filter(score %in% c(5, 6)),
             aes(x = date_, y = score, fill = sex),
             shape = 24, color = "white", size = 4, stroke = 1) +
  geom_point(data = photo_data %>% filter(score %in% c(4)),
             aes(x = date_, y = score, fill = sex),
             shape = 22, color = "white", size = 4, stroke = 1) +
  geom_point(data = photo_data %>% filter(score %in% c(2)),
             aes(x = date_, y = score, fill = sex),
             shape = 25, color = "white", size = 4, stroke = 1) +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(family = "Lato"),
        axis.ticks = element_blank(),
        axis.title.y = element_text(size = 11, colour = "black"),
        axis.text.y = element_text(size = 9, colour = "grey30"),
        # plot.margin = unit(c(0,5.5,5.5,0), "cm"),
        panel.spacing = unit(1, "cm"),
        # panel.border = element_rect(size = 0.25, colour = "grey70", fill = 'transparent'),
        panel.background = element_rect(fill = 'transparent', color = NA),
        plot.background = element_rect(fill = 'transparent', color = NA),
        strip.background = element_rect(fill = 'grey80', color = NA),
        strip.text = element_text(size = 13, colour = "grey40"),
        panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 9,
                                   angle = 45,
                                   hjust = 1,
                                   vjust = 1), colour = "grey30") +
  scale_x_date(date_labels = "%B", expand = c(0.01, 0.01),
               date_breaks = "2 month",
               limits = c(as.Date("2021-07-01"), as.Date("2022-05-01"))) +
  scale_y_continuous(limits = c(0.5, 7.5), breaks = c(1:7)) +
  ylab("Rufous band score\n\n\n") +
  annotate("text",
         x = min(found_dates_all_$date_) + (max(found_dates_all_$date_) - min(found_dates_all_$date_)) / 2,
         y = 0.5,
         label = "Nest initiation period",
         fontface = "italic",
         family = "Lato",
         size = 3,  # adjust size as needed
         hjust = 0.5,
         vjust = 0, colour = "grey40") +
  scale_fill_manual(values = c("#F28E2B", "#4E79A7")) +
  scale_color_manual(values = c("#F28E2B", "#4E79A7")) +
  # Add the double-headed arrow
  geom_segment(aes(x = min(found_dates_all_$date_), xend = max(found_dates_all_$date_), y = 1, yend = 1),
               arrow = arrow(ends = "both", type = "closed", length = unit(0.2, "cm")),
               linewidth = 0.7, color = "grey70")

fig2 <-
  (nesting_hist / sex_score_plot) +
   plot_layout(heights = c(1, 3)) & theme(plot.margin = margin(5.5, 5.5, 0, 0))

fig2

Figure 3 - posterior ridge plots for age model

#### t0: age_avg 
fit <- sex_age_avg_mod$fitted_model
dat <- sex_age_avg_mod$data

# --- counts of individuals (EDIT the age column name if needed) ---
# Assumes your original data has: sex, rings_comb, and an integer age column called `age`
counts_age_sex <- sex_age_avg_mod$data %>%
  distinct(sex, rings_comb, age) %>%
  # filter(age %in% c(1, 3, 8)) %>%
  count(sex, age, name = "n_ind")

# Season-specific scaling parameters (because date_J_scaled = scale(date_J) within season)
season_scale <- dat %>%
  group_by(season) %>%
  summarise(
    mu_dateJ = mean(date_J, na.rm = TRUE),
    sd_dateJ = sd(date_J, na.rm = TRUE),
    .groups = "drop"
  )

draws <- as_draws_df(fit)

# Grid: sex x age x season
grid <- expand_grid(
  sex    = c("F", "M"),
  age    = c(1:8),
  season = unique(dat$season)
) %>%
  left_join(season_scale, by = "season")

# Helper indices for speed/clarity
idx <- match(draws$.draw, draws$.draw)

# Back-transform t0 to date_J separately for each season
t0_post <- grid %>%
  crossing(.draw = draws$.draw) %>%
  mutate(
    # t0 on scaled scale (set age_dev = 0; plug "age" into avg_age)
    t0_scaled =
      draws$b_t0_Intercept[match(.draw, draws$.draw)] +
      if_else(sex == "M", draws$b_t0_sexM[match(.draw, draws$.draw)], 0) +
      0 * draws$b_t0_age_dev[match(.draw, draws$.draw)] +
      age * draws$b_t0_avg_age[match(.draw, draws$.draw)],

    # invert season-specific scaling: date_J = z * sd + mean
    t0_dateJ = t0_scaled * sd_dateJ + mu_dateJ,

    # map molt-year Julian (origin July 1) onto a dummy calendar year for shared x-axis labels
    origin = as.Date("2021-07-01"),
    date   = origin + round(t0_dateJ - 1),

    age_f  = factor(age, levels = c(1:8))
  ) %>% 
  filter(!(sex == "F" & age == 8))

# --- x-positions at the right-most "tip" (use a high quantile) ---
ann <- t0_post %>%
  mutate(age = as.integer(as.character(age_f))) %>%
  group_by(sex, age_f, age) %>%
  summarise(
    x_tip_age = as.Date(quantile(as.numeric(date), probs = 0.75, na.rm = TRUE),
                    origin = "1970-01-01"),
    x_tip_n = as.Date(quantile(as.numeric(date), probs = 0.995, na.rm = TRUE),
                    origin = "1970-01-01"),
    x_tip_n_ = as.Date("2022-05-01"),
    .groups = "drop"
  ) %>%
  left_join(counts_age_sex, by = c("sex", "age")) %>%
  mutate(
    lab_n = sprintf("paste(italic(n), ' = %s')", n_ind),
    lab_age = age,
    x_lab_age = x_tip_age + 3,
    x_lab_n = x_tip_n,
    x_lab_n_ = x_tip_n_
  ) %>% 
  filter(!(sex == "F" & age == 8)) %>% 
  mutate(x_lab_age = ifelse(sex == "F", x_lab_age - 2, x_lab_age)) %>% 
  mutate(x_lab_age = ifelse(sex == "F" & lab_age == 1, x_lab_age + 4, x_lab_age)) %>% 
  mutate(lab_age = ifelse(lab_age == 1, "age 1", lab_age)) %>% 
  mutate(x_lab_age = ifelse(lab_age == "age 1", x_lab_age - 17, x_lab_age) %>% as.Date())

# 2) annotation positions (one per facet)
sex_annot <- tibble::tibble(
  sex = c("F", "M"),
  date = as.Date(c("2022-01-18", "2022-01-18")),  # any dummy year; must match your molt-year axis mapping
  age_dev_f = factor(8, levels = c(1:8)),
  lab = c("\u2640", "\u2642")
)

plot_t0_age_avg <- 
  ggplot(t0_post, aes(x = date, y = age_f, fill = sex)) +
  geom_density_ridges(scale = 4, rel_min_height = 0.01, alpha = 0.7, color = "white") +
  geom_text(
    data = ann,
    aes(x = x_lab_age, y = age_f, label = lab_age),
    inherit.aes = FALSE,
    hjust = 0, vjust = -0.5,
    parse = FALSE,
    size = 3,
    color = "white"
  ) +
  geom_text(
    data = ann,
    aes(x = x_lab_n_ - 2, y = age_f, label = lab_n),
    inherit.aes = FALSE,
    hjust = 1, vjust = -0.5,
    parse = TRUE,
    size = 2,
    color = "grey40"
  ) +
  facet_wrap(sex ~ ., nrow = 2,
             labeller = as_labeller(c(F = "\u2640", M = "\u2642")),
             strip.position = "right") +
  scale_x_date(
    breaks = seq(from = as.Date("2021-12-21"), to = as.Date("2022-04-05"), by = "21 days"),
    date_labels = "%d %b",
    limits = c(as.Date("2021-12-21"), as.Date("2022-05-01"))
  ) +
  scale_y_discrete(expand = expansion(mult = c(0, 0.15))) +
  labs(
    x = expression("Date of molt inflection point (" * italic(t)[0] * ")"),
    y = "Between-individual age effect"
  ) +
  theme_minimal() +
  theme(
        panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey60"),
    text = element_text(family = "Lato"),
    legend.position = "none",
    axis.text.y = element_blank(),
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    axis.title.x = element_blank(),
    strip.placement   = "outside",
    strip.text.y.right = element_blank(),
    strip.background  = element_blank(),
    plot.margin = margin(5.5, 5.5, 5.5, 5.5)  # more right margin so text isn't cut off
  ) +
  scale_fill_manual(values = c("#F28E2B", "#4E79A7"))

#### baseline: age_avg 
fit <- sex_age_avg_mod$fitted_model
dat <- sex_age_avg_mod$data
draws <- as_draws_df(fit)

# --- counts of individuals ---
counts_age_sex <- dat %>%
  distinct(sex, rings_comb, age) %>%
  count(sex, age, name = "n_ind")

# Grid: sex x age (season not needed for baseline predictions)
grid_b <- expand_grid(
  sex = c("F", "M"),
  age = 1:8
)

# Posterior for baseline (population-level; age_dev = 0)
baseline_post <- grid_b %>%
  crossing(.draw = draws$.draw) %>%
  mutate(
    baseline =
      draws$b_baseline_Intercept[match(.draw, draws$.draw)] +
      if_else(sex == "M", draws$b_baseline_sexM[match(.draw, draws$.draw)], 0) +
      0 * draws$b_baseline_age_dev[match(.draw, draws$.draw)] +
      age * draws$b_baseline_avg_age[match(.draw, draws$.draw)],
    age_f = factor(age, levels = 1:8)
  ) %>%
  filter(!(sex == "F" & age == 8))

# --- annotation positions (right-tail "tips") ---
ann_b <- baseline_post %>%
  mutate(age = as.integer(as.character(age_f))) %>%
  group_by(sex, age_f, age) %>%
  summarise(
    x_tip_age = quantile(baseline, probs = 0.75, na.rm = TRUE),
    x_tip_n   = quantile(baseline, probs = 0.995, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  left_join(counts_age_sex, by = c("sex", "age")) %>%
  mutate(
    lab_n   = sprintf("paste(italic(n), ' = %s')", n_ind),
    lab_age = age,
    x_lab_age = x_tip_age + 0.05,   # small nudge right; adjust as needed
    x_lab_n   = x_tip_n
  ) %>%
  filter(!(sex == "F" & age == 8)) %>%
  mutate(
    # IMPORTANT: use if_else (keeps numeric). You used ifelse() before.
    x_lab_age = if_else(sex == "F", x_lab_age - 0.05, x_lab_age)
  )

# Sex symbols: place at a high baseline value per sex (so they sit to the right)
sex_annot <- baseline_post %>%
  group_by(sex) %>%
  summarise(
    baseline = quantile(baseline, probs = 0.90, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    age_dev_f = factor(8, levels = c(1:8)),
    lab = if_else(sex == "F", "\u2640", "\u2642")
  )

plot_baseline_age_avg <- 
  ggplot(baseline_post, aes(x = baseline, y = age_f, fill = sex)) +
  geom_density_ridges(scale = 4, rel_min_height = 0.01, alpha = 0.7, color = "white") +
  facet_wrap(sex ~ ., nrow = 2,
             labeller = as_labeller(c(F = "\u2640", M = "\u2642")),
             strip.position = "right")  +
  scale_x_continuous(#expand = expansion(mult = c(0.01, 0.20)), 
                     limits = c(3.8, 7), labels = c(4:7), breaks = c(4:7)) +
  labs(
    x = "Baseline score",
    y = "Between-individual age effect"
  ) +
  scale_y_discrete(expand = expansion(mult = c(0, 0.15))) +
  theme_minimal() +
  theme(
        panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"),
    text = element_text(family = "Lato"),
    legend.position = "none",
    axis.text.y = element_blank(),
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
    plot.margin = margin(5.5, 5.5, 5.5, 5.5),
    strip.placement   = "outside",
    strip.text.y.right = element_text(family = "Arial", size = 25, angle = 0),
  ) +
  scale_fill_manual(values = c("#F28E2B", "#4E79A7"))

#### t0: age_dev
fit <- sex_age_avg_mod$fitted_model
dat <- sex_age_avg_mod$data
draws <- as_draws_df(fit)

# --- sample size per sex (since age_dev levels are hypothetical) ---
counts_sex_age_dev <- dat %>%
  distinct(sex, rings_comb, round(age_dev)) %>%
  rename(age_dev = `round(age_dev)`) %>% 
  count(sex, age_dev, name = "n_ind")

# Season-specific scaling parameters (because date_J_scaled = scale(date_J) within season)
season_scale <- dat %>%
  group_by(season) %>%
  summarise(
    mu_dateJ = mean(date_J, na.rm = TRUE),
    sd_dateJ = sd(date_J, na.rm = TRUE),
    .groups = "drop"
  )

# Grid: sex x age_dev x season
avg_age_fix  <- 3
age_dev_vals <- c(-1, 0, 1)

grid_dev <- expand_grid(
  sex     = c("F", "M"),
  age_dev = age_dev_vals,
  season  = unique(dat$season)
) %>%
  left_join(season_scale, by = "season")

# Back-transform t0 to date_J separately for each season
t0_post_dev <- grid_dev %>%
  crossing(.draw = draws$.draw) %>%
  mutate(
    t0_scaled =
      draws$b_t0_Intercept[match(.draw, draws$.draw)] +
      if_else(sex == "M", draws$b_t0_sexM[match(.draw, draws$.draw)], 0) +
      age_dev * draws$b_t0_age_dev[match(.draw, draws$.draw)] +
      avg_age_fix * draws$b_t0_avg_age[match(.draw, draws$.draw)],

    t0_dateJ = t0_scaled * sd_dateJ + mu_dateJ,

    origin = as.Date("2021-07-01"),
    date   = origin + round(t0_dateJ - 1),

    age_dev_f = factor(age_dev, levels = c(-1, 0, 1))
  )

# --- x-positions for labels, plus labels ---
ann_dev <- t0_post_dev %>%
  group_by(sex, age_dev_f, age_dev) %>%
  summarise(
    x_tip_level = as.Date(quantile(as.numeric(date), probs = 0.87, na.rm = TRUE),
                          origin = "1970-01-01"),
    x_tip_n     = as.Date(quantile(as.numeric(date), probs = 0.999, na.rm = TRUE),
                          origin = "1970-01-01"),
    x_tip_n_ = as.Date("2022-05-01"),
    .groups = "drop"
  ) %>%
  left_join(counts_sex_age_dev) %>%
  mutate(
    lab_n     = sprintf("paste(italic(n), ' = %s')", n_ind),
    lab_level = as.character(age_dev),              # white label (the level itself)
    x_lab_level = x_tip_level,
    x_lab_n     = x_tip_n,
    x_lab_n_ = x_tip_n_
  ) %>%
  mutate(
    # keep as Date: use if_else, not ifelse
    x_lab_level = if_else(sex == "F", x_lab_level, x_lab_level)
  ) %>% 
  mutate(lab_level = ifelse(lab_level == -1, "Δage -1", lab_level)) %>% 
  mutate(x_lab_level = ifelse(lab_level == "Δage -1", x_lab_level - 19, x_lab_level))

# 2) annotation positions (one per facet)
sex_annot <- tibble::tibble(
  sex = c("F", "M"),
  date = as.Date(c("2022-01-18", "2022-01-08")),  # any dummy year; must match your molt-year axis mapping
  age_dev_f = factor(1, levels = c(-1, 0, 1)),
  lab = c("\u2640", "\u2642")
)

plot_t0_age_dev <- 
  ggplot(t0_post_dev, aes(x = date, y = age_dev_f, fill = sex)) +
  geom_density_ridges(scale = 4, rel_min_height = 0.01, alpha = 0.7, color = "white") +
  geom_text(
    data = ann_dev,
    aes(x = x_lab_level, y = age_dev_f, label = lab_level),
    inherit.aes = FALSE,
    hjust = 0, vjust = -0.5,
    parse = FALSE,
    size = 3,
    color = "white"
  ) +
  geom_text(
    data = ann_dev,
    aes(x = x_lab_n_ - 2, y = age_dev_f, label = lab_n),
    inherit.aes = FALSE,
    hjust = 1, vjust = -0.5,
    parse = TRUE,
    size = 2,
    color = "grey40"
  ) +
  facet_wrap(
  sex ~ ., nrow = 2,
             labeller = as_labeller(c(F = "\u2640", M = "\u2642")),
             strip.position = "right"
) +
  scale_x_date(
    breaks = seq(from = as.Date("2021-12-21"), to = as.Date("2022-04-05"), by = "21 days"),
    date_labels = "%d %b",
    limits = c(as.Date("2021-12-21"), as.Date("2022-05-01"))
  ) +
  scale_y_discrete(expand = expansion(mult = c(0, 0.2))) +
  labs(
    x = expression("Date of molt inflection point (" * italic(t)[0] * ")"),
    y = expression(atop(italic(age)[dev], scriptstyle("(within-individual age effect)")))
  ) +
  theme_minimal() +
  theme(
        panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey60"),
    text = element_text(family = "Lato"),
    legend.position = "none",
    axis.text.x = element_text(size = 9,
                               angle = 45,
                               hjust = 1,
                               vjust = 1, colour = "grey30"),
    axis.text.y = element_blank(),
    axis.title.y = element_blank(),
    plot.margin = margin(5.5, 5.5, 5.5, 5.5),
    strip.placement   = "outside",
    strip.text.y.right = element_blank(), 
    strip.background  = element_blank()
  ) +
  scale_fill_manual(values = c("#F28E2B", "#4E79A7"))

#### baseline: age_dev
fit <- sex_age_avg_mod$fitted_model
dat <- sex_age_avg_mod$data
draws <- as_draws_df(fit)

# --- sample size per sex x rounded age_dev (as in your t0 plot) ---
counts_sex_age_dev <- dat %>%
  distinct(sex, rings_comb, round(age_dev)) %>%
  rename(age_dev = `round(age_dev)`) %>%
  count(sex, age_dev, name = "n_ind")

# Grid: sex x age_dev
avg_age_fix  <- 3
age_dev_vals <- c(-1, 0, 1)

grid_dev <- expand_grid(
  sex     = c("F", "M"),
  age_dev = age_dev_vals
)

# Posterior for baseline (population-level; hold avg_age constant; vary age_dev)
baseline_post_dev <- grid_dev %>%
  crossing(.draw = draws$.draw) %>%
  mutate(
    baseline =
      draws$b_baseline_Intercept[match(.draw, draws$.draw)] +
      if_else(sex == "M", draws$b_baseline_sexM[match(.draw, draws$.draw)], 0) +
      age_dev * draws$b_baseline_age_dev[match(.draw, draws$.draw)] +
      avg_age_fix * draws$b_baseline_avg_age[match(.draw, draws$.draw)],

    age_dev_f = factor(age_dev, levels = c(-1, 0, 1))
  )

# --- x-positions for labels, plus labels ---
ann_dev <- baseline_post_dev %>%
  group_by(sex, age_dev_f, age_dev) %>%
  summarise(
    x_tip_level = quantile(baseline, probs = 0.87,  na.rm = TRUE),
    x_tip_n     = quantile(baseline, probs = 0.999, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  left_join(counts_sex_age_dev, by = c("sex", "age_dev")) %>%
  mutate(
    lab_n       = sprintf("paste(italic(n), ' = %s')", n_ind),
    lab_level   = as.character(age_dev),
    x_lab_level = x_tip_level + 0.08,   # <- adjust if you want more/less offset
    x_lab_n     = x_tip_n
  ) %>%
  mutate(
    x_lab_level = if_else(sex == "F", x_lab_level - 0.08, x_lab_level)
  )

# Sex symbols: place at a high baseline value per sex (so they sit to the right)
sex_annot <- baseline_post_dev %>%
  group_by(sex) %>%
  summarise(
    baseline = quantile(baseline, probs = 0.90, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    age_dev_f = factor(1, levels = c(-1, 0, 1)),
    lab = if_else(sex == "F", "\u2640", "\u2642")
  )

plot_baseline_age_dev <- 
  ggplot(baseline_post_dev, aes(x = baseline, y = age_dev_f, fill = sex)) +
  geom_density_ridges(scale = 4, rel_min_height = 0.01, alpha = 0.7, color = "white") +
  facet_wrap(sex ~ ., nrow = 2,
             labeller = as_labeller(c(F = "\u2640", M = "\u2642")),
             strip.position = "right")  +
  scale_x_continuous(limits = c(3.8, 7), labels = c(4:7), breaks = c(4:7)) +
  labs(
    x = "Baseline plumage score (<span style='font-family:NotoSans; font-style:italic;'>β</span>)",
    y = "Within-individual age effect"
  ) +
  scale_y_discrete(expand = expansion(mult = c(0, 0.15))) +
  theme_minimal() +
  theme(
        panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey90"),
    text = element_text(family = "Lato"),
    legend.position = "none",
    axis.text.y = element_blank(),
    plot.margin = margin(5.5, 5.5, 5.5, 5.5),
    strip.placement   = "outside",
    strip.text.y.right = element_text(family = "Arial", size = 25, angle = 0),
    axis.title.x = ggtext::element_markdown(family = "Lato"),
    axis.text.x = element_text(colour = "grey30")
  ) +
  scale_fill_manual(values = c("#F28E2B", "#4E79A7"))

#### predicted trajectories based on age class
fit_ <- sex_age_class_mod$fitted_model
dat_ <- sex_age_class_mod$data

dates_for_plot <- 
  data.frame(date_ = as.Date(c(as.Date("2021-07-01"):as.Date("2022-06-30"))), 
             date_J = c(1:365))

# Create new data for prediction
newdata_age_class_pop_prediction <- 
  expand.grid(
  date_J_scaled = seq(from = min(dat_$date_J_scaled),
               to = max(dat_$date_J_scaled),
               length.out = max(dat_$date_J) - min(dat_$date_J) + 1),
  sex = c("F", "M"),
  age_class = c("FY", "AFY"),
  rings_comb = NA,   # exclude group-level effects
  season = NA
) %>% 
  left_join(., data.frame(date_J = seq(from = min(dat_$date_J),
                                      to = max(dat_$date_J),
                                      length.out = max(dat_$date_J) - min(dat_$date_J) + 1),
                         date_J_scaled = seq(from = min(dat_$date_J_scaled),
                                             to = max(dat_$date_J_scaled),
                                             length.out = max(dat_$date_J) - min(dat_$date_J) + 1)),
            by = "date_J_scaled") %>% 
  left_join(., dates_for_plot, by = "date_J")

fitted_vals <- fitted(
  fit_,
  newdata = newdata_age_class_pop_prediction,
  re_formula = NA,  # exclude random effects
  summary = TRUE
)

newdata_age_class_pop_prediction$fit <- fitted_vals[, "Estimate"]
newdata_age_class_pop_prediction$lower <- fitted_vals[, "Q2.5"]
newdata_age_class_pop_prediction$upper <- fitted_vals[, "Q97.5"]

# Extract unique combinations of grouping variables
id_season_sex_age_class <- 
  fit_$data %>%
  distinct(rings_comb, season, sex, age_class)

# Expand grid for predictions
newdata_age_class <- expand_grid(
  date_J_scaled = seq(from = min(dat_$date_J_scaled),
               to = max(dat_$date_J_scaled),
               length.out = max(dat_$date_J) - min(dat_$date_J) + 1),
  id_season_sex_age_class
)

# Get predicted curves with uncertainty
newdata_age_class_ind_prediction <- 
  add_epred_draws(fit_, 
                  newdata = newdata_age_class, 
                  ndraws = 250) %>%   # You can increase `ndraws` for smoother curves
  left_join(., data.frame(date_J = seq(from = min(dat_$date_J),
                                      to = max(dat_$date_J),
                                      length.out = max(dat_$date_J) - min(dat_$date_J) + 1),
                         date_J_scaled = seq(from = min(dat_$date_J_scaled),
                                             to = max(dat_$date_J_scaled),
                                             length.out = max(dat_$date_J) - min(dat_$date_J) + 1)),
            by = "date_J_scaled") %>% 
  left_join(., dates_for_plot, by = "date_J")

age_class_pop_curves <-
  ggplot() +
  geom_rect(aes(xmin = min(found_dates_all_$date_), xmax = max(found_dates_all_$date_),
                ymin = -Inf, ymax = Inf),
            fill = "grey70", alpha = 0.2) +
  annotate("text",
         x = min(found_dates_all_$date_) + (max(found_dates_all_$date_) - min(found_dates_all_$date_)) / 2,
         y = 0.65,
         label = "Nest initiation period",
         fontface = "italic",
         family = "Lato",
         size = 2.5,  # adjust size as needed
         hjust = 0.5,
         vjust = 0, colour = "grey40") +
  # Add the double-headed arrow
  geom_segment(aes(x = min(found_dates_all_$date_), xend = max(found_dates_all_$date_), y = 0.55, yend = 0.55),
               arrow = arrow(ends = "both", type = "closed", length = unit(0.2, "cm")),
               linewidth = 0.7, color = "grey70") +
  geom_line(data = newdata_age_class_ind_prediction %>% filter(age_class == "FY"), 
            aes(x = date_, y = .epred, 
                group = interaction(rings_comb, .draw, season)), color =  "#E15759", alpha = 0.07) +
  geom_line(data = newdata_age_class_ind_prediction %>% filter(age_class == "AFY"), 
            aes(x = date_, y = .epred, 
                group = interaction(rings_comb, .draw, season)), color =  "#59A14F", alpha = 0.02) +
  geom_ribbon(
    data = newdata_age_class_pop_prediction %>% mutate(age_class = factor(age_class, levels = c("FY", "AFY")),
                                                       sex = factor(sex, levels = c("M", "F"))) %>% filter(sex == "M"),
    aes(x = date_, ymin = lower, ymax = upper, fill = age_class,
        group = interaction(sex, age_class)),
    alpha = 1,
    color = "white"
  ) +
  geom_ribbon(
    data = newdata_age_class_pop_prediction %>% mutate(age_class = factor(age_class, levels = c("FY", "AFY")),
                                                       sex = factor(sex, levels = c("M", "F"))) %>% filter(sex == "F"),
    aes(x = date_, ymin = lower, ymax = upper, fill = age_class,
        group = interaction(sex, age_class)),
    alpha = 1,
    color = "white"
  ) +
  scale_fill_manual(values =  c("#E15759", "#59A14F")) +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(family = "Lato"),
        axis.ticks = element_blank(),
        axis.title.y = element_text(size = 11, colour = "black"),
        axis.text.y = element_text(size = 9, colour = "grey30"),
        panel.spacing = unit(1, "cm"),
        panel.border = element_blank(),
        panel.background = element_rect(fill = 'transparent', color = NA),
        plot.background = element_rect(fill = 'transparent', color = NA), 
        strip.background = element_rect(fill = 'grey90', color = NA),
        strip.text = element_text(size = 13, colour = "grey40"),
        panel.grid.major.y = element_line(size = 0.25, lineend = "round", colour = "grey90"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_line(size = 0.25, lineend = "round", colour = "grey60"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(size = 9, 
                                   angle = 45, 
                                   hjust = 1, 
                                   vjust = 1, colour = "grey30"),
        plot.margin = unit(c(0, -15, 5.5, 0), "pt")) +
  scale_x_date(expand = c(0.01, 0.01),
               limits = c(as.Date("2021-07-20"), as.Date("2022-05-04")),
               breaks = seq(from = as.Date("2021-07-27"), to = as.Date("2022-04-12"), by = "21 days"),
               date_labels = "%d %b") +
  scale_y_continuous(limits = c(0.5, 7.5), breaks = seq(1, 7, 1)) +
  ylab("Rufous band score")

# --- desired legend anchor in data coords ---
x_target <- as.Date("2022-02-26")
y_target <- 6.75

# --- convert to npc (0-1) using your plot limits ---
x_limits <- c(as.Date("2021-07-20"), as.Date("2022-05-04"))
y_limits <- c(0.5, 7.5)

x_npc <- as.numeric(x_target - x_limits[1]) / as.numeric(x_limits[2] - x_limits[1])
y_npc <- (y_target - y_limits[1]) / (y_limits[2] - y_limits[1])

age_class_pop_curves <-
  age_class_pop_curves +
  scale_fill_manual(
    values = c(FY = "#E15759", AFY = "#59A14F"),
    breaks = c("FY", "AFY"),
    labels = c(FY = "First-year", AFY = "After-first-year"),
    name   = "Age class"
  ) +
  theme(
    legend.position = c(x_npc, y_npc),
    legend.justification = c(0.5, 0.5),  # center the legend at that point
    legend.direction = "horizontal",     # optional; remove if you want vertical
    legend.background = element_rect(fill = scales::alpha("white", 0.7), colour = NA),
    legend.title = element_text(size = 10),
    legend.text  = element_text(size = 9)
  ) +
  guides(
    fill = guide_legend(
      title.position = "top",
      title.hjust = 0.5
    )
  )

age_class_pop_curves <-
  age_class_pop_curves +
  theme(
    # transparent legend background
    legend.background = element_blank(),
    legend.box.background = element_blank(),

    # smaller legend text
    legend.title = element_text(size = 9),
    legend.text  = element_text(size = 8),

    # smaller legend keys
    legend.key.size = unit(8, "pt")
  ) +
  guides(
    fill = guide_legend(
      title.position = "top",
      title.hjust = 0,     # left-justify within the guide
      label.position = "right",
      keywidth  = unit(10, "pt"),
      keyheight = unit(8, "pt")
    )
  )

#### comboplot
left_col  <- plot_baseline_age_avg / plot_baseline_age_dev
right_col <- plot_t0_age_avg       / plot_t0_age_dev

combo_age_plot <- (left_col | right_col) +
  plot_layout(widths = c(0.92, 1.08))

final_plot <-
  age_class_pop_curves /
  combo_age_plot +
  plot_layout(heights = c(0.75, 2))

final_plot

Auxiliary analyses

Sensitivity analyses for the breeding phenology model

These two conservative variants of the breeding phenology dataset and model are presented here as critiques of the main breeding phenology analysis. dat_breeding_mod2 keeps only birds that were photographed after their last recorded breeding activity, while dat_breeding_mod3 excludes birds whose final breeding record was earlier than the seasonal median.

# subset to include only birds that were photographed after their last breeding activity
dat_breeding2 <- 
  dat_breeding %>%
  group_by(rings_comb, season) %>%
  mutate(
    max_date = max(date_J, na.rm = TRUE)
  ) %>%
  filter(last_breeding_date_J < max_date)

# wrangle data and subset to dat_breeding2
dat_breeding_mod2 <- dat_breeding2 %>%
  # subset to columns of interest
  select(score, date_J, rings_comb, sex, season, 
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum) %>%
  # reduce to distinct data
  distinct() %>% 
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>% 
  ungroup() %>%
  # drop NAs
  drop_na(score, date_J, sex, rings_comb, season, 
          last_breeding_date_J_snum, 
          first_breeding_date_J_snum, 
          last_nest_first_obs_J_snum) %>% 
  # convert to dataframe
  as.data.frame()

# wrangle data and subset to individuals that had late last breeding dates
dat_breeding_mod3 <- dat_breeding %>%
  # remove individuals that had early last breeding dates (i.e., likely dispersed and were not seen again)
  filter(last_breeding_date_J_snum > 0) %>% 
  # subset to columns of interest
  select(score, date_J, rings_comb, sex, season, 
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum) %>%
  # reduce to distinct data
  distinct() %>% 
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>% 
  ungroup() %>%
  # drop NAs
  drop_na(score, date_J, sex, rings_comb, season, 
          last_breeding_date_J_snum, 
          first_breeding_date_J_snum, 
          last_nest_first_obs_J_snum) %>% 
  # convert to dataframe
  as.data.frame()

fit_mod_sex_all_3_breed2 <- brm(
  formula = mod_sex_all_3_breed,
  data = dat_breeding_mod2,
  family = gaussian(),
  prior = moult_priors,
  control = list(adapt_delta = 0.999, max_treedepth = 15),
  chains = 4,
  cores = 4,
  iter = 10000,
  warmup = 5000,
  backend = "cmdstanr",
  seed = 123
)

fit_mod_sex_all_3_breed3 <- brm(
  formula = mod_sex_all_3_breed,
  data = dat_breeding_mod3,
  family = gaussian(),
  prior = moult_priors,
  control = list(adapt_delta = 0.999, max_treedepth = 15),
  chains = 4,
  cores = 4,
  iter = 10000,
  warmup = 5000,
  backend = "cmdstanr",
  seed = 123
)

Sensitivity analysis for the continuous age model

This auxiliary model repeats the continuous-age analysis after excluding age-1 observations, so that the estimated within- and between-individual age effects can be checked without the first-year transition.

# calculate average age in the dataset of individuals of known age, excluding age 1
avg_ages2 <- 
  dat_breeding %>%
  distinct(rings_comb, age) %>% 
  filter(!is.na(age)) %>% 
  filter(age > 1) %>% 
  group_by(rings_comb) %>% 
  summarise(avg_age = mean(age, na.rm = TRUE))

# wrangle data
dat_age2_mod <- dat_breeding %>%
  # subset to columns of interest
  select(score, date_J, rings_comb, sex, season, age) %>%
  # reduce to distinct data
  distinct() %>% 
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>%
  # center age
  filter(!is.na(age)) %>%
  filter(age > 1) %>% 
  left_join(., avg_ages2, by = "rings_comb") %>% 
  # age deviance
  ungroup() %>%
  mutate(age_dev = age - avg_age) %>% 
  # drop NAs
  drop_na(score, date_J, sex, rings_comb, 
          season, avg_age, age_dev) %>% 
  # convert to dataframe
  as.data.frame()

fit_mod_sex_age2_avg <- brm(
  formula = mod_sex_age_avg,
  data = dat_age2_mod,
  family = gaussian(),
  prior = moult_priors,
  control = list(adapt_delta = 0.999, max_treedepth = 15),
  chains = 4,
  cores = 4,
  iter = 10000,
  warmup = 5000,
  backend = "cmdstanr",
  save_pars = save_pars(all = TRUE),
  seed = 123
)

Do migratory strategies differ in breeding phenology?

# wrangle data
dat_migratory_status_breed <- 
  dat_breeding %>%
  # subset to columns of interest and reduce to distinct data
  distinct(score, date_J, rings_comb, 
         sex, season, migratory_status, 
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum) %>%
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>% 
  ungroup() %>%
  # drop NAs
  drop_na(score, date_J_scaled, rings_comb, 
         sex, season, migratory_status, 
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum) %>% 
  # convert to dataframe and reduce
  as.data.frame() %>% 
  distinct(rings_comb, migratory_status, 
           last_nest_first_obs_J_snum, 
           first_breeding_date_J_snum, 
           last_breeding_date_J_snum)

# Weakly-informative priors suitable for season-standardized outcomes
priors_gauss <- c(
  prior(normal(0, 2), class = "b"),          # fixed effects (incl. migratory_status)
  prior(normal(0, 5), class = "Intercept"),  # intercept
  prior(exponential(1), class = "sd"),       # random effect SDs
  prior(exponential(1), class = "sigma")     # residual SD
)

# 1) last_nest_first_obs_J_snum ~ migratory_status + (1 | rings_comb)
mod_mig_lay_brms <- brm(
  last_nest_first_obs_J_snum ~ migratory_status + (1 | rings_comb),
  data   = dat_migratory_status_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123
)

# 2) first_breeding_date_J_snum ~ migratory_status + (1 | rings_comb)
mod_mig_first_brms <- brm(
  first_breeding_date_J_snum ~ migratory_status + (1 | rings_comb),
  data   = dat_migratory_status_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123
)

# 3) last_breeding_date_J_snum ~ migratory_status + (1 | rings_comb)
mod_mig_last_brms <- brm(
  last_breeding_date_J_snum ~ migratory_status + (1 | rings_comb),
  data   = dat_migratory_status_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123
)

Fixed effects with 95% CrI

migratory status is not associated with any of the three breeding phenology traits

fixef(mod_mig_first_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% round(2)

                  Estimate Est.Error  Q2.5   Q50 Q97.5
Intercept            -2.08      2.95 -7.83 -2.06  3.72
migratory_statusR     0.25      1.95 -3.55  0.24  4.12

fixef(mod_mig_last_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% round(2)

                  Estimate Est.Error  Q2.5   Q50 Q97.5
Intercept             2.34      3.04 -3.55  2.34  8.37
migratory_statusR    -0.38      1.93 -4.17 -0.38  3.42

fixef(mod_mig_lay_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% round(2)

                  Estimate Est.Error  Q2.5  Q50 Q97.5
Intercept             2.67      2.83 -2.88 2.66  8.25
migratory_statusR     0.39      1.92 -3.42 0.39  4.18

Is age associated with breeding phenology?

# average ages including first year
avg_ages <- 
  dat_breeding %>%
  distinct(rings_comb, age) %>% 
  filter(!is.na(age)) %>% 
  group_by(rings_comb) %>% 
  summarise(avg_age = mean(age, na.rm = TRUE))

# average ages excluding first year
avg_ages2 <- 
  dat_breeding %>%
  distinct(rings_comb, age) %>% 
  filter(!is.na(age)) %>% 
  filter(age > 1) %>% 
  group_by(rings_comb) %>% 
  summarise(avg_age = mean(age, na.rm = TRUE))

# wrangle data
dat_age_breed <- dat_breeding %>%
  # subset to columns of interest
  select(score, date_J, rings_comb, sex, season, age, 
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum, 
         last_breeding_date_J_s1num, 
         first_breeding_date_J_s1num,
         last_nest_first_obs_J_s1num) %>%
  # reduce to distinct data
  distinct() %>% 
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>%
  # center age
  filter(!is.na(age)) %>% 
  left_join(., avg_ages, by = "rings_comb") %>% 
  # age deviation
  ungroup() %>%
  mutate(age_dev = age - avg_age) %>% 
  # mutate(age_centered = age - mean(age, na.rm = TRUE)) %>% 
  # drop NAs
  drop_na(score, date_J, sex, rings_comb, 
          season, avg_age, age_dev, 
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum, 
         last_breeding_date_J_s1num, 
         first_breeding_date_J_s1num,
         last_nest_first_obs_J_s1num) %>% 
  # convert to dataframe
  as.data.frame() %>% 
  select(rings_comb, avg_age, sex, age_dev, 
           last_breeding_date_J_snum, 
           first_breeding_date_J_snum,
           last_nest_first_obs_J_snum, 
           last_breeding_date_J_s1num, 
           first_breeding_date_J_s1num,
           last_nest_first_obs_J_s1num) %>% 
  distinct()

# wrangle a version of the data that excludes first year
dat_age2_breed <- dat_breeding %>%
  # subset to columns of interest
  select(score, date_J, rings_comb, sex, season, age, 
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum, 
         last_breeding_date_J_s1num, 
         first_breeding_date_J_s1num,
         last_nest_first_obs_J_s1num) %>%
  # reduce to distinct data
  distinct() %>% 
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>%
  # center age
  filter(!is.na(age)) %>%
  filter(age > 1) %>% 
  left_join(., avg_ages2, by = "rings_comb") %>% 
  # age deviance
  ungroup() %>%
  mutate(age_dev = age - avg_age) %>% 
  # drop NAs
  drop_na(score, date_J, sex, rings_comb, 
          season, avg_age, age_dev, 
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum, 
         last_breeding_date_J_s1num, 
         first_breeding_date_J_s1num,
         last_nest_first_obs_J_s1num) %>% 
  # convert to dataframe
  as.data.frame() %>% 
  select(rings_comb, sex, avg_age, age_dev, 
           last_breeding_date_J_snum, 
           first_breeding_date_J_snum,
           last_nest_first_obs_J_snum, 
           last_breeding_date_J_s1num, 
           first_breeding_date_J_s1num,
           last_nest_first_obs_J_s1num) %>% 
  distinct()

# Weakly-informative priors suitable for season-standardized outcomes
priors_gauss <- c(
  prior(normal(0, 2), class = "b"),          # fixed effects (incl. migratory_status)
  prior(normal(0, 5), class = "Intercept"),  # intercept
  prior(exponential(1), class = "sd"),       # random effect SDs
  prior(exponential(1), class = "sigma")     # residual SD
)

# 1) last_nest_first_obs_J_snum 
mod_age_lay_brms <- brm(
  last_nest_first_obs_J_snum ~ age_dev + avg_age + (1 | rings_comb),
  data   = dat_age_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

# 2) first_breeding_date_J_snum 
mod_age_first_brms <- brm(
  first_breeding_date_J_snum ~ age_dev + avg_age + (1 | rings_comb),
  data   = dat_age_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

# 3) last_breeding_date_J_snum
mod_age_last_brms <- brm(
  last_breeding_date_J_snum ~ age_dev + avg_age + (1 | rings_comb),
  data   = dat_age_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

# 1) last_nest_first_obs_J_snum 
mod_age2_lay_brms <- brm(
  last_nest_first_obs_J_snum ~ age_dev + avg_age + (1 | rings_comb),
  data   = dat_age2_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123
)

# 2) first_breeding_date_J_snum 
mod_age2_first_brms <- brm(
  first_breeding_date_J_snum ~ age_dev + avg_age + (1 | rings_comb),
  data   = dat_age2_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123
)

# 3) last_breeding_date_J_snum
mod_age2_last_brms <- brm(
  last_breeding_date_J_snum ~ age_dev + avg_age + (1 | rings_comb),
  data   = dat_age2_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123
)

Fixed effects with 95% CrI

Continuous age is not associated with breeding phenology

fixef(mod_age_first_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% 
  round(2)

          Estimate Est.Error   Q2.5  Q50 Q97.5
Intercept     0.71      7.02 -13.10 0.70 14.33
age_dev       0.10      1.86  -3.50 0.09  3.79
avg_age       0.69      1.59  -2.39 0.71  3.80

fixef(mod_age_last_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% 
  round(2)

          Estimate Est.Error  Q2.5   Q50 Q97.5
Intercept     7.22      7.50 -7.44  7.23 21.90
age_dev      -0.78      1.89 -4.51 -0.77  2.92
avg_age      -1.80      1.70 -5.13 -1.82  1.53

fixef(mod_age_lay_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% 
  round(2)

          Estimate Est.Error  Q2.5   Q50 Q97.5
Intercept     5.47      6.74 -7.72  5.48 18.50
age_dev      -0.34      1.83 -3.95 -0.34  3.25
avg_age      -0.18      1.54 -3.17 -0.17  2.84

fixef(mod_age2_first_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% 
  round(2)

          Estimate Est.Error   Q2.5   Q50 Q97.5
Intercept    -4.25      7.50 -18.87 -4.30 10.34
age_dev       0.04      1.89  -3.68  0.04  3.70
avg_age       1.64      1.65  -1.59  1.65  4.86

fixef(mod_age2_last_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% 
  round(2)

          Estimate Est.Error   Q2.5   Q50 Q97.5
Intercept     5.72      8.05 -10.28  5.86 21.32
age_dev      -0.77      1.90  -4.47 -0.77  2.93
avg_age      -1.48      1.79  -4.97 -1.50  2.04

fixef(mod_age2_lay_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% 
  round(2)

          Estimate Est.Error  Q2.5   Q50 Q97.5
Intercept     8.51      7.50 -6.20  8.71 23.04
age_dev      -0.29      1.86 -3.87 -0.31  3.37
avg_age      -0.78      1.65 -4.01 -0.79  2.49

Do males and females differ in breeding phenology?

dat_sex_breed <- dat_breeding %>%
  # subset to columns of interest
  select(score, date_J, rings_comb, 
         sex, season,
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum) %>%
  # reduce to distinct data
  distinct() %>% 
  # scale the date_J variable
  group_by(season) %>% 
  mutate(date_J_scaled = as.vector(scale(date_J))) %>% 
  ungroup() %>%
  # drop NAs
  drop_na(score, date_J_scaled, rings_comb, 
         sex, season,
         last_breeding_date_J_snum, 
         first_breeding_date_J_snum,
         last_nest_first_obs_J_snum) %>% 
  # convert to dataframe
  as.data.frame() %>% 
  distinct(rings_comb, sex, 
           last_nest_first_obs_J_snum, 
           first_breeding_date_J_snum, 
           last_breeding_date_J_snum)

# Weakly-informative priors suitable for season-standardized outcomes
priors_gauss <- c(
  prior(normal(0, 2), class = "b"),          # fixed effects (incl. migratory_status)
  prior(normal(0, 5), class = "Intercept"),  # intercept
  prior(exponential(1), class = "sd"),       # random effect SDs
  prior(exponential(1), class = "sigma")     # residual SD
)

# 1) last_nest_first_obs_J_snum ~ sex + (1 | rings_comb)
mod_sex_lay_brms <- brm(
  last_nest_first_obs_J_snum ~ sex + (1 | rings_comb),
  data   = dat_sex_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123
)

# 2) first_breeding_date_J_snum ~ sex + (1 | rings_comb)
mod_sex_first_brms <- brm(
  first_breeding_date_J_snum ~ sex + (1 | rings_comb),
  data   = dat_sex_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123
)

# 3) last_breeding_date_J_snum ~ sex + (1 | rings_comb)
mod_sex_last_brms <- brm(
  last_breeding_date_J_snum ~ sex + (1 | rings_comb),
  data   = dat_sex_breed,
  family = gaussian(),
  prior  = priors_gauss,
  chains = 4, cores = 4, iter = 4000, warmup = 1000,
  seed   = 123
)

mod_sex_lay_brms <- readRDS(here("out/mod_mig_lay_brms.rds"))
mod_sex_first_brms <- readRDS(here("out/mod_mig_first_brms.rds"))
mod_sex_last_brms <- readRDS(here("out/mod_mig_last_brms.rds"))

# Summaries (posterior medians/means + CrIs)
summary(mod_sex_lay_brms$fitted_model)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: last_nest_first_obs_J_snum ~ migratory_status + (1 | rings_comb) 
   Data: dat_migratory_status_breed (Number of observations: 75) 
  Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
         total post-warmup draws = 12000

Multilevel Hyperparameters:
~rings_comb (Number of levels: 35) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.99      0.97     0.02     3.54 1.00     9900     5994

Regression Coefficients:
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept             2.67      2.83    -2.88     8.25 1.00    17631     8692
migratory_statusR     0.39      1.92    -3.42     4.18 1.00    16993     8600

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    23.23      1.61    20.33    26.65 1.00    15816     8797

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

summary(mod_sex_first_brms$fitted_model)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: first_breeding_date_J_snum ~ migratory_status + (1 | rings_comb) 
   Data: dat_migratory_status_breed (Number of observations: 75) 
  Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
         total post-warmup draws = 12000

Multilevel Hyperparameters:
~rings_comb (Number of levels: 35) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     1.36      1.60     0.03     6.03 1.00     4521     3289

Regression Coefficients:
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept            -2.08      2.95    -7.83     3.72 1.00    13054     7736
migratory_statusR     0.25      1.95    -3.55     4.12 1.00    13736     7834

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    24.84      1.71    21.71    28.41 1.00     8557     4492

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

summary(mod_sex_last_brms$fitted_model)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: last_breeding_date_J_snum ~ migratory_status + (1 | rings_comb) 
   Data: dat_migratory_status_breed (Number of observations: 75) 
  Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
         total post-warmup draws = 12000

Multilevel Hyperparameters:
~rings_comb (Number of levels: 35) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     1.02      1.04     0.03     3.86 1.00    11108     7349

Regression Coefficients:
                  Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept             2.34      3.04    -3.55     8.37 1.00    18525     8927
migratory_statusR    -0.38      1.93    -4.17     3.42 1.00    19059     9387

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    26.93      1.78    23.69    30.70 1.00    18467     8882

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

# Random-effect SDs and correlations
VarCorr(mod_sex_first_brms$fitted_model)

$rings_comb
$rings_comb$sd
          Estimate Est.Error       Q2.5    Q97.5
Intercept 1.362425  1.599855 0.03031744 6.032592


$residual__
$residual__$sd
 Estimate Est.Error     Q2.5   Q97.5
 24.83942  1.711846 21.70594 28.4091

VarCorr(mod_sex_last_brms$fitted_model)

$rings_comb
$rings_comb$sd
          Estimate Est.Error       Q2.5    Q97.5
Intercept   1.0216  1.037058 0.02619513 3.862986


$residual__
$residual__$sd
 Estimate Est.Error     Q2.5    Q97.5
 26.92988  1.784677 23.68803 30.70043

VarCorr(mod_sex_lay_brms$fitted_model)

$rings_comb
$rings_comb$sd
           Estimate Est.Error       Q2.5    Q97.5
Intercept 0.9854944 0.9735003 0.02296337 3.535091


$residual__
$residual__$sd
 Estimate Est.Error     Q2.5    Q97.5
 23.23398  1.614583 20.33116 26.65414

Fixed effects with 95% Cr

Males and female exhibit the same breeding phenologies

fixef(mod_sex_first_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% round(2)

                  Estimate Est.Error  Q2.5   Q50 Q97.5
Intercept            -2.08      2.95 -7.83 -2.06  3.72
migratory_statusR     0.25      1.95 -3.55  0.24  4.12

fixef(mod_sex_last_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% round(2)

                  Estimate Est.Error  Q2.5   Q50 Q97.5
Intercept             2.34      3.04 -3.55  2.34  8.37
migratory_statusR    -0.38      1.93 -4.17 -0.38  3.42

fixef(mod_sex_lay_brms$fitted_model, probs = c(0.025, 0.5, 0.975)) %>% round(2)

                  Estimate Est.Error  Q2.5  Q50 Q97.5
Intercept             2.67      2.83 -2.88 2.66  8.25
migratory_statusR     0.39      1.92 -3.42 0.39  4.18

Is sex composition independent of age class?

Yes

sex_age_class_mod <- readRDS(here("out/sex_age_class_mod.rds"))
sex_age <- 
  sex_age_class_mod$data %>% 
  distinct(rings_comb, sex, age_class) %>% 
  group_by(rings_comb) %>% 
  arrange(desc(age_class)) %>% 
  slice(1)

tab_sex_age <- table(
  sex_age$sex,
  sex_age$age_class
)

fisher.test(tab_sex_age)

    Fisher's Exact Test for Count Data

data:  tab_sex_age
p-value = 0.4192
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  0.3313901 29.8213437
sample estimates:
odds ratio 
  2.427811 

Is sex composition independent of migratory status?

Yes

sex_migratory_status_mod <- readRDS(here("out/sex_migratory_status_mod.rds"))

sex_mig <- 
  sex_migratory_status_mod$data %>% 
  select(rings_comb, sex, migratory_status) %>% 
  distinct()

tab_sex_mig <- table(
  sex_mig$sex,
  sex_mig$migratory_status
)

fisher.test(tab_sex_mig)

    Fisher's Exact Test for Count Data

data:  tab_sex_mig
p-value = 0.1169
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.03495678 1.69930660
sample estimates:
odds ratio 
 0.2705294