What the code does:

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

Next steps

Load packages

## vector with package names
x <- c( "pbapply", "parallel", "ggplot2", "warbleR", "Rraven", "viridis", "readxl", "rptR", "MCMCglmm", "MuMIn", "corrplot", "lme4", "grid", "gridExtra")

aa <- lapply(x, function(y) {
  
  # check if installed, if not then install 
  if (!y %in% installed.packages()[,"Package"]) 
    install.packages(y) 

  # load package
  try(require(y, character.only = T), silent = T)
})

Functions and parameters

#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 <- 100000

# 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=10e-2), G=list(G1=list(V=1,n=10e-2)) )
    # 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")
} 
  
}

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

Repeatability

We estimated repeatability using linear mixed models with MCMC sampling (MCMCglmm() function)

# 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, c(rep_grid$outliers[x], 1- rep_grid$outliers[x]))
  
        if (!rep_grid$only.low.outliers[x])
    X <- X[X$cal.spl >= outlier_quant[1] & X$cal.spl <= outlier_quant[2],] else  X <- X[X$cal.spl >= outlier_quant[1],]
    
    # quantlie for each max quantile
    quant <- quantile(X$cal.spl, probs = 1 - rep_grid$max_quantile[x])
    
    # subset
   X <- X[X$cal.spl >= quant, ]
  })
  
  quant_subet <- do.call(rbind, quant_subet)
  
  # frequentist
    # out <- rpt(cal.spl ~ (1 | ID), grname = "ID", data = quant_subet, datatype = "Gaussian", npermut = 0, nboot = 100)

  # bayesian
    out <- rpt.mcmcLMM(y = quant_subet$cal.spl, 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.R[1], out$CI.R[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

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=.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)

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), ]

kable(head(repts_df))

	max_quantile	outliers	only.low.outliers	repeatability	lowCI	hiCI	range
116	0.4	0.10	FALSE	0.5790952	0.2912715	0.8124867	0.5212152
120	0.8	0.10	FALSE	0.5760627	0.2933315	0.8183214	0.5249900
86	0.7	0.01	FALSE	0.5742588	0.2970074	0.8256426	0.5286352
52	0.6	0.10	TRUE	0.5679999	0.2918631	0.8131767	0.5213137
6	0.4	0.00	TRUE	0.5633166	0.2903133	0.8141992	0.5238859
106	0.5	0.05	FALSE	0.5549716	0.2926263	0.8167255	0.5240992

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.4 upper quartile after excluding 0.2 lower tail outliers:

# 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, c(repts_df$outliers[1], 1))
    X <- X[X$cal.spl >= outlier_quant[1] & X$cal.spl <= outlier_quant[2],]
    
    # quantlie for each max quantile (0.6 was selected due to high repeatability)
    quant <- quantile(X$cal.spl, probs = repts_df$max_quantile[1])
    
    # subset
   X <- X[X$cal.spl >= quant, ]
  })

rm_outlier_amp <- do.call(rbind, rm_outlier_amp_l)

Hypothesis testing analysis

We use bayesian MCMC 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.

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)\]

# 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 ~ 1", "cal.spl ~ period")

period_mods <- pblapply(period_formulas, function(x){
  
  replicate(n = 3, MCMCglmm(as.formula(x), random = ~ ID + songtype, data = period_data, verbose = FALSE, nitt = itrns,  start = list(QUASI = FALSE)), simplify = FALSE)
  
})

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

Model selection

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

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

names(mod_list) <- gsub("cal.spl ~ ", "", period_formulas)

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

model_selection

	(Intercept)	period	family	df	logLik	DIC	delta	weight
period	98.38209	+	(NA)	5	-6192.009	12389.99	0.00000	1.000000e+00
1	99.28286	NA	(NA)	4	-6223.692	12452.40	62.40826	2.806839e-14

Best model summary

# 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

	post.mean	l-95% CI	u-95% CI	eff.samp
(Intercept)	98.382089	90.2424713	106.17904	4700
periodmorning	1.232936	0.9223919	1.52236	4700

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, 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

# 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

SPL vs condition

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) and lifting power

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)\]
Lifting power and body size as predictors: \[SPL \sim + lifting\ power + body\ size + (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)\]
Null model with no predictor (intercept only): \[SPL \sim + 1 + (1 | individual) + (1 | song\ type)\]

Correlation between body size parameters included in the PCA

caps <- read_excel("/home/m/Dropbox/LBH data/Additional data files/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

Variance explained by each PC

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

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("/home/m/Dropbox/LBH data/Morphology-condition/other files/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)
})

# 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 ~ 1", "cal.spl ~ PC1", "cal.spl ~ lift.weight", "cal.spl ~ lift.weight + PC1", "cal.spl ~ lift.weight *  PC1")


size_mods <- pblapply(size_formulas, 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 ~ ", "", size_formulas)
saveRDS(list(size_formulas = size_formulas, size_mods = size_mods), "./output/mcmcglmm_size_models.RDS")

Model selection

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

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

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

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

model_selection

	(Intercept)	PC1	lift.weight	lift.weight:PC1	family	df	logLik	DIC	delta	weight
PC1	101.09190	1.0046808	NA	NA	(NA)	5	-19465.59	38951.93	0.0000000	0.2097915
lift.weight * PC1	77.00884	10.1703441	1.523415	-0.6213227	(NA)	7	-19465.47	38951.96	0.0295211	0.2067176
lift.weight + PC1	84.96609	0.7948061	1.007591	NA	(NA)	6	-19465.58	38952.07	0.1427185	0.1953426
lift.weight	84.09664	NA	1.075396	NA	(NA)	5	-19465.63	38952.08	0.1527718	0.1943631
1	101.30179	NA	NA	NA	(NA)	4	-19465.59	38952.09	0.1587259	0.1937853

Best model summary

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

# 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)	84.9660852	48.977501	122.674261	9700.000
lift.weight	1.0075913	-1.218165	3.364722	9700.000
PC1	0.7948061	-2.060258	3.951657	9391.704

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, 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)) +
    geom_jitter(position = position_jitter(width = 0.2, height = 0.2), aes(size = count)) +
  labs(x = "Body size (PC1)", y = "Sound pressure level (dB)", color = "Period", size = "Sample size") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) + 
  # 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)))

ggplot(agg_size, aes(x = lift.weight, y = cal.spl)) +
    geom_jitter(position = position_jitter(width = 0.2, height = 0.2), aes(size = count)) +
  labs(x = "Lifiting power (g)", y = "Sound pressure level (dB)", color = "Period", size = "Sample size") +
  scale_color_viridis_d(begin = 0.1, end = 0.8, alpha = 0.7) +
  theme_classic(base_size = 24) + 
  # 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)))

Each dot represents the average for a single individual

Effect size posterior distribution

# 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)

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)\]

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 ~ Treatment  + (Treatment|ID), data = playback_dat, REML = FALSE)

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

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

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

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

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)

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

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

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

ggplot(agg_playback, aes(x = Treatment, y = cal.spl, color = Treatment)) +
  geom_point(size = 2, show.legend = FALSE) +
   geom_errorbar(width=.05, aes(ymin = cal.spl - sd, ymax = cal.spl + 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)

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)

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)\]

# 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 ~ Treatment  + (Treatment|ID), data = coord_data, REML = FALSE)

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

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

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)

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

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

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

ggplot(agg_coord, aes(x = Treatment, y = cal.spl, color = Treatment)) +
  geom_point(size = 2, show.legend = FALSE) +
   geom_errorbar(width=.05, aes(ymin = cal.spl - sd, ymax = cal.spl + 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)

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

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)\]

# 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[dat_interaction$year < 2021] <- dat_interaction$cal.spl[dat_interaction$year < 2021] + mean(dat_interaction$cal.spl[dat_interaction$year == 2021]) - mean(dat_interaction$cal.spl[dat_interaction$year < 2021])

dat_aggresive_inter <- rbind(dat_interaction, dat_chase)

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

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

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

# 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 ~ ", "", interaction_formulas)

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

Model selection

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

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

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

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

model_selection

	(Intercept)	inter_treatment	inter_type	inter_treatment:inter_type	family	df	logLik	DIC	delta	weight
inter_treatment * inter_type	100.40099	+	+	+	(NA)	6	-3020.352	6062.358	0.00000	1.000000e+00
inter_treatment + inter_type	99.66253	+	+	NA	(NA)	5	-3057.029	6134.672	72.31389	1.982618e-16
inter_treatment	100.08454	+	NA	NA	(NA)	4	-3072.151	6163.933	101.57488	8.775953e-23
1	101.15672	NA	NA	NA	(NA)	3	-3158.309	6335.201	272.84253	5.662346e-60

Best model summary

# 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

	post.mean	l-95% CI	u-95% CI	eff.samp
(Intercept)	100.400990	98.4147789	102.329394	9700
inter_treatment2_After	1.265453	0.7989197	1.747272	9700
inter_type2_Chase	1.004035	0.3499950	1.652210	9700
inter_treatment2_After:inter_type2_Chase	3.121944	2.4057612	3.805274	9700

# 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 ~ 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)) +
  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)

Effect size posterior distribution

# 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"))

# plot posteriors
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)

mcmc1	mcmc2	p
Before-Perch interaction	After-Perch interaction	0
Before-Chase	After-Chase	0

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

SPL and song structure

Song structure measured as spectrographic consistency and gap duration

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)\]

agg_amp_inter <- aggregate(cal.spl ~ 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 -  
aggregate(cal.spl ~ ID + inter_type, data = dat_aggresive_inter[dat_aggresive_inter$inter_treatment == "1_Before" & dat_aggresive_inter$extsn != "mp3",], mean)$cal.spl  

# 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], ]

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), 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), 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")

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 ~ ", "", song_consistency_formulas)

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

Model selection

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

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

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

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

model_selection

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

Best model summary

# 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

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)

# 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

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)\]

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 ~ ", "", gap_formulas)

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

Model selection

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

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

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

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

model_selection

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

Best model summary

# 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

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

# 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))

# 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) +
 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

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:

SPL and period of the day

mcmc_diagnostics(period_mods)

## [1] "model: Null"

## [1] "model: period"

SPL and condition

mcmc_diagnostics(size_mods)

## [1] "model: Null"

## [1] "model: PC1"

## [1] "model: lift.weight"

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

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

SPL and interactions

mcmc_diagnostics(inter_mods)

## [1] "model: Null"

## [1] "model: inter_treatment"

## [1] "model: inter_treatment + inter_type"

## [1] "model: inter_treatment * inter_type"

SPL and spectrographic consistency

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"

SPL and gap duration

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"

Session information

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=pt_PT.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=pt_PT.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=pt_PT.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       
## 
## attached base packages:
## [1] grid      parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] gridExtra_2.3      lme4_1.1-26        corrplot_0.84      MuMIn_1.43.17     
##  [5] MCMCglmm_2.32      ape_5.4-1          coda_0.19-4        Matrix_1.2-18     
##  [9] rptR_0.9.22        readxl_1.3.1       viridis_0.5.1      viridisLite_0.3.0 
## [13] Rraven_1.0.12      warbleR_1.1.27     NatureSounds_1.0.3 knitr_1.31        
## [17] seewave_2.1.6      tuneR_1.3.3        ggplot2_3.3.3      pbapply_1.4-3     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6        lattice_0.20-41   corpcor_1.6.9     fftw_1.0-6       
##  [5] digest_0.6.27     R6_2.5.0          cellranger_1.1.0  signal_0.7-6     
##  [9] stats4_4.0.3      evaluate_0.14     highr_0.8         pillar_1.4.6     
## [13] rlang_0.4.10      cubature_2.0.4.1  minqa_1.2.4       nloptr_1.2.2.2   
## [17] rmarkdown_2.4     labeling_0.4.2    splines_4.0.3     statmod_1.4.35   
## [21] stringr_1.4.0     RCurl_1.98-1.3    munsell_0.5.0     proxy_0.4-25     
## [25] compiler_4.0.3    xfun_0.22         pkgconfig_2.0.3   htmltools_0.5.1.1
## [29] tidyselect_1.1.0  tibble_3.0.4      tensorA_0.36.2    dtw_1.22-3       
## [33] crayon_1.4.1      dplyr_1.0.2       withr_2.4.1       MASS_7.3-53      
## [37] bitops_1.0-6      nlme_3.1-149      gtable_0.3.0      lifecycle_0.2.0  
## [41] magrittr_2.0.1    scales_1.1.1      stringi_1.5.3     farver_2.1.0     
## [45] ellipsis_0.3.1    generics_0.1.0    vctrs_0.3.6       boot_1.3-25      
## [49] rjson_0.2.20      tools_4.0.3       glue_1.4.2        purrr_0.3.4      
## [53] yaml_2.2.1        colorspace_2.0-0

Sound pressure level stats

Long-billed hermit song amplitude

Marcelo Araya-Salas, PhD & Melanie Talavera

20-04-2021

What the code does:

Next steps

Load packages

Functions and parameters

Repeatability

Hypothesis testing analysis

SPL vs day period (morning / afternoon)

Model selection

Best model summary

SPL vs condition

Correlation between body size parameters included in the PCA

Variance explained by each PC

Model selection

Best model summary

SPL vs playback

SPL vs coordination

SPL vs aggresive interaction

Model selection

Best model summary

SPL and song structure

Song spectrographic consistency

Model selection

Best model summary

Song gap duration

Model selection

Best model summary

Diagnostic plots for MCMCglmm models

SPL and period of the day

SPL and condition

SPL and interactions

SPL and spectrographic consistency

SPL and gap duration