Stats

Frog female response to male and female calls

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

Purpose

• Evaluate the effect of different male and female vocal stimuli on the response of female frogs

Load packages

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

# install/ load packages
sketchy::load_packages(packages = c("Rraven", "warbleR", "ggplot2",
    "brms", "brmsish", "viridis", "emmeans", "cowplot", "loo", "rstan"))

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

# to create several posterior predictive check plots out of a brms fit
custom_ppc <- function(fit, group = NULL, ndraws = 30) {
  by_group  <- if (!is.null(group)){
    TRUE 
    } else
    FALSE
  
  if (by_group)
    by_group  <-  if (any(names(fit$data) == group)){
      TRUE
      } else
    FALSE
  
  if (by_group)
    by_group <-
      if (is.character(fit$data[, group]) |
          is.factor(fit$data[, group])){
        TRUE
        } else
    FALSE
  
  
  if (by_group) {
    ppc_dens <- pp_check(fit,
                         ndraws = ndraws,
                         type = 'dens_overlay_grouped',
                         group = group)
    
    pp_mean <- pp_check(
      fit,
      type = "stat_grouped",
      stat = "mean",
      group = group,
      ndraws = ndraws
    )  + theme_classic()
    
    pp_scat <- pp_check(fit,
                        type = "scatter_avg",
                        # group = group,
                        ndraws = ndraws)
  } else {
    ppc_dens <- pp_check(fit,
                         ndraws = ndraws,
                         type = 'dens_overlay')
    
    pp_mean <- pp_check(fit,
                        type = "stat",
                        stat = "mean",
                        ndraws = ndraws) + theme_classic()
    
    pp_scat <-  pp_check(fit,
                         type = "scatter_avg",
                         ndraws = ndraws)
  }
  
  pp_stat2 <- pp_check(fit, type = "stat_2d", ndraws = ndraws)
  
  pp_plot_list <-
    list(ppc_dens, pp_mean, pp_scat,  pp_stat2)
  
  pp_plot_list[c(1, 3:4)] <-
    lapply(pp_plot_list[c(1, 3:4)], function(x)
      x  + scale_color_viridis_d(
        begin = 0.3,
        end = 0.8,
        alpha = 0.5,
        option = "mako",
      ) + scale_fill_viridis_d(
        begin = 0.3,
        end = 0.8,
        alpha = 0.5,
        option = "mako"
      ) + theme_classic())
  
  
  ppc_plot <- plot_grid(plotlist = pp_plot_list, ncol = 2)
  
  print(ppc_plot)
}

1 Read data

dat <- imp_raven(path = "/home/m/Dropbox/Projects/frog_playback_experiment/data/raw/Lpa_female playbacks_selection tables-20240917T211317Z-001/Lpa_female playbacks_selection tables/",
    all.data = TRUE, warbler.format = TRUE, name.from.file = TRUE,
    ext.case = "upper")

# fix names
dat$Name[grep("cour", dat$Name)] <- "PB_courtship"

2 Stats

2.1 Latency

# get response times
resp_time_list <- lapply(unique(dat$sound.files), function(x) {

    # print(x)
    X <- dat[dat$sound.files == x, ]
    X <- X[order(X$start), ]
    resp_indices <- which(X$Name == "fem_resp")

    X_stimuli <- X[X$Name != "fem_resp", ]

    if (length(resp_indices) > 0) {
        stimuli.latency <- NULL
        stimuli <- NULL
        stimuly.index <- NULL

        for (i in resp_indices) {

            start <- X[i, "start"]
            time_diffs <- start - X_stimuli$end
            time_diffs <- time_diffs[time_diffs >= 0]

            if (i < min(which(X$Name != "fem_resp"))) {
                stimuli[length(stimuli) + 1] <- NA
                stimuli.latency[length(stimuli.latency) + 1] <- start
                stimuly.index[length(stimuly.index) + 1] <- NA
            } else {
                stimuli[length(stimuli) + 1] <- X_stimuli$Name[which.min(time_diffs)]
                stimuli.latency[length(stimuli.latency) + 1] <- min(time_diffs)
                stimuly.index[length(stimuly.index) + 1] <- X_stimuli$selec[which.min(time_diffs)]
            }
        }

        out_df <- data.frame(sound.files = x, stimuli.latency = stimuli.latency,
            stimuli = stimuli, stimuly.index = stimuly.index)
    } else out_df <- data.frame(sound.files = x, stimuli.latency = NA,
        stimuli = NA, stimuly.index = NA)

    # remove those that responded to the same stimuli
    out_df <- out_df[!duplicated(out_df$stimuly.index) | is.na(out_df$stimuly.index),
        ]

    return(out_df)
})

resp_time_df <- do.call(rbind, resp_time_list)

resp_time_df <- resp_time_df[!is.na(resp_time_df$stimuli), ]

write.csv(resp_time_df, "./data/processed/resp_time.csv")

3 Plots

dat$sing.event <- dat$sound.files
dat$indiv <- ifelse(dat$Name == "fem_resp", "fem_resp", "stimuli")

agg <- aggregate(indiv ~ sing.event, data = dat, FUN = function(x) length(unique(x)))

dat_plot <- dat[!dat$sing.event %in% agg$sing.event[agg$indiv == 1],
    ]

ggs <- warbleR::plot_coordination(dat_plot, img = FALSE)

ggs <- lapply(ggs, function(x) x + theme_classic(base_size = 20) +
    theme(legend.position = "none"))

ggs

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

[[19]]

[[20]]

[[21]]

[[22]]

[[23]]

[[24]]

[[25]]

[[26]]

[[27]]

[[28]]

[[29]]

[[30]]

[[31]]

[[32]]

[[33]]

[[34]]

[[35]]

[[36]]

[[37]]

[[38]]

[[39]]

[[40]]

[[41]]

[[42]]

[[43]]

[[44]]

[[45]]

[[46]]

[[47]]

[[48]]

[[49]]

3.0.1 Add metadata

metadata <- read.csv("./data/raw/female_playbacks_datasheet.csv")
metadata$file_name <- tolower(metadata$file_name)

resp_time_df$experiment <- sapply(tolower(resp_time_df$sound.files),
    function(x) {

        metadata$stimuli[metadata$file_name == x]
    })

resp_time_df$experiment <- gsub("+", "|", resp_time_df$experiment,
    fixed = TRUE)

resp_time_df$location <- sapply(tolower(resp_time_df$sound.files),
    function(x) {

        metadata$location[metadata$file_name == x]
    })


resp_time_df$ID <- sapply(tolower(resp_time_df$sound.files), function(x) {

    metadata$female_ID[metadata$file_name == x]
})

# correct by stimuli level
resp_time_df$norm.latency <- ifelse(resp_time_df$experiment == "female call",
    resp_time_df$stimuli.latency * 24/6, resp_time_df$stimuli.latency)

# remove outlier
resp_time_df <- resp_time_df[resp_time_df$norm.latency < 100, ]

3.0.2 Explore data

aggregate(norm.latency ~ experiment, data = resp_time_df, FUN = mean)

experiment	norm.latency
female call	27.7992
female call | male courtship	25.8250
male courtship	25.0395
male trill	29.6180
male trill | female call	28.8122

ggplot(resp_time_df,
aes(
x = experiment,
y = stimuli.latency,
color = experiment,
fill = experiment
)) +
# add half-violin from {ggdist} package
ggdist::stat_halfeye(
# fill = fill_color,
alpha = 0.5,
# custom bandwidth
adjust = .5,
# adjust height
width = .6,
.width = 0,
# move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(# fill = fill_color,
width = .15,
# remove outliers
outlier.shape = NA) +
# add justified jitter from the {gghalves} package
gghalves::geom_half_point(
# color = fill_color,
# draw jitter on the left
side = "l",
# control range of jitter
range_scale = .4,
# add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
scale_color_viridis_d(option = "G", begin = 0.3, end = 0.8) +
scale_fill_viridis_d(option = "G", begin = 0.3, end = 0.8, alpha = 0.6) +
theme_classic(base_size = 20) +
theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(title = "Raw latency from previous stimuli by experiment", y = "Latency from previous stimuli", x = "Experiment type")

ggplot(resp_time_df,
aes(
x = experiment,
y = norm.latency,
color = experiment,
fill = experiment
)) +
# add half-violin from {ggdist} package
ggdist::stat_halfeye(
# fill = fill_color,
alpha = 0.5,
# custom bandwidth
adjust = .5,
# adjust height
width = .6,
.width = 0,
# move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(# fill = fill_color,
width = .15,
# remove outliers
outlier.shape = NA) +
# add justified jitter from the {gghalves} package
gghalves::geom_half_point(
# color = fill_color,
# draw jitter on the left
side = "l",
# control range of jitter
range_scale = .4,
# add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
scale_color_viridis_d(option = "G", begin = 0.3, end = 0.8) +
scale_fill_viridis_d(option = "G", begin = 0.3,
end = 0.8,
alpha = 0.6) +
theme_classic(base_size = 20) +
    # ylim(c(0, 75)) +
theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(title = "Standardized latency from previous stimuly by experiment", y = "Standardized latency", x = "Experiment type")

3.0.3 Regression models

3.0.3.1 Raw latency

Main model: \[ \text{log(latency + 1)} \sim \text{experiment} + \text{location} + (1 \mid \text{sound.files}) + (1 \mid \text{individual}) \] Interaction model:

\[ \text{log(latency + 1)} \sim \text{experiment} * \text{location} + (1 \mid \text{sound.files}) + (1 \mid \text{individual}) \]

mod <- brm(formula = log(stimuli.latency + 1) ~ experiment + location +
    (1 | sound.files) + (1 | ID), iter = iter, thin = 1, data = resp_time_df,
    family = gaussian(), silent = 2, chains = chains, cores = chains,
    control = list(adapt_delta = 0.99, max_treedepth = 15), backend = "cmdstanr",
    file_refit = "always", file = "./data/processed/raw_latency_by_experiment")

mod <- add_criterion(mod, criterion = c("loo"))

int_mod <- brm(formula = log(stimuli.latency + 1) ~ experiment * location +
    (1 | sound.files) + (1 | ID), iter = iter, thin = 1, data = resp_time_df,
    family = gaussian(), silent = 2, chains = chains, cores = chains,
    control = list(adapt_delta = 0.99, max_treedepth = 15), backend = "cmdstanr",
    file_refit = "always", file = "./data/processed/raw_latency_by_experiment_interaction")

int_mod <- add_criterion(int_mod, criterion = c("loo"))

null_mod <- brm(formula = log(stimuli.latency + 1) ~ 1 + (1 | sound.files) +
    (1 | ID), iter = iter, thin = 1, data = resp_time_df, family = gaussian(),
    silent = 2, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), backend = "cmdstanr", file_refit = "always",
    file = "./data/processed/null_raw_latency_by_experiment")

null_mod <- add_criterion(null_mod, criterion = c("loo"))

3.0.3.1.1 Results

3.0.3.1.1.1 Model performance vs null model

mod <- readRDS("./data/processed/raw_latency_by_experiment.rds")

int_mod <- readRDS("./data/processed/raw_latency_by_experiment_interaction.rds")

null_mod <- readRDS("./data/processed/null_raw_latency_by_experiment.rds")

loo_diff <- loo::loo_compare(mod, int_mod, null_mod)

loo_diff

         elpd_diff se_diff
int_mod    0.0       0.0  
mod       -3.5       2.9  
null_mod -14.7       4.1  

3.0.3.1.1.2 Model fit

extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7],
    trace.palette = my.viridis, highlight = TRUE, remove.intercepts = TRUE,
    print.name = FALSE)

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	Intercept-student_t(3, 2.8, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5)	log(stimuli.latency + 1) ~ experiment + location + (1 | sound.files) + (1 | ID)	10000	4	1	5000	0 (0%)	0	11241.1	11175.3	1206488579

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_experimentfemalecall|malecourtship	1.086	0.772	1.411	1	11528.5	12097.6
b_experimentmalecourtship	1.144	0.825	1.470	1	11856.3	13120.2
b_experimentmaletrill	1.458	1.133	1.798	1	11760.3	11829.1
b_experimentmaletrill|femalecall	1.253	0.938	1.584	1	11295.1	11727.8
b_locationField	-0.553	-0.928	-0.216	1	11241.1	11175.3

# Generate EMMs for the 'experiment' variable
emms_experiment <- emmeans(mod, ~experiment)

# Get pairwise contrasts between levels of 'experiment'
contrast_experiment <- contrast(emms_experiment, method = "pairwise")

# Display the contrasts
summ_contrasts <- summary(contrast_experiment)

names(summ_contrasts) <- c("Contrasts", "Estimate", "l-95% CI", "u-95% CI")


# print estimates
summ_contrasts <- html_format_coef_table(summ_contrasts, fill = viridis(10)[7],
    highlight = TRUE)

summ_contrasts

	Contrasts	Estimate	l-95% CI	u-95% CI
1	female call - female call | male courtship	-1.083	-1.406	-0.768
2	female call - male courtship	-1.143	-1.474	-0.831
3	female call - male trill	-1.456	-1.786	-1.122
4	female call - male trill | female call	-1.251	-1.566	-0.922
5	female call | male courtship - male courtship	-0.057	-0.403	0.290
6	female call | male courtship - male trill	-0.371	-0.740	-0.030
7	female call | male courtship - male trill | female call	-0.166	-0.512	0.175
8	male courtship - male trill	-0.314	-0.657	0.050
9	male courtship - male trill | female call	-0.109	-0.447	0.245
10	male trill - male trill | female call	0.203	-0.151	0.568

3.0.3.1.1.3 Posterior predictive checks

custom_ppc(fit = mod)

3.0.3.2 Standardized latency

Main model: \[ \text{log(standardized latency + 1)} \sim \text{experiment} + \text{location} + (1 \mid \text{sound.files}) + (1 \mid \text{individual}) \] Interaction model:

\[ \text{log(standardized latency + 1)} \sim \text{experiment} * \text{location} + (1 \mid \text{sound.files}) + (1 \mid \text{individual}) \]

mod <- brm(formula = log(norm.latency + 1) ~ experiment + location +
    (1 | sound.files) + (1 | ID), iter = iter, thin = 1, data = resp_time_df,
    family = gaussian(), silent = 2, chains = chains, cores = chains,
    control = list(adapt_delta = 0.99, max_treedepth = 15), backend = "cmdstanr",
    file_refit = "always", file = "./data/processed/norm_latency_by_experiment")

mod <- add_criterion(mod, criterion = c("loo"))

int_mod <- brm(formula = log(norm.latency + 1) ~ experiment * location +
    (1 | sound.files) + (1 | ID), iter = iter, thin = 1, data = resp_time_df,
    family = gaussian(), silent = 2, chains = chains, cores = chains,
    control = list(adapt_delta = 0.99, max_treedepth = 15), backend = "cmdstanr",
    file_refit = "always", file = "./data/processed/norm_latency_by_experiment_interaction")

int_mod <- add_criterion(int_mod, criterion = c("loo"))

null_mod <- brm(formula = log(norm.latency + 1) ~ 1 + (1 | sound.files) +
    (1 | ID), iter = iter, thin = 1, data = resp_time_df, family = gaussian(),
    silent = 2, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), backend = "cmdstanr", file_refit = "always",
    file = "./data/processed/null_norm_latency_by_experiment")

null_mod <- add_criterion(null_mod, criterion = c("loo"))

3.0.3.2.1 Results

3.0.3.2.1.1 Model performance vs null model

mod <- readRDS("./data/processed/norm_latency_by_experiment.rds")

int_mod <- readRDS("./data/processed/norm_latency_by_experiment_interaction.rds")

null_mod <- readRDS("./data/processed/null_norm_latency_by_experiment.rds")

loo_diff <- loo::loo_compare(mod, int_mod, null_mod)

loo_diff

         elpd_diff se_diff
int_mod   0.0       0.0   
mod      -2.6       2.8   
null_mod -2.7       3.8   

3.0.3.2.1.2 Model fit

extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7],
    trace.palette = my.viridis, highlight = TRUE, remove.intercepts = TRUE,
    print.name = FALSE)

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	Intercept-student_t(3, 3.3, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5)	log(norm.latency + 1) ~ experiment + location + (1 | sound.files) + (1 | ID)	10000	4	1	5000	0 (0%)	0	10534.4	10549.4	1188285862

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_experimentfemalecall|malecourtship	-0.135	-0.460	0.195	1	13378.5	13480.6
b_experimentmalecourtship	-0.078	-0.406	0.261	1	13628.6	13093.7
b_experimentmaletrill	0.234	-0.097	0.576	1	12946.8	13707.1
b_experimentmaletrill|femalecall	0.031	-0.293	0.366	1	13092.7	12604.4
b_locationField	-0.552	-0.948	-0.206	1	10534.4	10549.4

# Generate EMMs for the 'experiment' variable
emms_experiment <- emmeans(mod, ~experiment)

# Get pairwise contrasts between levels of 'experiment'
contrast_experiment <- contrast(emms_experiment, method = "pairwise")

# Display the contrasts
summ_contrasts <- summary(contrast_experiment)

names(summ_contrasts) <- c("Contrasts", "Estimate", "l-95% CI", "u-95% CI")


# print estimates
summ_contrasts <- html_format_coef_table(summ_contrasts, fill = viridis(10)[7],
    highlight = TRUE)

summ_contrasts

	Contrasts	Estimate	l-95% CI	u-95% CI
1	female call - female call | male courtship	0.135	-0.200	0.453
2	female call - male courtship	0.080	-0.253	0.412
3	female call - male trill	-0.234	-0.573	0.099
4	female call - male trill | female call	-0.029	-0.359	0.299
5	female call | male courtship - male courtship	-0.057	-0.418	0.309
6	female call | male courtship - male trill	-0.370	-0.734	-0.001
7	female call | male courtship - male trill | female call	-0.164	-0.523	0.204
8	male courtship - male trill	-0.312	-0.693	0.057
9	male courtship - male trill | female call	-0.109	-0.489	0.248
10	male trill - male trill | female call	0.204	-0.178	0.569

3.0.3.2.1.3 Posterior predictive checks

custom_ppc(fit = mod)

3.1 Vocal output

# get response times
vocal_output_list <- lapply(unique(dat$sound.files), function(x) {

    # print(x)
    X <- dat[dat$sound.files == x, ]
    X <- X[order(X$start), ]
    counts <- sum(X$Name == "fem_resp")
    out_df <- data.frame(sound.files = x, counts = counts)

    return(out_df)
})

vocal_output_df <- do.call(rbind, vocal_output_list)

3.1.1 Add metadata

metadata <- read.csv("./data/raw/female_playbacks_datasheet.csv")
metadata$file_name <- tolower(metadata$file_name)

vocal_output_df$experiment <- sapply(tolower(vocal_output_df$sound.files),
    function(x) {

        metadata$stimuli[metadata$file_name == x]
    })

vocal_output_df$experiment <- gsub("+", "|", vocal_output_df$experiment,
    fixed = TRUE)

vocal_output_df$location <- sapply(tolower(vocal_output_df$sound.files),
    function(x) {

        metadata$location[metadata$file_name == x]
    })

vocal_output_df$ID <- sapply(tolower(vocal_output_df$sound.files),
    function(x) {

        metadata$female_ID[metadata$file_name == x]
    })

# correct by stimuli level
vocal_output_df$norm.counts <- ifelse(vocal_output_df$experiment ==
    "female call", vocal_output_df$counts * 6/24, vocal_output_df$counts)

3.1.2 Explore data

aggregate(counts ~ experiment, data = vocal_output_df, FUN = sum)

experiment	counts
female call	97
female call | male courtship	70
male courtship	71
male trill	72
male trill | female call	70

ggplot(vocal_output_df,
aes(
x = experiment,
y = counts,
color = experiment,
fill = experiment
)) +
# add half-violin from {ggdist} package
ggdist::stat_halfeye(
# fill = fill_color,
alpha = 0.5,
# custom bandwidth
adjust = .5,
# adjust height
width = .6,
.width = 0,
# move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(# fill = fill_color,
width = .15,
# remove outliers
outlier.shape = NA) +
# add justified jitter from the {gghalves} package
gghalves::geom_half_point(
# color = fill_color,
# draw jitter on the left
side = "l",
# control range of jitter
range_scale = .4,
# add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
scale_color_viridis_d(option = "G", begin = 0.3, end = 0.8) +
scale_fill_viridis_d(option = "G", begin = 0.3,
end = 0.8,
alpha = 0.6) +
theme_classic(base_size = 20) +
    # ylim(c(0, 25)) +
theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(title = "Call count by experiment", y = "Call counts", x = "Experiment type")

ggplot(vocal_output_df,
aes(
x = experiment,
y = norm.counts,
color = experiment,
fill = experiment
)) +
# add half-violin from {ggdist} package
ggdist::stat_halfeye(
# fill = fill_color,
alpha = 0.5,
# custom bandwidth
adjust = .5,
# adjust height
width = .6,
.width = 0,
# move geom to the cright
justification = -.2,
point_colour = NA
) +
geom_boxplot(# fill = fill_color,
width = .15,
# remove outliers
outlier.shape = NA) +
# add justified jitter from the {gghalves} package
gghalves::geom_half_point(
# color = fill_color,
# draw jitter on the left
side = "l",
# control range of jitter
range_scale = .4,
# add some transparency
alpha = .5,
transformation = ggplot2::position_jitter(height = 0)
) +
scale_color_viridis_d(option = "G", begin = 0.3, end = 0.8) +
scale_fill_viridis_d(option = "G", begin = 0.3,
end = 0.8,
alpha = 0.6) +
theme_classic(base_size = 20) +
    # ylim(c(0, 25)) +
theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(title = "Standardized call count by experiment", y = "Call counts", x = "Experiment type")

3.1.3 Regression models

3.1.3.1 Raw counts

Main model: \[ \text{log(call count + 1)} \sim \text{experiment} + \text{location} + (1 \mid \text{individual}) \] Interaction model:

\[ \text{log(call count + 1)} \sim \text{experiment} * \text{location} + (1 \mid \text{individual}) \]

mod <- brm(formula = log(counts + 1) ~ experiment + location + (1 |
    ID), iter = iter, thin = 1, data = vocal_output_df, family = gaussian(),
    silent = 2, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), backend = "cmdstanr", file_refit = "always",
    file = "./data/processed/counts_by_experiment")

mod <- add_criterion(mod, criterion = c("loo"))

int_mod <- brm(formula = log(counts + 1) ~ experiment * location +
    (1 | ID), iter = iter, thin = 1, data = vocal_output_df, family = gaussian(),
    silent = 2, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), backend = "cmdstanr", file_refit = "always",
    file = "./data/processed/counts_by_experiment_interaction")

int_mod <- add_criterion(int_mod, criterion = c("loo"))

null_mod <- brm(formula = log(counts + 1) ~ 1 + (1 | ID), iter = iter,
    thin = 1, data = vocal_output_df, family = gaussian(), silent = 2,
    chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), backend = "cmdstanr", file_refit = "always",
    file = "./data/processed/null_counts_by_experiment")

null_mod <- add_criterion(null_mod, criterion = c("loo"))

3.1.3.2 Results

3.1.3.2.1 Model performance vs null model

mod <- readRDS("./data/processed/counts_by_experiment.rds")

int_mod <- readRDS("./data/processed/counts_by_experiment_interaction.rds")

null_mod <- readRDS("./data/processed/norm_counts_by_experiment.rds")

loo_diff <- loo::loo_compare(mod, int_mod, null_mod)

loo_diff

         elpd_diff se_diff
null_mod  0.0       0.0   
mod      -8.2       6.2   
int_mod  -9.1       5.0   

3.1.3.2.2 Model fit

extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7],
    trace.palette = my.viridis, highlight = TRUE, remove.intercepts = TRUE,
    print.name = FALSE)

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	Intercept-student_t(3, 1.9, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5)	log(counts + 1) ~ experiment + location + (1 | ID)	10000	4	1	5000	0 (0%)	0	6216.62	8477.45	38094587

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_experimentfemalecall|malecourtship	-0.170	-0.551	0.218	1	11602.49	13280.44
b_experimentmalecourtship	-0.280	-0.659	0.104	1	11634.38	12521.73
b_experimentmaletrill	-0.277	-0.652	0.111	1	12139.96	12565.89
b_experimentmaletrill|femalecall	-0.211	-0.592	0.176	1	11426.60	12515.15
b_locationField	0.741	-0.196	1.673	1.001	6216.62	8477.45

# Generate EMMs for the 'experiment' variable
emms_experiment <- emmeans(mod, ~experiment)

# Get pairwise contrasts between levels of 'experiment'
contrast_experiment <- contrast(emms_experiment, method = "pairwise")

# Display the contrasts
summ_contrasts <- summary(contrast_experiment)

names(summ_contrasts) <- c("Contrasts", "Estimate", "l-95% CI", "u-95% CI")


# print estimates
summ_contrasts <- html_format_coef_table(summ_contrasts, fill = viridis(10)[7],
    highlight = TRUE)

summ_contrasts

	Contrasts	Estimate	l-95% CI	u-95% CI
1	female call - female call | male courtship	0.172	-0.220	0.548
2	female call - male courtship	0.281	-0.096	0.664
3	female call - male trill	0.278	-0.109	0.653
4	female call - male trill | female call	0.212	-0.165	0.599
5	female call | male courtship - male courtship	0.109	-0.284	0.486
6	female call | male courtship - male trill	0.107	-0.276	0.481
7	female call | male courtship - male trill | female call	0.041	-0.359	0.415
8	male courtship - male trill	-0.002	-0.378	0.385
9	male courtship - male trill | female call	-0.069	-0.468	0.303
10	male trill - male trill | female call	-0.065	-0.451	0.308

3.1.3.2.3 Posterior predictive checks

custom_ppc(fit = mod)

3.1.3.3 Standardized counts

Main model: \[ \text{log(standardized call count + 1)} \sim \text{experiment} + \text{location} + (1 \mid \text{individual}) \] Interaction model:

\[ \text{log(standardized call count + 1)} \sim \text{experiment} * \text{location} + (1 \mid \text{individual}) \]

mod <- brm(formula = log(norm.counts + 1) ~ experiment + location +
    (1 | ID), iter = iter, thin = 1, data = vocal_output_df, family = gaussian(),
    silent = 2, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), backend = "cmdstanr", file_refit = "always",
    file = "./data/processed/norm_counts_by_experiment")

mod <- add_criterion(mod, criterion = c("loo"))

int_mod <- brm(formula = log(norm.counts + 1) ~ experiment * location +
    (1 | ID), iter = iter, thin = 1, data = vocal_output_df, family = gaussian(),
    silent = 2, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), backend = "cmdstanr", file_refit = "always",
    file = "./data/processed/norm_counts_by_experiment_interaction")

int_mod <- add_criterion(int_mod, criterion = c("loo"))

null_mod <- brm(formula = log(norm.counts + 1) ~ 1 + (1 | ID), iter = iter,
    thin = 1, data = vocal_output_df, family = gaussian(), silent = 2,
    chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), backend = "cmdstanr", file_refit = "always",
    file = "./data/processed/null_norm_counts_by_experiment")

null_mod <- add_criterion(null_mod, criterion = c("loo"))

3.1.3.3.1 Results

3.1.3.3.1.1 Model performance vs null model

mod <- readRDS("./data/processed/norm_counts_by_experiment.rds")

int_mod <- readRDS("./data/processed/norm_counts_by_experiment_interaction.rds")

null_mod <- readRDS("./data/processed/null_norm_counts_by_experiment.rds")

loo_diff <- loo::loo_compare(mod, int_mod, null_mod)

loo_diff

         elpd_diff se_diff
mod        0.0       0.0  
int_mod   -3.6       1.9  
null_mod -12.0       5.1  

3.1.3.3.1.2 Model fit

extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7],
    trace.palette = my.viridis, highlight = TRUE, remove.intercepts = TRUE,
    print.name = FALSE)

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	Intercept-student_t(3, 1.6, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5)	log(norm.counts + 1) ~ experiment + location + (1 | ID)	10000	4	1	5000	0 (0%)	0	5935.94	7298.32	141855414

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_experimentfemalecall|malecourtship	0.814	0.481	1.144	1	11105.29	12691.42
b_experimentmalecourtship	0.703	0.374	1.035	1	10321.71	12350.45
b_experimentmaletrill	0.707	0.371	1.038	1.001	10649.02	12886.07
b_experimentmaletrill|femalecall	0.773	0.440	1.104	1.001	9948.04	12641.66
b_locationField	0.675	-0.239	1.596	1.001	5935.94	7298.32

# Generate EMMs for the 'experiment' variable
emms_experiment <- emmeans(mod, ~experiment)

# Get pairwise contrasts between levels of 'experiment'
contrast_experiment <- contrast(emms_experiment, method = "pairwise")

# Display the contrasts
summ_contrasts <- summary(contrast_experiment)

names(summ_contrasts) <- c("Contrasts", "Estimate", "l-95% CI", "u-95% CI")


# print estimates
summ_contrasts <- html_format_coef_table(summ_contrasts, fill = viridis(10)[7],
    highlight = TRUE)

summ_contrasts

	Contrasts	Estimate	l-95% CI	u-95% CI
1	female call - female call | male courtship	-0.815	-1.140	-0.478
2	female call - male courtship	-0.703	-1.027	-0.369
3	female call - male trill	-0.707	-1.031	-0.365
4	female call - male trill | female call	-0.771	-1.111	-0.449
5	female call | male courtship - male courtship	0.110	-0.221	0.438
6	female call | male courtship - male trill	0.106	-0.224	0.429
7	female call | male courtship - male trill | female call	0.040	-0.295	0.376
8	male courtship - male trill	-0.006	-0.330	0.322
9	male courtship - male trill | female call	-0.070	-0.389	0.267
10	male trill - male trill | female call	-0.064	-0.399	0.259

3.1.3.3.1.3 Posterior predictive checks

custom_ppc(fit = mod)

Takeaways

---

Session information

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

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

other attached packages:
 [1] rstan_2.32.6        StanHeaders_2.32.10 loo_2.8.0          
 [4] cowplot_1.1.3       emmeans_1.10.3      viridis_0.6.5      
 [7] viridisLite_0.4.2   brmsish_1.0.0       brms_2.22.0        
[10] Rcpp_1.0.13         ggplot2_3.5.1       warbleR_1.1.32     
[13] NatureSounds_1.0.4  knitr_1.48          seewave_2.2.3      
[16] tuneR_1.4.7         Rraven_1.0.13      

loaded via a namespace (and not attached):
  [1] bitops_1.0-9         pbapply_1.7-2        gridExtra_2.3       
  [4] formatR_1.14         inline_0.3.19        remotes_2.5.0       
  [7] testthat_3.2.1.1     sandwich_3.1-0       rlang_1.1.4         
 [10] magrittr_2.0.3       multcomp_1.4-25      matrixStats_1.4.1   
 [13] compiler_4.4.1       reshape2_1.4.4       systemfonts_1.1.0   
 [16] vctrs_0.6.5          stringr_1.5.1        pkgconfig_2.0.3     
 [19] arrayhelpers_1.1-0   crayon_1.5.3         fastmap_1.2.0       
 [22] backports_1.5.0      labeling_0.4.3       utf8_1.2.4          
 [25] cmdstanr_0.8.1.9000  rmarkdown_2.28       ps_1.8.0            
 [28] purrr_1.0.2          xfun_0.48            jsonlite_1.8.9      
 [31] gghalves_0.1.4       tidybayes_3.0.7      parallel_4.4.1      
 [34] R6_2.5.1             stringi_1.8.4        brio_1.1.5          
 [37] estimability_1.5.1   zoo_1.8-12           bayesplot_1.11.1    
 [40] Matrix_1.7-0         splines_4.4.1        tidyselect_1.2.1    
 [43] rstudioapi_0.16.0    abind_1.4-8          yaml_2.3.10         
 [46] dtw_1.23-1           codetools_0.2-20     processx_3.8.4      
 [49] curl_5.2.3           pkgbuild_1.4.4       plyr_1.8.9          
 [52] lattice_0.22-6       tibble_3.2.1         withr_3.0.1         
 [55] bridgesampling_1.1-2 posterior_1.6.0      coda_0.19-4.1       
 [58] evaluate_1.0.1       signal_1.8-1         survival_3.7-0      
 [61] sketchy_1.0.3        proxy_0.4-27         RcppParallel_5.1.9  
 [64] ggdist_3.3.2         xml2_1.3.6           pillar_1.9.0        
 [67] tensorA_0.36.2.1     packrat_0.9.2        stats4_4.4.1        
 [70] checkmate_2.3.2      distributional_0.5.0 generics_0.1.3      
 [73] RCurl_1.98-1.16      rstantools_2.4.0     munsell_0.5.1       
 [76] scales_1.3.0         xtable_1.8-4         glue_1.8.0          
 [79] tools_4.4.1          xaringanExtra_0.8.0  mvtnorm_1.3-1       
 [82] grid_4.4.1           tidyr_1.3.1          ape_5.8             
 [85] QuickJSR_1.4.0       colorspace_2.1-1     nlme_3.1-165        
 [88] cli_3.6.3            kableExtra_1.4.0     fansi_1.0.6         
 [91] svUnit_1.0.6         svglite_2.1.3        Brobdingnag_1.2-9   
 [94] dplyr_1.1.4          V8_4.4.2             gtable_0.3.5        
 [97] fftw_1.0-9           digest_0.6.37        TH.data_1.1-2       
[100] farver_2.1.2         rjson_0.2.23         htmlwidgets_1.6.4   
[103] htmltools_0.5.8.1    lifecycle_1.0.4      MASS_7.3-61