1.1. Managing sound files: ss_find_files and ss_assess_files

Before getting into the details of the workflow, let’s talk about sound files, file naming, and the expected file structure. At this point, we assume you’ve successfully gone out and collected soundscape recordings for a certain period of time at one or multiple sites of interest. These files are now stored somewhere on your device, external hard drive, or HPC.

To manage these sound files and extract the required information from them, the soundscapeR workflow expects files to be named using the following format:

SITENAME_YYYYMMDD_HHMMSS.wav

Moreover, the soundscapeR package expects the sound files to be stored in directories following one of the following structures:

Figure 4: A schematic overview of the expected file structure

Now, we can tell soundscapeR where to look for these files using the ss_find_files function. This function will automatically detect all the soundfiles in a parent directory and all subdirectories, and save the location to each detected sound file in a list that maintains the directory architecture shown in Fig 4.

Let’s take a look at the input parameters of the ss_find_files function:

As we can see, the ss_find_files function only takes a single input argument: the name of a directory containing .wav files, or a parent directory containing named directories with .wav files.

Let’s give this a try:

  # Find the file locations

filelocs <- ss_find_files(parent_directory = output_dir)
names(filelocs) <- c("Andre", "Mascote_A1", "Mascote_A2", "Mascote_B1", "Mascote_B2")

  # Show us what the fileloc output looks like

summary(filelocs)

##            Length Class  Mode     
## Andre      1440   -none- character
## Mascote_A1 1440   -none- character
## Mascote_A2 1440   -none- character
## Mascote_B1 1440   -none- character
## Mascote_B2 1440   -none- character

head(filelocs[[1]])

## [1] "G:/soundscapeR_paper_data/Andre/G10_Andre_20151010_155500Z.wav"
## [2] "G:/soundscapeR_paper_data/Andre/G10_Andre_20151010_160000Z.wav"
## [3] "G:/soundscapeR_paper_data/Andre/G10_Andre_20151010_160500Z.wav"
## [4] "G:/soundscapeR_paper_data/Andre/G10_Andre_20151010_161000Z.wav"
## [5] "G:/soundscapeR_paper_data/Andre/G10_Andre_20151010_161500Z.wav"
## [6] "G:/soundscapeR_paper_data/Andre/G10_Andre_20151010_162000Z.wav"

The ss_find_files function created a list of length 5 (one for each soundscape plot in our case study). Each list element contains a 1,440 character strings which gives the file path to each 1-minute sound file for each soundscape.

However, it is possible that something went wrong during the sound recording stage, or when we transferred the files from our acoustic sensor to our storage device, leading to missing files or gaps in the soundscape recordings. This is where the ss_assess_files function comes into play. The ss_assess_files function will automatically detect the sampling regime of the study, and based on this, will check whether: (i) any files are missing; and (ii) all soundscapes in the filelocs list had the same sampling regime.

Let’s give this a try:

filelocs <- ss_assess_files(file_locs = filelocs, 
                            full_days = FALSE)

## Sampling regime for study: 1 minute in 5 minutes 
## Number of files per day: 288 files 
## The data still contains partially sampled days - be careful with downstream analysis...

As we can see, the ss_assess_files function successfully detected that our sound was sampled for 1 minute every 5 minutes, resulting in 288 files per day. Since no other error message was thrown, it appears no files were missing.

Let’s take a look at what would happen if files were missing:

  # Simulate a soundscape with missing data

filelocs_missing <- filelocs
filelocs_missing[[1]] <- filelocs_missing[[1]][-seq(1, length(filelocs_missing[[1]]), 5)]

  # Run the ss_assess_files function

filelocs_missing <- ss_assess_files(file_locs = filelocs_missing, 
                                    full_days = FALSE)

## Irregular timeinterval detected for:  Andre 
## Based on the expected sampling regime, there are missing files...

## Error in ss_assess_files(file_locs = filelocs_missing, full_days = FALSE): Irregular timeintervals detected - check files

As you can see, the ss_assess_files function generates an error message if the files in fileloc deviate from the expected sampling regime.

Finally, using the ss_assess_files function, we can also remove any sound files from the fileloc object that are not part of a full soundscape sample (24h period). We do this because, otherwise, we would have an uneven number of samples per time of day. To do so, we use the full_days = TRUE argument.

Let’s give this a try:

  # Create a soundscape with sound files that are not part of a complete soundscape sample

filelocs_incomplete <- filelocs
filelocs_incomplete[[1]] <- filelocs[[1]][5:length(filelocs[[1]])]

  # Run ss_assess_files again

filelocs_incomplete_fixed <- ss_assess_files(file_locs = filelocs_incomplete, 
                            full_days = TRUE)

## Sampling regime for study: 1 minute in 5 minutes 
## Number of files per day: 288 files 
## The file_locs argument has been subsetted to contain only full sampling days

It seems like there was an incomplete soundscape sample, ss_assess_files removed some soundfiles so only complete soundscape samples (24h period) were retained.

Let’s have a look at the before and after:

# Before

length(filelocs_incomplete[[1]])

## [1] 1436

# After

length(filelocs_incomplete_fixed[[1]])

## [1] 1152

The ss_assess_files function has detected an incomplete soundscape sample and removed sound files until the ‘previous’ complete sample was found - in this case 1,152 files, or 4 full sampling days (288 * 4 = 1,152).

Great, now that we know how these two file management functions work, we can proceed!

1.2. Calculate ecoacoustic indices: The ss_index_calc function

# ss_index_calc(fileloc = "filepath_here",
#               outputloc = "filepath_here",
#               samplerate = 44100,
#               window = 256,
#               parallel = FALSE)

  # 1. Download the metajam package 

# install.packages("metajam")

  # 2. Set folder location to save data

setwd("G:/soundscapeR_case_study")

dir.create(paste0(getwd(), "/case_study_data"))

output_dir <- paste0(getwd(), "/case_study_data")

  # 3. Download data and metadata

# metajam::download_d1_data(data_url = "https://knb.ecoinformatics.org/knb/d1/mn/v2/object/urn%3Auuid%3Acfccc678-27f3-473a-b34a-b4b5d4371e6b", path = ".")
# 
#   # 4. Unzip
# 
# to_unzip <- list.files(output_dir, full.names = TRUE)
# 
# sapply(to_unzip, function(x) unzip(zipfile = x,exdir = output_dir, overwrite = TRUE))

  # Specify the folder locations

folder_locs <- list.dirs(path = output_dir, recursive = FALSE)

island_metadata <- data.frame(plot = c("Andre", "Mascote_A1", "Mascote_A2", "Mascote_B1", "Mascote_B2"), 
                              first_day = c("2015-10-10", rep("2015-10-16", 4)), 
                              lat = c(-1.58462, -1.64506, -1.6489, -1.64406, -1.65936), 
                              lon = c(-59.87211, -59.82035, -59.83297, -59.84817, -59.83546))

island_metadata$folderloc <- folder_locs

1.3. Create a soundscape: The ss_create function

balbina <- vector("list")

for (i in 1:nrow(island_metadata)){
  
  balbina[[i]] <- ss_create(fileloc = island_metadata$folderloc[i], 
                            samplerate = 44100, 
                            window = 256, 
                            index = "CVR", 
                            date = island_metadata$first_day[i], 
                            lat = island_metadata$lat[i], 
                            lon = island_metadata$lon[i], 
                            method = "IsoData", 
                            output = "incidence_freq")
  
}
                            

names(balbina) <- c("Andre", "Mascote_A1", "Mascote_A2", "Mascote_B1", "Mascote_B2")

1.4. Introducing the soundscape object

As we previously mentioned, the ss_create function requires us to provide some additional metadata information, as well as several parameters on how the data should be processed. Using the ss_create function, we are creating a new type of data object which stores all the relevant information, and will form the basis on which all subsequent functions in the workflow are performed!

Let’s take a look at this newly created object for Andre island:

# Let's see what class this object is:

summary(balbina[["Andre"]])

##     Length      Class       Mode 
##          1 soundscape         S4

This new data object is an ‘S4’ object of the class soundscape.

# Let's see what sort of information this object holds

balbina[["Andre"]]

## 
## 1.  Soundscape metadata 
## 
##     Sampling point metadata:  
## 
##     First day of recording:  2015-10-10 
##     Latitude of sampling point:  -1.58462 
##     Longitude of sampling point:  -59.87211 
##     Time zone of sampling point:  America/Manaus 
##     Sunrise time at sampling point:  05:43:36 
##     Sunset time at sampling point:  17:51:47 
## 
##     Acoustic index metadata:  
## 
##     Path to raw sound files:  G:/soundscapeR_paper_data/Andre 
##     Spectral index used:  CVR 
##     Sampling rate of the recording:  44100  Hz 
##     Window size used in FFT:  256  samples 
##     Frequency resolution:  172.2656  Hz 
##     Temporal resolution:  0.005804989  ms 
## 
##     Data frame binarization metadata:  
## 
##     Used binarization algorithm:  IsoData 
##     Binarization threshold:  0.16 
## 
##     Aggregated data frame metadata:  
## 
##     Output format:  incidence_freq 
##     Data frame frequency range:  172 - 22050 Hz 
##     Data frame time range:  00:00:00 - 23:55:00 
## 
## 2.  Soundscape data 
## 
##     Merged data frame data:  
## 
## Columns 1 to 5 and rows 1 to 5 displayed 
## 
##         15:55:00   16:00:00   16:05:00   16:10:00   16:15:00
## 22050 0.19408888 0.18264127 0.18909356 0.18763659 0.20761786
## 21877 0.45009887 0.48069518 0.49037361 0.48579457 0.51306067
## 21705 0.49172651 0.50431887 0.50005203 0.50993860 0.51930482
## 21533 0.31439276 0.34571756 0.31616193 0.33884900 0.32802581
## 21360 0.01966906 0.02580914 0.01821209 0.02435217 0.01842023
## 
##     Binarized data frame data:  
## 
## Columns 1 to 5 and rows 1 to 5 displayed 
## 
##       15:55:00 16:00:00 16:05:00 16:10:00 16:15:00
## 22050        1        1        1        1        1
## 21877        1        1        1        1        1
## 21705        1        1        1        1        1
## 21533        1        1        1        1        1
## 21360        0        0        0        0        0
## 
##     Aggregated data frame data:  
## 
##     Aggregated data frame: 
## 
## Columns 1 to 5 and rows 1 to 5 displayed 
## 
##       00:00:00 00:05:00 00:10:00 00:15:00 00:20:00
## 22050      0.6        1      0.6      0.6        1
## 21877      1.0        1      1.0      1.0        1
## 21705      1.0        1      1.0      1.0        1
## 21533      1.0        1      1.0      1.0        1
## 21360      0.0        0      0.0      0.0        0
## 
##     Aggregated data frame per time: 
## 
## First list element displayed:  00:00:00 
## 
## Columns 1 to 5 and rows 1 to 5 displayed 
## 
##       00:00:00 00:00:00.1 00:00:00.2 00:00:00.3 00:00:00.4
## 22050        0          1          0          1          1
## 21877        1          1          1          1          1
## 21705        1          1          1          1          1
## 21533        1          1          1          1          1
## 21360        0          0          0          0          0
## 
##     Number of soundscape samples per time (sampling effort): 
## 
## List elements 1 to 5 displayed 
## 
## $`00:00:00`
## [1] 5
## 
## $`00:05:00`
## [1] 5
## 
## $`00:10:00`
## [1] 5
## 
## $`00:15:00`
## [1] 5
## 
## $`00:20:00`
## [1] 5

As we can see, this objects holds two types of information:

• Metadata: Information regarding the data collection and which parameters where used in the different steps of preparing the soundscape object
  
  

• Soundscape information: Various intermediate data objects that are created as the ss_create function works its way through the index concatenation, binarization and computation of OSU abundance steps.
  
  

To access the information stored in this object, we will use the ‘@’ symbol.

Let’s take a look at which types of data are stored in the object so far:

# Let's check what sort of data collection metadata is stored in the object

print(paste0("First day of data collection: ", balbina[["Andre"]]@first_day))

## [1] "First day of data collection: 2015-10-10"

print(paste0("Latitude at data collection site: ", balbina[["Andre"]]@lat))

## [1] "Latitude at data collection site: -1.58462"

print(paste0("Longitude at data collection site: ", balbina[["Andre"]]@lon))

## [1] "Longitude at data collection site: -59.87211"

print(paste0("Time zone at data collection site: ", balbina[["Andre"]]@tz))

## [1] "Time zone at data collection site: America/Manaus"

print(paste0("Sunrise at time of data collection: ", balbina[["Andre"]]@sunrise))

## [1] "Sunrise at time of data collection: 2015-10-10 05:43:36"

print(paste0("Sunset at time of data collection: ", balbina[["Andre"]]@sunset))

## [1] "Sunset at time of data collection: 2015-10-10 17:51:47"

The ss_create function has automatically calculated a set of important ecological variables, such as sunrise and sunset times, and timezones, based on the first day of recording and geographical coordinates.

Let’s continue looking at the data stored in the soundscape object:

# Let's check what sort of metadata the object has stored regarding past data processing steps

print(paste0("Where are the raw sound files located: ", balbina[["Andre"]]@fileloc))

## [1] "Where are the raw sound files located: G:/soundscapeR_paper_data/Andre"

print(paste0("What acoustic index are we using: ", balbina[["Andre"]]@index, " index"))

## [1] "What acoustic index are we using: CVR index"

print(paste0("What was the samplerate used to collect the data: ", balbina[["Andre"]]@samplerate, " Hz"))

## [1] "What was the samplerate used to collect the data: 44100 Hz"

print(paste0("What was the window length used during the FFT: ", balbina[["Andre"]]@window, " samples"))

## [1] "What was the window length used during the FFT: 256 samples"

The soundscape object has recorded where our raw data files are stored, which acoustic index we’re working with, what sampling rate was used during data collection, and what window length was used during the acoustic index calculation.

Let’s take a look at some of the data frames we created by running the ss_create function:

head(balbina[["Andre"]]@merged_df)[,1:5]

##         15:55:00   16:00:00   16:05:00   16:10:00   16:15:00
## 22050 0.19408888 0.18264127 0.18909356 0.18763659 0.20761786
## 21877 0.45009887 0.48069518 0.49037361 0.48579457 0.51306067
## 21705 0.49172651 0.50431887 0.50005203 0.50993860 0.51930482
## 21533 0.31439276 0.34571756 0.31616193 0.33884900 0.32802581
## 21360 0.01966906 0.02580914 0.01821209 0.02435217 0.01842023
## 21188 0.01207202 0.01800395 0.01186388 0.01592257 0.01321678

As we previously mentioned, the first step performed by the ss_create function is the chronological concatenation of the CVR-index files for a site into a time-frequency-index value data frame. This data frame is stored in the ‘@merged_df’ slot, and contains the time-of-recording as column names, the frequency bins as row names, and the CVR-index for each time-frequency pair as values. Each column contains the spectral index values of a single sound file.

Let’s inspect this data frame a little closer:

# How many columns does the data frame contain?

paste0("The data frame contains: ", ncol(balbina[["Andre"]]@merged_df), " columns")

## [1] "The data frame contains: 1440 columns"

# What are the column names?

head(colnames(balbina[["Andre"]]@merged_df))

## [1] "15:55:00" "16:00:00" "16:05:00" "16:10:00" "16:15:00" "16:20:00"

tail(colnames(balbina[["Andre"]]@merged_df))

## [1] "15:25:00" "15:30:00" "15:35:00" "15:40:00" "15:45:00" "15:50:00"

As we said, the number of columns equals the number of sound files collected during the acoustic survey - in this case, 1440 sound files. The name of each column corresponds to the time of day at which the recording was collected.

Next, let’s take a look at the rows:

# How many rows does the data frame contain?

paste0("The data frame contains: ", nrow(balbina[["Andre"]]@merged_df), " rows")

## [1] "The data frame contains: 128 rows"

# What do these row names look like?

  # The first five names

paste0("The first five rownames: ", paste0(rownames(balbina[["Andre"]]@merged_df)[1:5], collapse = ", "))

## [1] "The first five rownames: 22050, 21877, 21705, 21533, 21360"

  # The last five names

paste0("The last five rownames: ", paste0(rownames(balbina[["Andre"]]@merged_df)[123:128], collapse = ", "))

## [1] "The last five rownames: 1033, 861, 689, 516, 344, 172"

The data frame contains 128 rows, each corresponding to a unique frequency bin. The frequency bins range from 0 - 22,050 Hz, and are of approximately 172 Hz width.

Now, let’s inspect the CVR-index values:

# What is the minimum CVR-index value in our data frame?

paste0("The minimum CVR-value in our data frame is: ", min(balbina[["Andre"]]@merged_df))

## [1] "The minimum CVR-value in our data frame is: 0"

# What is the maximum CVR-index value in our data frame?

paste0("The max CVR-value in our data frame is: ", max(balbina[["Andre"]]@merged_df))

## [1] "The max CVR-value in our data frame is: 0.672286398168384"

As we can see, in our dataset, the CVR-index values range between 0 - 0.67. Remember, CVR-index values capture the proportion of cells in each noise-reduced frequency bin of a sound file that exceeds a 3-dB amplitude threshold. As such, the values can technically range between 0-1.

After the CVR-indices for a soundscape recording of a site have been chronologically concatenated, instead of using the raw CVR-values that we computed for every OSU in each 24h sample of the soundscape, the ss_create function determines a threshold value for each site, and converts the OSU’s CVR-index values into a binary detection (1) / non-detection (0) variable per 24h sample based on this threshold. This step is aimed at detecting the presence of sound for every OSU in each 24h sample of acoustic trait space while removing low-amplitude or transient sounds, which potentially have a non-biological origin, from the data. In this way, we hope to capture the acoustic structure of the soundscape while removing background noise. 

We can see that the soundscape objects produced by the ss_create function has some new metadata related to this binarization step:

# Which threshold algorithm was used for binarization?

paste0("The ", balbina[["Andre"]]@binarization_method, " method was used for binarization")

## [1] "The IsoData method was used for binarization"

# Which threshold was used for binarization?

paste0("The threshold used for binarization was: ", balbina[["Andre"]]@threshold)

## [1] "The threshold used for binarization was: 0.16"

Additionally, a binarized_df data frame was added to the object.

Let’s inspect this data frame:

# What are the dimensions of the binarized dataframe?

dim(balbina[["Andre"]]@binarized_df)

## [1]  128 1440

# What are the unique values contained in this data frame?

unique(unlist(balbina[["Andre"]]@binarized_df))

## [1] 1 0

As expected, the new binarized data frame contains the same number of rows and columns as the merged_df data frame. However, whereas previously we had CVR-index values ranging between anywhere between 0-1, due to the binarization step, the values are now strictly 0 or 1.

Finally, after binarization, the ss_create function calculates the relative abundance of OSUs across all 24h samples of acoustic trait space for a site. For Andre island, we have five 24h samples (5 recording days) per plot.

To do this, the function uses an incidence-based approach. Previously, the ss_create function computed the detection (1) / non-detection (0) of each OSU in each 24h soundscape sample. Next, per OSU, the function takes the mean of this binary variable across all 24h soundscape samples to get the relative frequency by which each OSU was detected. To avoid confusion between the frequency of OSU detection and the sound frequency (in Hz), we’ll refer to this OSU importance value as the ‘relative abundance’.

The soundscape object contains some information related to this step. Let’s investigate the soundscape object for Andre island to see what is new.

What are the dimensions of the new ‘aggregated_df’ data frame?

# First, let's check out the new 'aggregated_df' data frame

paste0("The aggregated_df data frame has ", 
       nrow(balbina[["Andre"]]@aggregated_df), 
       " rows")

## [1] "The aggregated_df data frame has 128 rows"

paste0("The aggregated_df data frame has ", 
       ncol(balbina[["Andre"]]@aggregated_df), 
       " columns")

## [1] "The aggregated_df data frame has 288 columns"

The number of rows in the ‘aggregated_df’ data frame that was added to the soundscape object is still 128, one for each frequency bins resulting from the Fast Fourier Transformation. However, the number of columns in the ‘aggregated_df’ data frame is 288 - that is five times less than the ‘binarized_df’ data frame! Well, since we averaged the detection/non-detection values across five 24h samples of the acoustic trait space, this makes sense…

Now, what range of values do the OSU relative abundances take?

paste0("The relative abundance of OSUs can take the following values: ", 
       paste(sort(unique(unlist(balbina[["Andre"]]@aggregated_df))), collapse = ","))

## [1] "The relative abundance of OSUs can take the following values: 0,0.2,0.4,0.6,0.8,1"

The relative abundance can take 6 values: 0, 0.2, 0.4, 0.6, 0.8 and 1. This makes sense. If an OSU was detected in 0 out of 5 soundscape samples, its relative abundance is 0. If an OSU was detected in 5 out of 5 soundscape samples, its relative abundance is 1. If it was detected somewhere in between, say 3 out of 5 samples, its relative abundance would be 0.6.

Great, we’ve got a soundscape object that contains a dataframe with our a unit of diversity measurement (OSUs) with an associated importance value (the relative abundance). We now have all the ingredients to start the next part of the workflow: exploration and diversity quantification of a single soundscape object.

The remaining functions in the soundscapeR package can be divided into two types:

• Phase 2: Functions for exploring and visualizing the diversity of a single soundscape
  
  

• Phase 3: Functions for visualizing and contrasting the diversity of multiple soundscapes
  
  

2.1. The ss_heatmap function

The ss_heatmap function allows us to visualize the distribution and relative abundance of OSUs in the 24h acoustic trait space. This function is a true workhorse - it is highly flexible, allowing us to subset the soundscape by specific time or frequency coordinates, portray the soundscape in either Cartesian or polar coordinates, annotate the soundscape with the sunset and sunrise time, and many more. Below, we’ll provide an overview of ss_heatmap’s abilities.

Let’s take a look at the soundscape for Andre island. We will start with the ss_heatmap function in it’s most basic form and add progressively more advanced plotting arguments.

ss_heatmap(soundscape_obj = balbina[["Andre"]])

ss_heatmap(soundscape_obj = balbina[["Andre"]],
           type = "polar")

# Regular heatmap

regular_heatmap_annotated <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]],
             type = "regular", 
             annotate = TRUE)

# Polar heatmap

polar_heatmap_annotated <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]],
             type = "polar", 
             annotate = TRUE)

# Combine

library(patchwork)

annotated_heatmaps <- regular_heatmap_annotated + polar_heatmap_annotated

annotated_heatmaps

ss_heatmap(soundscape_obj = balbina[["Andre"]],
           timeinterval = "4 hours", 
           freqinterval = 1000)

## Warning: Removed 542 rows containing missing values (`geom_tile()`).

ss_heatmap(soundscape_obj = balbina[["Andre"]],
           mintime = "06:00:00", 
           maxtime = "18:00:00" , 
           minfreq = 0 , 
           maxfreq = 11000)

## Warning: Removed 27855 rows containing missing values (`geom_tile()`).

# Try the palettes one by one

color_1 <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]],
           palette = "A")

color_2 <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]],
           palette = "B")

color_3 <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]],
           palette = "C")

color_4 <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]],
           palette = "D")

color_5 <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]],
           palette = "E")


# Combine into one plot

library(patchwork)

all_colors <- color_1 + color_2 + color_3 + color_4 + color_5

all_colors

## Warning: Removed 542 rows containing missing values (`geom_tile()`).
## Removed 542 rows containing missing values (`geom_tile()`).
## Removed 542 rows containing missing values (`geom_tile()`).
## Removed 542 rows containing missing values (`geom_tile()`).
## Removed 542 rows containing missing values (`geom_tile()`).

# Try the direction options

direction_1 <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]], 
             palette = "magma", 
             direction = 1)

direction_2 <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]], 
             palette = "magma", 
             direction = -1)

# Combine into one plot

library(patchwork)

all_directions <- direction_1 + direction_2

all_directions

## Warning: Removed 542 rows containing missing values (`geom_tile()`).
## Removed 542 rows containing missing values (`geom_tile()`).

# Try the direction options

zero_black_off <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]], 
             palette = "viridis", 
             zero.black = FALSE)

zero_black_on <- 
  ss_heatmap(soundscape_obj = balbina[["Andre"]], 
             palette = "viridis", 
             zero.black = TRUE)

# Combine into one plot

library(patchwork)

zero_black_all <- zero_black_off + zero_black_on

zero_black_all

## Warning: Removed 542 rows containing missing values (`geom_tile()`).
## Removed 542 rows containing missing values (`geom_tile()`).

2.2. A note on Hill numbers

Now that we know how to visually explore the diversity of our soundscape, let’s get to quantifying how much diversity we’re actually dealing with using the various soundscape metrics available in the soundscapeR package.

To calculate the soundscape metrics, we use the analytical framework of Hill numbers. In brief, Hill numbers provide a unified statistical framework to measure biological diversity in all its facets. The framework is highly robust and flexible, allowing us to quantify different dimensions of diversity (e.g. taxonomic, functional, phylogenetic and soundscape diversity). Doing so, we can answer common scientific questions by measuring, estimating, partitioning, and comparing different diversity types using a common analytical framework. Although a plethora of indices has been proposed to measure diversity, there is a growing consensus that Hill numbers are the most appropriate framework to separate system diversity into its various components.

Want to know more about Hill numbers? For a deep dive, visit this website and this GitHub tutorial!

For an overview of the benefits of Hill numbers over conventional diversity indices, expand the section below. If you’re familiar with Hill numbers, feel free to skip over this part.

Now that we know why we use Hill numbers, let’s have a look at how we can compute our soundscape metrics.

2.3. The ss_diversity and ss_evenness functions

The ss_diversity function is used to compute the soundscape richness and diversity values. The user can modulate the importance of common or rare OSUs on the diversity values using the q-parameter (order of diversity).

The ss_evenness function is used to compute the soundscape evenness. In Luypaert et al. (2022), we mentioned the soundscape evenness was calculated as 2D (soundscape diversiy with q = 2)/0D (soundscape richness with q = 0), following Jost et al. (2010). However, based on recent publications, a more proper way to calculate the soundscape evenness is: (2D - 1) / (0D - 1). Therefore, from here on out, we will only use the latter equation for the calculation of the soundscape evenness.

Both functions allow the soundscape metrics to be computed at a range of different scales and resolutions. For instance, the user can specify custom time-frequency limits, calculate the soundscape metrics for various built-in diurnal-phase presets (dawn, dusk, day, night), or track the soundscape metrics at each unique time of day.

# Soundscape richness

ss_diversity(soundscape_obj = balbina[["Andre"]], 
         qvalue = 0, 
         output = "raw")

## [1] 5577

# Soundscape diversity at q=1

ss_diversity(soundscape_obj = balbina[["Andre"]],
         qvalue = 1, 
        output = "raw")

## [1] 4532.245

# Soundscape diversity at q=2

ss_diversity(soundscape_obj = balbina[["Andre"]],
         qvalue = 2, 
        output = "raw")

## [1] 3931.819

# Soundscape richness

ss_diversity(soundscape_obj = balbina[["Andre"]],
         qvalue = 0, 
         output = "percentage")

## [1] 15.12858

# Soundscape diversity at q=1

ss_diversity(soundscape_obj = balbina[["Andre"]],
         qvalue = 1, 
        output = "percentage")

## [1] 12.2945

# Soundscape diversity at q=2

ss_diversity(soundscape_obj = balbina[["Andre"]],
         qvalue = 2, 
        output = "percentage")

## [1] 10.66574

# Soundscape evenness 

ss_evenness(soundscape_obj = balbina[["Andre"]])

## [1] 0.7049532

# Soundscape richness

  # subset = 'day'

ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             output = "percentage",
             subset = "day")

## [1] 9.364298

  # subset = 'night'

ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             output = "percentage",
             subset = "night")

## [1] 21.05524

  # subset = 'dawn' 

ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             output = "percentage",
             subset = "dawn")

## [1] 22.43924

  # subset = 'dusk'

ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             output = "percentage",
             subset = "dusk")

## [1] 8.767361

# Soundscape evenness

  # subset = 'day'

ss_evenness(soundscape_obj = balbina[["Andre"]], 
             subset = "day")

## [1] 0.6612279

  # subset = 'night'

ss_evenness(soundscape_obj = balbina[["Andre"]], 
             subset = "night")

## [1] 0.7251704

  # subset = 'dawn' 

ss_evenness(soundscape_obj = balbina[["Andre"]], 
             subset = "dawn")

## [1] 0.7056203

  # subset = 'dusk'

ss_evenness(soundscape_obj = balbina[["Andre"]], 
             subset = "dusk")

## [1] 0.6265451

# Soundscape richness

ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             subset = "dawn", 
             dawnstart = 3600, 
             dawnend = 3600)

## [1] 28.1901

ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             subset = "dusk",
             duskstart = 3600, 
             duskend = 3600)

## [1] 11.6862

# Soundscape richness

head(ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             output = "percentage", 
             subset = "tod"), n = 24)

##    soundscape_div time_of_day
## 1        22.65625    00:00:00
## 2        17.18750    00:05:00
## 3        19.53125    00:10:00
## 4        20.31250    00:15:00
## 5        17.18750    00:20:00
## 6        21.09375    00:25:00
## 7        21.09375    00:30:00
## 8        20.31250    00:35:00
## 9        22.65625    00:40:00
## 10       21.87500    00:45:00
## 11       19.53125    00:50:00
## 12       20.31250    00:55:00
## 13       21.09375    01:00:00
## 14       21.09375    01:05:00
## 15       25.00000    01:10:00
## 16       22.65625    01:15:00
## 17       17.18750    01:20:00
## 18       21.09375    01:25:00
## 19       21.09375    01:30:00
## 20       21.09375    01:35:00
## 21       21.09375    01:40:00
## 22       21.09375    01:45:00
## 23       22.65625    01:50:00
## 24       22.65625    01:55:00

ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             output = "percentage",
             mintime = "11:00:00", 
             maxtime = "13:00:00")

## [1] 6.9375

# Soundscape richness 

  # 0 - 10,000 Hz

ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             output = "percentage",
             minfreq = 0, 
             maxfreq = 10000)

## [1] 22.62931

  # 10,000 - 20,000 Hz

ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             output = "percentage",
             minfreq = 10000, 
             maxfreq = 20000)

## [1] 3.861351

ss_diversity(soundscape_obj = balbina[["Andre"]], 
             qvalue = 0, 
             output = "percentage",
             freqseq = TRUE, 
             nbins = 20)

##    soundscape_div    freq_interval
## 1       5.0347222      0 - 1033 Hz
## 2       8.6805556   1033 - 2067 Hz
## 3      36.4583333   2067 - 3273 Hz
## 4      21.9907407   3273 - 4306 Hz
## 5      32.5396825   4306 - 5512 Hz
## 6      31.7708333   5512 - 6546 Hz
## 7      12.9629630   6546 - 7579 Hz
## 8      16.7658730   7579 - 8785 Hz
## 9      36.4583333   8785 - 9819 Hz
## 10      6.1921296  9819 - 10852 Hz
## 11     20.1967593 11025 - 12058 Hz
## 12      6.6550926 12058 - 13092 Hz
## 13      2.3809524 13092 - 14298 Hz
## 14      0.6944444 14298 - 15331 Hz
## 15      0.4960317 15331 - 16537 Hz
## 16      0.0000000 16537 - 17571 Hz
## 17      0.0000000 17571 - 18604 Hz
## 18      0.0000000 18604 - 19810 Hz
## 19      0.0000000 19810 - 20844 Hz
## 20     50.0000000 20844 - 21877 Hz

2.3. The ss_diversity_plot function

The ss_diversity_plot function produces plots showing the variation in soundscape richness and diversity metrics by time-of-day and frequency subsets. Like the ss_diversity function, the ss_diversity_plot function takes three basic arguments: ‘soundscape_obj’, ‘qvalue’ and ‘output’ (see above).

Additionally, the ss_diversity_plot function can be used to create four types of plots by specifying the graphtype argument.

Let’s have a look at these plotting options one by one.

# Temporal patterns in the soundscape richness

ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "total", 
                  output = "percentage", 
                  smooth = FALSE)

## Warning: Removed 1 rows containing missing values (`geom_text()`).

# movavg = 6

plot_1 <- ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "total", 
                  output = "percentage", 
                  smooth = TRUE, 
                  movavg = 6, 
                  timeinterval = "4 hours") + 
  
  ggplot2::scale_y_continuous(limits = c(0, 50), 
                              expand = c(0,0))

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

# movavg = 12


plot_2 <- ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "total", 
                  output = "percentage", 
                  smooth = TRUE, 
                  movavg = 12, 
                  timeinterval = "4 hours")+ 
  
  ggplot2::scale_y_continuous(limits = c(0, 50), 
                              expand = c(0,0))

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

# movavg = 24


plot_3 <- ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "total", 
                  output = "percentage", 
                  smooth = TRUE, 
                  movavg = 24, 
                  timeinterval = "4 hours")+ 
  
  ggplot2::scale_y_continuous(limits = c(0, 50), 
                              expand = c(0,0))

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

# Combine plots 

library(patchwork)

plot_total <- plot_1 + plot_2 + plot_3

plot_total

## Warning: Removed 3 rows containing non-finite values (`stat_align()`).

## Warning: Removed 1 rows containing missing values (`geom_text()`).
## Removed 1 rows containing missing values (`geom_text()`).
## Removed 1 rows containing missing values (`geom_text()`).

  # Below 2,000 Hz

plot_1 <- ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "total", 
                  output = "percentage", 
                  smooth = TRUE, 
                  movavg = 24,  
                  maxfreq = 2000, 
                  timeinterval = "4 hours")+ 
  
  ggplot2::scale_y_continuous(limits = c(0, 55), 
                              expand = c(0,0))

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

  # Between 2,000 - 8,000 Hz

plot_2 <- ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "total", 
                  output = "percentage", 
                  smooth = TRUE, 
                  movavg = 24, 
                  minfreq = 2000, 
                  maxfreq = 8000, 
                  timeinterval = "4 hours")+ 
  
  ggplot2::scale_y_continuous(limits = c(0, 55), 
                              expand = c(0,0))

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

# Above 8,000 Hz

plot_3 <- ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "total", 
                  output = "percentage", 
                  smooth = TRUE, 
                  movavg =24, 
                  minfreq = 8000, 
                  timeinterval = "4 hours")+ 
  
  ggplot2::scale_y_continuous(limits = c(0, 55), 
                              expand = c(0,0))

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

# Combine the plots

plot_total <- plot_1 + plot_2 + plot_3

plot_total

## Warning: Removed 1 rows containing missing values (`geom_text()`).
## Removed 1 rows containing missing values (`geom_text()`).
## Removed 1 rows containing missing values (`geom_text()`).

ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "frequency", 
                  output = "percentage", 
                  smooth = TRUE, 
                  movavg = 24, 
                  maxfreq = 20000)

ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "frequency", 
                  output = "percentage", 
                  smooth = TRUE, 
                  movavg = 24, 
                  maxfreq = 20000, 
                  nbins = 5)

ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "normfreq", 
                  output = "percentage", 
                  smooth = TRUE, 
                  movavg = 24, 
                  maxfreq = 20000)

ss_diversity_plot(soundscape_obj = balbina[["Andre"]], 
                  qvalue = 0, 
                  graphtype = "linefreq", 
                  output = "percentage", 
                  smooth = TRUE, 
                  movavg = 24, 
                  maxfreq = 20000, 
                  nbins = 5)

3.1. The ss_compare function

The ss_compare function works in a very similar way to the ss_heatmap function and requires many of the same arguments. However, instead of portraying the presence and relative abundance of OSUs for a single soundscape, the ss_compare function is used to contrast OSUs between two soundscapes. To do so, the function creates a differential heatmap that visually illustrates the differences between two soundscapes.

Like the ss_heatmap function, ss_compare offers the option to subset the soundscape by specific time or frequency coordinates, portray the soundscape in either Cartesian or polar coordinates, and more.

Let’s have a look at the mandatory input argument.

The remaining arguments are the same as ss_heatmap, except for the annotate argument, which cannot be specified. This is because, when contrasting two different soundscapes, the sunrise and sunset times may differ.

Let’s have a look at how this works by comparing the soundscape of ‘Andre’ and ‘Mascote_A1’:

soundscapeR::ss_compare(soundscape_obj_A = balbina[["Andre"]],
                        soundscape_obj_B = balbina[["Mascote_A1"]], 
                        maxfreq = 20000)

Based on the differential heatmap produced by ss_compare, we can see that the OSU composition between Andre island and Mascote_A1 differs quite a bit. For instance, we can see that Mascote_A1 contains more sounds between 12,000-14,000 Hz at night time (redder values indicate OSUs were more common in soundscape_obj_B, or Mascote_A1). Additionally, Mascote_A1 clearly contains more sound in the daytime, showing three distinct bands at approx. 2,000 Hz, 6,000 Hz, and 8,000-10,000 Hz. Conversely, we can see that Andre island contains sound in places where little sound sound is present for Mascote_A1 (bluer values indicate OSUs were more common in soundscape_obj_A, or Andre island). For instance, the vertical peak around dawn is distinct in Andre, but not as clear is Mascote A1. Moreover, Andre contains some sounds between 3,000-4,000 Hz and 5,000-6,000 Hz at nighttime, which are not present for Mascote_A1.

Clearly, these differential heatmaps are a useful tool to quickly compare the OSU composition between two soundscapes.

3.2. The ss_pcoa function

If we want to go beyond heatmaps, and assess how similar or different multiple soundscapes actually are from one another, we can use the ss_pcoa function. The ss_pcoa function can be used to create principle coordinate plots, using a dimensionality reduction approach on the OSU incidence data of multiple soundscapes to plot the distances between these soundscapes in a two-dimensional space. To do this, the function calculates the Bray-Curtis dissimilarities between the soundscapes.

Let’s have a look at the mandatory arguments:

In it’s most basic form, the ss_pcoa function takes a single argument: a list of soundscape object to compare.

Let’s have a look at what this looks like:

ss_pcoa(soundscape_list =  balbina)

Great, this already reveals a lot! Based on this principle coordinate plot, we can see that the first PCo explains 52.9% of the variation between the soundscapes, and separates the plots on Mascote island from the plots on Andre island - clearly these is a difference in the soundscapes between these two islands. We can also see that a second PCo explains 22.4% of the variation and separates out plot Mascote_A2 from the other plots on Mascote island.

Now, let’s have a look at the other plotting parameters:

Let’s start of by taking a look at the grouping argument:

ss_pcoa(soundscape_list = balbina, 
        grouping = c("Andre", rep("Mascote", 4)))

The grouping argument is used to indicate which soundscapes in the balbina argument belong together, and colors them accordingly. This is already a little clearer than the previous plot.

Next, let’s take a look at the screeplot argument:

ss_pcoa(soundscape_list = balbina, 
        grouping = c("Andre", rep("Mascote", 4)), 
        screeplot = TRUE)

When we indicate screeplot = TRUE, the ss_pcoa function produces a second plot showing the cumulative percentage of the total variation that is explained as additional PCo’s are added. We can see that, at 4 PCo’s, all the variation between soundscapes is captured.

3.3. The ss_divpart function

In addition to quantifying the soundscape richness, diversity, and evenness components, the workflow proposed in this manuscript can be used to decompose the regional metacommunity diversity (γ-diversity) into its local diversities (α-diversity) and a community turnover component (β-diversity).

To do this, we use the following multiplicative framework provided by Hill numbers:

Here, N refers to the total number of sub-systems (soundscapes), j refers to each individual sub-system, and wj represents the relative weight given to each sub-system in the system. If all soundscapes are weighted equally, wj equals 1/N. The alpha diversity is the Hill number of the averaged basic sums of the soundscapes (Eqn. 2). The gamma diversity is computed by taking the average of the relative abundance of each OSU across the soundscapes in the system and calculating the Hill number of the pooled system (Eqn. 3). The beta diversity captures the degree of heterogeneity in the OSU composition across sites (Eqn. 4). It ranges from 1 to N and quantifies the relationship between the regional and local diversity, that is, how many times more diverse is the whole system in the effective number of OSUs compared to the sub-systems on average. The beta diversity can also be seen as the effective number of completely distinct soundscapes in the system.

Now that we know the theory, let’s take a look at how we can do this in practice. In its most basic form, the ss_divpart function takes the following arguments:

Now, we can use our island system to see how the overall system soundscape diversity (gamma) differs from the local soundscape diversity (alpha), and what the soundscape turnover (beta diversity) looks like.

Let’s try the code to decompose the soundscape richness into its various components:

ss_divpart(soundscape_list = balbina,
           qvalue = 0)

## # A tibble: 1 × 7
##   levels     q alpha_l1 gamma    N1    N2 beta_l1
##    <dbl> <dbl>    <dbl> <dbl> <dbl> <dbl>   <dbl>
## 1      2     0     31.4  57.2     5     1    1.82

We can see that the system’s alpha soundscape richness is 31.4%, the average richness of the local soundscape. However, overall, the consolidated soundscape richness of the system (gamma richness) is 57.2% - almost twice as high. A lower alpha richness at each location, but higher gamma richness when taking all locations together, suggests there might be some OSU turnover going on between our soundscapes - they are not identical. Indeed, this is confirmed by the soundscape turnover (beta), which is 1.82. In theory, this value can range from 1 (all soundscapes completely identical) to 5 (all soundscapes completely distinct). Our value here suggests that there is low to moderate turnover in the OSU composition between the soundscapes in our system.

But not all soundscapes in our study system are equal - some are more similar to others. For example, we have four soundscapes from Mascote island, but only one soundscape from Andre island. We can incorporate this ‘grouping’ or ‘hierarchical structure’ into our our diversity partitioning calculations.

Let’s have a look at the next argument of the ss_divpart function:

As we can see, we can supply the ss_divpart function with a table indicating the hierarchy between samples.

Let’s give this a try:

hierarchy_table <- as.matrix(
  data.frame(Sample = c("Andre", "Mascote_A1", "Mascote_A2", "Mascote_B1", "Mascote_B2"), 
             grouping_1 = c("Andre", "Mascote", "Mascote", "Mascote", "Mascote")))

ss_divpart(soundscape_list = balbina, 
           qvalue = 0, 
           hier_table = hierarchy_table)

## # A tibble: 1 × 10
##   levels     q alpha_l1 alpha_l2 gamma    N1    N2    N3 beta_l1 beta_l2
##    <dbl> <dbl>    <dbl>    <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>   <dbl>
## 1      3     0     31.4     35.5  57.2     5     2     1    1.13    1.61

As you can see, by defining the hierarchical relationships between the plots on our islands, the output has changed a little. We now obtain two alpha richness values, one calculated at the sample (or plot) level, and one calculated at the level of the grouping we specified (islands).

We also get two soundscape turnover values (beta), one for each grouping level:

1. β1 = γ-islands (35.5%) / α-plots (31.4%) = 1.13
2. β2 = γ-total (57.2%) / α-islands (35.5%) = 1.61

By using the hier_table argument, we can further decompose the soundscape metrics based on the hierarchical relationship between samples (plots).

Yet, because beta diversity ranges between 1-N, it is not independent of the number of soundscapes in the system, and can thus not be used as a measure of similarity directly. Instead, to compare the relative compositional difference between soundscapes across multiple systems with a different number of soundscapes, some simple transformations can be performed on the beta diversity to remove the dependence on the number of soundscapes. To do this, we will use the next function.

Note: We can also partition the soundscape metrics using subsets of the soundscapes using the minfreq, maxfreq, mintime and maxtime arguments we’ve seen before.

3.4. The ss_pairdis function

The framework of Hill numbers also allows us to define several measures of similarity between soundscapes in the wider system. Let’s take a look at how these are calculated:

Equations (5) and (6) are measures of overlap between soundscapes. The local or Sørensen-type overlap (CqN) quantifies the effective average proportion of a soundscape’s OSUs which are shared across all soundscapes. It captures the overlap between soundscapes from the sub-system’s perspective. For N soundscapes each having S equally common OSUs and sharing A OSUs between them, this function reduces to CqN=A/S. The regional or Jaccard-type overlap(UqN) quantifies the effective proportion of shared OSUs in a pooled assemblage of soundscapes, and thus captures the overlap between soundscapes from a regional perspective. Assume N soundscapes in a region with S unique and equally abundant OSUs. Here, R OSUs are shared between all soundscapes and the remaining OSUs (S-R) are distributed evenly among N soundscapes. In this scenario, Eqn. 8 reduces to UqN = R/S.

Equations (7) and (8) are measures of turnover in OSUs between soundscapes. The local or Sørensen-type turnover complement (VqN) quantifies the normalised OSU turnover rate with respect to the average soundscape. It measures the proportion of a typical soundscape which changes as one goes from one soundscape to the next. The regional or Jaccard-type turnover complement (SqN) quantifies the proportion of the regional soundscape diversity contained in the average assemblage and is a measure of regional homogeneity.

For all of the aforementioned similarity indices, the ss_pairdis function transforms the values into metrics of dissimilarity by taking their one-complement (e.g. 1 – XqN). Unlike the beta soundscape turnover, these dissimilarity indices range from 0-1, where 0 means the soundscapes are completely identical, and 1 indicates the soundscapes are completely unique.

If you want to learn more about the computation of similarity metrics using the framework of Hill numbers, have a look at the hilldiv GitHub tutorial here and this paper.

To use the ss_pairdis function, we need to supply the same parameters as the ss_divpart function: balbina, qvalue and hier_table (optional).

Let’s have a look:

ss_pairdis(soundscape_list = balbina, 
           qvalue = 0)

## $Level_1_beta
##               Andre Mascote_A1 Mascote_A2 Mascote_B1 Mascote_B2
## Andre            NA         NA         NA         NA         NA
## Mascote_A1 1.608346         NA         NA         NA         NA
## Mascote_A2 1.606310   1.409356         NA         NA         NA
## Mascote_B1 1.614435   1.312740   1.374302         NA         NA
## Mascote_B2 1.623355   1.264852   1.367437   1.274029         NA
## 
## $L1_CqN
##                Andre Mascote_A1 Mascote_A2 Mascote_B1 Mascote_B2
## Andre             NA         NA         NA         NA         NA
## Mascote_A1 0.6083460         NA         NA         NA         NA
## Mascote_A2 0.6063103  0.4093560         NA         NA         NA
## Mascote_B1 0.6144346  0.3127399  0.3743021         NA         NA
## Mascote_B2 0.6233550  0.2648524  0.3674368  0.2740286         NA
## 
## $L1_UqN
##                Andre Mascote_A1 Mascote_A2 Mascote_B1 Mascote_B2
## Andre             NA         NA         NA         NA         NA
## Mascote_A1 0.7564865         NA         NA         NA         NA
## Mascote_A2 0.7549106  0.5809122         NA         NA         NA
## Mascote_B1 0.7611762  0.4764689  0.5447159         NA         NA
## Mascote_B2 0.7679836  0.4187878  0.5374095  0.4301766         NA
## 
## $L1_VqN
##                Andre Mascote_A1 Mascote_A2 Mascote_B1 Mascote_B2
## Andre             NA         NA         NA         NA         NA
## Mascote_A1 0.6083460         NA         NA         NA         NA
## Mascote_A2 0.6063103  0.4093560         NA         NA         NA
## Mascote_B1 0.6144346  0.3127399  0.3743021         NA         NA
## Mascote_B2 0.6233550  0.2648524  0.3674368  0.2740286         NA
## 
## $L1_SqN
##                Andre Mascote_A1 Mascote_A2 Mascote_B1 Mascote_B2
## Andre             NA         NA         NA         NA         NA
## Mascote_A1 0.7564865         NA         NA         NA         NA
## Mascote_A2 0.7549106  0.5809122         NA         NA         NA
## Mascote_B1 0.7611762  0.4764689  0.5447159         NA         NA
## Mascote_B2 0.7679836  0.4187878  0.5374095  0.4301766         NA

Indeed, based on this output, we can see that the pairwise dissimilarity between the plots in our case study is greatest between the soundscape of Andre island and all other plots. Furthermore, within Mascote island, we can see that plot Mascote_A2 is most dissimilar from the other plots on the island - a pattern which we previously observed used the ss_pcoa function.