Song loudness analysis

Long-billed hermit song amplitude

Author

Marcelo Araya-Salas

Published

February 15, 2024

1 What the code does:

Calculates all stats for song sound pressure level (i.e. amplitude) covariates

2 Next steps

Make analysis for vegetation density
Change randome slope analysis with absolute differences in before and after treatments (both in interactions and coordinations) -See if eliminating interaction videos allows a cleaner analysis (leave only calibrated songs)

2.0.0.1 Load packages

Code

## vector with package names
x <- c("pbapply", "parallel", "ggplot2", "warbleR", "Rraven", "viridis",
    "readxl", "rptR", "brms", "MuMIn", "corrplot", "lme4", "grid",
    "gridExtra", "dynaSpec", github = "maRce10/brmsish", "coda", "MCMCglmm",
    "knitr")

source("https://raw.githubusercontent.com/maRce10/sketchy/main/R/load_packages.R")

load_packages(x)

2.0.0.2 Functions and parameters

Code

# functions and parameters
knitr::opts_knit$set(root.dir = normalizePath(".."))

knitr::opts_chunk$set(dpi = 50, fig.width = 12, warning = FALSE, message = FALSE)

# ggplot2 theme theme_set(theme_classic(base_size = 20))

cut_path <- "./data/raw/cuts"

treatments <- c("Calibration", "Regular_singing", "Coordination",
    "After_chase", "Before_playback", "After_playback", "Before_interaction",
    "After_interaction", "Before_noise", "After_noise")

# iterations for MCMCglmm models
itrns <- 1e+05

# functions from https://rdrr.io/rforge/rptR/src/R/rpt.mcmcLMM.R
rpt.mcmcLMM <- function(y, groups, CI = 0.95, prior = NULL, verbose = FALSE,
    ...) {
    # initial checks
    if (length(y) != length(groups))
        stop("y and group are of unequal length")
    # preparation
    groups <- factor(groups)
    if (is.null(prior))
        prior <- list(R = list(V = 1, n = 0.1), G = list(G1 = list(V = 1,
            n = 0.1)))
    # point estimation according to model 8 and equation 9
    mod <- MCMCglmm(y ~ 1, random = ~groups, family = "gaussian",
        data = data.frame(y = y, groups = groups), prior = prior,
        verbose = verbose, ...)
    var.a <- mod$VCV[, "groups"]
    var.e <- mod$VCV[, "units"]
    postR <- var.a/(var.a + var.e)
    # point estimate
    R <- posterior.mode(postR)
    # credibility interval estimation from paterior distribution
    CI.R <- coda::HPDinterval(postR, CI)[1, ]
    se <- sd(postR)
    # 'significance test'
    P <- NA
    res = list(call = match.call(), datatype = "Gaussian", method = "LMM.MCMC",
        CI = CI, R = R, CI.R = CI.R, se = se, P = P)
    # class(res) <- 'rpt'
    return(res)
}


## print Gelman-Rubin convergence statistics, plots traces and
## autocorrelations
mcmc_diagnostics <- function(rep_mods_list) {

    for (w in 1:length(rep_mods_list)) {
        mod_name <- names(rep_mods_list)[w]

        if (mod_name == "1")
            mod_name <- "Null"
        print(paste("model:", mod_name))

        Y <- lapply(rep_mods_list[[w]], "[[", "Sol")

        ## add global plots and gelman test gelman_diagnostic
        gel_diag <- as.data.frame(gelman.diag(mcmc.list(Y))$psrf)

        # add estimate as column
        gel_diag$estimate <- rownames(gel_diag)

        # reorder columns
        gel_diag <- gel_diag[, c(3, 1, 2)]

        par(mfrow = c(1, 1))

        # plot table
        grid.newpage()
        grid.draw(tableGrob(gel_diag, rows = NULL, theme = ttheme_default(base_size = 25)))

        par(mfrow = c(1, 4))

        traceplot(Y, col = adjustcolor(c("yellow", "blue", "red"),
            alpha.f = 0.6))

        autocorr.plot(x = Y[[1]], auto.layout = FALSE, lwd = 4, col = "red")
    }

}

Code

amp <- read.csv("./output/calibrated_amplitude_all_songs.csv")
amp$Treatment[amp$Treatment == "Regular_sining"] <- "Regular_singing"

sound_meter_calib <- readxl::read_excel("./data/raw/calibrate_sound_meter_based_on_new_sound_meter.xlsx")

# calibrate SPL based on new sound meter
mean_spl_diff <- mean(sound_meter_calib$nuevo[1:19] - sound_meter_calib$viejo[1:19])

amp$cal.spl <- amp$cal.spl + mean_spl_diff

# extrapolate SPL at 1m using the inverse square law function
amp$cal.spl.1m <- amp$cal.spl + 20 * log10(40/100)

3 Repeatability

We estimated repeatability using linear mixed models:

Code

# keep only morning recs, only regular singing and only from
# 2021
amp_morn <- amp[amp$period == "morning" & amp$year == 2021 & amp$Treatment ==
    "Regular_singing", ]

# keep only males recorded at least twice
count_sf <- table(amp_morn$ID[!duplicated(amp_morn$org.sound.file)])

amp_morn_rep <- amp_morn[amp_morn$ID %in% names(count_sf)[count_sf >
    1], ]

max_quantile <- c(0, 1, seq(10, 90, by = 10))/100
outliers <- c(0, 1, 2, 5, 10, 20)/100

rep_grid <- expand.grid(max_quantile = max_quantile, outliers = outliers,
    only.low.outliers = c(TRUE, FALSE))

repts <- pblapply(1:nrow(rep_grid), cl = 6, function(x) {

    quant_subet <- lapply(unique(amp_morn_rep$org.sound.file), function(y) {

        X <- amp_morn_rep[amp_morn_rep$org.sound.file == y, ]

        # remove outliers
        outlier_quant <- quantile(X$cal.spl.1m, c(rep_grid$outliers[x],
            1 - rep_grid$outliers[x]))

        if (!rep_grid$only.low.outliers[x])
            X <- X[X$cal.spl.1m >= outlier_quant[1] & X$cal.spl.1m <=
                outlier_quant[2], ] else X <- X[X$cal.spl.1m >= outlier_quant[1], ]

        # quantlie for each max quantile
        quant <- quantile(X$cal.spl.1m, probs = 1 - rep_grid$max_quantile[x])

        # subset
        X <- X[X$cal.spl.1m >= quant, ]
    })

    quant_subet <- do.call(rbind, quant_subet)

    # frequentist
    out <- rpt(cal.spl.1m ~ (1 | ID), grname = "ID", data = quant_subet,
        datatype = "Gaussian", npermut = 0, nboot = 100)

    # bayesian out <- rpt.mcmcLMM(y = quant_subet$cal.spl.1m,
    # groups = quant_subet$ID, nitt = itrns)

    out_df <- data.frame(rep_grid$max_quantile[x], rep_grid$outliers[x],
        rep_grid$only.low.outliers[x], out$R, out$CI_emp[1], out$CI_emp[2])

    return(out_df)
})

repts_df <- do.call(rbind, repts)

names(repts_df) <- c("max_quantile", "outliers", "only.low.outliers",
    "repeatability", "lowCI", "hiCI")

# write.csv(repts_df, './output/repeatability_optimization.csv',
# row.names = FALSE)

Different subsets of data for upper quantiles (x axis)
Pannels show the proportion of outliers removed

Code

repts_df <- read.csv("./output/repeatability_optimization.csv")

pd <- position_dodge(width = 0.1)
ggplot(data = repts_df, aes(x = 1 - max_quantile, y = repeatability,
    color = only.low.outliers, group = only.low.outliers)) + geom_hline(yintercept = 0.5,
    col = adjustcolor("red", alpha.f = 0.5)) + geom_point(size = 2,
    position = pd) + geom_errorbar(width = 0.05, aes(ymin = lowCI,
    ymax = hiCI), position = pd) + scale_color_viridis(discrete = TRUE,
    begin = 0.2, end = 0.8, alpha = 0.7) + geom_line(position = pd) +
    labs(y = "Repeatability", x = "Upper quantile used") + ylim(c(0,
    1)) + xlim(c(1, 0)) + facet_wrap(~outliers, scales = "fixed") +
    theme_classic(base_size = 24)

Code

repts_df$range <- repts_df$hiCI - repts_df$lowCI

# repts_df[order(repts_df$repeatability), c('max_quantile',
# 'outliers', 'only.low.outliers', 'repeatability', 'range')]

repts_df <- repts_df[order(repts_df$repeatability, decreasing = TRUE),
    ]

head(repts_df, 23)

	max_quantile	outliers	only.low.outliers	repeatability	lowCI	hiCI	range
130	0.7	0.20	FALSE	0.5168506	0.1817020	0.7240848	0.5423828
116	0.4	0.10	FALSE	0.5166605	0.1858151	0.7234608	0.5376457
132	0.9	0.20	FALSE	0.5159134	0.1080931	0.7197176	0.6116245
115	0.3	0.10	FALSE	0.5157500	0.2128287	0.7303455	0.5175168
104	0.3	0.05	FALSE	0.5156010	0.1834339	0.7242149	0.5407810
119	0.7	0.10	FALSE	0.5153219	0.2520312	0.7077252	0.4556940
131	0.8	0.20	FALSE	0.5152836	0.1372292	0.6673814	0.5301522
127	0.4	0.20	FALSE	0.5151673	0.1778775	0.7069726	0.5290951
118	0.6	0.10	FALSE	0.5151024	0.2096563	0.7167910	0.5071348
120	0.8	0.10	FALSE	0.5150953	0.1912322	0.6905403	0.4993081
129	0.6	0.20	FALSE	0.5150763	0.1747855	0.7476081	0.5728226
107	0.6	0.05	FALSE	0.5150239	0.1327933	0.6867971	0.5540038
125	0.2	0.20	FALSE	0.5148523	0.1525131	0.7216677	0.5691546
117	0.5	0.10	FALSE	0.5146927	0.1537053	0.7293507	0.5756453
128	0.5	0.20	FALSE	0.5143596	0.2197766	0.7512893	0.5315128
121	0.9	0.10	FALSE	0.5142971	0.2000092	0.6690270	0.4690178
105	0.4	0.05	FALSE	0.5141632	0.2041509	0.7079257	0.5037748
108	0.7	0.05	FALSE	0.5139307	0.2218817	0.7142121	0.4923304
114	0.2	0.10	FALSE	0.5139231	0.2163404	0.7265155	0.5101751
106	0.5	0.05	FALSE	0.5139198	0.1483935	0.7195049	0.5711114
126	0.3	0.20	FALSE	0.5137413	0.1497227	0.6850615	0.5353388
109	0.8	0.05	FALSE	0.5136567	0.1744174	0.7230653	0.5486479
39	0.4	0.05	TRUE	0.5133548	0.1709177	0.7022720	0.5313543

Code

# keep only morning recs, only regular singing and only from
# 2021
amp_morn <- amp[amp$period == "morning" & amp$year == 2021 & amp$Treatment ==
    "Regular_singing", ]
# keep only males recorded at least twice
count_sf <- table(amp_morn$ID[!duplicated(amp_morn$org.sound.file)])

amp_morn_rep <- amp_morn[amp_morn$ID %in% names(count_sf)[count_sf >
    1], ]

amp_morn_rep <- amp_morn

quant_subet <- lapply(unique(amp_morn_rep$org.sound.file), function(y) {
    X <- amp_morn_rep[amp_morn_rep$org.sound.file == y, ]

    # remove outliers
    outlier_quant <- quantile(X$cal.spl.1m, c(0.1, 1 - 0.1))

    X <- X[X$cal.spl.1m >= outlier_quant[1], ]

    # quantlie for each max quantile
    quant <- quantile(X$cal.spl.1m, probs = 1 - 0.2)

    # subset
    X <- X[X$cal.spl.1m >= quant, ]
})

quant_subet <- do.call(rbind, quant_subet)

ggplot(quant_subet, aes(x = factor(ID), y = cal.spl.1m, width = 2)) +
    geom_violin(fill = viridis(10, alpha = 0.5)[4]) + theme_classic(base_size = 30) +
    labs(x = "Individual", y = "Sound pressure level (dB)") + theme(axis.text.x = element_text(angle = 90))

Code

# ggsave('./output/individual_differences_in_spl.jpeg', width =
# 14, height = 7)

Excluding lowest values increases repeatability
Excluding outliers doesn’t have a strong effect but removing 0.1 outliers produced the highest repeatability
Removing lower tail outliers or both works similarly (we excluded only low tail outliers to keep more data)

Subsetting the 0.2 upper quartile after excluding below 0.1 or above 0.9 outliers:

Code

# compose variable to remove low values and outliers based on
# repeatabiliy
amp$osf.treat <- paste(amp$org.sound.file, amp$Treatment, sep = "-")

rm_outlier_amp_l <- lapply(unique(amp$osf.treat), function(y) {
    X <- amp[amp$osf.treat == y, ]

    # remove outliers
    outlier_quant <- quantile(X$cal.spl.1m, c(0.1, 0.9))
    X <- X[X$cal.spl.1m >= outlier_quant[1] & X$cal.spl.1m <= outlier_quant[2],
        ]

    # quantlie for each max quantile (0.6 was selected due to
    # high repeatability)
    quant <- quantile(X$cal.spl.1m, probs = 0.2)

    # subset
    X <- X[X$cal.spl.1m >= quant, ]
})

rm_outlier_amp <- do.call(rbind, rm_outlier_amp_l)

4 Final data set

Code

id.count <- as.data.frame(table(rm_outlier_amp$ID))

names(id.count) <- c("ID", "n")

Total number of observations: 16153 (16110 excluding unidentified individuals)
Median and range of observations per individual: 178 (5-2242)

Code

id.count

ID	n
0	43
9	5
36	9
54	7
108	8
124	24
146	12
176	19
178	8
179	8
384	209
390	151
391	175
395	123
397	1507
398	449
400	919
402	85
403	474
413	881
415	977
416	705
417	1706
419	2242
422	607
423	744
424	1427
432	414
433	53
435	45
437	442
438	981
444	181
446	77
448	332
449	104

5 Hypothesis testing analysis

We use bayesian generalized linear models with individual ID and song type as random factors (random intercepts) except indicated. An intercept-only (null) model is also included in the analyses. A loop is used to run these models. Each model is replicated three times with starting values sampled from a Z-distribution (“start” argument in MCMCglmm()). Diagnostic plots for MCMC model performance are shown at the end of this report.

6 SPL vs day period (morning / afternoon)

Only for individuals recorded in both periods

Period of the day as predictor: \[SPL \sim + period + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

Code

iter <- 1000
chains <- 4
priors <- c(prior(normal(0, 2), class = "b"))

# subset data
period_data <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" &
    rm_outlier_amp$year == 2021 & rm_outlier_amp$ID %in% unique(rm_outlier_amp$ID[rm_outlier_amp$period ==
    "afternoon"]), ]

# model SPL by body size
period_formulas <- c("cal.spl.1m ~ 1 + (1|ID) + (1|songtype)", "cal.spl.1m ~ period + (1|ID) + (1|songtype)")

period_mods <- pblapply(period_formulas, function(x) {

    mod <- brm(as.formula(x), data = period_data, iter = iter, chains = chains,
        cores = chains)

    return(mod)
})

names(period_mods) <- gsub("cal.spl.1m ~ ", "", period_formulas)
saveRDS(list(period_formulas = period_formulas, period_mods = period_mods),
    "./output/brms_period_models.RDS")

6.1 Model selection

Code

attach(readRDS("./output/brms_period_models.RDS"))

loo_compare(add_criterion(period_mods[[1]], "waic"), add_criterion(period_mods[[2]],
    "waic"), criterion = "waic", model_names = period_formulas)

                                            elpd_diff se_diff
cal.spl.1m ~ period + (1|ID) + (1|songtype)   0.0       0.0  
cal.spl.1m ~ 1 + (1|ID) + (1|songtype)      -11.0       5.0

6.2 Best model summary

Code

# fixed effects with HPD intervals
best_mod_period <- period_mods[[which(names(period_mods) == row.names(model_selection)[1])]][[1]]

sm <- as.data.frame(summary(best_mod_period)$solutions[, -5])

sm

Code

dat_period <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" &
    rm_outlier_amp$ID %in% unique(rm_outlier_amp$ID[rm_outlier_amp$period ==
        "afternoon"]), ]

dat_period$period <- factor(dat_period$period, levels = c("morning",
    "afternoon"))

ggplot(dat_period, aes(x = as.factor(ID), y = cal.spl.1m, color = period,
    fill = period)) + geom_violin(position = pd) + labs(x = "Individual",
    y = "Sound pressure level (dB)", fill = "Period") + scale_color_viridis_d(begin = 0.1,
    end = 0.8, alpha = 0.6) + scale_fill_viridis_d(begin = 0.1, end = 0.8,
    alpha = 0.6) + guides(color = FALSE) + theme_classic(base_size = 24)

Effect size posterior distribution

Code

# simplify names
colnames(best_mod_period$Sol) <- c("Intercept", "Morning vs afternoon")

# stack posteriors
Y <- stack(as.data.frame(best_mod_period$Sol[, -1]))
Y$ind <- "Morning vs afternoon"

# plot posteriors
ggplot(Y, aes(y = values, x = ind)) + geom_hline(yintercept = 0, col = "red",
    lty = 2) + geom_violin(color = "gray40", fill = viridis(n = 1,
    alpha = 0.25, begin = 0.4)) + labs(y = "Effect size posterior distribution",
    x = "Predictor") + theme_classic(base_size = 24) + theme(axis.text.x = element_text(angle = 45,
    vjust = 0.9, hjust = 1))

Morning songs have significantly higher SPL than afternoon songs within individuals

7 SPL vs condition and habitat structure

Only for individuals recorded in the morning (afternoon recordings were excluded)
Condition parameters: body size (PC1 from a PCA on total culmen, flattened wing length, body mass and central rectrix), lifting power and vegetation density

Body size as predictor: \[SPL \sim + body\ size + (1 | individual) + (1 | song\ type)\]
Lifting power as predictor: \[SPL \sim + lifting\ power + (1 | individual) + (1 | song\ type)\]
Vegetation density as predictor: \[SPL \sim + vegetation\ density + (1 | individual) + (1 | song\ type)\]
Body size and vegetation density as predictors: \[SPL \sim + body\ size + vegetation\ density + (1 | individual) + (1 | song\ type)\]
Lifting power and vegetation density as predictors: \[SPL \sim + lifting\ power + vegetation\ density + (1 | individual) + (1 | song\ type)\]
Lifting power and body size as predictors: \[SPL \sim + lifting\ power + body\ size + (1 | individual) + (1 | song\ type)\]
Lifting power, body size and vegetation density as predictors: \[SPL \sim + lifting\ power + body\ size + vegetation\ density + (1 | individual) + (1 | song\ type)\]
Interaction of lifting power and body size as predictor: \[SPL \sim + lifting\ power * body\ size + (1 | individual) + (1 | song\ type)\]
Interaction of lifting power, body size and vegetation density as predictor: \[SPL \sim + lifting\ power * body\ size * vegetation\ density + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

7.0.1 Correlation between body size parameters included in the PCA

Code

caps <- read_excel("./data/raw/LBH captures data.xlsx")

# same variables as in paper on spatial memory and body size
size_vars <- caps[caps$`Bird ID` %in% amp$ID, c("Bird ID", "Total culmen",
    "Flattened wing length", "Weight", "Central rectriz")]

# replace 419 central rectriz with mean
size_vars$`Central rectriz`[size_vars$`Bird ID` == 419] <- mean(size_vars$`Central rectriz`,
    na.rm = TRUE)

cor(size_vars[, c("Total culmen", "Flattened wing length", "Weight",
    "Central rectriz")], use = "pairwise.complete.obs")

                      Total culmen Flattened wing length     Weight
Total culmen           1.000000000           -0.13969106 0.14426785
Flattened wing length -0.139691059            1.00000000 0.08646515
Weight                 0.144267849            0.08646515 1.00000000
Central rectriz        0.007536892            0.34107800 0.10235729
                      Central rectriz
Total culmen              0.007536892
Flattened wing length     0.341077996
Weight                    0.102357292
Central rectriz           1.000000000

7.0.2 Variance explained by each PC

Code

mean_size <- aggregate(. ~ `Bird ID`, data = size_vars, FUN = mean, na.action = "na.omit")

pca <- prcomp(mean_size[, c("Total culmen", "Flattened wing length", "Weight", "Central rectriz")], scale. = TRUE)

summary(pca)

Importance of components:
                          PC1    PC2    PC3    PC4
Standard deviation     1.2488 1.0841 0.8726 0.7098
Proportion of Variance 0.3899 0.2938 0.1903 0.1259
Cumulative Proportion  0.3899 0.6837 0.8741 1.0000

Code

mean_size$PC1 <- pca$x[, 1]


rm_outlier_amp$PC1 <- sapply(1:nrow(rm_outlier_amp), function(x) mean_size$PC1[mean_size$`Bird ID` == rm_outlier_amp$ID[x]][1])

rm_outlier_amp$weight <- sapply(1:nrow(rm_outlier_amp), function(x) mean_size$Weight[mean_size$`Bird ID` == rm_outlier_amp$ID[x]][1])


### add weight lifted
lift <- read_excel("./data/raw/Uplift power experiment.xlsx", sheet = "Results")

rm_outlier_amp$lift.weight <- sapply(1:nrow(rm_outlier_amp), function(x){
  
  if(rm_outlier_amp$ID[x] == 0) weight <- NA else {
    
    if (rm_outlier_amp$year[x] == 2019)
      sub_lift <- lift[grep("\\.2019\\.", lift$Video), ] else
              sub_lift <- lift[grep("\\.2021\\.", lift$Video), ]

    
    weights <- sub_lift$Weight[sub_lift$ID == rm_outlier_amp$ID[x]]
   
    if (length(weights) > 0)
    # get the mean of the 2 highest flights
     weight <- mean(sort(weights, decreasing = TRUE)[1:2]) else
    weight <- NA
  }  
  
  return(weight)
})


# read vegetation density
veg <- read_excel(path = "./data/raw/Vegetation_density_lbh.xlsx")

veg$mean_count <- apply(veg[, c("N", "S", "W", "E")], 1, mean, na.rm = TRUE) 

# divide by 16 (max number of lines on pole)
veg$veg_density <- 1 - veg$mean_count / 16

agg_veg <- aggregate(veg_density ~ ID, veg, mean)

rm_outlier_amp$veg_density <- sapply(1:nrow(rm_outlier_amp), function(x) agg_veg$veg_density[agg_veg$ID == as.numeric(rm_outlier_amp$ID[x])][1])

Code

# subset data
body_size_data <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" &
    rm_outlier_amp$year == 2021 & !is.na(rm_outlier_amp$PC1) & !is.na(rm_outlier_amp$lift.weight) &
    rm_outlier_amp$period == "morning", ]

# model SPL by body condition
size_formulas <- c("cal.spl.1m ~ 1", "cal.spl.1m ~ PC1", "cal.spl.1m ~ lift.weight",
    "cal.spl.1m ~ lift.weight + PC1", "cal.spl.1m ~ lift.weight *  PC1",
    "cal.spl.1m ~ veg_density", "cal.spl.1m ~ PC1 + veg_density",
    "cal.spl.1m ~ lift.weight + veg_density", "cal.spl.1m ~ lift.weight + PC1 + veg_density",
    "cal.spl.1m ~ lift.weight *  PC1 * veg_density")


size_mods <- pblapply(size_formulas, cl = 8, function(x) {

    replicate(n = 3, MCMCglmm(as.formula(x), random = ~ID + songtype,
        data = body_size_data, verbose = FALSE, nitt = itrns, start = list(QUASI = FALSE)),
        simplify = FALSE)

})

names(size_mods) <- gsub("cal.spl.1m ~ ", "", size_formulas)
saveRDS(list(size_formulas = size_formulas, size_mods = size_mods),
    "./output/mcmcglmm_size_models.RDS")

# how predictors covary body_size_sub <-
# body_size_data[!duplicated(body_size_data$lift.weight), ]
# cor(body_size_sub[, c('lift.weight', 'veg_density', 'PC1',
# 'weight')]) md <- MCMCglmm(veg_density ~ lift.weight + PC1 +
# weight, random = ~ lek, data = body_size_sub, verbose = FALSE,
# nitt = itrns, start = list(QUASI = FALSE)) summary(md)

7.1 Model selection

Code

attach(readRDS("./output/mcmcglmm_size_models.RDS"))

mod_list <- lapply(size_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl.1m ~ ", "", size_formulas)

model_selection <- model.sel(mod_list, rank = "DIC")

model_selection <- as.data.frame(model_selection)

model_selection[, grep("(Intercept)|family|df|DIC|logLik", names(model_selection),
    invert = TRUE)]

	PC1	lift.weight	lift.weight:PC1	veg_density	lift.weight:veg_density	PC1:veg_density	lift.weight:PC1:veg_density	delta	weight
cal.spl ~ lift.weight * PC1	10.1358837	1.5025819	-0.6198906	NA	NA	NA	NA	0.0000000	0.10479150
cal.spl ~ lift.weight * PC1 * veg_density	-68.8026969	-2.3173798	4.3269399	-107.57043	7.443306	170.527	-10.71843	0.0215217	0.10366990
cal.spl ~ lift.weight + PC1	0.8329336	0.9772124	NA	NA	NA	NA	NA	0.0347903	0.10298440
cal.spl ~ PC1	1.0205350	NA	NA	NA	NA	NA	NA	0.0718532	0.10109352
cal.spl ~ lift.weight + PC1 + veg_density	1.2250997	1.0562670	NA	13.66817	NA	NA	NA	0.0743268	0.10096856
cal.spl ~ 1	NA	NA	NA	NA	NA	NA	NA	0.1086323	0.09925144
cal.spl ~ lift.weight + veg_density	NA	1.1660739	NA	11.13103	NA	NA	NA	0.1092599	0.09922030
cal.spl ~ lift.weight	NA	1.0661066	NA	NA	NA	NA	NA	0.1407326	0.09767116
cal.spl ~ PC1 + veg_density	1.4331449	NA	NA	13.29240	NA	NA	NA	0.1504457	0.09719796
cal.spl ~ veg_density	NA	NA	NA	10.58292	NA	NA	NA	0.2354962	0.09315125

7.2 Best model summary

None of the models provided a better fit compared to the null model

Code

# fixed effects with HPD intervals
best_mod_size <- size_mods[[which(names(size_mods) == "lift.weight + PC1")]][[1]]

sm_size <- as.data.frame(summary(best_mod_size)$solutions[, -5])

sm_size

	post.mean	l-95% CI	u-95% CI	eff.samp
(Intercept)	85.4496757	50.174598	123.015641	9164.555
lift.weight	0.9772124	-1.198775	3.346698	9213.649
PC1	0.8329336	-2.169439	3.872608	9870.742

Code

body_size_data <- rm_outlier_amp[rm_outlier_amp$Treatment == "Regular_singing" &
    rm_outlier_amp$year == 2021 & !is.na(rm_outlier_amp$PC1) & !is.na(rm_outlier_amp$lift.weight) &
    rm_outlier_amp$period == "morning", ]

agg_size <- aggregate(cbind(cal.spl.1m, lift.weight, PC1) ~ ID, data = body_size_data,
    mean)
agg_size$count <- aggregate(PC1 ~ ID, data = body_size_data, length)$PC1


ggplot(agg_size, aes(x = PC1, y = cal.spl.1m)) + geom_jitter(position = position_jitter(width = 0.2,
    height = 0.2), aes(size = count), color = viridis(10)[1]) + labs(x = "Body size (PC1)",
    y = "Sound pressure level (dB)", color = "Period", size = "Sample\n size") +
    # scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7)
    # +
theme_classic(base_size = 30)

Code

# geom_abline(slope = c(ES_continious = sm_size[2, 1],
# ES_continious_plus_ES_interaction= sm_size[2, 1] + sm_size[4,
# 1]), intercept = c(intercept = sm_size[1, 1],
# intercept_plus_ES_categorical = sm_size[1, 1] + sm_size[3,
# 1]), color = viridis(2, begin = 0.2, end = 0.8, alpha =
# 0.6)[2:1], size = 1) + guides(colour =
# guide_legend(override.aes = list(size=8)))

# ggsave('./output/body_size_vs_SPL.jpeg', width = 11, height =
# 7)

ggplot(agg_size, aes(x = lift.weight, y = cal.spl.1m)) + geom_jitter(position = position_jitter(width = 0.2,
    height = 0.2), aes(size = count), color = viridis(10)[1]) + labs(x = "Lifiting power (g)",
    y = "Sound pressure level (dB)", color = "Period", size = "Sample\n size") +
    # scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7)
    # +
theme_classic(base_size = 30)

Code

# geom_abline(slope = c(ES_continious = sm_size[2, 1],
# ES_continious_plus_ES_interaction= sm_size[2, 1] + sm_size[4,
# 1]), intercept = c(intercept = sm_size[1, 1],
# intercept_plus_ES_categorical = sm_size[1, 1] + sm_size[3,
# 1]), color = viridis(2, begin = 0.2, end = 0.8, alpha =
# 0.6)[2:1], size = 1) + guides(colour =
# guide_legend(override.aes = list(size=8)))

# ggsave('./output/lifting_power_vs_SPL.jpeg', width = 11,
# height = 7)

Each dot represents the average for a single individual

Effect size posterior distribution

Code

# simplify names
colnames(best_mod_size$Sol) <- c("Intercept", "lifted weight", "Body size (PC1)")

# stack posteriors
Y <- stack(as.data.frame(best_mod_size$Sol[, -1]))

# plot posteriors
ggplot(Y, aes(y = values, x = ind)) + geom_hline(yintercept = 0, col = "red",
    lty = 2) + geom_violin(color = "gray40", fill = viridis(n = 1,
    alpha = 0.25, begin = 0.4)) + labs(y = "Effect size posterior distribution",
    x = "Predictor") + theme_classic(base_size = 24) + theme(axis.text.x = element_text(angle = 45,
    vjust = 1, hjust = 1))

Song PSL is not affected by condition (either body size or lifting power)

8 SPL vs playback

Using data from 2019 and 2021
Recordings from 2019 were not calibrated and were center around SPL mean from 2021
Comparing before playback vs after playback (during regular singing)
Including condition parameters as covariates
No songtype random intercept was included

Playback as predictor and random slope: \[SPL \sim + playback + (playback | individual)\]
Playback and body size interaction as predictor and random intercept: \[SPL \sim + playback * body\ size + (1 | individual) \]
Playback and lifting power interaction as predictor and random intercept: \[SPL \sim + playback * lifting\ power + (1 | individual) \]
Playback as predictor and only random intercept: \[SPL \sim + playback + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (random intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

Code

playback_dat <- rm_outlier_amp[rm_outlier_amp$Treatment %in% c("Before_playback",
    "After_playback"), ]

# label 2019 recordings
playback_dat$ID <- ifelse(grepl("2019", playback_dat$org.sound.file),
    paste(playback_dat$ID, "(2019)"), playback_dat$ID)

# remove indiviudal recorded in 2019 and 2021
playback_dat <- playback_dat[playback_dat$ID != "397 (2019)", ]

random_slope_mod <- lmer(cal.spl.1m ~ Treatment + (Treatment | ID),
    data = playback_dat, REML = FALSE)

interaction_size_mod <- lmer(cal.spl.1m ~ PC1 + Treatment + PC1:Treatment +
    (1 | ID), data = playback_dat, REML = FALSE)

interaction_lift_mod <- lmer(cal.spl.1m ~ PC1 + lift.weight + PC1:lift.weight +
    (1 | ID), data = playback_dat, REML = FALSE)

random_intercept_only_mod <- lmer(cal.spl.1m ~ Treatment + (1 | ID) +
    (1 | songtype), data = playback_dat, REML = FALSE)

null_mod <- lmer(cal.spl.1m ~ 1 + (1 | ID) + (1 | songtype), data = playback_dat,
    REML = FALSE)

Code

d1 <- (-2 * logLik(random_intercept_only_mod)) + (2 * logLik(random_slope_mod))

pval_rand_slope_vs_rand_intercept <- 0.5 * pchisq(d1[1], 0, lower.tail = FALSE) +
    0.5 * pchisq(d1[1], 1, lower.tail = FALSE)

d2 <- (-2 * logLik(null_mod)) + (2 * logLik(random_intercept_only_mod))

pval_rand_intercept_vs_null <- 0.5 * pchisq(d2[1], 0, lower.tail = FALSE) +
    0.5 * pchisq(d2[1], 1, lower.tail = FALSE)

d3 <- (-2 * logLik(null_mod)) + (2 * logLik(random_slope_mod))

pval_rand_slope_vs_null <- 0.5 * pchisq(d3[1], 0, lower.tail = FALSE) +
    0.5 * pchisq(d3[1], 1, lower.tail = FALSE)

d4 <- (-2 * logLik(interaction_size_mod)) + (2 * logLik(random_slope_mod))

pval_rand_slope_vs_interaction_size <- 0.5 * pchisq(d4[1], 0, lower.tail = FALSE) +
    0.5 * pchisq(d4[1], 1, lower.tail = FALSE)

d5 <- (-2 * logLik(interaction_lift_mod)) + (2 * logLik(random_slope_mod))

pval_rand_slope_vs_interaction_lift <- 0.5 * pchisq(d5[1], 0, lower.tail = FALSE) +
    0.5 * pchisq(d5[1], 1, lower.tail = FALSE)

Code

agg_playback <- aggregate(cal.spl.1m ~ Treatment + ID, data = playback_dat,
    mean)

agg_playback$sd <- aggregate(cal.spl.1m ~ Treatment + ID, data = playback_dat,
    sd)$cal.spl.1m

agg_playback$Treatment <- factor(gsub("_playback", "", agg_playback$Treatment),
    levels = c("Before", "After"))

ggplot(agg_playback, aes(x = Treatment, y = cal.spl.1m, color = Treatment)) +
    geom_point(size = 2, show.legend = FALSE) + geom_errorbar(width = 0.05,
    aes(ymin = cal.spl.1m - sd, ymax = cal.spl.1m + sd), show.legend = FALSE) +
    labs(x = "Playback treatment", y = "Sound pressure level (dB)",
        color = "Period") + scale_color_viridis_d(begin = 0.1, end = 0.8,
    alpha = 1) + facet_wrap(~ID, nrow = 3, scales = "free_y") + theme_classic(base_size = 24)

Code

sub_agg <- agg_playback[agg_playback$ID %in% unique(agg_playback$ID)[1:8],
    ]
sub_agg$ID <- substr(sub_agg$ID, 0, 3)

ggplot(sub_agg, aes(x = Treatment, y = cal.spl.1m, color = Treatment)) +
    geom_point(size = 4) + geom_errorbar(width = 0.3, aes(ymin = cal.spl.1m -
    sd, ymax = cal.spl.1m + sd), show.legend = FALSE, lwd = 2) + labs(x = "Playback treatment",
    y = "Sound pressure level (dB)", color = "Period") + scale_color_viridis_d(begin = 0.1,
    end = 0.8, alpha = 1) + facet_wrap(~ID, nrow = 2, scales = "free_y") +
    theme_classic(base_size = 30) + theme(axis.text.x = element_blank())

Code

# ggsave('./output/SPL_by_playback.jpeg', width = 14, height =
# 7)

The random slope model is significantly better than the random intercept model (p = < 0.0001), the interaction with body size model (p = < 0.0001), the interaction with lifting power model (p = < 0.0001) and the null model (p = < 0.0001)

The playback stimuli elicits a significant change in SPL but the direction of the response varies between individuals
Direction of the response doesn’t seem to be affected by individual condition (neither body size nor lifting power)

9 SPL vs coordination

Random slope models, using data from 2019 and 2021
Only for individuals in which coordination was observed when recorded
SPL of songs sang during coordinated singing were compared to regular singing songs from the same recording session
No song type random intercept was included

Coordination as predictor and random slope: \[SPL \sim + coordination + (coordination | individual)\]
Coordination as predictor and random intercept only: \[SPL \sim + coordination + (1 | individual)\]
Null model with no predictor (random intercept only): \[SPL \sim + 1 + (1 | individual)\]

Code

# fixed effects with HPD intervals
coord_data <- rm_outlier_amp[rm_outlier_amp$org.sound.file %in% rm_outlier_amp$org.sound.file[rm_outlier_amp$Treatment ==
    "Coordination"], ]

coord_data <- coord_data[coord_data$Treatment %in% c("Regular_singing",
    "Coordination"), ]

random_slope_mod <- lmer(cal.spl.1m ~ Treatment + (Treatment | ID),
    data = coord_data, REML = FALSE)

random_intercept_only_mod <- lmer(cal.spl.1m ~ Treatment + (1 | ID),
    data = coord_data, REML = FALSE)

null_mod <- lmer(cal.spl.1m ~ 1 + (1 | ID), data = coord_data, REML = FALSE)

Code

d1 <- (-2 * logLik(random_intercept_only_mod)) + (2 * logLik(random_slope_mod))

random_slope_vs_random_intercept <- 0.5 * pchisq(d1[1], 0, lower.tail = FALSE) +
    0.5 * pchisq(d1[1], 1, lower.tail = FALSE)

d2 <- (-2 * logLik(null_mod)) + (2 * logLik(random_intercept_only_mod))

random_intercept_vs_null <- 0.5 * pchisq(d2[1], 0, lower.tail = FALSE) +
    0.5 * pchisq(d2[1], 1, lower.tail = FALSE)

d3 <- (-2 * logLik(null_mod)) + (2 * logLik(random_slope_mod))

random_slope_vs_null <- 0.5 * pchisq(d3[1], 0, lower.tail = FALSE) +
    0.5 * pchisq(d3[1], 1, lower.tail = FALSE)

Code

agg_coord <- aggregate(cal.spl.1m ~ Treatment + ID, data = coord_data,
    mean)

agg_coord$sd <- aggregate(cal.spl.1m ~ Treatment + ID, data = coord_data,
    sd)$cal.spl.1m

agg_coord$Treatment <- factor(ifelse(agg_coord$Treatment == "Coordination",
    "Coordinated", "Solo"), levels = c("Solo", "Coordinated"))

ggplot(agg_coord, aes(x = Treatment, y = cal.spl.1m, color = Treatment)) +
    geom_point(size = 2, show.legend = FALSE) + geom_errorbar(width = 0.05,
    aes(ymin = cal.spl.1m - sd, ymax = cal.spl.1m + sd), show.legend = FALSE) +
    labs(x = "coord treatment", y = "Sound pressure level (dB)", color = "Period") +
    scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 1) + facet_wrap(~ID,
    nrow = 2, scales = "free_y") + theme_classic(base_size = 24)

Code

sub_agg <- agg_coord[agg_coord$ID %in% unique(agg_coord$ID)[c(1:4,
    11:8)], ]
sub_agg$ID <- substr(sub_agg$ID, 0, 3)

ggplot(sub_agg, aes(x = Treatment, y = cal.spl.1m, color = Treatment)) +
    geom_point(size = 4, show.legend = TRUE) + geom_errorbar(width = 0.3,
    aes(ymin = cal.spl.1m - sd, ymax = cal.spl.1m + sd), show.legend = FALSE,
    lwd = 2) + labs(x = "Singing type", y = "Sound pressure level (dB)",
    color = "Singing type") + scale_color_viridis_d(begin = 0.1, end = 0.8,
    alpha = 1) + facet_wrap(~ID, nrow = 2, scales = "free_y") + theme_classic(base_size = 30) +
    theme(axis.text.x = element_blank())

Code

# ggsave('./output/SPL_by_coordination.jpeg', width = 14, height
# = 7)

The random slope model is significantly better than the random intercept model (p = < 0.0001) and the null model (p = < 0.0001)

The playback stimuli elicits a significant change in SPL

10 SPL vs aggresive interaction

2 types of interactions: perch interaction and chases
SPL of songs before (i.e. regular singing) and after interactions were compared
Interactions encoded in 2 ways:
- 1. factor clumping perch interactions and chases in a single category (2 levels: before and after interaction)
- 1. Interaction type (chase of perch interaction)

Interaction as predictor: \[SPL \sim + interaction + (1 | individual)\]
Interaction and interaction type as predictors: \[SPL \sim + interaction + interaction\ type + (1 | individual)\]
Statistical interaction between “Interaction” and “interaction type” as predictor: \[SPL \sim + interaction * interaction\ type + (1 | individual)\]
Null model with no predictor (intercept only): \[SPL \sim + 1 + (1 | individual)\]

Interactions from videos were omitted

Code

# subset data
dat_chase <- rm_outlier_amp[rm_outlier_amp$org.sound.file %in% rm_outlier_amp$org.sound.file[rm_outlier_amp$Treatment ==
    "After_chase"], ]

dat_chase <- dat_chase[dat_chase$Treatment %in% c("Regular_singing",
    "After_chase"), ]

dat_chase$inter_treatment <- ifelse(grepl("After", dat_chase$Treatment),
    "2_After", "1_Before")

dat_chase$inter_type <- "2_Chase"

dat_interaction <- rm_outlier_amp[rm_outlier_amp$Treatment %in% c("After_interaction",
    "Before_interaction"), ]


dat_interaction$inter_treatment <- ifelse(grepl("After", dat_interaction$Treatment),
    "2_After", "1_Before")
dat_interaction$inter_type <- "1_Perch_interaction"


# fix range of non-calibrated recs
dat_interaction$cal.spl.1m[dat_interaction$year < 2021] <- dat_interaction$cal.spl.1m[dat_interaction$year <
    2021] + mean(dat_interaction$cal.spl.1m[dat_interaction$year ==
    2021]) - mean(dat_interaction$cal.spl.1m[dat_interaction$year <
    2021])

dat_aggresive_inter <- rbind(dat_interaction, dat_chase)

# fix low value one
dat_aggresive_inter$cal.spl.1m[dat_aggresive_inter$sound.files ==
    dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl.1m)]] <- dat_aggresive_inter$cal.spl.1m[dat_aggresive_inter$sound.files ==
    dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl.1m)]] +
    mean(dat_aggresive_inter$cal.spl.1m[dat_aggresive_inter$sound.files !=
        dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl.1m)]]) -
    mean(dat_aggresive_inter$cal.spl.1m[dat_aggresive_inter$sound.files ==
        dat_aggresive_inter$sound.files[which.min(dat_aggresive_inter$cal.spl.1m)]])

# table(dat_aggresive_inter$ID,
# dat_aggresive_inter$inter_treatment,
# dat_aggresive_inter$inter_type)

Code

# model SPL by interaction
interaction_formulas <- c("cal.spl.1m ~ 1", "cal.spl.1m ~ inter_treatment",
    "cal.spl.1m ~ inter_treatment + inter_type", "cal.spl.1m ~ inter_treatment * inter_type")

# remove video data
dat_aggresive_inter <- dat_aggresive_inter[grep("2013|2014", dat_aggresive_inter$sound.files,
    invert = TRUE), ]

# dat_aggresive_inter$inter_treatment <-
# factor(dat_aggresive_inter$inter_treatment, levels =
# c('After', 'Before'))

# dat_aggresive_inter$inter_type <-
# factor(dat_aggresive_inter$inter_type, levels = c('Perch
# interaction', 'Chase'))

inter_mods <- pblapply(interaction_formulas, function(x) {

    replicate(n = 3, MCMCglmm(as.formula(x), random = ~ID, data = dat_aggresive_inter,
        verbose = FALSE, nitt = itrns, start = list(QUASI = FALSE)),
        simplify = FALSE)

})

names(inter_mods) <- gsub("cal.spl.1m ~ ", "", interaction_formulas)

saveRDS(list(interaction_formulas = interaction_formulas, inter_mods = inter_mods),
    "./output/mcmcglmm_interaction_models.RDS")

10.1 Model selection

Code

attach(readRDS("./output/mcmcglmm_interaction_models.RDS"))

mod_list <- lapply(inter_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl.1m ~ ", "", interaction_formulas)

model_selection <- model.sel(mod_list, rank = "DIC")

model_selection

	(Intercept)	inter_treatment	inter_type	inter_treatment:inter_type	df	logLik	DIC	delta	weight
cal.spl ~ inter_treatment * inter_type	100.7272	+	+	+	6	-2678.677	5369.257	0.00000	9.999479e-01
cal.spl ~ inter_treatment + inter_type	100.1180	+	+	NA	5	-2688.952	5388.980	19.72329	5.213387e-05
cal.spl ~ inter_treatment	101.0589	+	NA	NA	4	-2708.553	5427.038	57.78129	2.837462e-13
cal.spl ~ 1	101.9855	NA	NA	NA	3	-2802.094	5613.167	243.91042	1.085180e-53

10.2 Best model summary

Code

# fixed effects with HPD intervals
best_mod_interaction <- inter_mods[[which(names(inter_mods) == row.names(model_selection)[1])]][[1]]

sm_interaction <- as.data.frame(summary(best_mod_interaction)$solutions[,
    -5])

sm_interaction

Code

# subset data dat_aggresive_inter$inter_treatment <-
# factor(dat_aggresive_inter$inter_treatment, levels =
# c('1_Before', '2_After'))

dat_aggresive_inter$year <- factor(dat_aggresive_inter$year, levels = c("2013",
    "2014", "2019", "2021"))

# dat_aggresive_inter$inter_type <- gsub('_', ' ',
# dat_aggresive_inter$inter_type)

# dat_aggresive_inter$inter_type <-
# factor(dat_aggresive_inter$inter_type, levels =
# c('1_Perch_interaction', '2_Chase'))

agg_aggresive_inter <- aggregate(cal.spl.1m ~ ID + inter_type + inter_treatment +
    year, data = dat_aggresive_inter, mean)

agg_aggresive_inter$ID <- as.character(agg_aggresive_inter$ID)

ggplot(dat_aggresive_inter, aes(x = inter_treatment, y = cal.spl.1m)) +
    geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25,
        begin = 0.4)) + geom_point(size = 4, data = agg_aggresive_inter,
    aes(group = ID, color = year)) + geom_line(data = agg_aggresive_inter,
    aes(group = ID, color = year)) + labs(x = "Context", y = "Sound pressure level (dB)") +
    facet_wrap(~inter_type) + scale_color_viridis_d(begin = 0.1, end = 0.8,
    alpha = 0.7) + theme_classic(base_size = 24)

Code

agg_dat_agg <- aggregate(cal.spl.1m ~ inter_treatment + inter_type,
    data = dat_aggresive_inter, mean)
agg_dat_agg$sd <- aggregate(cal.spl.1m ~ inter_treatment + inter_type,
    data = dat_aggresive_inter, sd)$cal.spl.1m

agg_dat_agg$inter_treatment <- factor(substr(agg_dat_agg$inter_treatment,
    3, 100), levels = substr(agg_dat_agg$inter_treatment, 3, 100)[1:2])

ggplot(agg_dat_agg, aes(x = inter_treatment, y = cal.spl.1m, color = inter_treatment)) +
    geom_point(size = 4, show.legend = FALSE) + geom_errorbar(width = 0.3,
    aes(ymin = cal.spl.1m - sd, ymax = cal.spl.1m + sd), show.legend = FALSE,
    lwd = 2) + labs(x = "Period", y = "Sound pressure level (dB)") +
    scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 1) + facet_wrap(~inter_type,
    nrow = 1, scales = "free_y", labeller = labeller(inter_type = c(`1_Perch_interaction` = "Perch interaction",
        `2_Chase` = "Chase"))) + theme_classic(base_size = 30)

Code

# ggsave('./output/interaction_vs_SPL.jpeg', width = 11, height
# = 5.8)

Effect size posterior distribution

Code

# extract mcmcs
mcmcs <- best_mod_interaction$Sol

# make empty list
mcmc_l <- list()

## perch interaction before playback
mcmc_l[[length(mcmc_l) + 1]] <- data.frame(spl = mcmcs[, "(Intercept)"],
    treatment = "Before", type = "Perch interaction")

## perch interaction after playback
mcmc_l[[length(mcmc_l) + 1]] <- data.frame(spl = mcmcs[, "(Intercept)"] +
    mcmcs[, "inter_treatment2_After"], treatment = "After", type = "Perch interaction")

# chase before playback
mcmc_l[[length(mcmc_l) + 1]] <- data.frame(spl = mcmcs[, "(Intercept)"] +
    mcmcs[, "inter_type2_Chase"], treatment = "Before", type = "Chase")

# chase after playback
mcmc_l[[length(mcmc_l) + 1]] <- data.frame(spl = mcmcs[, "(Intercept)"] +
    mcmcs[, "inter_type2_Chase"] + mcmcs[, "inter_treatment2_After:inter_type2_Chase"],
    treatment = "After", type = "Chase")

mcmc_df <- do.call(rbind, mcmc_l)
colnames(mcmc_df)[1] <- "spl"

mcmc_df$treatment_type <- paste(mcmc_df$treatment, mcmc_df$type, sep = "-")

# p-values
combs <- combn(x = unique(mcmc_df$treatment_type), m = 2)[, c(1, 6)]


pvals_l <- lapply(1:ncol(combs), function(x) {

    mcmc1 <- mcmc_df$spl[mcmc_df$treatment_type == combs[1, x]]
    mcmc2 <- mcmc_df$spl[mcmc_df$treatment_type == combs[2, x]]
    p <- if (mean(mcmc2) > mean(mcmc1))
        sum(mcmc1 > mcmc2)/length(mcmc1) else sum(mcmc2 > mcmc1)/length(mcmc1)

    out <- data.frame(mcmc1 = combs[1, x], mcmc2 = combs[2, x], p = p)

    return(out)
})

pvals <- do.call(rbind, pvals_l)

mcmc_df$treatment <- factor(mcmc_df$treatment, levels = c("Before",
    "After"))
mcmc_df$type <- factor(mcmc_df$type, levels = c("Perch interaction",
    "Chase"))

agg <- aggregate(spl ~ type + treatment, mcmc_df, mean)

# plot posteriors
ggplot(mcmc_df, aes(y = spl, x = treatment)) + geom_point(data = agg) +
    # geom_line(data = agg, stat = 'identity') +
geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25,
    begin = 0.4)) + labs(y = "SPL posterior distribution", x = "Predictor") +
    theme_classic(base_size = 24) + facet_wrap(~type) + theme(axis.text.x = element_text(angle = 45,
    vjust = 0.9, hjust = 1))

kable(pvals)

Male-male aggressive interactions result in a significant SPL increase
Chase effect cannot be compare against perch interactions as not all recordings were calibrated

11 SPL and song structure

Song structure measured as spectrographic consistency and gap duration

11.1 Song spectrographic consistency

spectrographic consistency: mean spectrographic cross-correlation of all songs in a singing bout
Recordings from interactions were used as they showed the highest increase in SPL
Spectrographic consistency was measured for all songs within contexts (before & after interaction, clumping chase and perch interactions)
Spectrographic consistency was compared between the contexts

interaction as predictor and spectrographic consistency as response: \[consistency \sim + interaction + (1 | individual) + (1 | song\ type)\]
interaction and SPL as predictors and spectrographic consistency as response: \[consistency \sim + interaction + SPL + (1 | individual) + (1 | song\ type)\]
interaction between “interaction” and SPL as predictor and spectrographic consistency as response: \[consistency \sim + interaction * SPL + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (intercept only): \[consistency \sim + 1 + (1 | individual) + (1 | song\ type)\]

Code

agg_amp_inter <- aggregate(cal.spl.1m ~ ID + inter_type, data = dat_aggresive_inter[dat_aggresive_inter$inter_treatment ==
    "2_After" & dat_aggresive_inter$extsn != "mp3", ], mean)

agg_amp_inter$amp.change <- agg_amp_inter$cal.spl.1m - aggregate(cal.spl.1m ~
    ID + inter_type, data = dat_aggresive_inter[dat_aggresive_inter$inter_treatment ==
    "1_Before" & dat_aggresive_inter$extsn != "mp3", ], mean)$cal.spl.1m

# sort by change
agg_amp_inter <- agg_amp_inter[order(agg_amp_inter$amp.change, decreasing = TRUE),
    ]

# subset data
dat_song_consistency <- dat_aggresive_inter[dat_aggresive_inter$ID %in%
    agg_amp_inter$ID[agg_amp_inter$amp.change >= 0.5], ]

Code

dat_song_consistency$ID.treatment <- paste(dat_song_consistency$ID,
    dat_song_consistency$inter_treatment, sep = "-")

dat_song_consistency_l <- lapply(unique(dat_song_consistency$ID.treatment),
    function(x) {

        X <- dat_song_consistency[dat_song_consistency$ID.treatment ==
            x, ]

        print(x)
        if (nrow(X) > 1) {
            # measure cross-correlation
            xcrr <- (cross_correlation(X, path = "./data/raw/cuts",
                pb = FALSE, wl = 150))

            # measure gaps
            gps <- gaps(X, pb = FALSE)$gaps
            gps <- gps[gps < 1]
            mean.gaps <- mean(gps, na.rm = TRUE)

            # mean xcorr
            if (!is.na(xcrr[1])) {
                xcrr <- mean(xcrr[lower.tri(xcrr)])
            }

            out_df <- data.frame(sound.files = x, ID = X$ID[1], Treatment = X$Treatment[1],
                xcrr = xcrr, mean.gaps = mean.gaps, mean.spl = mean(X$cal.spl.1m),
                songtype = X$songtype[1], treatment = X$inter_treatment[1])
        } else out_df <- data.frame(sound.files = x, ID = X$ID[1], Treatment = X$Treatment[1],
            xcrr = NA, mean.gaps = NA, mean.spl = mean(X$cal.spl.1m),
            songtype = X$songtype[1], treatment = X$inter_treatment[1])

        return(out_df)

    })

dat_song_consistency <- do.call(rbind, dat_song_consistency_l)

dat_song_consistency <- dat_song_consistency[!is.na(dat_song_consistency$xcrr),
    ]

saveRDS(dat_song_consistency, "./output/dat_song_spectrographic_consistency_and_gaps.RDS")

Code

dat_song_consistency <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")


# model SPL by body size
song_consistency_formulas <- c("xcrr ~ 1", "xcrr ~ treatment", "xcrr ~ treatment + mean.spl",
    "xcrr ~ treatment * mean.spl", "xcrr ~mean.spl")

song_consistency_mods <- pblapply(song_consistency_formulas, function(x) {

    replicate(n = 3, MCMCglmm(as.formula(x), random = ~ID + songtype,
        data = dat_song_consistency, verbose = FALSE, nitt = itrns,
        start = list(QUASI = FALSE)), simplify = FALSE)

})

names(song_consistency_mods) <- gsub("cal.spl.1m ~ ", "", song_consistency_formulas)

saveRDS(list(song_consistency_formulas = song_consistency_formulas,
    song_consistency_mods = song_consistency_mods), "./output/mcmcglmm_song_consistency_models.RDS")

11.1.1 Model selection

Code

attach(readRDS("./output/mcmcglmm_song_consistency_models.RDS"))

mod_list <- lapply(song_consistency_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl.1m ~ ", "", song_consistency_formulas)

model_selection <- model.sel(mod_list, rank = "DIC")

model_selection

	(Intercept)	treatment	mean.spl	mean.spl:treatment	df	logLik	DIC	delta	weight
xcrr ~ treatment	0.7255519	+	NA	NA	5	22.10251	-39.67280	0.0000000	0.3600905
xcrr ~ 1	0.7419545	NA	NA	NA	4	21.48487	-39.15445	0.5183422	0.2778786
xcrr ~ treatment + mean.spl	0.7170405	+	-0.0063874	NA	6	21.35235	-38.20927	1.4635309	0.1732248
xcrr ~mean.spl	0.7409975	NA	-0.0049058	NA	5	20.71608	-37.77691	1.8958865	0.1395485
xcrr ~ treatment * mean.spl	0.7167435	+	-0.0072560	+	7	20.40215	-35.69422	3.9785794	0.0492577

11.1.2 Best model summary

Code

# summary
best_mod_song_consistency <- song_consistency_mods[[which(names(song_consistency_mods) ==
    row.names(model_selection)[1])]][[1]]

sm_song_consistency <- as.data.frame(summary(best_mod_song_consistency)$solutions)

sm_song_consistency[, -5]

	post.mean	l-95% CI	u-95% CI	eff.samp
(Intercept)	0.7255519	0.6577871	0.7976582	10109.16
treatmentBefore	0.0321469	-0.0364210	0.0982015	9700.00

Code

dat_song_consistency <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")

# subset data
dat_song_consistency$ID <- as.character(dat_song_consistency$ID)

ggplot(dat_song_consistency, aes(x = mean.spl, y = xcrr)) + geom_point(size = 4,
    color = viridis(5, alpha = 0.5)[3]) + labs(y = "Mean cross-correlation",
    x = "Sound pressure level (dB)") + scale_color_viridis_d(begin = 0.1,
    end = 0.8, alpha = 0.7) + theme_classic(base_size = 24)

Code

# stack posteriors
Y <- stack(as.data.frame(best_mod_song_consistency$Sol))
# Y$ind <- 'mean.spl'

# plot posteriors
ggplot(Y, aes(y = values, x = ind)) + geom_hline(yintercept = 0, col = "red",
    lty = 2) + geom_violin(color = "gray40", fill = viridis(n = 1,
    alpha = 0.25, begin = 0.4)) + labs(y = "Effect size posterior distribution",
    x = "Predictor") + theme_classic(base_size = 24) + theme(axis.text.x = element_text(angle = 45,
    vjust = 0.9, hjust = 1))

Aggressive interaction but not SPL affects the spectrographic consistency of songs

11.2 Song gap duration

-Gaps: silent time intervals between songs - Mean SPL has been mean-centered to facilitate interpretation

interaction as predictor and gap duration as response: \[gap\ duration \sim + interaction + (1 | individual) + (1 | song\ type)\]
interaction and SPL as predictors and gap duration as response: \[gap\ duration \sim + interaction + SPL + (1 | individual) + (1 | song\ type)\]
interaction between “interaction” and SPL as predictor and gap duration as response: \[gap\ duration \sim + interaction * SPL + (1 | individual) + (1 | song\ type)\]
Null model with no predictor (intercept only): \[gap\ duration \sim + 1 + (1 | individual) + (1 | song\ type)\]

Code

dat_gaps <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")

dat_gaps$mean.spl <- mean(dat_gaps$mean.spl) - dat_gaps$mean.spl

# model SPL by body size
gap_formulas <- c("mean.gaps ~ 1", "mean.gaps ~ treatment", "mean.gaps ~ treatment + mean.spl",
    "mean.gaps ~ treatment * mean.spl", "mean.gaps ~ mean.spl")

gap_mods <- pblapply(gap_formulas, function(x) {

    replicate(n = 3, MCMCglmm(as.formula(x), random = ~ID + songtype,
        data = dat_gaps, verbose = FALSE, nitt = itrns, start = list(QUASI = FALSE)),
        simplify = FALSE)

})

names(gap_mods) <- gsub("cal.spl.1m ~ ", "", gap_formulas)

saveRDS(list(gap_formulas = gap_formulas, gap_mods = gap_mods), "./output/mcmcglmm_gap_models.RDS")

11.2.1 Model selection

Code

attach(readRDS("./output/mcmcglmm_gap_models.RDS"))

mod_list <- lapply(gap_mods, "[[", 1)

names(mod_list) <- gsub("cal.spl.1m ~ ", "", gap_formulas)

model_selection <- model.sel(mod_list, rank = "DIC")

model_selection

	(Intercept)	treatment	mean.spl	mean.spl:treatment	df	logLik	DIC	delta	weight
mean.gaps ~ treatment * mean.spl	0.6678307	+	-0.0092550	+	7	25.85278	-46.42077	0.000000	0.71994310
mean.gaps ~ treatment + mean.spl	0.6729941	+	-0.0048448	NA	6	23.75377	-43.12839	3.292379	0.13879285
mean.gaps ~ treatment	0.6797477	+	NA	NA	5	23.11416	-42.36029	4.060472	0.09453179
mean.gaps ~ 1	0.7056103	NA	NA	NA	4	21.21979	-39.80056	6.620202	0.02628691
mean.gaps ~ mean.spl	0.7053243	NA	-0.0035014	NA	5	21.25474	-39.29793	7.122835	0.02044534

11.2.2 Best model summary

Code

# summary
best_mod_gap <- gap_mods[[which(names(gap_mods) == row.names(model_selection)[1])]][[1]]

sm_gap <- as.data.frame(summary(best_mod_gap)$solutions[, -5])

sm_gap

	post.mean	l-95% CI	u-95% CI	eff.samp
(Intercept)	0.6678307	0.6203597	0.7130255	9700.000
treatmentBefore	0.0636878	0.0114525	0.1166435	9700.000
mean.spl	-0.0092550	-0.0168684	-0.0014809	9313.085
treatmentBefore:mean.spl	0.0087411	-0.0012228	0.0179838	10105.701

Code

dat_gaps <- readRDS("./output/dat_song_spectrographic_consistency_and_gaps.RDS")

# subset data
dat_gaps$ID <- as.character(dat_gaps$ID)

ggplot(dat_gaps, aes(x = mean.gaps, y = xcrr)) + geom_point(size = 4,
    color = viridis(5, alpha = 0.5)[3]) + labs(y = "Mean gap duration (s)",
    x = "Sound pressure level (dB)") + scale_color_viridis_d(begin = 0.1,
    end = 0.8, alpha = 0.7) + theme_classic(base_size = 24)

Effect size posterior distribution

Code

# stack posteriors
Y <- stack(as.data.frame(best_mod_gap$Sol[, -1]))

# plot posteriors
ggplot(Y, aes(y = values, x = ind)) + geom_hline(yintercept = 0, col = "red",
    lty = 2) + geom_violin(color = "gray40", fill = viridis(n = 1,
    alpha = 0.25, begin = 0.4)) + labs(y = "Effect size posterior distribution",
    x = "Predictor") + theme_classic(base_size = 24) + theme(axis.text.x = element_text(angle = 45,
    vjust = 0.9, hjust = 1))

# zoom in
Y <- stack(as.data.frame(best_mod_gap$Sol[, -c(1, 2)]))
ggplot(Y, aes(y = values, x = ind)) + geom_hline(yintercept = 0, col = "red",
    lty = 2) + geom_violin(color = "gray40", fill = viridis(n = 1,
    alpha = 0.25, begin = 0.4)) + labs(y = "Effect size posterior distribution",
    x = "Predictor") + theme_classic(base_size = 24) + theme(axis.text.x = element_text(angle = 45,
    vjust = 0.9, hjust = 1))

Code

# extract mcmcs
mcmcs <- best_mod_gap$Sol

# make empty list
mcmc_l <- list()

## after interaction 0 spl
mcmc_l[[length(mcmc_l) + 1]] <- data.frame(spl = mcmcs[, "(Intercept)"],
    treatment = "After", type = "0 spl")

## before interaction 0 spl
mcmc_l[[length(mcmc_l) + 1]] <- data.frame(spl = mcmcs[, "(Intercept)"] +
    mcmcs[, "treatmentBefore"], treatment = "Before", type = "0 spl")


## before interaction 1 spl
mcmc_l[[length(mcmc_l) + 1]] <- data.frame(spl = mcmcs[, "(Intercept)"] +
    mcmcs[, "treatmentBefore"] + mcmcs[, "mean.spl"], treatment = "Before",
    type = "1 spl")

## after interaction 1 spl
mcmc_l[[length(mcmc_l) + 1]] <- data.frame(spl = mcmcs[, "(Intercept)"] +
    mcmcs[, "mean.spl"], treatment = "After", type = "1 spl")

mcmc_df <- do.call(rbind, mcmc_l)
colnames(mcmc_df)[1] <- "spl"

mcmc_df$treatment_type <- paste(mcmc_df$treatment, mcmc_df$type, sep = "-")

# p-values
combs <- combn(x = unique(mcmc_df$treatment_type), m = 2)


pvals_l <- lapply(1:ncol(combs), function(x) {

    mcmc1 <- mcmc_df$spl[mcmc_df$treatment_type == combs[1, x]]
    mcmc2 <- mcmc_df$spl[mcmc_df$treatment_type == combs[2, x]]
    p <- if (mean(mcmc2) > mean(mcmc1))
        sum(mcmc1 > mcmc2)/length(mcmc1) else sum(mcmc2 > mcmc1)/length(mcmc1)

    out <- data.frame(mcmc1 = combs[1, x], mcmc2 = combs[2, x], p = p)

    return(out)
})

pvals <- do.call(rbind, pvals_l)

mcmc_df$treatment <- factor(mcmc_df$treatment, levels = c("Before",
    "After"))
mcmc_df$type <- factor(mcmc_df$type, levels = c("Perch interaction",
    "Chase"))

# plot posteriors
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + ggplot(mcmc_df,
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + aes(y
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + =
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + spl,
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + x
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + =
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + treatment))
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + +
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + #
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + geom_hline(yintercept
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + =
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + 0,
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + col
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + =
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + 'red',
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + lty
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + =
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + 2)
ggplot(mcmc_df, aes(y = spl, x = treatment)) + # geom_hline(yintercept = 0, col = 'red', lty = 2) + +
geom_violin(color = "gray40", fill = viridis(n = 1, alpha = 0.25,
    begin = 0.4)) + labs(y = "SPL posterior distribution", x = "Predictor") +
    theme_classic(base_size = 24) + facet_wrap(~type) + theme(axis.text.x = element_text(angle = 45,
    vjust = 0.9, hjust = 1))

kable(pvals)

No detectable effect on gap duration

11.3 Diagnostic plots for MCMCglmm models

Include Gelman/Rubin’s convergence diagnostic, MCMC chain trace (all tree replicates in a single plot: yellow, blue and red colors) and autocorrelation plots:

11.3.1 SPL and period of the day

Code

mcmc_diagnostics(period_mods)

11.3.2 SPL and condition

Code

mcmc_diagnostics(size_mods)

[1] "model: Null"

[1] "model: PC1"

[1] "model: lift.weight"

[1] "model: lift.weight + PC1"

[1] "model: lift.weight *  PC1"

[1] "model: veg_density"

[1] "model: PC1 + veg_density"

[1] "model: lift.weight + veg_density"

[1] "model: lift.weight + PC1 + veg_density"

[1] "model: lift.weight *  PC1 * veg_density"

11.3.3 SPL and interactions

Code

mcmc_diagnostics(inter_mods)

[1] "model: Null"

[1] "model: inter_treatment"

[1] "model: inter_treatment + inter_type"

[1] "model: inter_treatment * inter_type"

11.3.4 SPL and spectrographic consistency

Code

mcmc_diagnostics(song_consistency_mods)

[1] "model: xcrr ~ 1"

[1] "model: xcrr ~ treatment"

[1] "model: xcrr ~ treatment + mean.spl"

[1] "model: xcrr ~ treatment * mean.spl"

[1] "model: xcrr ~mean.spl"

11.3.5 SPL and gap duration

Code

mcmc_diagnostics(gap_mods)

[1] "model: mean.gaps ~ 1"

[1] "model: mean.gaps ~ treatment"

[1] "model: mean.gaps ~ treatment + mean.spl"

[1] "model: mean.gaps ~ treatment * mean.spl"

[1] "model: mean.gaps ~ mean.spl"

Code

wv <- readWave("./data/raw/cuts/397_SUR_03-Feb_2021_08.06_10.57_A56-11.wav")

spectro(cutw(wv, from = 1.2, to = 1.6, output = "Wave"), wl = 0, scale = FALSE,
    ovlp = 95, flim = c(1, 10))

scrolling_spectro(wave = wv, wl = 300, t.display = 1.7, pal = viridis,
    ovlp = 90, grid = FALSE, flim = c(1, 9), colbg = "black", width = 1000,
    height = 500, collevels = seq(-50, 0, 5), res = 120, file.name = "397_SUR_03-Feb_2021_08.06_10.57_A56-11_v2.mp4")

scrolling_spectro(wave = wv, wl = 300, speed = 0.5, t.display = 1.7,
    pal = viridis, ovlp = 90, grid = FALSE, flim = c(1, 9), colbg = "black",
    width = 1000, height = 500, collevels = seq(-50, 0, 5), res = 120,
    file.name = "397_SUR_03-Feb_2021_08.06_10.57_A56-11_v3.mp4")

scrolling_spectro(wave = wv, wl = 300, t.display = 1.7, ovlp = 90,
    pal = magma, grid = FALSE, flim = c(1, 10), width = 1000, height = 500,
    res = 120, collevels = seq(-40, 0, 5), file.name = "397_SUR_03-Feb_2021_08.06_10.57_A56-11_v4.mp4",
    colbg = "black", speed = 0.2, axis.type = "none", loop = 3)


data("Phae.long4")

scrolling_spectro(wave = Phae.long4, wl = 300, t.display = 1.7, ovlp = 90,
    pal = magma, grid = FALSE, flim = c(1, 10), width = 1000, height = 500,
    res = 120, collevels = seq(-50, 0, 5), file.name = "no_axis_short.mp4",
    colbg = "black", speed = 1, axis.type = "none", loop = 3)

scrolling_spectro(wave = Phae.long4, wl = 300, t.display = 1.7, ovlp = 90,
    pal = viridis, grid = FALSE, flim = c(1, 10), width = 1000, height = 500,
    res = 120, collevels = seq(-50, 0, 5), file.name = "no_axis_viridis2.mp4",
    colbg = "black", speed = 1, axis.type = "none", loop = 3)


scrolling_spectro(wave = wv, wl = 150, t.display = 1.2, ovlp = 95,
    pal = magma, grid = FALSE, flim = c(1, 10), width = 1000, height = 500,
    res = 120, collevels = seq(-40, 0, 5), file.name = "3397_SUR_03-Feb_2021_08.06_10.57_A56-11_v5.mp4",
    colbg = "black", speed = 0.5, axis.type = "none", loop = 3)


scrolling_spectro(wave = Phae.long4, wl = 300, t.display = 1.7, ovlp = 90,
    pal = magma, grid = FALSE, flim = c(1, 10), width = 1000, height = 500,
    res = 120, collevels = seq(-40, 0, 5), file.name = "Phae.long4.mp4",
    colbg = "black", speed = 1, axis.type = "minimal", loop = 3)

Phae.long4.2 <- cutw(Phae.long4, from = 0, to = 2.8, output = "Wave")

scrolling_spectro(wave = Phae.long4.2, wl = 300, t.display = 1.7,
    ovlp = 90, pal = viridis, grid = FALSE, flim = c(1, 10), width = 1000,
    height = 500, res = 120, collevels = seq(-40, 0, 5), file.name = "Phae.long4_viridis.mp4",
    colbg = "black", speed = 1, axis.type = "minimal", loop = 5)

Session information

R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3 
LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3;  LAPACK version 3.9.0

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_BR.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] grid      parallel  stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] MCMCglmm_2.35      ape_5.7-1          coda_0.19-4.1      brmsish_1.0.0     
 [5] dynaSpec_1.0.1     gridExtra_2.3      lme4_1.1-35.1      Matrix_1.6-5      
 [9] corrplot_0.92      MuMIn_1.47.5       brms_2.20.4        Rcpp_1.0.12       
[13] rptR_0.9.22        readxl_1.4.3       viridis_0.6.5      viridisLite_0.4.2 
[17] Rraven_1.0.14      warbleR_1.1.30     NatureSounds_1.0.4 knitr_1.45        
[21] seewave_2.2.3      tuneR_1.4.6        ggplot2_3.4.4      pbapply_1.7-2     

loaded via a namespace (and not attached):
  [1] tensorA_0.36.2.1     rstudioapi_0.15.0    jsonlite_1.8.8      
  [4] magrittr_2.0.3       TH.data_1.1-2        estimability_1.4.1  
  [7] farver_2.1.1         nloptr_2.0.3         rmarkdown_2.25      
 [10] vctrs_0.6.5          minqa_1.2.6          RCurl_1.98-1.14     
 [13] base64enc_0.1-3      htmltools_0.5.7      distributional_0.3.2
 [16] signal_1.8-0         cellranger_1.1.0     StanHeaders_2.32.5  
 [19] htmlwidgets_1.6.4    plyr_1.8.9           testthat_3.2.1      
 [22] sandwich_3.1-0       emmeans_1.10.0       zoo_1.8-12          
 [25] igraph_1.6.0         mime_0.12            lifecycle_1.0.4     
 [28] pkgconfig_2.0.3      colourpicker_1.3.0   R6_2.5.1            
 [31] fastmap_1.1.1        shiny_1.8.0          digest_0.6.34       
 [34] colorspace_2.1-0     crosstalk_1.2.1      labeling_0.4.3      
 [37] fansi_1.0.6          abind_1.4-5          compiler_4.3.1      
 [40] proxy_0.4-27         withr_3.0.0          backports_1.4.1     
 [43] inline_0.3.19        shinystan_2.6.0      QuickJSR_1.1.3      
 [46] pkgbuild_1.4.3       MASS_7.3-60.0.1      corpcor_1.6.10      
 [49] rjson_0.2.21         gtools_3.9.5         loo_2.6.0           
 [52] tools_4.3.1          httpuv_1.6.14        fftw_1.0-8          
 [55] threejs_0.3.3        glue_1.7.0           nlme_3.1-164        
 [58] promises_1.2.1       checkmate_2.3.1      reshape2_1.4.4      
 [61] generics_0.1.3       gtable_0.3.4         xml2_1.3.6          
 [64] utf8_1.2.4           pillar_1.9.0         ggdist_3.3.1        
 [67] markdown_1.12        stringr_1.5.1        posterior_1.5.0     
 [70] later_1.3.2          splines_4.3.1        dplyr_1.1.4         
 [73] lattice_0.22-5       survival_3.5-7       tidyselect_1.2.0    
 [76] miniUI_0.1.1.1       svglite_2.1.3        stats4_4.3.1        
 [79] xfun_0.41            bridgesampling_1.1-2 brio_1.1.4          
 [82] matrixStats_1.2.0    DT_0.31              rstan_2.32.5        
 [85] stringi_1.8.3        yaml_2.3.8           boot_1.3-28.1       
 [88] kableExtra_1.4.0     evaluate_0.23        codetools_0.2-19    
 [91] dtw_1.23-1           tibble_3.2.1         cli_3.6.2           
 [94] RcppParallel_5.1.7   systemfonts_1.0.5    shinythemes_1.2.0   
 [97] xtable_1.8-4         munsell_0.5.0        rstantools_2.4.0    
[100] ellipsis_0.3.2       cubature_2.1.0       dygraphs_1.1.1.6    
[103] bayesplot_1.11.0     Brobdingnag_1.2-9    bitops_1.0-7        
[106] mvtnorm_1.2-4        scales_1.3.0         xts_0.13.2          
[109] rlang_1.1.3          cowplot_1.1.3        multcomp_1.4-25     
[112] formatR_1.14         shinyjs_2.1.0