Social learning is often diagnosed by mapping the geographic variation of behavior. Behavioral variation at a small geographical scale that shows both sharp differences among localities and consistency within localities is indicative of social learning of local traditions. This pattern translates into a pretty straightforward statistical hypothesis: the behavior is more similar within than between groups (although absence of this pattern doesn’t necessarily imply a lack of learning!). In other words, if there is social learning going on, we can expect a group level signature.

This post is about how to test this pattern in vocal signals (i.e. testing vocal learning) using pairwise similarities derived from spectrographic cross-correlation. In fact, that’s exactly what we did in our recent paper on the co-ocurrence of social learning in vocal and visual signals in the long billed hermit. So the post reproduces the code I used for the acoustic analysis in that paper.

Male long-billed hermits have a vocal repertoire consisting of a single song type. Different song types are found at different sites, leks and even within leks: plot of chunk song 1

Location of three study sites in northeastern Costa Rica. Spectrograms of the 10 observed song types grouped by lek are also shown. Maps of leks at La Selva Biological Station are shown in greater detail in the lower left map. The two song neighbourhoods at lek SUR are shown in lower right map (similar song neighbourhoods were found at the other leks with 2 song types—SJA, TR2, and HC1—but are not shown); polygons represent lekking male territories and the coloured borders delineate the song neighbourhoods corresponding to the coloured borders of the spectrograms at the far right.

All annotations and acoustic data used on that paper were made available on Dryad. We just need to download the extended selection table (R object including acoustic data + annotations) from Dryad and unzip the file as follows (it could take a while):

# set temporary working directory
setwd(tempdir())

url <- "https://datadryad.org/bitstream/handle/10255/dryad.216487/extended%20selection%20table%20LBH%20songs.zip?sequence=2"

download.file(url = url, destfile = "lbh_est.zip", mode="wb")

unzip("lbh_est.zip")

Now we can read the file and take a look at the data:

library(warbleR)

lbh_est <- readRDS("extended selection table LBH songs.RDS")

lbh_est
## object of class 'extended_selection_table' 
##  contains a selection table data frame with 7646 rows and 14 columns: 
##                  sound.files selec start     end comm bottom.freq top.freq
## 1 0.HC1.2013.6.9.10.08.WAV_1     1   0.1 0.25607    A      1.9493   10.747
## 2 0.HC1.2013.6.9.10.08.WAV_2     1   0.1 0.25607    A      1.9493   10.747
## 3 0.HC1.2013.6.9.10.08.WAV_3     1   0.1 0.25607    A      1.9493   10.747
## 4 0.HC1.2013.6.9.10.08.WAV_4     1   0.1 0.25607    A      1.9493   10.747
## 5 0.HC1.2013.6.9.10.08.WAV_5     1   0.1 0.25607    A      1.9493   10.747
## 6 0.HC1.2013.6.9.10.08.WAV_6     1   0.1 0.25607    A      1.9493   10.747
##              file.name tailored song.type Lek Bird.ID old.selec
## 1 0.HC1.2013.6.9.10.08        y         A HC1  Male17   1+AC0-1
## 2 0.HC1.2013.6.9.10.08        y         A HC1  Male17   1+AC0-1
## 3 0.HC1.2013.6.9.10.08        y         A HC1  Male17   1+AC0-1
## 4 0.HC1.2013.6.9.10.08        y         A HC1  Male17   1+AC0-1
## 5 0.HC1.2013.6.9.10.08        y         A HC1  Male17   1+AC0-1
## 6 0.HC1.2013.6.9.10.08        y         A HC1  Male17   1+AC0-1
##   lek.song.type
## 1         HC1-A
## 2         HC1-A
## 3         HC1-A
## 4         HC1-A
## 5         HC1-A
## 6         HC1-A
## ... and 7640 more rows 
## 7646 wave objects (as attributes): 
## [1] "0.HC1.2013.6.9.10.08.WAV_1" "0.HC1.2013.6.9.10.08.WAV_2"
## [3] "0.HC1.2013.6.9.10.08.WAV_3" "0.HC1.2013.6.9.10.08.WAV_4"
## [5] "0.HC1.2013.6.9.10.08.WAV_5" "0.HC1.2013.6.9.10.08.WAV_6"
## ... and 7640 more 
## and a data frame (check.results) generated by checkres() (as attribute) 
## the selection table was created by element (see 'class_extended_selection_table')

The data contains several songs per individual, and several individuals per lek. This is a pretty big data set, so it takes a while to run the analysis. To speed it up a bit (and avoid pseudoreplication!), we will keep only the highest signal-to-noise ratio song for each individual:

# set warbleR global options 
warbleR_options(wl = 200)

# measure SNR
lbh_est <- sig2noise(lbh_est, mar = 0.1)

# subset ext. sel. tab.
sub_lbh_est <- lbh_est[ave(x = lbh_est$SNR, 
    paste0(lbh_est$Lek, lbh_est$Bird.ID), 
    FUN = function(x) rank(x, ties.method = "first")) == 1, ]

Now we can run the cross-correlation analysis as follows:

xc_mat <- x_corr(sub_lbh_est)

The output is a similarity matrix with dimensions 60 x 60 (for simplicity only the first 20 columns/rows are shown):

0.HC1.2013.6.9.10.08.WAV_7-1 0.HC1.2013.6.9.7.20.WAV_73-1 0.HC1.2013.6.9.7.31.WAV_218-1 0.HC1.2013.6.9.7.50.WAV_334-1 0.HC1.2013.6.9.8.09.WAV_605-1 0.HC1.2013.6.9.8.18.WAV_635-1 0.HC1.2013.6.9.8.28.WAV_815-1 0.HC1.2013.6.9.8.49.WAV_835-1 0.HC1.2013.6.9.8.53.WAV_935-1 0.HC1.2013.6.9.9.58.WAV_1016-1 0.TR2.2013.4.16.8.53.WAV_1161-1 0.TR2.2013.4.16.8.59.WAV_1197-1 0.TR2.2013.4.16.9.25.WAV_1281-1 0.TR2.2013.4.17.8.49.WAV_1301-1 0.TR2.2013.4.17.9.00.WAV_1501-1 0.TR2.2013.4.17.9.02.WAV_1528-1 0.TR2.2013.4.17.9.05.WAV_1639-1 0.TR2.2013.4.17.9.09.WAV_1657-1 0.TR2.2013.4.17.9.18.WAV_1722-1 108.LOC.2013.3.26.8.15.WAV_1810-1
0.HC1.2013.6.9.10.08.WAV_7-1 1.00000 0.63182 0.63959 0.55934 0.68051 0.66121 0.64135 0.69259 0.66024 0.71202 0.69887 0.51545 0.38529 0.38876 0.39564 0.39147 0.48834 0.51111 0.52910 0.37781
0.HC1.2013.6.9.7.20.WAV_73-1 0.63182 1.00000 0.75495 0.54539 0.66660 0.66511 0.61887 0.69640 0.66803 0.63830 0.60416 0.54407 0.41529 0.45166 0.48237 0.42941 0.42745 0.51071 0.46030 0.46680
0.HC1.2013.6.9.7.31.WAV_218-1 0.63959 0.75495 1.00000 0.51396 0.64127 0.67599 0.61162 0.71097 0.69376 0.65245 0.66341 0.49186 0.35839 0.38409 0.42478 0.40135 0.44090 0.52388 0.46107 0.41861
0.HC1.2013.6.9.7.50.WAV_334-1 0.55934 0.54539 0.51396 1.00000 0.67408 0.62532 0.64077 0.65590 0.58549 0.52716 0.54930 0.44475 0.30083 0.28660 0.41001 0.37448 0.49695 0.57942 0.50946 0.19398
0.HC1.2013.6.9.8.09.WAV_605-1 0.68051 0.66660 0.64127 0.67408 1.00000 0.73936 0.75308 0.69222 0.64462 0.66669 0.62863 0.53840 0.42151 0.42051 0.50448 0.48591 0.55599 0.60296 0.52239 0.35289
0.HC1.2013.6.9.8.18.WAV_635-1 0.66121 0.66511 0.67599 0.62532 0.73936 1.00000 0.71308 0.69201 0.66034 0.66862 0.64829 0.48065 0.29665 0.31890 0.39544 0.35343 0.44982 0.55360 0.48631 0.28779
0.HC1.2013.6.9.8.28.WAV_815-1 0.64135 0.61887 0.61162 0.64077 0.75308 0.71308 1.00000 0.66206 0.59760 0.61746 0.55406 0.46961 0.39846 0.39799 0.47109 0.43293 0.51701 0.55178 0.44557 0.37665
0.HC1.2013.6.9.8.49.WAV_835-1 0.69259 0.69640 0.71097 0.65590 0.69222 0.69201 0.66206 1.00000 0.79225 0.66612 0.60694 0.49104 0.35715 0.36894 0.43598 0.40826 0.51970 0.54033 0.51483 0.33840
0.HC1.2013.6.9.8.53.WAV_935-1 0.66024 0.66803 0.69376 0.58549 0.64462 0.66034 0.59760 0.79225 1.00000 0.67312 0.60852 0.46700 0.32351 0.32463 0.41088 0.37981 0.48276 0.53686 0.53730 0.30993
0.HC1.2013.6.9.9.58.WAV_1016-1 0.71202 0.63830 0.65245 0.52716 0.66669 0.66862 0.61746 0.66612 0.67312 1.00000 0.63694 0.48499 0.34605 0.35287 0.40963 0.36159 0.47811 0.52213 0.52346 0.36947
0.TR2.2013.4.16.8.53.WAV_1161-1 0.69887 0.60416 0.66341 0.54930 0.62863 0.64829 0.55406 0.60694 0.60852 0.63694 1.00000 0.47882 0.30888 0.30329 0.36293 0.35559 0.43950 0.55284 0.51489 0.32043
0.TR2.2013.4.16.8.59.WAV_1197-1 0.51545 0.54407 0.49186 0.44475 0.53840 0.48065 0.46961 0.49104 0.46700 0.48499 0.47882 1.00000 0.53392 0.54120 0.54787 0.47744 0.47772 0.53889 0.57507 0.42549
0.TR2.2013.4.16.9.25.WAV_1281-1 0.38529 0.41529 0.35839 0.30083 0.42151 0.29665 0.39846 0.35715 0.32351 0.34605 0.30888 0.53392 1.00000 0.75587 0.63558 0.56687 0.51914 0.43910 0.48098 0.50214
0.TR2.2013.4.17.8.49.WAV_1301-1 0.38876 0.45166 0.38409 0.28660 0.42051 0.31890 0.39799 0.36894 0.32463 0.35287 0.30329 0.54120 0.75587 1.00000 0.65570 0.55677 0.46035 0.40279 0.44048 0.53632
0.TR2.2013.4.17.9.00.WAV_1501-1 0.39564 0.48237 0.42478 0.41001 0.50448 0.39544 0.47109 0.43598 0.41088 0.40963 0.36293 0.54787 0.63558 0.65570 1.00000 0.56429 0.55023 0.52848 0.45378 0.52557
0.TR2.2013.4.17.9.02.WAV_1528-1 0.39147 0.42941 0.40135 0.37448 0.48591 0.35343 0.43293 0.40826 0.37981 0.36159 0.35559 0.47744 0.56687 0.55677 0.56429 1.00000 0.50411 0.44731 0.40129 0.41843
0.TR2.2013.4.17.9.05.WAV_1639-1 0.48834 0.42745 0.44090 0.49695 0.55599 0.44982 0.51701 0.51970 0.48276 0.47811 0.43950 0.47772 0.51914 0.46035 0.55023 0.50411 1.00000 0.63954 0.53222 0.27489
0.TR2.2013.4.17.9.09.WAV_1657-1 0.51111 0.51071 0.52388 0.57942 0.60296 0.55360 0.55178 0.54033 0.53686 0.52213 0.55284 0.53889 0.43910 0.40279 0.52848 0.44731 0.63954 1.00000 0.59617 0.25242
0.TR2.2013.4.17.9.18.WAV_1722-1 0.52910 0.46030 0.46107 0.50946 0.52239 0.48631 0.44557 0.51483 0.53730 0.52346 0.51489 0.57507 0.48098 0.44048 0.45378 0.40129 0.53222 0.59617 1.00000 0.23123
108.LOC.2013.3.26.8.15.WAV_1810-1 0.37781 0.46680 0.41861 0.19398 0.35289 0.28779 0.37665 0.33840 0.30993 0.36947 0.32043 0.42549 0.50214 0.53632 0.52557 0.41843 0.27489 0.25242 0.23123 1.00000

We will need a second matrix representing lek membership. It has to be a pairwise matrix in which 0 denotes pairs of individuals that belong to the same lek and 1 pairs that belong to different leks. The following function creates this type of matrix:

#function to create group membership binary matrix
bi_mats <- function(X, labels) {
  
  # create empty matrix to store memebership matrix
  mat <- matrix(nrow = ncol(X), ncol = ncol(X))
 
  # add labels to row and col names
  rownames(mat) <- colnames(mat) <- labels
  
  # add 0 if same lek and 1 if else 
  out <- lapply(1:(length(labels) - 1), function(i){
  sapply((i + 1):length(labels), function(j) 
    if (labels[i] == labels[j]) 0 else 1)  
    })

  # add to mat
  mat[lower.tri(mat)] <- unlist(out)

  # retunr as distance matrix
  return(as.dist(mat))
  }

The function takes as arguments the cross-correlation similarity matrix and a label indicating lek membership:

# create lek membership from column names
lbls <- sapply(strsplit(colnames(xc_mat), ".", fixed = TRUE), "[[", 2)

# create binary matrix
lek_bi_mat <- bi_mats(xc_mat, lbls)

# look at the first 15 rows/cols
as.matrix(lek_bi_mat)[1:15, 1:15]
##     HC1 HC1 HC1 HC1 HC1 HC1 HC1 HC1 HC1 HC1 TR2 TR2 TR2 TR2 TR2
## HC1   0   0   0   0   0   0   0   0   0   0   1   1   1   1   1
## HC1   0   0   0   0   0   0   0   0   0   0   1   1   1   1   1
## HC1   0   0   0   0   0   0   0   0   0   0   1   1   1   1   1
## HC1   0   0   0   0   0   0   0   0   0   0   1   1   1   1   1
## HC1   0   0   0   0   0   0   0   0   0   0   1   1   1   1   1
## HC1   0   0   0   0   0   0   0   0   0   0   1   1   1   1   1
## HC1   0   0   0   0   0   0   0   0   0   0   1   1   1   1   1
## HC1   0   0   0   0   0   0   0   0   0   0   1   1   1   1   1
## HC1   0   0   0   0   0   0   0   0   0   0   1   1   1   1   1
## HC1   0   0   0   0   0   0   0   0   0   0   1   1   1   1   1
## TR2   1   1   1   1   1   1   1   1   1   1   0   0   0   0   0
## TR2   1   1   1   1   1   1   1   1   1   1   0   0   0   0   0
## TR2   1   1   1   1   1   1   1   1   1   1   0   0   0   0   0
## TR2   1   1   1   1   1   1   1   1   1   1   0   0   0   0   0
## TR2   1   1   1   1   1   1   1   1   1   1   0   0   0   0   0

The 2 matrices are then input into a Mantel test to evaluate if acoustic similarity is higher within than between leks. Note that the cross-correlation matrix is transformed into a distance matrix by subtracting it from 1:

# install vegan if necessary
library(vegan)

# convert xcorr mat to distance
xc_dist <- as.dist(1 - xc_mat)

# run mantel test
mantel(xc_dist, lek_bi_mat, permutations = 10000)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = xc_dist, ydis = lek_bi_mat, permutations = 10000) 
## 
## Mantel statistic r: 0.38 
##       Significance: 1e-04 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0299 0.0403 0.0507 0.0620 
## Permutation: free
## Number of permutations: 10000

That’s a pretty solid association between lek membership and acoustic similarity, r = 0.37. So there is a lek level acoustic signature.

Using other metrics

The same test can be done using other acoustic structure metrics. For instance dynamic time warping returns a distance matrix which can be directly input into the mantel test:

# measure DTW distances
dtw_dist <- df_DTW(sub_lbh_est, img = FALSE)

# run mantel test
mantel(as.dist(dtw_dist), lek_bi_mat, permutations = 10000)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = as.dist(dtw_dist), ydis = lek_bi_mat, permutations = 10000) 
## 
## Mantel statistic r: 0.278 
##       Significance: 1e-04 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0307 0.0404 0.0491 0.0593 
## Permutation: free
## Number of permutations: 10000

We can also use non-pairwise metrics of acoustic structure like spectrographic parameters or descriptors of cepstral coefficients. However, this would need an extra step for converting those metrics into a distance matrix. We can do that with the base R function dist():

# measure spectrographic parameters
sp <- specan(sub_lbh_est)

# create distance matrix
dist_sp <- dist(sp[, -4:-1])

# run mantel test
mantel(dist_sp, lek_bi_mat, permutations = 10000)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = dist_sp, ydis = lek_bi_mat, permutations = 10000) 
## 
## Mantel statistic r: 0.0482 
##       Significance: 0.0357 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0311 0.0425 0.0532 0.0649 
## Permutation: free
## Number of permutations: 10000
# measure descriptors of mel cepstral coefficients
mfcc <- mfcc_stats(sub_lbh_est, bp = c(1, 11), nbands = 20)

# create distance matrix
dist_cc <- dist(mfcc[, -2:-1])

# run mantel test
mantel(dist_cc, lek_bi_mat, permutations = 10000)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = dist_cc, ydis = lek_bi_mat, permutations = 10000) 
## 
## Mantel statistic r: -0.0297 
##       Significance: 0.906 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0313 0.0418 0.0513 0.0630 
## Permutation: free
## Number of permutations: 10000

Mantel tests using dynamic time warping and spectrographic parameters were also able to detect lek level signatures. That was not the case for cepstral coefficients.

That’s it!

Session information 
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 19.04
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] vegan_2.5-5        lattice_0.20-38    permute_0.9-5     
## [4] kableExtra_1.1.0   warbleR_1.1.16     NatureSounds_1.0.1
## [7] seewave_2.1.3      tuneR_1.3.3        maps_3.3.0        
## 
## loaded via a namespace (and not attached):
##  [1] xfun_0.7          pbapply_1.4-0     splines_3.6.1    
##  [4] colorspace_1.4-1  htmltools_0.3.6   viridisLite_0.3.0
##  [7] mgcv_1.8-28       rlang_0.3.4       pracma_2.2.5     
## [10] pillar_1.4.0      glue_1.3.1        jpeg_0.1-8       
## [13] stringr_1.4.0     munsell_0.5.0     rvest_0.3.3      
## [16] evaluate_0.14     knitr_1.23        fftw_1.0-5       
## [19] parallel_3.6.1    highr_0.8         Rcpp_1.0.1       
## [22] readr_1.3.1       scales_1.0.0      webshot_0.5.1    
## [25] soundgen_1.4.0    Sim.DiffProc_4.3  Deriv_3.8.5      
## [28] rjson_0.2.20      hms_0.4.2         packrat_0.5.0    
## [31] digest_0.6.19     stringi_1.4.3     dtw_1.20-1       
## [34] grid_3.6.1        tools_3.6.1       bitops_1.0-6     
## [37] magrittr_1.5      RCurl_1.95-4.12   proxy_0.4-23     
## [40] tibble_2.1.1      cluster_2.1.0     crayon_1.3.4     
## [43] pkgconfig_2.0.2   Matrix_1.2-17     MASS_7.3-51.4    
## [46] xml2_1.2.0        rmarkdown_1.13    httr_1.4.0       
## [49] rstudioapi_0.10   iterators_1.0.10  R6_2.4.0         
## [52] nlme_3.1-140      signal_0.7-6      compiler_3.6.1