Statistical analysis

Emerald glassfrog mating activity

Author

Published

February 22, 2026

Source code and data found at https://github.com/maRce10/emerald_glassfrog_mating_activity

Code

# options to customize chunk outputs
knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE
 )

Purpose

Infer the effects of environmental variables on emerald glassfrog mating activity

Analysis flowchart

flowchart LR
  A[Read data] --> B(Format data) 
  B --> C(Define DAGs)
  C --> D(Statistical analysis)
  D --> E(Model summary) 

style A fill:#382A5433
style B fill:#395D9C33
style C fill:#3497A933
style D fill:#60CEAC33

Load packages

Code

# knitr is require for creating html/pdf/word reports
# formatR is used for soft-wrapping code

# install/ load packages
sketchy::load_packages(packages = c("ggplot2", "viridis", "dagitty", "ggdag", "lunar", "brms", "brmsish", "rstan", "ggraph", "ggparty"))

# set theme globally
  theme_set(theme_classic(base_size = 20))

source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")

1 Read data

Code

dat <- read.csv("./data/raw/males and pairs_all years_updated.csv")

dat$org.date <- dat$date

dat$date <- paste(dat$year,  dat$date, sep = "-")

# convert "8-Jun" to month-day format
dat$date <- as.Date(dat$date, format = "%Y-%d-%b")

#convert to "1970-01-02" format
dat$date <- as.Date(dat$date, origin = "1970-01-01")

## add moon
dat$moonlight <- lunar.illumination(dat$date, shift = -6)

Prepare data

Code

# scale continuous variables
dat$mean_temp_sc <- scale(dat$mean_temp)
dat$min_temp_sc <- scale(dat$min_temp)
dat$max_temp_sc <- scale(dat$max_temp)
dat$mean_rh_sc <- scale(dat$mean_rh)
dat$total_rainfall_sc <- scale(dat$total_rainfall)
dat$prev_rainfall_sc <- scale(dat$prev_rainfall)
dat$moonlight_sc <- scale(dat$moonlight)
dat$year <- as.factor(dat$year)
dat$num_males_sc <- scale(dat$num_males)

2 Covariance temperature

Code

cor(dat[, c("mean_temp", "max_temp", "min_temp", "mean_rh")])

           mean_temp   max_temp   min_temp    mean_rh
mean_temp  1.0000000  0.7305338  0.6893214 -0.6662954
max_temp   0.7305338  1.0000000  0.3513733 -0.5508679
min_temp   0.6893214  0.3513733  1.0000000 -0.2277868
mean_rh   -0.6662954 -0.5508679 -0.2277868  1.0000000

3 DAG + stats

3.1 Male activity

Code

dag_l <- dagify(
    total_rainfall_sc ~ evapotranspiration,
    prev_rainfall_sc ~ evapotranspiration,
    mean_temp_sc ~ total_rainfall_sc,
    mean_rh_sc ~ total_rainfall_sc,
    num_males ~ mean_rh_sc,
    num_males ~ moonlight_sc,
    mean_rh_sc ~ mean_temp_sc,
    mean_rh_sc ~ prev_rainfall_sc,
    mean_rh_sc ~ total_rainfall_sc,
    num_males ~ mean_temp_sc,
    num_males ~ prev_rainfall_sc,
    num_males ~ total_rainfall_sc,
    labels = c(
        num_males = "Active\nmales",
        mean_rh_sc = "Relative\nhumidity",
        total_rainfall_sc = "Rainfall",
        prev_rainfall_sc = "Previous\nrainfall",
        moonlight_sc = "Moonlight",
        mean_temp_sc = "Temperature",
        evapotranspiration = "Evapotrans-\npiration",
        latent = c("evapotrans-\npiration"),
        outcome = "num_males"
    )
)

set.seed(5)
tidy_dag <- tidy_dagitty(dag_l)
tidy_dag$data$type <- ifelse(is.na(tidy_dag$data$to), "outcome", "predictor")
tidy_dag$data$type[tidy_dag$data$name %in% c("evapotranspiration")] <- "latent"

ggplot(tidy_dag, aes(
    x = x,
    y = y,
    xend = xend,
    yend = yend
)) + scale_color_viridis_d(
    option = "G",
    begin = 0.2,
    end = 0.8,
    alpha = 0.5
) + geom_dag_edges_fan(
    edge_color = mako(10, alpha = 0.4)[2],
    arrow = grid::arrow(length = grid::unit(10, "pt"), type = "closed")
) + geom_dag_point(aes(color = type), size = 24) + 
    geom_dag_text(color = "black", aes(label = label, color = label)) +
    theme_dag() + 
    theme(legend.position = "bottom")

3.1.1 SEM model

Code

# ----------------------------
# Submodels with year random intercept
# ----------------------------

bf_total_rain <- bf(
  total_rainfall_sc ~ 1 + (1 | year)
)

bf_prev_rain <- bf(
  prev_rainfall_sc ~ 1 + (1 | year)
)

bf_temp <- bf(
  temp_sc ~ 1 + total_rainfall_sc + (1 | year)
)

bf_rh <- bf(
  mean_rh_sc ~ 1 +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    (1 | year)
)

bf_males <- bf(
  num_males ~ 
    mean_rh_sc +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    moonlight_sc +
    (1 | year),
  family = negbinomial()
)



# ----------------------------
# Weakly informative priors
# ----------------------------

priors <- c(

  # ----------------------------
  # Linear slopes (all z-scored)
  # ----------------------------

  prior(normal(0, 0.5), class = "b", resp = "meanrhsc"),
  prior(normal(0, 0.5), class = "b", resp = "tempsc"),
  prior(normal(0, 0.5), class = "b", resp = "nummales"),

    # Negative binomial dispersion

  prior(exponential(1), class = "shape", resp = "nummales")
)

dat$temp_sc <- dat$mean_temp_sc

# run a prior sampling model to check the priors
sem_prior_only <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    set_rescor(FALSE),
  data = dat,
  prior = priors,
  sample_prior = "only",
  iter = 2000,
  chains = 2,
  cores = 2,
  # backend = "cmdstanr",
  file = "male_sem_prior_only.RDS",
  file_refit = "on_change")

# check there are not flat priors
ps <- prior_summary(sem_prior_only)

any(grepl("flat", ps$prior))


# ----------------------------
# Fit full SEM
# ----------------------------

sem_model <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    set_rescor(FALSE),
  data = dat,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
    file = "./data/processed/male_sem_model.RDS",
    backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(sem_model, criterion = c("loo"))

dat$temp_sc <- dat$max_temp_sc

sem_model_max_temp <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    set_rescor(FALSE),
  data = dat,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
    file = "./data/processed/male_sem_model_max_temp.RDS",
    backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(sem_model_max_temp, criterion = c("loo"))


dat$temp_sc <- dat$min_temp_sc


sem_model_min_temp <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    set_rescor(FALSE),
  data = dat,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
    file = "./data/processed/male_sem_model_min_temp.RDS",
    backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(sem_model_min_temp, criterion = c("loo"))

3.1.2 Model selection

Code

mean_temp_males_mod <- readRDS("./data/processed/male_sem_model.RDS")
max_temp <- readRDS("./data/processed/male_sem_model_max_temp.RDS")
min_temp <- readRDS("./data/processed/male_sem_model_min_temp.RDS")


loo_diff <- loo::loo_compare(mean_temp_males_mod, max_temp, min_temp)

rows <- rownames(loo_diff)
loo_diff <- as.data.frame(loo_diff)

loo_diff$model <- rows

loo_diff

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic	model
mean_temp_males_mod	0.00000	0.000000	-1182.076	29.90885	38.20647	6.311699	2364.152	59.81770	mean_temp_males_mod
max_temp	-26.98459	8.477540	-1209.060	29.08594	37.45516	6.215772	2418.121	58.17189	max_temp
min_temp	-30.25239	9.519574	-1212.328	29.69509	38.43242	6.650247	2424.657	59.39019	min_temp

3.1.2.1 Best model

Code

extended_summary(
  mean_temp_males_mod,
  highlight = TRUE,
  trace.palette = mako,
  print.name = FALSE
)

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(totalrainfallsc = list(formula = total_rainfall_sc ~ 1 + (1 \| year), pforms = list(), pfix = list(), resp = “totalrainfallsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), prevrainfallsc = list(formula = prev_rainfall_sc ~ 1 + (1 \| year), pforms = list(), pfix = list(), resp = “prevrainfallsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), tempsc = list(formula = temp_sc ~ 1 + total_rainfall_sc + (1 \| year), pforms = list(), pfix = list(), resp = “tempsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), meanrhsc = list(formula = mean_rh_sc ~ 1 + temp_sc + prev_rainfall_sc + total_rainfall_sc + (1 \| year), pforms = list(), pfix = list(), resp = “meanrhsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nummales = list(formula = num_males ~ mean_rh_sc + temp_sc + prev_rainfall_sc + total_rainfall_sc + moonlight_sc + (1 \| year), pforms = list(), pfix = list(), family = list(family = “negbinomial”, link = “log”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “shape”), type = “int”, ybounds = c(0, Inf), closed = c(TRUE, NA), ad = c(“weights”, “subset”, “cens”, “trunc”, “rate”, “index”), specials = “sbi_log”, prior = function (dpar, link = “identity”, …) { if (dpar == “shape” && link == “identity”) { return(“inv_gamma(0.4, 0.3)”) } NULL }, link_shape = “log”), resp = “nummales”, mecor = TRUE))	()	b-normal(0, 0.5) b-normal(0, 0.5) b-normal(0, 0.5) Intercept-student_t(3, 0.1, 2.5) Intercept-student_t(3, 3.8, 2.5) Intercept-student_t(3, -0.4, 2.5) Intercept-student_t(3, 0, 2.5) Intercept-student_t(3, -0.5, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) shape-exponential(1) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5)	4000	4	1	2000	15 (0.0019%)	0	2330.063	2814.962	1242506459

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_totalrainfallsc_Intercept	0.022	-0.392	0.468	1.001	3781.668	3336.185
b_prevrainfallsc_Intercept	-0.019	-0.308	0.278	1	4177.449	3321.976
b_tempsc_Intercept	0.004	-0.603	0.639	1.004	2330.063	2814.962
b_meanrhsc_Intercept	0.015	-0.533	0.548	1	2471.338	3315.966
b_nummales_Intercept	3.727	3.189	4.234	1.002	2356.185	3036.242
b_tempsc_total_rainfall_sc	-0.444	-0.580	-0.305	1.001	12412.168	5850.807
b_meanrhsc_temp_sc	-0.584	-0.747	-0.418	1.001	8130.692	6256.263
b_meanrhsc_prev_rainfall_sc	-0.125	-0.260	0.007	1	11037.834	6062.898
b_meanrhsc_total_rainfall_sc	0.088	-0.054	0.232	1.001	9593.055	6373.230
b_nummales_mean_rh_sc	0.067	-0.020	0.154	1	7876.228	5932.116
b_nummales_temp_sc	-0.132	-0.224	-0.041	1	7104.219	5895.543
b_nummales_prev_rainfall_sc	-0.055	-0.118	0.009	1.001	9241.917	6099.592
b_nummales_total_rainfall_sc	0.003	-0.068	0.073	1.001	8856.459	5722.541
b_nummales_moonlight_sc	-0.021	-0.080	0.036	1	12094.706	5597.106

Code

es <- extended_summary(
  mean_temp_males_mod,
  highlight = TRUE,
  trace.palette = mako,
  remove.intercepts = TRUE,
  print.name = FALSE,
  beta.prefix = "b_nummales",
  trace = FALSE,
  fill = mako(10, alpha = 0.8)[7],
  gsub.pattern = c(
    "b_nummales_",
    "_sc",
    "mean_rh",
    "moonlight",
    "prev_rainfall",
    "temp",
    "total_rainfall"
  ),
  gsub.replacement = c(
    "",
    "",
    "Relative humidity",
    "Moonlight",
    "Previous rainfall",
    "Temperature",
    "Total rainfall"
  ),
  return = TRUE, 
  ylab = "Predictor"
)

es$plot

3.1.2.2 Evaluate “year” effect

Code

# ----------------------------
# Submodels with year random intercept
# ----------------------------

bf_total_rain <- bf(
  total_rainfall_sc ~ 1
)

bf_prev_rain <- bf(
  prev_rainfall_sc ~ 1
)

bf_temp <- bf(
  temp_sc ~ 1 + total_rainfall_sc
)

bf_rh <- bf(
  mean_rh_sc ~ 1 +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc
)

bf_males <- bf(
  num_males ~ 
    s(mean_rh_sc) +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    moonlight_sc,
  family = negbinomial()
)



# ----------------------------
# Weakly informative priors
# ----------------------------

priors <- c(

  # ----------------------------
  # Linear slopes (all z-scored)
  # ----------------------------

  prior(normal(0, 0.5), class = "b", resp = "meanrhsc"),
  prior(normal(0, 0.5), class = "b", resp = "tempsc"),
  prior(normal(0, 0.5), class = "b", resp = "nummales"),

  # ----------------------------
  # Smooth term (humidity → males)
  # ----------------------------

  prior(normal(0, 1), class = "sds", resp = "nummales"),

  # ----------------------------
  # Year random effects
  # ----------------------------

  # ----------------------------
  # Negative binomial dispersion
  # ----------------------------

  prior(exponential(1), class = "shape", resp = "nummales")
)

dat$temp_sc <- dat$mean_temp_sc



# ----------------------------
# Fit full SEM
# ----------------------------

sem_model_no_year <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    set_rescor(FALSE),
  data = dat,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
    file = "./data/processed/male_sem_model_no_year.RDS",
    backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(sem_model_no_year, criterion = c("loo"))

Code

no_year <- readRDS("./data/processed/male_sem_model_no_year.RDS")

loo_diff <- loo::loo_compare(mean_temp_males_mod, no_year)

rows <- rownames(loo_diff)
loo_diff <- as.data.frame(loo_diff)

loo_diff$model <- rows

loo_diff

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic	model
mean_temp_males_mod	0.00000	0.0000	-1182.076	29.90885	38.20647	6.311699	2364.152	59.8177	mean_temp_males_mod
no_year	-94.93883	11.4158	-1277.015	29.13650	26.99000	5.911663	2554.029	58.2730	no_year

Code

R2_full <- bayes_R2(mean_temp_males_mod, resp = "nummales",  summary = FALSE)

R2_fixed <- bayes_R2(
  mean_temp_males_mod,
  resp = "nummales",
  re_formula = NA,
   summary = FALSE
)
    
delta_R2 <- R2_full - R2_fixed

quantile(delta_R2, c(0.025, 0.5, 0.975))

     2.5%       50%     97.5% 
0.3976816 0.5806765 0.6840495

3.1.3 Parametized DAG

Code

edge_data <- tidy_dag$data
fixef_df <- as.data.frame(fixef(mean_temp_males_mod))
fixef_df$term <- rownames(fixef_df)

fixef_df <- fixef_df[grep(pattern = "Intercept", fixef_df$term, invert = TRUE), ]

split_terms <- strsplit(fixef_df$term, "_")

fixef_df$to <- sapply(split_terms, `[`, 1)
fixef_df$to <- gsub("sc", "_sc", fixef_df$to)
fixef_df$to <- gsub("mean", "mean_", fixef_df$to)
fixef_df$to <- gsub("num", "num_", fixef_df$to)
fixef_df$to <- gsub("temp_", "mean_temp_", fixef_df$to)

# split from "_" to end
fixef_df$from <- sapply(split_terms, function(x) paste(x[-1], collapse = "_"))
fixef_df$from <- gsub("temp_", "mean_temp_", fixef_df$from)

sem_estimates <- data.frame(
  from = fixef_df$from,
  to   = fixef_df$to,
  est  = fixef_df$Estimate,
  l95  = fixef_df$Q2.5,
  u95  = fixef_df$Q97.5,
  p    = 2 * pmin(
            pnorm(fixef_df$Estimate / fixef_df$Est.Error),
            1 - pnorm(fixef_df$Estimate / fixef_df$Est.Error)
         )
)


edge_data <- merge(edge_data, sem_estimates[, c("from", "to", "est", "l95", "u95")], by.x = c("name", "to"), by.y = c("from","to"), all.x = TRUE)

edge_data$sig <- ifelse(
  edge_data$l95 * edge_data$u95 > 0 & !is.na(edge_data$l95),
  "solid",
  "dashed"
)

edge_data$alpha <- ifelse(
  edge_data$l95 * edge_data$u95 > 0 & !is.na(edge_data$l95),
  1,
  0
)

edge_data$color <- ifelse(
  edge_data$l95 * edge_data$u95 > 0 & !is.na(edge_data$l95),
  "black",
  "gray"
)

edge_data$color[edge_data$name == "num_males"] <- "black"
edge_data$est <- ifelse(
!is.na(edge_data$l95),
  edge_data$est,
  0
)


tidy_dag_param <- tidy_dag

tidy_dag_param$data <- edge_data

ggplot(tidy_dag_param, aes(
    x = x,
    y = y,
    xend = xend,
    yend = yend,
    edge_linetype = sig,
    edge_alpha = alpha
)) + scale_color_viridis_d(
    option = "G",
    begin = 0.2,
    end = 0.8,
    alpha = 0.5
) + 
    geom_dag_edges_fan(aes(size = (est / max(abs(est), na.rm = TRUE))^2),
    arrow = grid::arrow(length = grid::unit(10, "pt"), type = "closed")
) +
    geom_dag_point(aes(color = type), 
                   size = 24) + 
    geom_dag_text(aes(label = label), color = c(edge_data$color, edge_data$color)[!duplicated(c(edge_data$name, edge_data$to))]) +
    scale_edge_alpha(range = c(0.2, 1), guide = NULL) +
    scale_edge_linetype_manual(values = c("solid" = "solid", "dashed" = "dashed"), guide = NULL) +
    theme_dag() + 
    theme(legend.position = "bottom") +
    guides(size = "none")

3.1.4 Graph raw data

Temperature

Code

temp_grid <- data.frame(
  mean_temp = seq(min(dat$mean_temp, na.rm = TRUE),
                  max(dat$mean_temp, na.rm = TRUE),
                  length.out = 100)
)
temp_grid$mean_rh_sc        <- 0
temp_grid$prev_rainfall_sc <- 0
temp_grid$total_rainfall_sc <- 0
temp_grid$moonlight_sc     <- 0


temp_center <- attr(dat$mean_temp_sc, "scaled:center")
temp_scale  <- attr(dat$mean_temp_sc, "scaled:scale")

temp_grid$temp_sc <- (temp_grid$mean_temp - temp_center) / temp_scale


pred <- fitted(
  mean_temp_males_mod,
  resp = "nummales",
  newdata = temp_grid,
  re_formula = NA,   # population-level effect
  summary = TRUE
)

temp_grid <- cbind(temp_grid, pred)


ggplot(data = dat, aes(x = mean_temp, y = num_males)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
  geom_line(data = temp_grid,
            aes(x = mean_temp, y = Estimate),
            linewidth = 1.2,
            color = mako(10, alpha = 0.8)[3]) +
  geom_ribbon(data = temp_grid,
              aes(x = mean_temp, y = Estimate, ymin = Q2.5, ymax = Q97.5),
              alpha = 0.1,
              fill = "black") +
       labs(x = "Mean temperature (°C)", y = "Number of males")

Code

# by year
temp_grid_year <- data.frame(
  mean_temp = seq(min(dat$mean_temp, na.rm = TRUE),
                  max(dat$mean_temp, na.rm = TRUE),
                  length.out = 100),
  year = unique(dat$year)
)
temp_grid_year$mean_rh_sc        <- 0
temp_grid_year$prev_rainfall_sc <- 0
temp_grid_year$total_rainfall_sc <- 0
temp_grid_year$moonlight_sc     <- 0


temp_center <- attr(dat$mean_temp_sc, "scaled:center")
temp_scale  <- attr(dat$mean_temp_sc, "scaled:scale")

temp_grid_year$temp_sc <- (temp_grid_year$mean_temp - temp_center) / temp_scale


pred <- fitted(
  mean_temp_males_mod,
  resp = "nummales",
  newdata = temp_grid_year,
  re_formula = NULL,   # <-- this is the key
  summary = TRUE
)

temp_grid_year <- cbind(temp_grid_year, pred)

ggplot(data = dat, aes(x = mean_temp, y = num_males)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
  geom_line(data = temp_grid_year,
            aes(x = mean_temp, y = Estimate),
            linewidth = 1.2,
            color = mako(10, alpha = 0.8)[3]) +
    facet_wrap(~year) +
  geom_ribbon(data = temp_grid_year,
              aes(x = mean_temp, y = Estimate, ymin = Q2.5, ymax = Q97.5),
              alpha = 0.1,
              fill = "black") +
    labs(x = "Mean temperature (°C)", y = "Number of males")

Previous rainfall

Code

ggplot(data = dat, aes(x = prev_rainfall, y = num_males)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
       labs(x = "Previous rainfall (mm)", y = "Number of males")

Total rainfall

Code

ggplot(data = dat, aes(x = total_rainfall, y = num_males)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
       labs(x = "Rainfall (mm)", y = "Number of males")

Moonlight

Code

ggplot(data = dat, aes(x = moonlight * 100, y = num_males)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
       labs(x = "Moonlight (%)", y = "Number of males")

Relative humidity

Code

ggplot(data = dat, aes(x = mean_rh, y = num_males)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
       labs(x = "Relative humidity (%)", y = "Number of males")

3.2 Female activity

Code

dag_l <- dagify(
    num_males_sc ~ mean_rh_sc,
    num_males_sc ~ moonlight_sc,
    num_males_sc ~ mean_temp_sc,
    num_males_sc ~ prev_rainfall_sc,
    num_males_sc ~ total_rainfall_sc,
    total_rainfall_sc ~ evapotranspiration,
    prev_rainfall_sc ~ evapotranspiration,
    mean_temp_sc ~ total_rainfall_sc,
    mean_rh_sc ~ total_rainfall_sc,
    num_pairs ~ mean_rh_sc,
    num_pairs ~ moonlight_sc,
    mean_rh_sc ~ mean_temp_sc,
    mean_rh_sc ~ prev_rainfall_sc,
    mean_rh_sc ~ total_rainfall_sc,
    num_pairs ~ mean_temp_sc,
    num_pairs ~ prev_rainfall_sc,
    num_pairs ~ total_rainfall_sc,
    num_pairs ~ num_males_sc,
    labels = c(
        num_pairs = "Active\nfemales",
        mean_rh_sc = "Relative\nhumidity",
        total_rainfall_sc = "Rainfall",
        prev_rainfall_sc = "Previous\nrainfall",
        moonlight_sc = "Moonlight",
        mean_temp_sc = "Temperature",
        num_males_sc = "Active\nmales",
        evapotranspiration = "Evapotrans-\npiration",
        latent = c("evapotrans-\npiration"),
        outcome = "num_pairs"
    )
)

set.seed(5)
tidy_dag_fem <- tidy_dagitty(dag_l)
tidy_dag_fem$data$type <- ifelse(is.na(tidy_dag_fem$data$to), "outcome", "predictor")
tidy_dag_fem$data$type[tidy_dag_fem$data$name %in% c("evapotranspiration")] <- "latent"

ggplot(tidy_dag_fem, aes(
    x = x,
    y = y,
    xend = xend,
    yend = yend
)) + scale_color_viridis_d(
    option = "G",
    begin = 0.2,
    end = 0.8,
    alpha = 0.5
) + geom_dag_edges_fan(
    edge_color = mako(10, alpha = 0.4)[2],
    arrow = grid::arrow(length = grid::unit(10, "pt"), type = "closed")
) + geom_dag_point(aes(color = type), size = 24) + 
    geom_dag_text(color = "black", aes(label = label, color = color)) +
    theme_dag() + 
    theme(legend.position = "bottom")

3.2.1 Full data set

3.2.1.1 SEM model

Code

# ----------------------------
# Submodels
# ----------------------------

bf_total_rain <- bf(
  total_rainfall_sc ~ 1 + (1 | year)
)

bf_prev_rain <- bf(
  prev_rainfall_sc ~ 1 + (1 | year)
)

bf_temp <- bf(
  temp_sc ~ 1 +
    total_rainfall_sc +
    (1 | year)
)

bf_rh <- bf(
  mean_rh_sc ~ 1 +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    (1 | year)
)

# mediator: active males (count)
bf_males <- bf(
  num_males ~ 
    mean_rh_sc +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    moonlight_sc +
    (1 | year),
  family = negbinomial()
)

# outcome: active females (count)
bf_pairs <- bf(
  num_pairs ~ 
    mean_rh_sc +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    moonlight_sc +
    num_males_sc +
    (1 | year),
  family = negbinomial()
)

# ----------------------------
# Weakly informative priors
# ----------------------------

priors <- c(

  # Linear slopes (all predictors z-scored)
  prior(normal(0, 0.5), class = "b", resp = "meanrhsc"),
  prior(normal(0, 0.5), class = "b", resp = "tempsc"),
  prior(normal(0, 0.5), class = "b", resp = "nummales"),
  prior(normal(0, 0.5), class = "b", resp = "numpairs"),

  # Negative binomial dispersion
  prior(exponential(1), class = "shape", resp = "nummales"),
  prior(exponential(1), class = "shape", resp = "numpairs")
)

# ----------------------------
# Prior predictive check
# ----------------------------
dat$temp_sc <- dat$mean_temp_sc

sem_prior_only <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    bf_pairs +
    set_rescor(FALSE),
  data = dat,
  prior = priors,
  sample_prior = "only",
  iter = 2000,
  chains = 2,
  cores = 2,
  backend = "cmdstanr",
    file = "female_sem_prior_only.RDS",
  file_refit = "on_change"
)

ps <- prior_summary(sem_prior_only)

any(grepl("flat", ps$prior))


# ----------------------------
# Fit full SEM
# ----------------------------

fem_sem_model <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    bf_pairs +
    set_rescor(FALSE),
  data = dat,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
  file = "./data/processed/female_sem_model.RDS",
  backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(fem_sem_model, criterion = c("loo"))

dat$temp_sc <- dat$max_temp_sc

fem_max_temp <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    bf_pairs +
    set_rescor(FALSE),
  data = dat,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
  file = "./data/processed/female_sem_model_max_temp.RDS",
  backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(fem_max_temp, criterion = c("loo"))

dat$temp_sc <- dat$min_temp_sc

fem_min_temp <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    bf_pairs +
    set_rescor(FALSE),
  data = dat,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
  file = "./data/processed/female_sem_model_min_temp.RDS",
  backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(fem_min_temp, criterion = c("loo"))

3.2.1.2 Model selection

Code

mean_temp_females_mod <- readRDS("./data/processed/female_sem_model.RDS")
max_temp <- readRDS("./data/processed/female_sem_model_max_temp.RDS")
min_temp <- readRDS("./data/processed/female_sem_model_min_temp.RDS")

loo_diff <- loo::loo_compare(mean_temp_females_mod, max_temp, min_temp)

rows <- rownames(loo_diff)
loo_diff <- as.data.frame(loo_diff)

loo_diff$model <- rows

loo_diff

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic	model
mean_temp_females_mod	0.00000	0.000000	-1344.815	33.34107	48.00978	7.811394	2689.630	66.68214	mean_temp_females_mod
max_temp	-27.08086	8.748104	-1371.896	32.82315	47.93526	8.279532	2743.792	65.64629	max_temp
min_temp	-30.37348	9.050054	-1375.189	33.05050	48.36581	8.172025	2750.377	66.10101	min_temp

3.2.1.3 Best model

Code

extended_summary(
  mean_temp_females_mod,
  highlight = TRUE, 
  trace.palette = mako, 
  print.name = FALSE
)

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(totalrainfallsc = list(formula = total_rainfall_sc ~ 1 + (1 \| year), pforms = list(), pfix = list(), resp = “totalrainfallsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), prevrainfallsc = list(formula = prev_rainfall_sc ~ 1 + (1 \| year), pforms = list(), pfix = list(), resp = “prevrainfallsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), tempsc = list(formula = temp_sc ~ 1 + total_rainfall_sc + (1 \| year), pforms = list(), pfix = list(), resp = “tempsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), meanrhsc = list(formula = mean_rh_sc ~ 1 + temp_sc + prev_rainfall_sc + total_rainfall_sc + (1 \| year), pforms = list(), pfix = list(), resp = “meanrhsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nummales = list(formula = num_males ~ mean_rh_sc + temp_sc + prev_rainfall_sc + total_rainfall_sc + moonlight_sc + (1 \| year), pforms = list(), pfix = list(), family = list(family = “negbinomial”, link = “log”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “shape”), type = “int”, ybounds = c(0, Inf), closed = c(TRUE, NA), ad = c(“weights”, “subset”, “cens”, “trunc”, “rate”, “index”), specials = “sbi_log”, prior = function (dpar, link = “identity”, …) { if (dpar == “shape” && link == “identity”) { return(“inv_gamma(0.4, 0.3)”) } NULL }, link_shape = “log”), resp = “nummales”, mecor = TRUE), numpairs = list(formula = num_pairs ~ mean_rh_sc + temp_sc + prev_rainfall_sc + total_rainfall_sc + moonlight_sc + num_males_sc + (1 \| year), pforms = list(), pfix = list(), family = list(family = “negbinomial”, link = “log”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “shape”), type = “int”, ybounds = c(0, Inf), closed = c(TRUE, NA), ad = c(“weights”, “subset”, “cens”, “trunc”, “rate”, “index”), specials = “sbi_log”, prior = function (dpar, link = “identity”, …) { if (dpar == “shape” && link == “identity”) { return(“inv_gamma(0.4, 0.3)”) } NULL }, link_shape = “log”), resp = “numpairs”, mecor = TRUE))	()	b-normal(0, 0.5) b-normal(0, 0.5) b-normal(0, 0.5) b-normal(0, 0.5) Intercept-student_t(3, 0.1, 2.5) Intercept-student_t(3, 3.7, 2.5) Intercept-student_t(3, 0.7, 2.5) Intercept-student_t(3, -0.4, 2.5) Intercept-student_t(3, 0.1, 2.5) Intercept-student_t(3, -0.5, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) shape-exponential(1) shape-exponential(1) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5)	4000	4	1	2000	2 (0.00025%)	0	2722.17	3228.553	115649959

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_totalrainfallsc_Intercept	-0.028	-0.531	0.464	1.001	3514.416	3245.842
b_prevrainfallsc_Intercept	-0.008	-0.280	0.261	1.002	5616.636	3961.449
b_tempsc_Intercept	0.004	-0.606	0.591	1.001	3169.978	3721.181
b_meanrhsc_Intercept	-0.015	-0.627	0.605	1.001	2959.416	3599.657
b_nummales_Intercept	3.733	3.185	4.273	1.001	2722.170	3228.553
b_numpairs_Intercept	0.808	-0.082	1.697	1.001	3050.907	3540.453
b_tempsc_total_rainfall_sc	-0.471	-0.621	-0.318	1.001	10837.245	5911.413
b_meanrhsc_temp_sc	-0.569	-0.734	-0.402	1	7514.632	6407.383
b_meanrhsc_prev_rainfall_sc	-0.114	-0.249	0.024	1.001	9047.581	6405.015
b_meanrhsc_total_rainfall_sc	0.087	-0.063	0.240	1	8710.609	6358.005
b_nummales_mean_rh_sc	0.068	-0.025	0.159	1	8102.973	5679.870
b_nummales_temp_sc	-0.143	-0.241	-0.046	1	7041.100	6299.009
b_nummales_prev_rainfall_sc	-0.054	-0.122	0.015	1	8955.078	6018.077
b_nummales_total_rainfall_sc	-0.001	-0.079	0.078	1.001	8389.213	5950.326
b_nummales_moonlight_sc	-0.015	-0.074	0.044	1.001	9619.170	6000.170
b_numpairs_mean_rh_sc	0.112	-0.134	0.359	1.001	8333.929	6476.236
b_numpairs_temp_sc	-0.187	-0.444	0.065	1	7227.806	5906.879
b_numpairs_prev_rainfall_sc	0.072	-0.091	0.243	1	8798.608	6266.832
b_numpairs_total_rainfall_sc	0.219	0.039	0.402	1.001	8757.315	5892.540
b_numpairs_moonlight_sc	-0.107	-0.261	0.044	1	9398.779	5895.316
b_numpairs_num_males_sc	0.748	0.424	1.076	1.001	6693.949	5080.851

Code

es <- extended_summary(
  mean_temp_females_mod,
  highlight = TRUE,
  trace.palette = mako,
  remove.intercepts = TRUE,
  print.name = FALSE,
  beta.prefix = "b_numpairs",
  trace = FALSE,
  fill = mako(10, alpha = 0.8)[7],
  gsub.pattern = c(
    "b_numpairs_",
    "_sc",
    "mean_rh",
    "moonlight",
    "prev_rainfall",
    "temp",
    "total_rainfall",
    "num_males"
  ),
  gsub.replacement = c(
    "",
    "",
    "Relative humidity",
    "Moonlight",
    "Previous rainfall",
    "Temperature",
    "Total rainfall",
    "Number of males"
  ),
  return = TRUE, 
  ylab = "Predictor"
)

es$plot

3.2.2 Excluding outlier

3.2.2.1 SEM model

Code

# ----------------------------
# Submodels
# ----------------------------

bf_total_rain <- bf(
  total_rainfall_sc ~ 1 + (1 | year)
)

bf_prev_rain <- bf(
  prev_rainfall_sc ~ 1 + (1 | year)
)

bf_temp <- bf(
  temp_sc ~ 1 +
    total_rainfall_sc +
    (1 | year)
)

bf_rh <- bf(
  mean_rh_sc ~ 1 +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    (1 | year)
)

# mediator: active males (count)
bf_males <- bf(
  num_males ~ 
    mean_rh_sc +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    moonlight_sc +
    (1 | year),
  family = negbinomial()
)

# outcome: active females (count)
bf_pairs <- bf(
  num_pairs ~ 
    mean_rh_sc +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    moonlight_sc +
    num_males_sc +
    (1 | year),
  family = negbinomial()
)

# ----------------------------
# Weakly informative priors
# ----------------------------

priors <- c(

  # Linear slopes (all predictors z-scored)
  prior(normal(0, 0.5), class = "b", resp = "meanrhsc"),
  prior(normal(0, 0.5), class = "b", resp = "tempsc"),
  prior(normal(0, 0.5), class = "b", resp = "nummales"),
  prior(normal(0, 0.5), class = "b", resp = "numpairs"),

  # Negative binomial dispersion
  prior(exponential(1), class = "shape", resp = "nummales"),
  prior(exponential(1), class = "shape", resp = "numpairs")
)


# exclude outlier
dat_no_outlier <- dat[dat$total_rainfall < 75, ]

# ----------------------------
# Prior predictive check
# ----------------------------
dat_no_outlier$temp_sc <- dat_no_outlier$mean_temp_sc

sem_prior_only <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    bf_pairs +
    set_rescor(FALSE),
  data = dat_no_outlier,
  prior = priors,
  sample_prior = "only",
  iter = 2000,
  chains = 2,
  cores = 2,
  backend = "cmdstanr",
    file = "female_sem_prior_only_no_outlier.RDS",
  file_refit = "on_change"
)

ps <- prior_summary(sem_prior_only)

any(grepl("flat", ps$prior))


# ----------------------------
# Fit full SEM
# ----------------------------

fem_sem_model <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    bf_pairs +
    set_rescor(FALSE),
  data = dat_no_outlier,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
  file = "./data/processed/female_sem_model_no_outlier.RDS",
  backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(fem_sem_model, criterion = c("loo"))

dat_no_outlier$temp_sc <- dat_no_outlier$max_temp_sc

fem_max_temp <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    bf_pairs +
    set_rescor(FALSE),
  data = dat_no_outlier,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
  file = "./data/processed/female_sem_model_max_temp_no_outlier.RDS",
  backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(fem_max_temp, criterion = c("loo"))

dat_no_outlier$temp_sc <- dat_no_outlier$min_temp_sc

fem_min_temp <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    bf_pairs +
    set_rescor(FALSE),
  data = dat_no_outlier,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
  file = "./data/processed/female_sem_model_min_temp_no_outlier.RDS",
  backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(fem_min_temp, criterion = c("loo"))

3.2.2.2 Model selection

Code

mean_temp_females_mod_no_outlier <- readRDS("./data/processed/female_sem_model_no_outlier.RDS")
max_temp_no_outlier <- readRDS("./data/processed/female_sem_model_max_temp_no_outlier.RDS")
min_temp_no_outlier <- readRDS("./data/processed/female_sem_model_min_temp_no_outlier.RDS")

loo_diff <- loo::loo_compare(mean_temp_females_mod_no_outlier, max_temp_no_outlier, min_temp_no_outlier)

rows <- rownames(loo_diff)
loo_diff <- as.data.frame(loo_diff)

loo_diff$model <- rows

loo_diff

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic	model
mean_temp_females_mod_no_outlier	0.00000	0.000000	-1310.734	26.81340	43.27293	4.460509	2621.468	53.62680	mean_temp_females_mod_no_outlier
max_temp_no_outlier	-26.89591	8.781583	-1337.630	26.01189	42.29543	4.274889	2675.260	52.02378	max_temp_no_outlier
min_temp_no_outlier	-29.64930	8.977484	-1340.383	26.13289	42.74375	4.237709	2680.767	52.26578	min_temp_no_outlier

3.2.2.3 Best model

Code

extended_summary(
  mean_temp_females_mod_no_outlier,
  highlight = TRUE, 
  trace.palette = mako, 
  print.name = FALSE
)

	formula	family	priors	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	list(totalrainfallsc = list(formula = total_rainfall_sc ~ 1 + (1 \| year), pforms = list(), pfix = list(), resp = “totalrainfallsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), prevrainfallsc = list(formula = prev_rainfall_sc ~ 1 + (1 \| year), pforms = list(), pfix = list(), resp = “prevrainfallsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), tempsc = list(formula = temp_sc ~ 1 + total_rainfall_sc + (1 \| year), pforms = list(), pfix = list(), resp = “tempsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), meanrhsc = list(formula = mean_rh_sc ~ 1 + temp_sc + prev_rainfall_sc + total_rainfall_sc + (1 \| year), pforms = list(), pfix = list(), resp = “meanrhsc”, family = list(family = “gaussian”, link = “identity”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “sigma”), type = “real”, ybounds = c(-Inf, Inf), closed = c(NA, NA), ad = c(“weights”, “subset”, “se”, “cens”, “trunc”, “mi”, “index”), normalized = c(“_time_hom”, “_time_het”, “_lagsar”, “_errorsar”, “_fcor”), specials = c(“residuals”, “rescor”)), mecor = TRUE), nummales = list(formula = num_males ~ mean_rh_sc + temp_sc + prev_rainfall_sc + total_rainfall_sc + moonlight_sc + (1 \| year), pforms = list(), pfix = list(), family = list(family = “negbinomial”, link = “log”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “shape”), type = “int”, ybounds = c(0, Inf), closed = c(TRUE, NA), ad = c(“weights”, “subset”, “cens”, “trunc”, “rate”, “index”), specials = “sbi_log”, prior = function (dpar, link = “identity”, …) { if (dpar == “shape” && link == “identity”) { return(“inv_gamma(0.4, 0.3)”) } NULL }, link_shape = “log”), resp = “nummales”, mecor = TRUE), numpairs = list(formula = num_pairs ~ mean_rh_sc + temp_sc + prev_rainfall_sc + total_rainfall_sc + moonlight_sc + num_males_sc + (1 \| year), pforms = list(), pfix = list(), family = list(family = “negbinomial”, link = “log”, linkfun = function (mu) link(mu, link = slink), linkinv = function (eta) inv_link(eta, link = slink), dpars = c(“mu”, “shape”), type = “int”, ybounds = c(0, Inf), closed = c(TRUE, NA), ad = c(“weights”, “subset”, “cens”, “trunc”, “rate”, “index”), specials = “sbi_log”, prior = function (dpar, link = “identity”, …) { if (dpar == “shape” && link == “identity”) { return(“inv_gamma(0.4, 0.3)”) } NULL }, link_shape = “log”), resp = “numpairs”, mecor = TRUE))	()	b-normal(0, 0.5) b-normal(0, 0.5) b-normal(0, 0.5) b-normal(0, 0.5) Intercept-student_t(3, 0, 2.5) Intercept-student_t(3, 3.8, 2.5) Intercept-student_t(3, 0.7, 2.5) Intercept-student_t(3, -0.4, 2.5) Intercept-student_t(3, 0.1, 2.5) Intercept-student_t(3, -0.5, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) sd-student_t(3, 0, 2.5) shape-exponential(1) shape-exponential(1) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5)	4000	4	1	2000	10 (0.0012%)	0	2513.228	2938.475	896167421

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_totalrainfallsc_Intercept	-0.088	-0.388	0.242	1.001	4374.462	3150.871
b_prevrainfallsc_Intercept	-0.001	-0.278	0.282	1.002	4881.716	3878.943
b_tempsc_Intercept	-0.018	-0.628	0.617	1.002	2613.281	3288.814
b_meanrhsc_Intercept	-0.013	-0.606	0.586	1.001	2551.544	2938.475
b_nummales_Intercept	3.736	3.190	4.273	1.003	2513.228	3132.459
b_numpairs_Intercept	0.811	-0.116	1.785	1.001	2736.804	3220.637
b_tempsc_total_rainfall_sc	-0.557	-0.736	-0.379	1	12578.392	5404.858
b_meanrhsc_temp_sc	-0.560	-0.735	-0.389	1	7996.921	6111.927
b_meanrhsc_prev_rainfall_sc	-0.116	-0.247	0.023	1.001	9901.264	5912.441
b_meanrhsc_total_rainfall_sc	0.123	-0.062	0.308	1	10276.103	6569.520
b_nummales_mean_rh_sc	0.065	-0.026	0.154	1.001	7744.325	4693.111
b_nummales_temp_sc	-0.139	-0.234	-0.043	1	7644.709	6677.752
b_nummales_prev_rainfall_sc	-0.054	-0.122	0.012	1	11302.672	5540.365
b_nummales_total_rainfall_sc	0.019	-0.073	0.113	1.001	11183.893	6330.711
b_nummales_moonlight_sc	-0.019	-0.080	0.043	1.001	14502.490	5881.696
b_numpairs_mean_rh_sc	0.105	-0.137	0.342	1.003	10176.743	5900.705
b_numpairs_temp_sc	-0.181	-0.432	0.073	1	8167.349	6235.935
b_numpairs_prev_rainfall_sc	0.070	-0.093	0.241	1.001	11606.833	6640.369
b_numpairs_total_rainfall_sc	0.262	0.048	0.481	1.001	11591.560	6388.269
b_numpairs_moonlight_sc	-0.118	-0.279	0.038	1	13629.260	5617.375
b_numpairs_num_males_sc	0.741	0.420	1.066	1	7156.220	5767.970

Code

es <- extended_summary(
  mean_temp_females_mod_no_outlier,
  highlight = TRUE,
  trace.palette = mako,
  remove.intercepts = TRUE,
  print.name = FALSE,
  beta.prefix = "b_numpairs",
  trace = FALSE,
  fill = mako(10, alpha = 0.8)[7],
  gsub.pattern = c(
    "b_numpairs_",
    "_sc",
    "mean_rh",
    "moonlight",
    "prev_rainfall",
    "temp",
    "total_rainfall",
    "num_males"
  ),
  gsub.replacement = c(
    "",
    "",
    "Relative humidity",
    "Moonlight",
    "Previous rainfall",
    "Temperature",
    "Total rainfall",
    "Number of males"
  ),
  return = TRUE, 
  ylab = "Predictor"
)

es$plot

3.2.2.4 Evaluate “year” effect

Code

# ----------------------------
# Submodels with year random intercept
# ----------------------------


bf_total_rain <- bf(
  total_rainfall_sc ~ 1 
)

bf_prev_rain <- bf(
  prev_rainfall_sc ~ 1 
)

bf_temp <- bf(
  temp_sc ~ 1 +
    total_rainfall_sc
)

bf_rh <- bf(
  mean_rh_sc ~ 1 +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc
)

# mediator: active males (count)
bf_males <- bf(
  num_males ~ 
    s(mean_rh_sc) +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    moonlight_sc,
  family = negbinomial()
)

# outcome: active females (count)
bf_pairs <- bf(
  num_pairs ~ 
    mean_rh_sc +
    temp_sc +
    prev_rainfall_sc +
    total_rainfall_sc +
    moonlight_sc +
    num_males_sc,
  family = negbinomial()
)

# ----------------------------
# Weakly informative priors
# ----------------------------

priors <- c(

  # Linear slopes (all predictors z-scored)
  prior(normal(0, 0.5), class = "b", resp = "meanrhsc"),
  prior(normal(0, 0.5), class = "b", resp = "tempsc"),
  prior(normal(0, 0.5), class = "b", resp = "nummales"),
  prior(normal(0, 0.5), class = "b", resp = "numpairs"),

  # Smooth terms (humidity effects)
  prior(normal(0, 1), class = "sds", resp = "nummales"),
 
  # Negative binomial dispersion
  prior(exponential(1), class = "shape", resp = "nummales"),
  prior(exponential(1), class = "shape", resp = "numpairs")
)

# Fit full SEM
fem_sem_model_no_year <- brm(
  bf_total_rain +
    bf_prev_rain +
    bf_temp +
    bf_rh +
    bf_males +
    bf_pairs +
    set_rescor(FALSE),
  data = dat,
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 4000,
  control = list(adapt_delta = 0.95),
  file = "./data/processed/female_sem_model_no_year.RDS",
  backend = "cmdstanr",
  file_refit = "on_change"
)

add_criterion(fem_sem_model_no_year, criterion = c("loo"))

Code

no_year <- readRDS("./data/processed/female_sem_model_no_year.RDS")

loo_diff <- loo::loo_compare(mean_temp_females_mod, no_year)

rows <- rownames(loo_diff)
loo_diff <- as.data.frame(loo_diff)

loo_diff$model <- rows

loo_diff

	elpd_diff	se_diff	elpd_loo	se_elpd_loo	p_loo	se_p_loo	looic	se_looic	model
mean_temp_females_mod	0.0000	0.00000	-1344.815	33.34107	48.00978	7.811394	2689.630	66.68214	mean_temp_females_mod
no_year	-110.4541	11.38276	-1455.269	34.02235	34.44033	7.684772	2910.539	68.04469	no_year

Code

R2_full <- bayes_R2(mean_temp_females_mod, resp = "numpairs",  summary = FALSE)

R2_fixed <- bayes_R2(
  mean_temp_females_mod,
  resp = "numpairs",
  re_formula = NA,
   summary = FALSE
)

delta_R2 <- R2_full - R2_fixed

quantile(delta_R2, c(0.025, 0.5, 0.975))

       2.5%         50%       97.5% 
-0.15196871 -0.03475038  0.08708981

3.2.3 Parametized DAG

Code

edge_data <- tidy_dag_fem$data
fixef_df <- as.data.frame(fixef(mean_temp_females_mod))

fixef_df$term <- rownames(fixef_df)

fixef_df <- fixef_df[grep(pattern = "Intercept", fixef_df$term, invert = TRUE), ]

split_terms <- strsplit(fixef_df$term, "_")

fixef_df$to <- sapply(split_terms, `[`, 1)
fixef_df$to <- gsub("sc", "_sc", fixef_df$to)
fixef_df$to <- gsub("mean", "mean_", fixef_df$to)
fixef_df$to <- gsub("num", "num_", fixef_df$to)
fixef_df$to <- gsub("temp_", "mean_temp_", fixef_df$to)
fixef_df$to <- gsub("males", "males_sc", fixef_df$to)

# split from "_" to end
fixef_df$from <- sapply(split_terms, function(x) paste(x[-1], collapse = "_"))
fixef_df$from <- gsub("temp_", "mean_temp_", fixef_df$from)
fixef_df$from <- gsub("males$", "males_sc", fixef_df$from)

sem_estimates <- data.frame(
  from = fixef_df$from,
  to   = fixef_df$to,
  est  = fixef_df$Estimate,
  l95  = fixef_df$Q2.5,
  u95  = fixef_df$Q97.5,
  p    = 2 * pmin(
            pnorm(fixef_df$Estimate / fixef_df$Est.Error),
            1 - pnorm(fixef_df$Estimate / fixef_df$Est.Error)
         )
)


edge_data <- merge(edge_data, sem_estimates[, c("from", "to", "est", "l95", "u95")], by.x = c("name", "to"), by.y = c("from","to"), all.x = TRUE)

edge_data$sig <- ifelse(
  edge_data$l95 * edge_data$u95 > 0 & !is.na(edge_data$l95),
  "solid",
  "dashed"
)

edge_data$alpha <- ifelse(
  edge_data$l95 * edge_data$u95 > 0 & !is.na(edge_data$l95),
  1,
  0
)

edge_data$color <- ifelse(
  edge_data$l95 * edge_data$u95 > 0 & !is.na(edge_data$l95),
  "black",
  "gray"
)

edge_data$color[edge_data$name == "num_pairs"] <- "black"
edge_data$est <- ifelse(
!is.na(edge_data$l95),
  edge_data$est,
  0
)


tidy_dag_param_fem <- tidy_dag_fem

tidy_dag_param_fem$data <- edge_data

ggplot(tidy_dag_param_fem, aes(
    x = x,
    y = y,
    xend = xend,
    yend = yend,
    edge_linetype = sig,
    edge_alpha = alpha
)) + scale_color_viridis_d(
    option = "G",
    begin = 0.2,
    end = 0.8,
    alpha = 0.5
) + 
    geom_dag_edges_fan(aes(size = (est / max(abs(est), na.rm = TRUE))^2),
    arrow = grid::arrow(length = grid::unit(10, "pt"), type = "closed")
) +
    geom_dag_point(aes(color = type), 
                   size = 24) + 
    geom_dag_text(aes(label = label), color = c(NA, c(edge_data$color, edge_data$color)[!duplicated(c(edge_data$name, edge_data$to))])) +
    scale_edge_alpha(range = c(0.2, 1), guide = NULL) +
    scale_edge_linetype_manual(values = c("solid" = "solid", "dashed" = "dashed"), guide = NULL) +
    theme_dag() + 
    theme(legend.position = "bottom") +
    guides(size = "none")

3.2.4 Graph raw data

Temperature

Code

temp_grid <- data.frame(
  mean_temp = seq(min(dat$mean_temp, na.rm = TRUE),
                  max(dat$mean_temp, na.rm = TRUE),
                  length.out = 100)
)
temp_grid$mean_rh_sc        <- 0
temp_grid$prev_rainfall_sc <- 0
temp_grid$total_rainfall_sc <- 0
temp_grid$moonlight_sc     <- 0
temp_grid$num_males_sc     <- 0


temp_center <- attr(dat$mean_temp_sc, "scaled:center")
temp_scale  <- attr(dat$mean_temp_sc, "scaled:scale")

temp_grid$temp_sc <- (temp_grid$mean_temp - temp_center) / temp_scale


pred <- fitted(
  mean_temp_females_mod,
  resp = "numpairs",
  newdata = temp_grid,
  re_formula = NA,   # population-level effect
  summary = TRUE
)

temp_grid <- cbind(temp_grid, pred)

ggplot(data = dat, aes(x = mean_temp, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3)  +
       labs(x = "Mean temperature (°C)", y = "Number of amplexus")

Code

ggplot(data = dat, aes(x = mean_temp, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
  geom_line(data = temp_grid,
            aes(x = mean_temp, y = Estimate),
            linewidth = 1.2,
            color = mako(10, alpha = 0.8)[3]) +
  geom_ribbon(data = temp_grid,
              aes(x = mean_temp, y = Estimate, ymin = Q2.5, ymax = Q97.5),
              alpha = 0.1,
              fill = "black") +
       labs(x = "Mean temperature (°C)", y = "Number of amplexus")

Code

# by year
temp_grid_year <- data.frame(
  mean_temp = seq(min(dat$mean_temp, na.rm = TRUE),
                  max(dat$mean_temp, na.rm = TRUE),
                  length.out = 100),
  year = unique(dat$year)
)
temp_grid_year$mean_rh_sc        <- 0
temp_grid_year$prev_rainfall_sc <- 0
temp_grid_year$total_rainfall_sc <- 0
temp_grid_year$moonlight_sc     <- 0
temp_grid_year$num_males_sc     <- 0

temp_center <- attr(dat$mean_temp_sc, "scaled:center")
temp_scale  <- attr(dat$mean_temp_sc, "scaled:scale")

temp_grid_year$temp_sc <- (temp_grid_year$mean_temp - temp_center) / temp_scale


pred <- fitted(
  mean_temp_females_mod,
  resp = "numpairs",
  newdata = temp_grid_year,
  re_formula = NULL,   # <-- this is the key
  summary = TRUE
)

temp_grid_year <- cbind(temp_grid_year, pred)

ggplot(data = dat, aes(x = mean_temp, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
    facet_wrap(~year) +
    labs(x = "Mean temperature (°C)", y = "Number of amplexus")

Code

ggplot(data = dat, aes(x = mean_temp, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
  geom_line(data = temp_grid_year,
            aes(x = mean_temp, y = Estimate),
            linewidth = 1.2,
            color = mako(10, alpha = 0.8)[3]) +
    facet_wrap(~year) +
  geom_ribbon(data = temp_grid_year,
              aes(x = mean_temp, y = Estimate, ymin = Q2.5, ymax = Q97.5),
              alpha = 0.1,
              fill = "black") +
    labs(x = "Mean temperature (°C)", y = "Number of amplexus")

Total rainfall

Code

temp_grid <- data.frame(
  total_rainfall = seq(min(dat$total_rainfall, na.rm = TRUE),
                  max(dat$total_rainfall, na.rm = TRUE),
                  length.out = 100)
)

temp_grid$temp_sc     <- 0
temp_grid$mean_rh_sc        <- 0
temp_grid$prev_rainfall_sc <- 0
temp_grid$moonlight_sc     <- 0
temp_grid$num_males_sc     <- 0


temp_center <- attr(dat$total_rainfall_sc, "scaled:center")
temp_scale  <- attr(dat$total_rainfall_sc, "scaled:scale")

temp_grid$total_rainfall_sc <- (temp_grid$total_rainfall - temp_center) / temp_scale


pred <- fitted(
  mean_temp_females_mod,
  resp = "numpairs",
  newdata = temp_grid,
  re_formula = NA,   # population-level effect
  summary = TRUE
)

temp_grid <- cbind(temp_grid, pred)


ggplot(data = dat, aes(x = total_rainfall, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
  geom_line(data = temp_grid,
            aes(x = total_rainfall, y = Estimate),
            linewidth = 1.2,
            color = mako(10, alpha = 0.8)[3]) +
  geom_ribbon(data = temp_grid,
              aes(x = total_rainfall, y = Estimate, ymin = Q2.5, ymax = Q97.5),
              alpha = 0.1,
              fill = "black") +
       labs(x = "Total rainfall (mm)", y = "Number of amplexus")

Code

# by year
temp_grid_year <- expand.grid(
  total_rainfall = seq(min(dat$total_rainfall, na.rm = TRUE),
                  max(dat$total_rainfall, na.rm = TRUE)),
                  length.out = 100,
  year = unique(dat$year))

temp_grid_year$temp_sc     <- 0
temp_grid_year$mean_rh_sc        <- 0
temp_grid_year$prev_rainfall_sc <- 0
temp_grid_year$moonlight_sc     <- 0
temp_grid_year$num_males_sc     <- 0

temp_center <- attr(dat$total_rainfall_sc, "scaled:center")
temp_scale  <- attr(dat$total_rainfall_sc, "scaled:scale")

temp_grid_year$total_rainfall_sc <- (temp_grid_year$total_rainfall - temp_center) / temp_scale


pred <- fitted(
  mean_temp_females_mod,
  resp = "numpairs",
  newdata = temp_grid_year,
  re_formula = NULL,   # <-- this is the key
  summary = TRUE
)

temp_grid_year <- cbind(temp_grid_year, pred)

ggplot(data = dat, aes(x = total_rainfall, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
  geom_line(data = temp_grid_year,
            aes(x = total_rainfall, y = Estimate),
            linewidth = 1.2,
            color = mako(10, alpha = 0.8)[3]) +
    facet_wrap(~year) +
  geom_ribbon(data = temp_grid_year,
              aes(x = total_rainfall, y = Estimate, ymin = Q2.5, ymax = Q97.5),
              alpha = 0.1,
              fill = "black") +
    labs(x = "Total rainfall (mm)", y = "Number of amplexus")

Total rainfall excluding outlier

Code

dat_no_outlier <- dat[dat$total_rainfall < 75 & !is.na(dat$num_pairs), ]

dat_no_outlier$total_rainfall_sc <- scale(dat_no_outlier$total_rainfall)

temp_grid <- data.frame(
  total_rainfall = seq(min(dat_no_outlier$total_rainfall, na.rm = TRUE),
                  max(dat_no_outlier$total_rainfall, na.rm = TRUE),
                  length.out = 100)
)

temp_grid$temp_sc     <- 0
temp_grid$mean_rh_sc        <- 0
temp_grid$prev_rainfall_sc <- 0
temp_grid$moonlight_sc     <- 0
temp_grid$num_males_sc     <- 0


temp_center <- attr(dat_no_outlier$total_rainfall_sc, "scaled:center")
temp_scale  <- attr(dat_no_outlier$total_rainfall_sc, "scaled:scale")

temp_grid$total_rainfall_sc <- (temp_grid$total_rainfall - temp_center) / temp_scale


pred <- fitted(
  mean_temp_females_mod_no_outlier,
  resp = "numpairs",
  newdata = temp_grid,
  re_formula = NA,   # population-level effect
  summary = TRUE
)

temp_grid <- cbind(temp_grid, pred)


ggplot(data = dat_no_outlier, aes(x = total_rainfall, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
  geom_line(data = temp_grid,
            aes(x = total_rainfall, y = Estimate),
            linewidth = 1.2,
            color = mako(10, alpha = 0.8)[3]) +
  geom_ribbon(data = temp_grid,
              aes(x = total_rainfall, y = Estimate, ymin = Q2.5, ymax = Q97.5),
              alpha = 0.1,
              fill = "black") +
       labs(x = "Total rainfall (mm)", y = "Number of amplexus")

Code

# by year
temp_grid_year <- expand.grid(
  total_rainfall = seq(min(dat_no_outlier$total_rainfall, na.rm = TRUE),
                  max(dat_no_outlier$total_rainfall, na.rm = TRUE)),
                  length.out = 100,
  year = unique(dat_no_outlier$year))

temp_grid_year$temp_sc     <- 0
temp_grid_year$mean_rh_sc        <- 0
temp_grid_year$prev_rainfall_sc <- 0
temp_grid_year$moonlight_sc     <- 0
temp_grid_year$num_males_sc     <- 0

temp_center <- attr(dat_no_outlier$total_rainfall_sc, "scaled:center")
temp_scale  <- attr(dat_no_outlier$total_rainfall_sc, "scaled:scale")

temp_grid_year$total_rainfall_sc <- (temp_grid_year$total_rainfall - temp_center) / temp_scale


pred <- fitted(
  mean_temp_females_mod_no_outlier,
  resp = "numpairs",
  newdata = temp_grid_year,
  re_formula = NULL,   # <-- this is the key
  summary = TRUE
)

temp_grid_year <- cbind(temp_grid_year, pred)

ggplot(data = dat_no_outlier, aes(x = total_rainfall, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
  geom_line(data = temp_grid_year,
            aes(x = total_rainfall, y = Estimate),
            linewidth = 1.2,
            color = mako(10, alpha = 0.8)[3]) +
    facet_wrap(~year) +
  geom_ribbon(data = temp_grid_year,
              aes(x = total_rainfall, y = Estimate, ymin = Q2.5, ymax = Q97.5),
              alpha = 0.1,
              fill = "black") +
    labs(x = "Total rainfall (mm)", y = "Number of amplexus")

Number of males

Code

temp_grid <- data.frame(
  num_males = seq(min(dat$num_males, na.rm = TRUE),
                  max(dat$num_males, na.rm = TRUE),
                  length.out = 100)
)

temp_grid$temp_sc     <- 0
temp_grid$mean_rh_sc        <- 0
temp_grid$prev_rainfall_sc <- 0
temp_grid$total_rainfall_sc <- 0
temp_grid$moonlight_sc     <- 0


temp_center <- attr(dat$num_males_sc, "scaled:center")
temp_scale  <- attr(dat$num_males_sc, "scaled:scale")

temp_grid$num_males_sc <- (temp_grid$num_males - temp_center) / temp_scale


pred <- fitted(
  mean_temp_females_mod,
  resp = "numpairs",
  newdata = temp_grid,
  re_formula = NA,   # population-level effect
  summary = TRUE
)

temp_grid <- cbind(temp_grid, pred)

ggplot(data = dat, aes(x = num_males, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
  geom_line(data = temp_grid,
            aes(x = num_males, y = Estimate),
            linewidth = 1.2,
            color = mako(10, alpha = 0.8)[3]) +
  geom_ribbon(data = temp_grid,
              aes(x = num_males, y = Estimate, ymin = Q2.5, ymax = Q97.5),
              alpha = 0.1,
              fill = "black") +
       labs(x = "Number of males", y = "Number of amplexus")

Code

# by year
temp_grid_year <- expand.grid(
  num_males = seq(min(dat$num_males, na.rm = TRUE),
                  max(dat$num_males, na.rm = TRUE)),
                  length.out = 100,
  year = unique(dat$year))

temp_grid_year$temp_sc     <- 0
temp_grid_year$mean_rh_sc        <- 0
temp_grid_year$prev_rainfall_sc <- 0
temp_grid_year$total_rainfall_sc <- 0
temp_grid_year$moonlight_sc     <- 0

temp_center <- attr(dat$num_males_sc, "scaled:center")
temp_scale  <- attr(dat$num_males_sc, "scaled:scale")

temp_grid_year$num_males_sc <- (temp_grid_year$num_males - temp_center) / temp_scale


pred <- fitted(
  mean_temp_females_mod,
  resp = "numpairs",
  newdata = temp_grid_year,
  re_formula = NULL,   # <-- this is the key
  summary = TRUE
)

temp_grid_year <- cbind(temp_grid_year, pred)

ggplot(data = dat, aes(x = num_males, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
  geom_line(data = temp_grid_year,
            aes(x = num_males, y = Estimate),
            linewidth = 1.2,
            color = mako(10, alpha = 0.8)[3]) +
    facet_wrap(~year) +
  geom_ribbon(data = temp_grid_year,
              aes(x = num_males, y = Estimate, ymin = Q2.5, ymax = Q97.5),
              alpha = 0.1,
              fill = "black") +
    labs(x = "Number of males", y = "Number of amplexus")

Previous rainfall

Code

ggplot(data = dat, aes(x = prev_rainfall, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
       labs(x = "Previous rainfall (mm)", y = "Number of amplexus")

Moonlight

Code

ggplot(data = dat, aes(x = moonlight * 100, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
       labs(x = "Moonlight (%)", y = "Number of amplexus")

Relative humidity

Code

ggplot(data = dat, aes(x = mean_rh, y = num_pairs)) +
    geom_point(alpha = 0.4, color = mako(6, alpha = 0.6)[3], size = 3) +
       labs(x = "Relative humidity (%)", y = "Number of amplexus")

Takeaways

Session information

Click to see

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.2 (2023-10-31)
 os       Ubuntu 22.04.2 LTS
 system   x86_64, linux-gnu
 ui       X11
 language es_ES:en
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Costa_Rica
 date     2026-02-22
 pandoc   3.8 @ /home/marce/anaconda3/bin/ (via rmarkdown)
 quarto   1.4.549 @ /usr/local/bin/quarto

─ Packages ───────────────────────────────────────────────────────────────────
 package        * version  date (UTC) lib source
 abind            1.4-8    2024-09-12 [1] CRAN (R 4.3.2)
 ape              5.8-1    2024-12-16 [1] CRAN (R 4.3.2)
 arrayhelpers     1.1-0    2020-02-04 [1] CRAN (R 4.3.2)
 backports        1.5.0    2024-05-23 [1] CRAN (R 4.3.2)
 bayesplot        1.15.0   2025-12-12 [1] CRAN (R 4.3.2)
 boot             1.3-28   2021-05-03 [4] CRAN (R 4.1.1)
 bridgesampling   1.2-1    2025-11-19 [1] CRAN (R 4.3.2)
 brms           * 2.23.0   2025-09-09 [1] CRAN (R 4.3.2)
 brmsish        * 1.0.0    2025-12-18 [1] CRANs (R 4.3.2)
 Brobdingnag      1.2-9    2022-10-19 [1] CRAN (R 4.3.2)
 cachem           1.1.0    2024-05-16 [1] CRAN (R 4.3.2)
 checkmate        2.3.4    2026-02-03 [1] CRAN (R 4.3.2)
 cli              3.6.5    2025-04-23 [1] CRAN (R 4.3.2)
 cmdstanr         0.8.1    2024-07-08 [1] https://stan-dev.r-universe.dev (R 4.3.2)
 coda             0.19-4.1 2024-01-31 [1] CRAN (R 4.3.2)
 codetools        0.2-18   2020-11-04 [4] CRAN (R 4.0.3)
 cowplot          1.2.0    2025-07-07 [1] CRAN (R 4.3.2)
 crayon           1.5.3    2024-06-20 [1] CRAN (R 4.3.2)
 curl             7.0.0    2025-08-19 [1] CRAN (R 4.3.2)
 dagitty        * 0.3-4    2023-12-07 [1] CRAN (R 4.3.2)
 devtools         2.4.6    2025-10-03 [1] CRAN (R 4.3.2)
 digest           0.6.39   2025-11-19 [1] CRAN (R 4.3.2)
 distributional   0.6.0    2026-01-14 [1] CRAN (R 4.3.2)
 dplyr            1.1.4    2023-11-17 [1] CRAN (R 4.3.2)
 ellipsis         0.3.2    2021-04-29 [3] CRAN (R 4.1.1)
 emmeans          2.0.1    2025-12-16 [1] CRAN (R 4.3.2)
 estimability     1.5.1    2024-05-12 [1] CRAN (R 4.3.2)
 evaluate         1.0.5    2025-08-27 [1] CRAN (R 4.3.2)
 farver           2.1.2    2024-05-13 [1] CRAN (R 4.3.2)
 fastmap          1.2.0    2024-05-15 [1] CRAN (R 4.3.2)
 Formula          1.2-5    2023-02-24 [1] CRAN (R 4.3.2)
 fs               1.6.6    2025-04-12 [1] CRAN (R 4.3.2)
 generics         0.1.4    2025-05-09 [1] CRAN (R 4.3.2)
 ggdag          * 0.2.13   2024-07-22 [1] CRAN (R 4.3.2)
 ggdist           3.3.3    2025-04-23 [1] CRAN (R 4.3.2)
 ggforce          0.5.0    2025-06-18 [1] CRAN (R 4.3.2)
 ggparty        * 1.0.0.1  2025-07-10 [1] CRAN (R 4.3.2)
 ggplot2        * 4.0.2    2026-02-03 [1] CRAN (R 4.3.2)
 ggraph         * 2.2.2    2025-08-24 [1] CRAN (R 4.3.2)
 ggrepel          0.9.6    2024-09-07 [1] CRAN (R 4.3.2)
 glue             1.8.0    2024-09-30 [1] CRAN (R 4.3.2)
 graphlayouts     1.2.2    2025-01-23 [1] CRAN (R 4.3.2)
 gridExtra        2.3      2017-09-09 [3] CRAN (R 4.0.1)
 gtable           0.3.6    2024-10-25 [1] CRAN (R 4.3.2)
 htmltools        0.5.9    2025-12-04 [1] CRAN (R 4.3.2)
 htmlwidgets      1.6.4    2023-12-06 [1] CRAN (R 4.3.2)
 igraph           2.2.1    2025-10-27 [1] CRAN (R 4.3.2)
 inline           0.3.21   2025-01-09 [1] CRAN (R 4.3.2)
 inum             1.0-5    2023-03-09 [1] CRAN (R 4.3.2)
 jsonlite         2.0.0    2025-03-27 [1] CRAN (R 4.3.2)
 kableExtra       1.4.0    2024-01-24 [1] CRAN (R 4.3.2)
 knitr            1.51     2025-12-20 [1] CRAN (R 4.3.2)
 labeling         0.4.3    2023-08-29 [1] CRAN (R 4.3.2)
 lattice          0.20-45  2021-09-22 [4] CRAN (R 4.1.1)
 libcoin        * 1.0-10   2023-09-27 [1] CRAN (R 4.3.2)
 lifecycle        1.0.5    2026-01-08 [1] CRAN (R 4.3.2)
 loo              2.9.0    2025-12-23 [1] CRAN (R 4.3.2)
 lunar          * 0.2-1    2022-06-13 [1] CRAN (R 4.3.2)
 magrittr         2.0.4    2025-09-12 [1] CRAN (R 4.3.2)
 MASS             7.3-55   2022-01-13 [4] CRAN (R 4.1.2)
 Matrix           1.6-5    2024-01-11 [1] CRAN (R 4.3.2)
 matrixStats      1.5.0    2025-01-07 [1] CRAN (R 4.3.2)
 memoise          2.0.1    2021-11-26 [3] CRAN (R 4.1.2)
 multcomp         1.4-29   2025-10-20 [1] CRAN (R 4.3.2)
 mvtnorm        * 1.3-3    2025-01-10 [1] CRAN (R 4.3.2)
 nlme             3.1-155  2022-01-13 [4] CRAN (R 4.1.2)
 otel             0.2.0    2025-08-29 [1] CRAN (R 4.3.2)
 packrat          0.9.3    2025-06-16 [1] CRAN (R 4.3.2)
 partykit       * 1.2-24   2025-05-02 [1] CRAN (R 4.3.2)
 pbapply          1.7-4    2025-07-20 [1] CRAN (R 4.3.2)
 pillar           1.11.1   2025-09-17 [1] CRAN (R 4.3.2)
 pkgbuild         1.4.8    2025-05-26 [1] CRAN (R 4.3.2)
 pkgconfig        2.0.3    2019-09-22 [3] CRAN (R 4.0.1)
 pkgload          1.4.1    2025-09-23 [1] CRAN (R 4.3.2)
 plyr             1.8.9    2023-10-02 [1] CRAN (R 4.3.2)
 polyclip         1.10-7   2024-07-23 [1] CRAN (R 4.3.2)
 posterior        1.6.1    2025-02-27 [1] CRAN (R 4.3.2)
 processx         3.8.6    2025-02-21 [1] CRAN (R 4.3.2)
 ps               1.9.1    2025-04-12 [1] CRAN (R 4.3.2)
 purrr            1.2.1    2026-01-09 [1] CRAN (R 4.3.2)
 QuickJSR         1.8.1    2025-09-20 [1] CRAN (R 4.3.2)
 R6               2.6.1    2025-02-15 [1] CRAN (R 4.3.2)
 RColorBrewer     1.1-3    2022-04-03 [1] CRAN (R 4.3.2)
 Rcpp           * 1.1.1    2026-01-10 [1] CRAN (R 4.3.2)
 RcppParallel     5.1.11-1 2025-08-27 [1] CRAN (R 4.3.2)
 remotes          2.5.0    2024-03-17 [1] CRAN (R 4.3.2)
 reshape2         1.4.5    2025-11-12 [1] CRAN (R 4.3.2)
 rlang            1.1.7    2026-01-09 [1] CRAN (R 4.3.2)
 rmarkdown        2.30     2025-09-28 [1] CRAN (R 4.3.2)
 rpart            4.1.16   2022-01-24 [4] CRAN (R 4.1.2)
 rstan          * 2.32.7   2025-03-10 [1] CRAN (R 4.3.2)
 rstantools       2.6.0    2026-01-10 [1] CRAN (R 4.3.2)
 rstudioapi       0.18.0   2026-01-16 [1] CRAN (R 4.3.2)
 S7               0.2.1    2025-11-14 [1] CRAN (R 4.3.2)
 sandwich         3.1-1    2024-09-15 [1] CRAN (R 4.3.2)
 scales           1.4.0    2025-04-24 [1] CRAN (R 4.3.2)
 sessioninfo      1.2.3    2025-02-05 [1] CRAN (R 4.3.2)
 sketchy          1.0.6    2025-12-04 [1] local (/home/marce/Dropbox/R_package_testing/sketchy)
 StanHeaders    * 2.32.10  2024-07-15 [1] CRAN (R 4.3.2)
 stringi          1.8.7    2025-03-27 [1] CRAN (R 4.3.2)
 stringr          1.6.0    2025-11-04 [1] CRAN (R 4.3.2)
 survival         3.2-13   2021-08-24 [4] CRAN (R 4.1.1)
 svglite          2.2.2    2025-10-21 [1] CRAN (R 4.3.2)
 svUnit           1.0.8    2025-08-26 [1] CRAN (R 4.3.2)
 systemfonts      1.3.1    2025-10-01 [1] CRAN (R 4.3.2)
 tensorA          0.36.2.1 2023-12-13 [1] CRAN (R 4.3.2)
 textshaping      1.0.4    2025-10-10 [1] CRAN (R 4.3.2)
 TH.data          1.1-5    2025-11-17 [1] CRAN (R 4.3.2)
 tibble           3.3.1    2026-01-11 [1] CRAN (R 4.3.2)
 tidybayes        3.0.7    2024-09-15 [1] CRAN (R 4.3.2)
 tidygraph        1.3.1    2024-01-30 [1] CRAN (R 4.3.2)
 tidyr            1.3.2    2025-12-19 [1] CRAN (R 4.3.2)
 tidyselect       1.2.1    2024-03-11 [1] CRAN (R 4.3.2)
 tweenr           2.0.3    2024-02-26 [1] CRAN (R 4.3.2)
 usethis          3.2.1    2025-09-06 [1] CRAN (R 4.3.2)
 V8               8.0.1    2025-10-10 [1] CRAN (R 4.3.2)
 vctrs            0.7.1    2026-01-23 [1] CRAN (R 4.3.2)
 viridis        * 0.6.5    2024-01-29 [1] CRAN (R 4.3.2)
 viridisLite    * 0.4.3    2026-02-04 [1] CRAN (R 4.3.2)
 withr            3.0.2    2024-10-28 [1] CRAN (R 4.3.2)
 xaringanExtra    0.8.0    2024-05-19 [1] CRAN (R 4.3.2)
 xfun             0.56     2026-01-18 [1] CRAN (R 4.3.2)
 xml2             1.5.2    2026-01-17 [1] CRAN (R 4.3.2)
 xtable           1.8-4    2019-04-21 [3] CRAN (R 4.0.1)
 yaml             2.3.12   2025-12-10 [1] CRAN (R 4.3.2)
 zoo              1.8-15   2025-12-15 [1] CRAN (R 4.3.2)

 [1] /home/marce/R/x86_64-pc-linux-gnu-library/4.3
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library
 * ── Packages attached to the search path.

──────────────────────────────────────────────────────────────────────────────