S3. Model fitting workflow

Location–scale and hierarchical Bayesian modeling pipeline

This chapter documents the model-fitting workflow used in the study. It includes Gaussian location–scale models for bloom phenology, as well as hierarchical Bayesian models for phenological mismatch and reproductive performance. Most code chunks are included as a reproducible workflow record and are not executed during rendering of the supplementary website.

How to read this chapter

This page is intended as a transparent workflow record rather than a results chapter. Executable code is shown primarily to document data preparation, model specification, and automation logic used in the analyses reported elsewhere in the website.

S3.1 Computational setup

Threading can accelerate computation by distributing sampling across CPU cores. This requires the cmdstanr package and a working CmdStan installation.

Show CmdStan installation and checks

options(mc.cores = parallel::detectCores(), brms.backend = "cmdstanr")

Threading can accelerate computation by distributing sampling across CPU cores. This requires the cmdstanr package and a working CmdStan installation.

cmdstanr::install_cmdstan()
cmdstanr::check_cmdstan_toolchain()
cmdstanr::cmdstan_version()

S3.2 Location–scale models of bloom phenology

These models were used to test whether both the mean timing and variability of bloom onset changed over time. Bloom onset was modeled as a Gaussian response, with calendar year included in both the location (3bc) and scale (3c3) submodels.

Variable definitions

delta5: bloom onset date (days since winter solstice) at the 5% threshold
delta15: bloom onset date (days since winter solstice) at the 15% threshold
delta25: bloom onset date (days since winter solstice) at the 25% threshold
TIME: calendar year, centered but not scaled

Data transformation

Bloom onset dates were converted to days elapsed since the previous winter solstice.

Show bloom data transformation

dat<-read.csv(here("..","data","bloom_dates.csv"))
dat <- dat %>%
  mutate(
    # Convert dates
    bloom_start_thr5_as_date = as.Date(bloom_start_thr5, format = "%d/%m/%Y"),
    bloom_start_thr15_as_date = as.Date(bloom_start_thr15, format = "%d/%m/%Y"),
    bloom_start_thr25_as_date = as.Date(bloom_start_thr25, format = "%d/%m/%Y"),
    winter_start_date = as.Date(paste0(season_year - 1, "-12-21")),

    # Calculate days since solstice
    delta5 = as.numeric(difftime(bloom_start_thr5_as_date, winter_start_date, units = "days")),
    delta15 = as.numeric(difftime(bloom_start_thr15_as_date, winter_start_date, units = "days")),
    delta25 = as.numeric(difftime(bloom_start_thr25_as_date, winter_start_date, units = "days")),

    # Factor year
    TIME = as.numeric(scale(season_year, center = TRUE, scale = FALSE)))%>%
  select( -starts_with("bloom_start_thr"), -winter_start_date)

Wrapper function for location–scale models

A wrapper function was used to fit one Gaussian location–scale model per bloom threshold, with TIME included in both the mean and residual-variance components.

Show location–scale model function

run_location_scale_models <- function(response_vars, data) {
  if (!dir.exists(here("Rdata"))) dir.create(here("Rdata"))
  
  purrr::walk(response_vars, function(var) {
    model_formula <- bf(
      as.formula(paste(var, "~ 1 + TIME")),
      sigma ~ 1 + TIME
    )
    
    file_path <- here("Rdata", paste0("model_", var, ".rds"))
    model_family <- gaussian()
    model_priors <- get_prior(model_formula, data = data, family = model_family)
    
    fit <- tryCatch({
      brm(
        formula = model_formula,
        data = data,
        prior = model_priors,
        family = model_family,
        iter = 6000, warmup = 2000, chains = 4, seed = 123, init = 0,
        control = list(adapt_delta = 0.99, max_treedepth = 15),
        file = file_path,
        file_refit = "on_change"
      )
    }, error = function(e) {
      cat("Error fitting model for:", var, "\n")
      cat("Error message:", e$message, "\n")
      return(NULL)
    })
    
    if (!is.null(fit)) cat("Saved model for:", var, "to", file_path, "\n")
  })
}

Model execution

Show location–scale model execution

# Bloom onset (delta variables)
delta_vars <- c("delta5", "delta15", "delta25")
run_location_scale_models(response_vars = delta_vars, data = dat)

S3.3 Hierarchical Bayesian models of reproductive performance

The primary breeding dataset was prepared separately for females and males prior to fitting Gaussian mismatch models and cumulative ordinal models for fledging success.

Preprocessing overview

Key preprocessing steps included:

conversion of focal predictors to numeric
encoding FLEDS as an ordered factor for cumulative ordinal models
centering, but not scaling, continuous predictors
conversion of YEAR and RING to factors for use as grouping variables

Variable definitions

X5MISMATCH: phenological mismatch (days) using the 5% bloom threshold
FLEDS: number of fledged chicks (ordinal response: 0–3)
TIME: calendar year, centered but not scaled
ZX5MISMATCH: centered phenological mismatch (5% threshold)
ABSMATCH: centered absolute mismatch (5% threshold)
AGE: age, centered
AFR: age at first reproduction, centered
LONGEVITY: age at last observation (proxy for lifespan), centered
YEAR: calendar year as a factor
RING: individual identifier as a factor

Show female data preparation

dat <- read.csv(here("..","data", "Lay_date_mismatches_250702.csv")) %>%
  subset(SEX == 1) %>%
  mutate(
    # First, ensure columns that will be used for scaling are numeric
    across(c(X5MISMATCH, AGE, AFR, LONGEVITY), as.numeric),
    
    # Second, create new scaled variables from existing columns
    TIME = scale(YEAR, center = TRUE, scale = FALSE)[,1],
    ZX5MISMATCH = scale(X5MISMATCH, center = TRUE, scale = FALSE)[,1],
    ABSMATCH = scale(abs(X5MISMATCH), center =TRUE, scale =FALSE)[,1],
    # Third, scale the other numeric predictors
    AGE = scale(AGE, center = TRUE, scale = FALSE)[,1],
    AFR = scale(AFR, center = TRUE, scale = FALSE)[,1],
    LONGEVITY = scale(LONGEVITY, center = TRUE, scale = FALSE)[,1],
    
    # Finally, convert grouping and outcome variables to factors
    YEAR = as.factor(YEAR),
    RING = as.factor(RING),
    FLEDS = factor(FLEDS, levels = 0:3, ordered = TRUE),
    
  )

Show male data preparation

dat0 <- read.csv(here("..","data", "Lay_date_mismatches_250702.csv")) %>%
  subset(SEX ==0) %>%
  mutate(
    across(c(X5MISMATCH, AGE, AFR, LONGEVITY), as.numeric),
    
    TIME = scale(YEAR, center = TRUE, scale = FALSE)[,1],
    ZX5MISMATCH = scale(X5MISMATCH, center = TRUE, scale = FALSE)[,1],
    
    AGE = scale(AGE, center = TRUE, scale = FALSE)[,1],
    AFR = scale(AFR, center = TRUE, scale = FALSE)[,1],
    LONGEVITY = scale(LONGEVITY, center = TRUE, scale = FALSE)[,1],
    
    YEAR = as.factor(YEAR),
    RING = as.factor(RING),
    FLEDS = factor(FLEDS, levels = 0:3, ordered = TRUE),
    
  )

Models syntax

Two model families were used:

Gaussian models for continuous responses such as phenological mismatch
cumulative ordinal probit models for ordered categorical responses such as number of fledglings

S3.4 Automated model fitting workflow

Candidate models were fitted systematically using a custom function, run_analysis_workflow(), which automated model specification, storage, and leave-one-out cross-validation.

What the workflow does

The automated workflow:

iterates across candidate fixed- and random-effects structures
supports both Gaussian and ordinal response families
optionally models auxiliary parameters (sigma for Gaussian models and disc for ordinal models)
skips models that have already been fitted
saves fitted model objects and corresponding LOO results
uses tryCatch() so that individual model failures do not halt the full pipeline

Fixed effects

Show fixed-effects definitions

FIXED_EFFECTS_ORDINAL_FULL  <- "ZX5MISMATCH + TIME + AGE + I(AGE^2) + AFR + LONGEVITY"

FIXED_EFFECTS_ORDINAL_QUADMATCH  <- "ZX5MISMATCH + I(ZX5MISMATCH^2) + TIME + AGE + I(AGE^2) + AFR + LONGEVITY"

FIXED_EFFECTS_ORDINAL_ABSMATCH  <- "ABSMATCH + TIME + AGE + I(AGE^2) + AFR + LONGEVITY"


FIXED_EFFECTS_GAUSSIAN      <- "AGE + I(AGE^2) + TIME + AFR + LONGEVITY"

Master workflow function

Show automated fitting workflow code

run_analysis_workflow <- function(analysis_name, response_var, fixed_effects, model_family, data, model_subset_key=NULL) {
  
  message(paste0("\n\n\n###--- Starting Analysis: '", analysis_name, "' ---###"))
  
  output_dir <- here("Rdata", analysis_name)
  dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  
  is_gaussian <- model_family$family == "gaussian"
  
  # --- Define Model Specifications based on analysis type ---
  model_specs <- if (is_gaussian) {
    tribble(
      ~name,                      ~mean_random,               ~aux_formula,
      "m1",            "(1|RING) + (1|YEAR)",      NA, 
      "m2",       "(1+AGE|RING) + (1|YEAR)",  NA,
      "m3",      "(1+TIME|RING) + (1|YEAR)", NA,
      "m4",       "(1|q|RING) + (1|YEAR)",    ~ 1 + (1|q|RING),
      "m5",      "(1+AGE|q|RING) + (1|YEAR)",    ~ 1 + (1|q|RING),
      "m6",     "(1+TIME|q|RING) + (1|YEAR)",   ~ 1 + (1|q|RING),
      "m7",    "(1|q|RING) + (1|YEAR)",    ~ 1 + AGE + TIME + (1|q|RING),
      "m8",   "(1+AGE|q|RING) + (1|YEAR)",    ~ 1 + AGE + TIME + (1|q|RING),
      "m9",  "(1+TIME|q|RING) + (1|YEAR)",   ~ 1 + AGE + TIME + (1|q|RING),
      "m10",      "(1|q|RING) + (1|YEAR)",    ~ 1 + AGE + I(AGE^2) + TIME + (1|q|RING),
      "m11",     "(1+AGE|q|RING) + (1|YEAR)",    ~ 1 + AGE + I(AGE^2) + TIME + (1|q|RING),
      "m12",    "(1+TIME|q|RING) + (1|YEAR)",   ~ 1 + AGE + I(AGE^2) + TIME + (1|q|RING)
    )
  } else {
    tribble(
      ~name,                          ~mean_random,                          ~aux_formula,
    #LINEAR MISMATCH  
      "m1",                "(1|RING) + (1|YEAR)",                 ~ 1,
      "m2",           "(1+AGE|RING) + (1|YEAR)",             ~ 1,
      "m3",          "(1+TIME|RING) + (1|YEAR)",            ~ 1,
      "m4",      "(1+ZX5MISMATCH|RING) + (1|YEAR)",     ~ 1,
      "m5",            "(1|q|RING) + (1|YEAR)",               ~ 1 + (1|q|RING),
      "m6",           "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + (1|q|RING),
      "m7",          "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + (1|q|RING),
      "m8",      "(1+ZX5MISMATCH|q|RING) + (1|YEAR)",   ~ 1 + (1|q|RING),
      "m9",         "(1|q|RING) + (1|YEAR)",               ~ 1 + ZX5MISMATCH + AGE + TIME + (1|q|RING),
      "m10",        "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + ZX5MISMATCH + AGE + TIME + (1|q|RING),
      "m11",       "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + ZX5MISMATCH + AGE + TIME + (1|q|RING),
      "m12",   "(1+ZX5MISMATCH|q|RING) + (1|YEAR)",   ~ 1 + ZX5MISMATCH + AGE + TIME + (1|q|RING),
      "m13",           "(1|q|RING) + (1|YEAR)",               ~ 1 + ZX5MISMATCH + AGE + I(AGE^2) + TIME + (1|q|RING),
      "m14",          "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + ZX5MISMATCH + AGE + I(AGE^2) + TIME + (1|q|RING),
      "m15",         "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + ZX5MISMATCH + AGE + I(AGE^2) + TIME + (1|q|RING),
      "m16",     "(1+ZX5MISMATCH|q|RING) + (1|YEAR)",   ~ 1 + ZX5MISMATCH + AGE + I(AGE^2) + TIME + (1|q|RING),
      
      
# QUAD MISMATCH      
      "Q1",                "(1|RING) + (1|YEAR)",                 ~ 1,
      "Q2",           "(1+AGE|RING) + (1|YEAR)",             ~ 1,
      "Q3",          "(1+TIME|RING) + (1|YEAR)",            ~ 1,
      "Q4",      "(1+ZX5MISMATCH|RING) + (1|YEAR)",     ~ 1,
      "Q5",            "(1|q|RING) + (1|YEAR)",               ~ 1 + (1|q|RING),
      "Q6",           "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + (1|q|RING),
      "Q7",          "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + (1|q|RING),
      "Q8",      "(1+ZX5MISMATCH|q|RING) + (1|YEAR)",   ~ 1 + (1|q|RING),
       "Q9",         "(1|q|RING) + (1|YEAR)",               ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) + AGE + TIME + (1|q|RING),
      "Q10",        "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + TIME + (1|q|RING),
      "Q11",       "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + TIME + (1|q|RING),
      "Q12",   "(1+ZX5MISMATCH|q|RING) + (1|YEAR)",   ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + TIME + (1|q|RING),
      "Q13",           "(1|q|RING) + (1|YEAR)",               ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + I(AGE^2) + TIME + (1|q|RING),
      "Q14",          "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + I(AGE^2) + TIME + (1|q|RING),
      "Q15",         "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + I(AGE^2) + TIME + (1|q|RING),
      "Q16",     "(1+ZX5MISMATCH|q|RING) + (1|YEAR)",   ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + I(AGE^2) + TIME + (1|q|RING),

       "Q17",         "(1|q|RING) + (1|YEAR)",               ~ 1 + ZX5MISMATCH +AGE + TIME + (1|q|RING),
      "Q18",        "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + ZX5MISMATCH +AGE + TIME + (1|q|RING),
      "Q19",       "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + ZX5MISMATCH  +AGE + TIME + (1|q|RING),
      "Q20",   "(1+ZX5MISMATCH|q|RING) + (1|YEAR)",   ~ 1 + ZX5MISMATCH  +AGE + TIME + (1|q|RING),
      "Q21",           "(1|q|RING) + (1|YEAR)",               ~ 1 + ZX5MISMATCH  +AGE + I(AGE^2) + TIME + (1|q|RING),
      "Q22",          "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + ZX5MISMATCH  +AGE + I(AGE^2) + TIME + (1|q|RING),
      "Q23",         "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + ZX5MISMATCH  +AGE + I(AGE^2) + TIME + (1|q|RING),
      "Q24",     "(1+ZX5MISMATCH|q|RING) + (1|YEAR)",   ~ 1 + ZX5MISMATCH  +AGE + I(AGE^2) + TIME + (1|q|RING),
      
      
      
      
#ABS MISMATCH      
      "ABS1",                "(1|RING) + (1|YEAR)",                 ~ 1,
      "ABS2",           "(1+AGE|RING) + (1|YEAR)",             ~ 1,
      "ABS3",          "(1+TIME|RING) + (1|YEAR)",            ~ 1,
      "ABS4",      "(1+ABSMATCH|RING) + (1|YEAR)",     ~ 1,
      "ABS5",            "(1|q|RING) + (1|YEAR)",               ~ 1 + (1|q|RING),
      "ABS6",           "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + (1|q|RING),
      "ABS7",          "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + (1|q|RING),
      "ABS8",      "(1+ABSMATCH|q|RING) + (1|YEAR)",   ~ 1 + (1|q|RING),
       "ABS9",         "(1|q|RING) + (1|YEAR)",               ~ 1 + ABSMATCH + AGE + TIME + (1|q|RING),
      "ABS10",        "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + ABSMATCH+ AGE + TIME + (1|q|RING),
      "ABS11",       "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + ABSMATCH + AGE + TIME + (1|q|RING),
      "ABS12",   "(1+ABSMATCH|q|RING) + (1|YEAR)",   ~ 1 + ABSMATCH+ AGE + TIME + (1|q|RING),
      "ABS13",           "(1|q|RING) + (1|YEAR)",               ~ 1 + ABSMATCH + AGE + I(AGE^2) + TIME + (1|q|RING),
      "ABS14",          "(1+AGE|q|RING) + (1|YEAR)",           ~ 1 + ABSMATCH + AGE + I(AGE^2) + TIME + (1|q|RING),
      "ABS15",         "(1+TIME|q|RING) + (1|YEAR)",          ~ 1 + ABSMATCH + AGE + I(AGE^2) + TIME + (1|q|RING),
      "ABS16",     "(1+ABSMATCH|q|RING) + (1|YEAR)",   ~ 1 + ABSMATCH+ AGE + I(AGE^2) + TIME + (1|q|RING),
      
      
      
      
      
      
      
      
      
      
      
      
    )
  }
  
if (!is.null(model_subset_key)) {
  message(paste0("Filtering models for subset: '", model_subset_key, "'"))
  
  if (model_subset_key == "linear") {
   # Keep only models that DON'T start with "Q_" or "ABS_"
   model_specs <- model_specs %>%
    filter(!grepl("^Q", name) & !grepl("^ABS", name))
  } else if (model_subset_key == "quad") {
   # Keep only models that DO start with "Q_"
   model_specs <- model_specs %>%
    filter(grepl("^Q", name))
  } else if (model_subset_key == "abs") {
   # Keep only models that DO start with "ABS_"
   model_specs <- model_specs %>%
    filter(grepl("^ABS", name))
  }
  
  # Safety check
  if (nrow(model_specs) == 0) {
   warning(paste0("No models found for subset key: '", model_subset_key, "'. Check your model names and key."))
   return(invisible(NULL)) # Stop if no models to run
  }
 }
  
  
  # --- Internal function to fit a single model ---
  fit_and_save_model <- function(name, mean_random, aux_formula, ...) {
    model_path <- file.path(output_dir, paste0("female_fit_", name, ".rds"))
    loo_path <- file.path(output_dir, paste0("female_loo_", name, ".rds"))
    
    if (file.exists(model_path)) {
      message(paste0("Skipping '", name, "': Model file already exists."))
      return(invisible(NULL))
    }
    
    message(paste0("\n--- Fitting model: ", name, " ---"))
    
    formula_str <- paste(response_var, "~", fixed_effects, "+", mean_random)
    
    # --- CORRECTED LOGIC ---
    # Use inherits() for a robust check if aux_formula is a real formula or NA.
    brm_formula <- if (!inherits(aux_formula, "formula")) {
      bf(formula_str, nl = FALSE)
    } else if (is_gaussian) {
      bf(formula_str, sigma = aux_formula, nl = FALSE)
    } else {
      bf(formula_str, disc = aux_formula, nl = FALSE)
    }
    
    model_priors <- get_prior(brm_formula, data = data, family = model_family)
    
    fit <- tryCatch({
      brm(
        formula = brm_formula, data = data, prior = model_priors, family = model_family,
        iter = 6000, warmup = 2000, chains = 4, seed = 123, init = 0,backend = "cmdstanr",threads = threading(12),
        control = list(adapt_delta = 0.99, max_treedepth = 15)
      )
    }, error = function(e) { message(paste0("ERROR fitting '", name, "': ", e$message)); NULL })
    
    if (!is.null(fit)) {
      saveRDS(fit, model_path)
      loo_result <- tryCatch(loo(fit), error = function(e) { message(paste0("ERROR with LOO: ", e$message)); NULL })
      if (!is.null(loo_result)) saveRDS(loo_result, loo_path)
    }
  }
  
  # --- Run all models and compare LOO results ---
  pwalk(model_specs, fit_and_save_model)
  

}

S3.5 Executing the workflow

The workflow was run separately for males (dat0) and females (dat) for two classes of response variables:

phenological mismatch as a Gaussian response
fledging success (FLEDS) as an ordinal cumulative probit response

Show Gaussian mismatch workflow for males

run_analysis_workflow(
  analysis_name = "pheno_mismatch",
  response_var = "X5MISMATCH",
  fixed_effects = FIXED_EFFECTS_GAUSSIAN,
  model_family = gaussian(),
  data = dat0
)

Here, the response variable (FLEDS) is an ordered categorical outcome (0, 1, 2, or 3 fledglings), so models were fitted using a cumulative probit family

Show ordinal fledging workflow for males

run_analysis_workflow(
  analysis_name = "fleds_full",
  response_var = "FLEDS",
  fixed_effects = FIXED_EFFECTS_ORDINAL_FULL,
  model_family = cumulative(link = "probit"),
  data = dat0,
  model_subset_key = "linear"
)

run_analysis_workflow(
 analysis_name = "fleds_full",
 response_var = "FLEDS",
 fixed_effects = FIXED_EFFECTS_ORDINAL_QUADMATCH,
 model_family = cumulative(link = "probit"),
 data = dat0,
 model_subset_key = "quad" 
)

run_analysis_workflow(
 analysis_name = "fleds_full",
 response_var = "FLEDS",
 fixed_effects = FIXED_EFFECTS_ORDINAL_ABSMATCH ,
 model_family = cumulative(link = "probit"),
 data = dat0,
 model_subset_key = "abs" 
)

Show Gaussian mismatch workflow for females

run_analysis_workflow(
  analysis_name = "pheno_mismatch",
  response_var = "X5MISMATCH",
  fixed_effects = FIXED_EFFECTS_GAUSSIAN,
  model_family = gaussian(),
  data = dat
)

Show ordinal fledging workflow for females

run_analysis_workflow(
  analysis_name = "fleds_full",
  response_var = "FLEDS",
  fixed_effects = FIXED_EFFECTS_ORDINAL_FULL,
  model_family = cumulative(link = "probit"),
  data = dat,
  model_subset_key = "linear"
)

run_analysis_workflow(
 analysis_name = "fleds_full",
 response_var = "FLEDS",
 fixed_effects = FIXED_EFFECTS_ORDINAL_QUADMATCH,
 model_family = cumulative(link = "probit"),
 data = dat,
 model_subset_key = "quad" 
)

run_analysis_workflow(
 analysis_name = "fleds_full",
 response_var = "FLEDS",
 fixed_effects = FIXED_EFFECTS_ORDINAL_ABSMATCH ,
 model_family = cumulative(link = "probit"),
 data = dat,
 model_subset_key = "abs" 
)

Reproducibility

Session information

R version 4.5.2 (2025-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
  LAPACK version 3.12.1

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

time zone: America/Edmonton
tzcode source: internal

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

other attached packages:
 [1] DT_0.34.0       loo_2.9.0       brms_2.23.0     Rcpp_1.1.1     
 [5] here_1.0.2      lubridate_1.9.4 forcats_1.0.1   stringr_1.6.0  
 [9] dplyr_1.1.4     purrr_1.2.1     readr_2.1.6     tidyr_1.3.2    
[13] tibble_3.3.1    ggplot2_4.0.1   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] gtable_0.3.6          tensorA_0.36.2.1      bslib_0.10.0         
 [4] xfun_0.56             htmlwidgets_1.6.4     lattice_0.22-7       
 [7] tzdb_0.5.0            crosstalk_1.2.2       vctrs_0.7.1          
[10] tools_4.5.2           generics_0.1.4        parallel_4.5.2       
[13] sandwich_3.1-1        pkgconfig_2.0.3       Matrix_1.7-4         
[16] checkmate_2.3.3       RColorBrewer_1.1-3    S7_0.2.1             
[19] distributional_0.6.0  RcppParallel_5.1.11-1 lifecycle_1.0.5      
[22] compiler_4.5.2        farver_2.1.2          Brobdingnag_1.2-9    
[25] codetools_0.2-20      sass_0.4.10           htmltools_0.5.9      
[28] bayesplot_1.15.0      yaml_2.3.12           jquerylib_0.1.4      
[31] pillar_1.11.1         MASS_7.3-65           cachem_1.1.0         
[34] bridgesampling_1.2-1  abind_1.4-8           multcomp_1.4-29      
[37] nlme_3.1-168          posterior_1.6.1       tidyselect_1.2.1     
[40] digest_0.6.39         mvtnorm_1.3-3         stringi_1.8.7        
[43] splines_4.5.2         rprojroot_2.1.1       fastmap_1.2.0        
[46] grid_4.5.2            cli_3.6.5             magrittr_2.0.4       
[49] survival_3.8-6        TH.data_1.1-5         withr_3.0.2          
[52] scales_1.4.0          backports_1.5.0       timechange_0.3.0     
[55] estimability_1.5.1    rmarkdown_2.30        matrixStats_1.5.0    
[58] emmeans_2.0.1         otel_0.2.0            zoo_1.8-15           
[61] hms_1.1.4             coda_0.19-4.1         evaluate_1.0.5       
[64] knitr_1.51            rstantools_2.6.0      rlang_1.1.7          
[67] xtable_1.8-4          glue_1.8.0            rstudioapi_0.18.0    
[70] jsonlite_2.0.0        R6_2.6.1

--- title: "S3. Model fitting workflow" subtitle: "Location–scale and hierarchical Bayesian modeling pipeline" unnumbered: true toc: true toc-depth: 3 --- ::: lead This chapter documents the model-fitting workflow used in the study. It includes Gaussian location–scale models for bloom phenology, as well as hierarchical Bayesian models for phenological mismatch and reproductive performance. Most code chunks are included as a reproducible workflow record and are not executed during rendering of the supplementary website. ::: ::: {.callout-note appearance="simple"} ## How to read this chapter This page is intended as a transparent workflow record rather than a results chapter. Executable code is shown primarily to document data preparation, model specification, and automation logic used in the analyses reported elsewhere in the website. ::: ```{r} #| label: setup-model-fitting #| include: false # For data manipulation (loads dplyr, tidyr, purrr, tibble, etc.) library(tidyverse) # For file path management library(here) # The core package for Bayesian hierarchical modeling library(brms) # For model comparison and cross-validation (LOO) library(loo) # For creating interactive HTML tables library(DT) ``` ## S3.1 Computational setup Threading can accelerate computation by distributing sampling across CPU cores. This requires the cmdstanr package and a working CmdStan installation. ```{r} #| label: cmdstan-installation #| eval: false #| code-fold: true #| code-summary: "Show CmdStan installation and checks" #| options(mc.cores = parallel::detectCores(), brms.backend = "cmdstanr") ``` ::: {.callout-tip appearance="simple" icon="false"} Threading can accelerate computation by distributing sampling across CPU cores. This requires the cmdstanr package and a working CmdStan installation. ```{r} #| code-fold: FALSE #| eval: false #| label: Set up CmdStan and confirm proper installation cmdstanr::install_cmdstan() cmdstanr::check_cmdstan_toolchain() cmdstanr::cmdstan_version() ``` ::: ## S3.2 Location–scale models of bloom phenology These models were used to test whether both the mean timing and variability of bloom onset changed over time. Bloom onset was modeled as a Gaussian response, with calendar year included in both the location (\u03bc) and scale (\u03c3) submodels. ::: {.callout-note title="Variable definitions"} - delta5: bloom onset date (days since winter solstice) at the 5% threshold - delta15: bloom onset date (days since winter solstice) at the 15% threshold - delta25: bloom onset date (days since winter solstice) at the 25% threshold - TIME: calendar year, centered but not scaled ::: ### Data transformation Bloom onset dates were converted to days elapsed since the previous winter solstice. ```{r} #| label: transform-bloom-data #| eval: false #| code-fold: true #| code-summary: "Show bloom data transformation" dat<-read.csv(here("..","data","bloom_dates.csv")) dat <- dat %>% mutate( # Convert dates bloom_start_thr5_as_date = as.Date(bloom_start_thr5, format = "%d/%m/%Y"), bloom_start_thr15_as_date = as.Date(bloom_start_thr15, format = "%d/%m/%Y"), bloom_start_thr25_as_date = as.Date(bloom_start_thr25, format = "%d/%m/%Y"), winter_start_date = as.Date(paste0(season_year - 1, "-12-21")), # Calculate days since solstice delta5 = as.numeric(difftime(bloom_start_thr5_as_date, winter_start_date, units = "days")), delta15 = as.numeric(difftime(bloom_start_thr15_as_date, winter_start_date, units = "days")), delta25 = as.numeric(difftime(bloom_start_thr25_as_date, winter_start_date, units = "days")), # Factor year TIME = as.numeric(scale(season_year, center = TRUE, scale = FALSE)))%>% select( -starts_with("bloom_start_thr"), -winter_start_date) ``` ### Wrapper function for location–scale models A wrapper function was used to fit one Gaussian location–scale model per bloom threshold, with TIME included in both the mean and residual-variance components. ```{r} #| label: bloom-location-scale-function #| eval: false #| code-fold: true #| code-summary: "Show location–scale model function" #| run_location_scale_models <- function(response_vars, data) { if (!dir.exists(here("Rdata"))) dir.create(here("Rdata")) purrr::walk(response_vars, function(var) { model_formula <- bf( as.formula(paste(var, "~ 1 + TIME")), sigma ~ 1 + TIME ) file_path <- here("Rdata", paste0("model_", var, ".rds")) model_family <- gaussian() model_priors <- get_prior(model_formula, data = data, family = model_family) fit <- tryCatch({ brm( formula = model_formula, data = data, prior = model_priors, family = model_family, iter = 6000, warmup = 2000, chains = 4, seed = 123, init = 0, control = list(adapt_delta = 0.99, max_treedepth = 15), file = file_path, file_refit = "on_change" ) }, error = function(e) { cat("Error fitting model for:", var, "\n") cat("Error message:", e$message, "\n") return(NULL) }) if (!is.null(fit)) cat("Saved model for:", var, "to", file_path, "\n") }) } ``` ### Model execution ```{r} #| label: run-bloom-location-scale-models #| eval: false #| code-fold: true #| code-summary: "Show location–scale model execution" # Bloom onset (delta variables) delta_vars <- c("delta5", "delta15", "delta25") run_location_scale_models(response_vars = delta_vars, data = dat) ``` ## S3.3 Hierarchical Bayesian models of reproductive performance The primary breeding dataset was prepared separately for females and males prior to fitting Gaussian mismatch models and cumulative ordinal models for fledging success. ### Preprocessing overview Key preprocessing steps included: - conversion of focal predictors to numeric - encoding FLEDS as an ordered factor for cumulative ordinal models - centering, but not scaling, continuous predictors - conversion of YEAR and RING to factors for use as grouping variables ::: {.callout-note title="Variable definitions"} - X5MISMATCH: phenological mismatch (days) using the 5% bloom threshold - FLEDS: number of fledged chicks (ordinal response: 0–3) - TIME: calendar year, centered but not scaled - ZX5MISMATCH: centered phenological mismatch (5% threshold) - ABSMATCH: centered absolute mismatch (5% threshold) - AGE: age, centered - AFR: age at first reproduction, centered - LONGEVITY: age at last observation (proxy for lifespan), centered - YEAR: calendar year as a factor - RING: individual identifier as a factor ::: ::: panel-tabset ## Female data ```{r} #| label: prepare-female-data #| eval: false #| code-fold: true #| code-summary: "Show female data preparation" dat <- read.csv(here("..","data", "Lay_date_mismatches_250702.csv")) %>% subset(SEX == 1) %>% mutate( # First, ensure columns that will be used for scaling are numeric across(c(X5MISMATCH, AGE, AFR, LONGEVITY), as.numeric), # Second, create new scaled variables from existing columns TIME = scale(YEAR, center = TRUE, scale = FALSE)[,1], ZX5MISMATCH = scale(X5MISMATCH, center = TRUE, scale = FALSE)[,1], ABSMATCH = scale(abs(X5MISMATCH), center =TRUE, scale =FALSE)[,1], # Third, scale the other numeric predictors AGE = scale(AGE, center = TRUE, scale = FALSE)[,1], AFR = scale(AFR, center = TRUE, scale = FALSE)[,1], LONGEVITY = scale(LONGEVITY, center = TRUE, scale = FALSE)[,1], # Finally, convert grouping and outcome variables to factors YEAR = as.factor(YEAR), RING = as.factor(RING), FLEDS = factor(FLEDS, levels = 0:3, ordered = TRUE), ) ``` ## Male data ```{r} #| label: prepare-male-data #| eval: false #| code-fold: true #| code-summary: "Show male data preparation" dat0 <- read.csv(here("..","data", "Lay_date_mismatches_250702.csv")) %>% subset(SEX ==0) %>% mutate( across(c(X5MISMATCH, AGE, AFR, LONGEVITY), as.numeric), TIME = scale(YEAR, center = TRUE, scale = FALSE)[,1], ZX5MISMATCH = scale(X5MISMATCH, center = TRUE, scale = FALSE)[,1], AGE = scale(AGE, center = TRUE, scale = FALSE)[,1], AFR = scale(AFR, center = TRUE, scale = FALSE)[,1], LONGEVITY = scale(LONGEVITY, center = TRUE, scale = FALSE)[,1], YEAR = as.factor(YEAR), RING = as.factor(RING), FLEDS = factor(FLEDS, levels = 0:3, ordered = TRUE), ) ``` ::: ### Models syntax Two model families were used: - Gaussian models for continuous responses such as phenological mismatch - cumulative ordinal probit models for ordered categorical responses such as number of fledglings ::: panel-tabset ## Gaussian models (Phenological mismatch) ```{r} #| label: Gaussian models #| echo: false dfg <- read.csv(here("summaries","gaussian.csv")) DT::datatable(dfg, options = list(pageLength = 5)) ``` ## Ordinal models (Number of fledglings) ```{r} #| label: Ordinal models #| echo: false dfo <- read.csv(here("summaries","ordinal.csv")) DT::datatable(dfo, options = list(pageLength = 5)) ``` ::: ## S3.4 Automated model fitting workflow Candidate models were fitted systematically using a custom function, run_analysis_workflow(), which automated model specification, storage, and leave-one-out cross-validation. ::: {.callout-note appearance="simple"} What the workflow does The automated workflow: - iterates across candidate fixed- and random-effects structures - supports both Gaussian and ordinal response families - optionally models auxiliary parameters (sigma for Gaussian models and disc for ordinal models) - skips models that have already been fitted - saves fitted model objects and corresponding LOO results - uses tryCatch() so that individual model failures do not halt the full pipeline ::: ### Fixed effects ```{r} #| label: fixed-effects-sets #| eval: false #| code-fold: true #| code-summary: "Show fixed-effects definitions" FIXED_EFFECTS_ORDINAL_FULL <- "ZX5MISMATCH + TIME + AGE + I(AGE^2) + AFR + LONGEVITY" FIXED_EFFECTS_ORDINAL_QUADMATCH <- "ZX5MISMATCH + I(ZX5MISMATCH^2) + TIME + AGE + I(AGE^2) + AFR + LONGEVITY" FIXED_EFFECTS_ORDINAL_ABSMATCH <- "ABSMATCH + TIME + AGE + I(AGE^2) + AFR + LONGEVITY" FIXED_EFFECTS_GAUSSIAN <- "AGE + I(AGE^2) + TIME + AFR + LONGEVITY" ``` ### Master workflow function ```{r} #| label: master-model-workflow #| eval: false #| code-fold: true #| code-summary: "Show automated fitting workflow code" run_analysis_workflow <- function(analysis_name, response_var, fixed_effects, model_family, data, model_subset_key=NULL) { message(paste0("\n\n\n###--- Starting Analysis: '", analysis_name, "' ---###")) output_dir <- here("Rdata", analysis_name) dir.create(output_dir, recursive = TRUE, showWarnings = FALSE) is_gaussian <- model_family$family == "gaussian" # --- Define Model Specifications based on analysis type --- model_specs <- if (is_gaussian) { tribble( ~name, ~mean_random, ~aux_formula, "m1", "(1|RING) + (1|YEAR)", NA, "m2", "(1+AGE|RING) + (1|YEAR)", NA, "m3", "(1+TIME|RING) + (1|YEAR)", NA, "m4", "(1|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "m5", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "m6", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "m7", "(1|q|RING) + (1|YEAR)", ~ 1 + AGE + TIME + (1|q|RING), "m8", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + AGE + TIME + (1|q|RING), "m9", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + AGE + TIME + (1|q|RING), "m10", "(1|q|RING) + (1|YEAR)", ~ 1 + AGE + I(AGE^2) + TIME + (1|q|RING), "m11", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + AGE + I(AGE^2) + TIME + (1|q|RING), "m12", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + AGE + I(AGE^2) + TIME + (1|q|RING) ) } else { tribble( ~name, ~mean_random, ~aux_formula, #LINEAR MISMATCH "m1", "(1|RING) + (1|YEAR)", ~ 1, "m2", "(1+AGE|RING) + (1|YEAR)", ~ 1, "m3", "(1+TIME|RING) + (1|YEAR)", ~ 1, "m4", "(1+ZX5MISMATCH|RING) + (1|YEAR)", ~ 1, "m5", "(1|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "m6", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "m7", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "m8", "(1+ZX5MISMATCH|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "m9", "(1|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + AGE + TIME + (1|q|RING), "m10", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + AGE + TIME + (1|q|RING), "m11", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + AGE + TIME + (1|q|RING), "m12", "(1+ZX5MISMATCH|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + AGE + TIME + (1|q|RING), "m13", "(1|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + AGE + I(AGE^2) + TIME + (1|q|RING), "m14", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + AGE + I(AGE^2) + TIME + (1|q|RING), "m15", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + AGE + I(AGE^2) + TIME + (1|q|RING), "m16", "(1+ZX5MISMATCH|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + AGE + I(AGE^2) + TIME + (1|q|RING), # QUAD MISMATCH "Q1", "(1|RING) + (1|YEAR)", ~ 1, "Q2", "(1+AGE|RING) + (1|YEAR)", ~ 1, "Q3", "(1+TIME|RING) + (1|YEAR)", ~ 1, "Q4", "(1+ZX5MISMATCH|RING) + (1|YEAR)", ~ 1, "Q5", "(1|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "Q6", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "Q7", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "Q8", "(1+ZX5MISMATCH|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "Q9", "(1|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) + AGE + TIME + (1|q|RING), "Q10", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + TIME + (1|q|RING), "Q11", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + TIME + (1|q|RING), "Q12", "(1+ZX5MISMATCH|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + TIME + (1|q|RING), "Q13", "(1|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + I(AGE^2) + TIME + (1|q|RING), "Q14", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + I(AGE^2) + TIME + (1|q|RING), "Q15", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + I(AGE^2) + TIME + (1|q|RING), "Q16", "(1+ZX5MISMATCH|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH + I(ZX5MISMATCH^2) +AGE + I(AGE^2) + TIME + (1|q|RING), "Q17", "(1|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH +AGE + TIME + (1|q|RING), "Q18", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH +AGE + TIME + (1|q|RING), "Q19", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH +AGE + TIME + (1|q|RING), "Q20", "(1+ZX5MISMATCH|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH +AGE + TIME + (1|q|RING), "Q21", "(1|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH +AGE + I(AGE^2) + TIME + (1|q|RING), "Q22", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH +AGE + I(AGE^2) + TIME + (1|q|RING), "Q23", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH +AGE + I(AGE^2) + TIME + (1|q|RING), "Q24", "(1+ZX5MISMATCH|q|RING) + (1|YEAR)", ~ 1 + ZX5MISMATCH +AGE + I(AGE^2) + TIME + (1|q|RING), #ABS MISMATCH "ABS1", "(1|RING) + (1|YEAR)", ~ 1, "ABS2", "(1+AGE|RING) + (1|YEAR)", ~ 1, "ABS3", "(1+TIME|RING) + (1|YEAR)", ~ 1, "ABS4", "(1+ABSMATCH|RING) + (1|YEAR)", ~ 1, "ABS5", "(1|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "ABS6", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "ABS7", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "ABS8", "(1+ABSMATCH|q|RING) + (1|YEAR)", ~ 1 + (1|q|RING), "ABS9", "(1|q|RING) + (1|YEAR)", ~ 1 + ABSMATCH + AGE + TIME + (1|q|RING), "ABS10", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + ABSMATCH+ AGE + TIME + (1|q|RING), "ABS11", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + ABSMATCH + AGE + TIME + (1|q|RING), "ABS12", "(1+ABSMATCH|q|RING) + (1|YEAR)", ~ 1 + ABSMATCH+ AGE + TIME + (1|q|RING), "ABS13", "(1|q|RING) + (1|YEAR)", ~ 1 + ABSMATCH + AGE + I(AGE^2) + TIME + (1|q|RING), "ABS14", "(1+AGE|q|RING) + (1|YEAR)", ~ 1 + ABSMATCH + AGE + I(AGE^2) + TIME + (1|q|RING), "ABS15", "(1+TIME|q|RING) + (1|YEAR)", ~ 1 + ABSMATCH + AGE + I(AGE^2) + TIME + (1|q|RING), "ABS16", "(1+ABSMATCH|q|RING) + (1|YEAR)", ~ 1 + ABSMATCH+ AGE + I(AGE^2) + TIME + (1|q|RING), ) } if (!is.null(model_subset_key)) { message(paste0("Filtering models for subset: '", model_subset_key, "'")) if (model_subset_key == "linear") { # Keep only models that DON'T start with "Q_" or "ABS_" model_specs <- model_specs %>% filter(!grepl("^Q", name) & !grepl("^ABS", name)) } else if (model_subset_key == "quad") { # Keep only models that DO start with "Q_" model_specs <- model_specs %>% filter(grepl("^Q", name)) } else if (model_subset_key == "abs") { # Keep only models that DO start with "ABS_" model_specs <- model_specs %>% filter(grepl("^ABS", name)) } # Safety check if (nrow(model_specs) == 0) { warning(paste0("No models found for subset key: '", model_subset_key, "'. Check your model names and key.")) return(invisible(NULL)) # Stop if no models to run } } # --- Internal function to fit a single model --- fit_and_save_model <- function(name, mean_random, aux_formula, ...) { model_path <- file.path(output_dir, paste0("female_fit_", name, ".rds")) loo_path <- file.path(output_dir, paste0("female_loo_", name, ".rds")) if (file.exists(model_path)) { message(paste0("Skipping '", name, "': Model file already exists.")) return(invisible(NULL)) } message(paste0("\n--- Fitting model: ", name, " ---")) formula_str <- paste(response_var, "~", fixed_effects, "+", mean_random) # --- CORRECTED LOGIC --- # Use inherits() for a robust check if aux_formula is a real formula or NA. brm_formula <- if (!inherits(aux_formula, "formula")) { bf(formula_str, nl = FALSE) } else if (is_gaussian) { bf(formula_str, sigma = aux_formula, nl = FALSE) } else { bf(formula_str, disc = aux_formula, nl = FALSE) } model_priors <- get_prior(brm_formula, data = data, family = model_family) fit <- tryCatch({ brm( formula = brm_formula, data = data, prior = model_priors, family = model_family, iter = 6000, warmup = 2000, chains = 4, seed = 123, init = 0,backend = "cmdstanr",threads = threading(12), control = list(adapt_delta = 0.99, max_treedepth = 15) ) }, error = function(e) { message(paste0("ERROR fitting '", name, "': ", e$message)); NULL }) if (!is.null(fit)) { saveRDS(fit, model_path) loo_result <- tryCatch(loo(fit), error = function(e) { message(paste0("ERROR with LOO: ", e$message)); NULL }) if (!is.null(loo_result)) saveRDS(loo_result, loo_path) } } # --- Run all models and compare LOO results --- pwalk(model_specs, fit_and_save_model) } ``` ## S3.5 Executing the workflow The workflow was run separately for males (dat0) and females (dat) for two classes of response variables: - phenological mismatch as a Gaussian response - fledging success (FLEDS) as an ordinal cumulative probit response ::::: panel-tabset ## Analyses for males ::: panel-tabset ## Phenological mismatch Analysis ```{r} #| label: run-male-mismatch-models #| eval: false #| code-fold: true #| code-summary: "Show Gaussian mismatch workflow for males" run_analysis_workflow( analysis_name = "pheno_mismatch", response_var = "X5MISMATCH", fixed_effects = FIXED_EFFECTS_GAUSSIAN, model_family = gaussian(), data = dat0 ) ``` ## Number of fledglings Analysis Here, the response variable (`FLEDS`) is an ordered categorical outcome (0, 1, 2, or 3 fledglings), so models were fitted using a cumulative probit family ```{r} #| label: run-male-fledging-models #| eval: false #| code-fold: true #| code-summary: "Show ordinal fledging workflow for males" run_analysis_workflow( analysis_name = "fleds_full", response_var = "FLEDS", fixed_effects = FIXED_EFFECTS_ORDINAL_FULL, model_family = cumulative(link = "probit"), data = dat0, model_subset_key = "linear" ) run_analysis_workflow( analysis_name = "fleds_full", response_var = "FLEDS", fixed_effects = FIXED_EFFECTS_ORDINAL_QUADMATCH, model_family = cumulative(link = "probit"), data = dat0, model_subset_key = "quad" ) run_analysis_workflow( analysis_name = "fleds_full", response_var = "FLEDS", fixed_effects = FIXED_EFFECTS_ORDINAL_ABSMATCH , model_family = cumulative(link = "probit"), data = dat0, model_subset_key = "abs" ) ``` ::: ## Analyses for females ::: panel-tabset ## Phenological mismatch Analysis ```{r} #| label: run-female-mismatch-models #| eval: false #| code-fold: true #| code-summary: "Show Gaussian mismatch workflow for females" run_analysis_workflow( analysis_name = "pheno_mismatch", response_var = "X5MISMATCH", fixed_effects = FIXED_EFFECTS_GAUSSIAN, model_family = gaussian(), data = dat ) ``` ## Number of fledglings Analysis ```{r} #| label: run-female-fledging-models #| eval: false #| code-fold: true #| code-summary: "Show ordinal fledging workflow for females" run_analysis_workflow( analysis_name = "fleds_full", response_var = "FLEDS", fixed_effects = FIXED_EFFECTS_ORDINAL_FULL, model_family = cumulative(link = "probit"), data = dat, model_subset_key = "linear" ) run_analysis_workflow( analysis_name = "fleds_full", response_var = "FLEDS", fixed_effects = FIXED_EFFECTS_ORDINAL_QUADMATCH, model_family = cumulative(link = "probit"), data = dat, model_subset_key = "quad" ) run_analysis_workflow( analysis_name = "fleds_full", response_var = "FLEDS", fixed_effects = FIXED_EFFECTS_ORDINAL_ABSMATCH , model_family = cumulative(link = "probit"), data = dat, model_subset_key = "abs" ) ``` ::: ::::: ## Reproducibility ::: {.callout-tip appearance="simple" icon="false"} ## Session information ```{r} #| echo: false #| warning: false #| message: false utils::sessionInfo() ``` :::