Evaluating the effects of tracking devices on survival, breeding success, behaviour, and condition of a small, partially migratory shorebird

Supplementary Material A: Journal of Avian Biology 10.1002/jav.03490

Authors
Affiliations

Department of Ornithology, Max Planck Institute for Biological Intelligence, Eberhard-Gwinner Straße 7, 82319 Seewiesen, Germany

Emma Williams

Fauna Science Team, Department of Conservation Christchurch Office, Christchurch Mail Centre, Private Bag 4715, Christchurch, 8140, New Zealand

Ailsa McGilvary-Howard

Kaikōura Banded Dotterel Project, 1 Maui Street, Kaikōura, New Zealand

Ted Howard

Kaikōura Banded Dotterel Project, 1 Maui Street, Kaikōura, New Zealand

Tony Habraken

Port Waikato Banded Dotterel Project, Port Waikato, New Zealand

Colin O`Donnell

Fauna Science Team, Department of Conservation Christchurch Office, Christchurch Mail Centre, Private Bag 4715, Christchurch, 8140, New Zealand

Clemens Küpper

Behavioural Genetics and Evolutionary Ecology, Max Planck Institute for Biological Intelligence, Eberhard-Gwinner Straße 5, 82319 Seewiesen, Germany

Bart Kempenaers

Department of Ornithology, Max Planck Institute for Biological Intelligence, Eberhard-Gwinner Straße 7, 82319 Seewiesen, Germany

Published

July 2, 2025

Prerequisites

R packages

# Set a CRAN mirror before installing packages
options(repos = c(CRAN = "https://cloud.r-project.org"))

# a vector of all the packages needed in the project
packages_required_in_project <-
  c("apis", "BaSTA", "brms", "broom.mixed", "colorspace", "corrplot", "data.table", 
    "effects", "elevatr", "geosphere", "ggblend", "ggeffects", "ggmap",
    "ggnewscale", #"ggsn", 
    "giscoR", "glue", "googlesheets4", "gt", "gtsummary", 
    "here", "hms", "lme4", "lubridate", "mapview", "mgcv", "multcomp", "MuMIn", "nnet", 
    "partR2", "patchwork", "RColorBrewer", "readxl", "remotes", "RMark", "rnaturalearth", 
    "rptR", "scales", "sf", "showtext", "smatr", "sp", "stringr", "tidybayes", 
    "tidyterra", "tidyverse", "adehabitatLT", "magrittr", "MASS", "leaflet", "shiny", 
    "leaflet.extras", "ggrepel", "ggpmisc", "grid", "knitr", "kableExtra")

# of the required packages, check if some need to be installed
new.packages <- 
  packages_required_in_project[!(packages_required_in_project %in% 
                                   installed.packages()[,"Package"])]

# install all packages that are not locally available
if(length(new.packages)) install.packages(new.packages)

# load all the packages into the current R session
lapply(packages_required_in_project, require, character.only = TRUE)

Custom functions

# function that does the opposite of "%in%"
`%!in%` = Negate(`%in%`)

# scaled mass index calculation from Peig and Green (2009)
# https://doi.org/10.1111/j.1365-2435.2010.01751.x
scaledMassIndex <-
  function(x, y, x.0 = mean(x)) {
    logM.ols <- lm(log(y) ~ log(x))
    logM.rob <- rlm(log(y) ~ log(x), method = "M")
    b.msa.ols <- coef(sma(log(y) ~ log(x)))[2]
    b.msa.rob <- coef(sma(log(y) ~ log(x), robust = T))[2]
    SMI.ols <- y * (x.0 / x) ^ b.msa.ols
    SMI.rob <- y * (x.0 / x) ^ b.msa.rob
    res <- data.frame(SMI.ols, SMI.rob, x, y)
    pred.DT <-
      data.table(x = seq(min(x), max(x), length = 100)) %>%
      .[, y.ols := predict(logM.ols, newdata = .) %>% exp] %>%
      .[, y.rob := predict(logM.rob, newdata = .) %>% exp]
    attr(res, "b.msa") <- c(ols = b.msa.ols, rob = b.msa.rob)
    return(res)
  }

Body Condition Dynamics

specify attachment weights

harness = 0.29
colour_band = 0.08
metal_band = 0.16
PTT = 1.8
PTT2 = 2.0
GPS = 1.2

PCA of structural traits

body_traits <- 
  read.csv(here("data/body_traits_season.csv"))

check most relevent structural traits with PCA

body_traits_pca <- 
  body_traits %>% 
  group_by(ring) %>% 
  arrange(measure) %>% 
  slice(1) %>% 
  ungroup() %>% 
  dplyr::select(culmen, tarsus, wing) %>%
  na.omit() %>% 
  princomp()

# check the PCA results
summary(body_traits_pca)
Importance of components:
                          Comp.1     Comp.2     Comp.3
Standard deviation     3.0375002 0.92702190 0.70667715
Proportion of Variance 0.8716353 0.08118619 0.04717851
Cumulative Proportion  0.8716353 0.95282149 1.00000000
biplot(body_traits_pca, cex = 0.7)

loadings(body_traits_pca)

Loadings:
       Comp.1 Comp.2 Comp.3
culmen         0.232  0.973
tarsus         0.973 -0.232
wing    1.000              

               Comp.1 Comp.2 Comp.3
SS loadings     1.000  1.000  1.000
Proportion Var  0.333  0.333  0.333
Cumulative Var  0.333  0.667  1.000
# bind PC scores to the original dataframe
body_traits_pca_scores <-
  body_traits %>% 
  group_by(ring) %>% 
  arrange(measure) %>% 
  slice(1) %>% 
  bind_cols(., body_traits_pca$scores[, 1], body_traits_pca$scores[, 2], body_traits_pca$scores[, 3]) %>% 
  rename(structure_pc1 = `...13`,
         structure_pc2 = `...14`,
         structure_pc3 = `...15`) %>% 
  na.omit() %>%
  ungroup()

check correlations of traits with weight to determine best trait to use for SMI

visualize the relationship between the principle components and body mass

visualize the relationship between the wing length and body mass

visualize the relationship between the tarsus length and body mass

Scaled Mass Index wrangle

average all repeated measures of structural traits for each bird (i.e., variation within individuals for these traits is assumed to be observer/instrument error)

body_traits_avg <- 
  body_traits %>% 
  group_by(ring) %>% 
  mutate(culmen = mean(culmen, na.rm = TRUE), 
         tarsus = mean(tarsus, na.rm = TRUE), 
         wing = mean(wing, na.rm = TRUE))

determine which structural trait to use for the scaled mass index calculation

see text in Peig and Green (2009) after Eq. 2: “…We recommend the use of the single ‘L’ variable which has the strongest correlation with M on a log-log scale, since this is likely to be the L that best explains that fraction of mass associated with structural size” Conclude to use wing for SMI transformation because it has the largest r and based on the PCA investigation above, wing did just as good as PC1 when explaining variation in weight

cor.test(log(body_traits_avg$tarsus), log(body_traits_avg$weight)) # r = 0.130

    Pearson's product-moment correlation

data:  log(body_traits_avg$tarsus) and log(body_traits_avg$weight)
t = 1.2261, df = 87, p-value = 0.2235
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.08010378  0.32963699
sample estimates:
      cor 
0.1303271 
cor.test(log(body_traits_avg$wing), log(body_traits_avg$weight)) # r = 0.226

    Pearson's product-moment correlation

data:  log(body_traits_avg$wing) and log(body_traits_avg$weight)
t = 2.1671, df = 87, p-value = 0.03296
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.01894248 0.41500412
sample estimates:
      cor 
0.2263065 
cor.test(log(body_traits_avg$culmen), log(body_traits_avg$weight)) # r = 0.090

    Pearson's product-moment correlation

data:  log(body_traits_avg$culmen) and log(body_traits_avg$weight)
t = 0.84253, df = 87, p-value = 0.4018
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1205525  0.2927350
sample estimates:
       cor 
0.08996281 

Use scaledMassIndex() function [defined in chunk above] to derive the SMI for each body mass measure (see Eq 2 of Peig and Green (2009))

smi_ <- 
  scaledMassIndex(x = body_traits_avg %>% 
                    ungroup() %>% 
                    dplyr::select(weight, wing) %>% 
                    na.omit() %>% 
                    pull(wing), 
                  y = body_traits_avg %>% 
                    ungroup() %>% 
                    dplyr::select(weight, wing) %>% 
                    na.omit() %>% 
                    pull(weight))

body_traits_smi <-
  body_traits_avg  %>% 
  filter(!is.na(weight)) %>% 
  bind_cols(., smi_) %>% 
  dplyr::select(ring, tag_type, tag_type_, treatment_group, measure, SMI.ols, weight, day_of_year, date_diff) %>% 
  rename(smi_wing = SMI.ols)

inspect relationship between raw body mass and scaled mass index (lines connect repeated measures)

Scaled Mass Index sequential measures model

sample size summary

# A tibble: 3 × 2
  tag_type     n
  <chr>    <int>
1 GPS         10
2 PTT         10
3 control     49

summary of attachment mass

summary of attachment weights
expressed as relative percentage of scaled mass index or raw body mass
1.2 g GPS (%) 1.8-2.0 g PTT (%) only leg bands (%)
raw body mass
only tracking device mean 2.0 3.0 NA
median 2.0 3.0 NA
min 1.9 2.8 NA
max 2.1 3.5 NA
sd 0.0 0.1 NA
all attachments mean 2.8 3.8 0.8
median 2.8 3.8 0.8
min 2.6 3.6 0.7
max 3.0 4.3 0.9
sd 0.1 0.2 0.0
scaled mass index
only tracking device mean 2.0 3.1 NA
median 2.0 3.0 NA
min 1.8 2.7 NA
max 2.3 3.6 NA
sd 0.1 0.3 NA
all attachments mean 2.9 3.9 0.8
median 2.9 3.8 0.8
min 2.6 3.4 0.6
max 3.2 4.4 0.9
sd 0.2 0.3 0.0

sample size summary

# A tibble: 3 × 2
  tag_type     n
  <chr>    <int>
1 GPS         10
2 PTT         10
3 control     49

distribution of time differences between repeated measures

hist(as.numeric(body_traits_smi$date_diff))

body_traits_smi %>% filter(!is.na(date_diff)) %>% pull(date_diff) %>% min()
[1] 2
body_traits_smi %>% filter(!is.na(date_diff)) %>% pull(date_diff) %>% max()
[1] 1097

modelling

smi_wing_lmer_DOY_season <- lme4::lmer(
  smi_wing ~ tag_type * measure + scale(day_of_year) + (1 | ring),
  data = body_traits_smi, REML = TRUE
)

contrast interactive and non-interactive model

smi_wing_lmer_DOY_season_interaction <- 
  lmer(smi_wing ~ tag_type * measure + scale(day_of_year) + (1 | ring), 
       data = body_traits_smi, REML = FALSE)

smi_wing_lmer_DOY_season_noninteraction <- 
  lmer(smi_wing ~ tag_type + measure + scale(day_of_year) + (1 | ring), 
       data = body_traits_smi, REML = FALSE)
Likelihood Ratio Test for Model Comparison
term npar AIC BIC logLik minus2logL statistic df p.value
smi_wing_lmer_DOY_season_noninteraction 8 514.554 534.464 -249.277 498.554 NA NA NA
smi_wing_lmer_DOY_season_interaction 12 521.687 551.551 -248.843 497.687 0.868 4 0.929

statistics

# set seed for reproducibility
set.seed(1234)

# Derive confidence intervals of effect sizes from parametric bootstrapping
tidy_smi_wing_lmer_DOY_season <-
  tidy(smi_wing_lmer_DOY_season, conf.int = TRUE, conf.method = "boot", nsim = 1000)

# run rptR to obtain repeatabilities of random effects
rpt_smi_wing_lmer_DOY_season <-
  rpt(smi_wing ~ tag_type * measure + scale(day_of_year) + (1 | ring),
      grname = c("ring", "Fixed"),
      data = body_traits_smi,
      datatype = "Gaussian",
      nboot = 1000, npermut = 1000, ratio = TRUE,
      adjusted = TRUE, ncores = 4, parallel = TRUE)

# run partR2 on each model to obtain marginal R2, parameter estimates, and beta
# weights
R2m_smi_wing_lmer_DOY_season <-
  partR2(smi_wing_lmer_DOY_season,
         partvars = c("tag_type",
                      "measure",
                      "scale(day_of_year)"),
         R2_type = "marginal",
         nboot = 1000,
         CI = 0.95,
         max_level = 1)

R2c_smi_wing_lmer_DOY_season <-
  partR2(smi_wing_lmer_DOY_season,
         partvars = c("tag_type",
                      "measure",
                      "scale(day_of_year)"),
         R2_type = "conditional",
         nboot = 1000,
         CI = 0.95,
         max_level = 1)

stats_smi_wing_lmer_DOY_season <-
  list(mod = smi_wing_lmer_DOY_season,
       tidy = tidy_smi_wing_lmer_DOY_season,
       rptR = rpt_smi_wing_lmer_DOY_season,
       partR2m = R2m_smi_wing_lmer_DOY_season,
       partR2c = R2c_smi_wing_lmer_DOY_season,
       data = body_traits_smi)
# saveRDS(stats_smi_wing_lmer_DOY_season, here("out/smi_date_results.rds"))
Scaled mass index (SMI) dynamics Mean estimate 95% confidence interval
Fixed effects 𝛽 (standardized; units = ĝ)
Intercept (control group 1st capture) 59.03 [57.72, 60.31]
Julian date of measure 0.02 [-0.19, 0.24]
GPS group 1st capture −1.28 [-4.66, 1.74]
PTT group 1st capture 0.19 [-3.31, 3.64]
Control group 2nd capture (within season) −2.28 [-4.42, 0.15]
Control group 2nd capture (between seasons) −3.13 [-7.72, 0.99]
GPS group*2nd capture (within season) 0.18 [-4.88, 5.27]
PTT group*2nd capture (within season) 1.89 [-4, 7.26]
GPS group*2nd capture (between seasons) 0.66 [-7.41, 9.16]
PTT group*2nd capture (between seasons) 1.09 [-3.87, 6.69]
Partitioned 𝑹²
Total Marginal 𝑹² 0.05 [0.04, 0.22]
Tracking device type 0.00 [0, 0.17]
Measurement sequence 0.00 [0, 0.17]
Julian date of measure 0.00 [0, 0.17]
Total Conditional 𝑹² 0.73 [0.5, 0.91]
Random effects 𝜎²
Individual 14.95 [6.53, 23.37]
Residual 5.96 [2.4, 12.02]
Adjusted repeatability 𝑟
Individual 0.71 [0.42, 0.89]
Residual 0.29 [0.11, 0.58]
Sample sizes 𝑛
Individuals 69
Observations 89
gtsave(smi_wing_lmer_table, here("tabs_figs/table_1.rtf"))

Breeding status and outcome

breed_fates <-
  read.csv(here("data/breed.csv")) 

Breeding data summary

number of distinct nests
n_nests
117
, , n_seasons = 1

      tag_type
parent GPS PTT
  RBBR   0   0
  RBBY   1   0
  RBGB   0   0
  RBOG   0   0
  RBOY   0   0
  RBRL   0   0
  RBWR   0   0
  RLYL   1   0
  RLYR   1   0
  RRBR   0   0
  RRGO   0   0
  RRLR   1   0
  RROW   0   1
  RRRG   0   1
  RRYY   0   0
  RWBG   0   0
  RWBY   0   0
  RWLB   0   0
  RWLG   0   0

, , n_seasons = 2

      tag_type
parent GPS PTT
  RBBR   1   0
  RBBY   0   0
  RBGB   0   1
  RBOG   0   1
  RBOY   0   1
  RBRL   0   1
  RBWR   0   1
  RLYL   0   0
  RLYR   0   0
  RRBR   1   0
  RRGO   0   1
  RRLR   0   0
  RROW   0   0
  RRRG   0   0
  RRYY   1   0
  RWBG   1   0
  RWBY   1   0
  RWLB   1   0
  RWLG   0   1
   nest_id parent tag_type season  state       date
1 2021_N22   RBBR      GPS   2021   fail 2021-10-09
2 2022_N05   RBBR      GPS   2022   fail 2022-09-01
3 2022_N15   RBBR      GPS   2022   fail 2022-09-19
4 2022_N39   RBBR  control   2022 fledge 2022-11-12
5 2023_N02   RBBR  control   2023   fail 2023-08-28
6 2023_N10   RBBR  control   2023   fail 2023-09-17
7 2023_N30   RBBR  control   2023   fail 2023-10-31
8 2023_N44   RBBR  control   2023   fail 2023-12-07
proportion of nests associated with each fate
state n_nests
aband 0.0256410
brood 0.1709402
fail 0.7435897
fledge 0.0598291
number of distinct nests for each fate
state n_nests
aband 3
brood 20
fail 87
fledge 7
number of distinct nests for each treatment
tag_type n_nests
GPS 22
PTT 26
control 88
number of distinct birds for each treatment
tag_type n_ind
GPS 10
PTT 9
control 42
number of distinct birds for each treatment
tag_type n_ind
GPS 10
PTT 9
control 42

check if breeding fates of subsequent seasons with the treatment are different than the deployment season

# Determine first treatment season per parent and tag_type
first_treatment_season <- breed_fates %>%
  filter(tag_type %in% c("GPS", "PTT")) %>%
  group_by(parent, tag_type) %>%
  summarise(first_treatment = min(season), .groups = "drop")

# Join back to main data and classify treatment stage
breed_fates_ <- breed_fates %>%
  left_join(first_treatment_season, by = c("parent", "tag_type")) %>%
  mutate(
    treatment_stage = case_when(
      tag_type == "control" ~ "control",
      season == first_treatment ~ "first_treatment_season",
      season > first_treatment ~ "subsequent_treatment_season",
      TRUE ~ NA_character_
    )
  ) %>% 
  mutate(treatment_stage = relevel(factor(treatment_stage), 
                                   ref = "control"))
# set seed for reproducibility
set.seed(1234)

breed_brm2 <- brm(
  state ~ treatment_stage + (1 | season) + (1 | parent),
  data = breed_fates_,
  family = categorical(link = "logit"),
  chains = 4,
  iter = 4000,
  control = list(adapt_delta = 0.99, max_treedepth = 15)
)

# consolidate results
stats_breed_brm_ <-
  list(mod = breed_brm2,
       data = breed_fates_)

modelling

# set seed for reproducibility
set.seed(1234)

# Fit the mixed-effects multinomial logistic regression model
breed_brm <- brm(
  state ~ tag_type + (1 | season) + (1 | parent),
  data = breed_fates,
  family = categorical(link = "logit"),
  chains = 4,
  iter = 4000,
  control = list(adapt_delta = 0.99, max_treedepth = 15)  
)

# consolidate results
stats_breed_brm <-
  list(mod = breed_brm,
       data = breed_fates)

Apparent survival

RMark data preparation

ch <- read.csv(here("data/ch.csv"))

# Specify the states and create the dataframe for RMark
ch_processed <- process.data(ch, model = "Multistrata", groups = NULL)

ch_ddl <- make.design.data(ch_processed)

# Define parameter formulas
S.model <- list(formula = ~ stratum)  # Survival depends on the state
p.model <- list(formula = ~ stratum)  # Detection probability depends on the state
Psi.model.stratum <- list(formula = ~ -1 + stratum)  
Psi.model.stratumto <- list(formula = ~ -1 + stratum:tostratum) 
Psi.model.stratumtotime <- list(formula = ~ -1 + stratum : tostratum : time)  # Transition probabilities depend on the state and destination state

RMark modelling

model.list <- create.model.list("Multistrata") 
surv.results <- mark.wrapper(model.list, data = ch_processed, ddl = ch_ddl,
                             threads = 4, brief = TRUE, delete = TRUE)

Figure 2 - multi-panelled plot of all results

Behavioural Observations

data summary

summary of behavioural observations
tagged n_obs
untagged 76
tagged 148

energy budget estimation

Table 5: Morrier and McNeil (1991). Time-activity budgets of Wilson's and Semipalmated Plovers in a tropical environment
avg_prop_Feeding 12.0340219
avg_prop_Preening 1.1922180
avg_prop_Resting 35.7261638
avg_prop_Flight 1.1573963
avg_prop_Walking 0.6282573
avg_prop_Fight 0.0553076
avg_prop_Ground_Display 0.1514821
avg_prop_Alert 2.0594441
avg_prop_Running 0.1839188
avg_prop_Nocturnal_Resting 46.8075087

modelling

inc_brood_model <- 
  glmer(inc_brood ~ tagged + sin_time + cos_time + (1 | bird), 
        data = behav, family = binomial)
Fixed Effects of the GLMM on incubation & brooding behaviour
effect term estimate std.error statistic p.value
fixed (Intercept) 43.355 15.721 2.758 0.006
fixed taggedY 0.030 0.488 0.061 0.952
fixed sin_time -3.004 1.254 -2.396 0.017
fixed cos_time -6.464 2.217 -2.916 0.004
walk_model <- 
  glmer(walk ~ tagged + sin_time + cos_time + (1 | bird), 
        data = behav, family = binomial)
Fixed Effects of the GLMM on walking behaviour
effect term estimate std.error statistic p.value
fixed (Intercept) -20.247 22.829 -0.887 0.375
fixed taggedY 0.691 0.783 0.883 0.377
fixed sin_time 1.586 1.963 0.808 0.419
fixed cos_time 3.336 3.164 1.054 0.292
feed_model <- 
  glmer(feed ~ tagged + sin_time + cos_time + (1 | bird), 
        data = behav, family = binomial)
Fixed Effects of the GLMM on feeding behaviour
effect term estimate std.error statistic p.value
fixed (Intercept) 14.710 12.465 1.180 0.238
fixed taggedY -0.398 0.301 -1.324 0.186
fixed sin_time -1.373 1.023 -1.342 0.180
fixed cos_time -2.066 1.746 -1.183 0.237
fly_model <- 
  glmer(fly ~ tagged + sin_time + cos_time + (1 | bird), 
        data = behav, family = binomial)
Fixed Effects of the GLMM on flying behaviour
effect term estimate std.error statistic p.value
fixed (Intercept) 1.664 17.726 0.094 0.925
fixed taggedY 0.946 0.738 1.282 0.200
fixed sin_time 0.107 1.491 0.072 0.943
fixed cos_time -0.931 2.473 -0.376 0.707

Movements and tag operation

summary of deployment duration (weeks)
tag_type mean max min
GPS 53.92857 60.71429 50.5714286
PTT 42.05357 57.14286 0.1428571
summary of fixes obtained per week
tag_type GPS PTT
mean 0.6310789 11.3666695
sd 0.1959389 9.4714222
min 0.389518 2.061186
max 0.9462716 35.3164070
summary of total number of fixes acquired
tag_type total_points
GPS 104
PTT 5544
summary of operational duration (weeks)
tag_type mean sd min max
GPS 26.99767 11.04417 13.98582 39.71011
PTT 44.57389 30.31189 3.88126 91.51782

plot of tag data collected over time

Supplementary Analysis

Breeding status and outcome (pairs)

summary of pair breeding fate data

number of mates and tag-type treatments over the three-year period

breed_fates_wide %>% 
  filter(parent_2 != "unkn") %>% 
  group_by(parent_1) %>% 
  summarise(n_mates = n_distinct(parent_2),
            n_combos = n_distinct(tag_type_combo)) %>% 
  arrange(desc(n_mates))
# A tibble: 21 × 3
   parent_1 n_mates n_combos
   <chr>      <int>    <int>
 1 RBBY           2        1
 2 RBOG           2        1
 3 RBOY           2        1
 4 RBGB           1        1
 5 RBRL           1        1
 6 RBWR           1        1
 7 RLYL           1        2
 8 RRBB           1        1
 9 RRBG           1        1
10 RRBR           1        1
# ℹ 11 more rows

number distinct nesting attempts with different treatment combinations

breed_fates_wide %>% 
  group_by(tag_type_combo) %>% 
  summarise(n_nests = n_distinct(nest_id))
# A tibble: 5 × 2
  tag_type_combo  n_nests
  <chr>             <int>
1 GPS-GPS               5
2 GPS-control          17
3 PTT-PTT               8
4 PTT-control          18
5 control-control      69

proportion of nests of different fates

breed_fates_wide %>%
  dplyr::select(nest_id, state) %>% 
  distinct() %>%
  group_by(state) %>% 
  summarise(n_nests = n()/117)
# A tibble: 4 × 2
  state  n_nests
  <chr>    <dbl>
1 aband   0.0256
2 brood   0.171 
3 fail    0.744 
4 fledge  0.0598

modelling

breed_pair_brm_ <- brm(
  state ~ tag_type_combo + (1 | season),
  data = breed_fates_wide,
  family = categorical(link = "logit"),
  chains = 4,
  iter = 4000,
  control = list(adapt_delta = 0.99, max_treedepth = 15)
)

stats_breed_pair_brm <-
  list(mod = breed_pair_brm_,
       data = breed_fates_wide)

supplementary figure 1