README

# African elephants address one another with individually specific calls

We investigated the hypothesis that elephants address individual members
of their family group with name-like calls. We recorded contact,
greeting, and caregiving rumbles from wild African elephants in Samburu
& Buffalo Springs National Reserves, northern Kenya and Amboseli
National Park, southern Kenya, noting when possible the identity of the
caller and the identity of the receiver. We measured a suite of acoustic
features on each call and found that calls were specific to individual
receivers and receiver identity could be predicted from call structure
at better than chance levels. We also played back calls to the elephants
and found that elephants responded more quickly and vocalized more in
response to calls that were originally addressed to them compared to
calls from the same caller that were originally addressed to someone
else. These results provide the first evidence that elephants address
one another with individually specific labels.

Description of the data and file structure

“20231225_observational_data_for_upload_revision1_n469.csv” This
spreadsheet contains the associated observational behavioral data for
each call used in the analysis. Each row represents a single call, and
columns indicate the identity of the caller, the identity of the
receiver, the behavioral context, the certainty with which each of these
values were known, and other relevant information, including the file
name of the sound file in which the call occurs, the start and end times
of the call in that sound file, the file name of the Raven Pro selection
table used to annotate the sound file, and the number of the selection
box in that table corresponding to the call in question. The sound files
and Raven selection tables themselves are not included in this upload.
The full set of variables in this spreadsheet are defined in the
spreadsheet named
“20240123_variable_definitions_observational_data_for_upload_revision1.csv”.
The derived acoustic features used in the analysis are not included in
this spreadsheet. Rather, those measurements can be calculated using the
code provided
(“20240104_vocal_labeling_code_for_upload_revision1.R”) and the
raw acoustic measurements in the .RDS file named
“20231225_acoustic_contours_for_upload_revision1_n469.rds”. Values
of “NA” in this spreadsheet indicate that the information in question
was unknown for a given observation.

“20240123_variable_definitions_observational_data_for_upload_revision1.csv”
This spreadsheet defines each of the variables in
“20231225_observational_data_for_upload_revision1_n469.csv”

“20231225_acoustic_contours_for_upload_revision1_n469.rds” This
RDS file is an R list object which contains the raw acoustic contours
measured on each call. The code provided
(“20240104_vocal_labeling_code_for_upload_revision1.R”) uses these
raw acoustic contours to calculate a suite of derived acoustic features
which are then used for subsequent statistical analyses. The list is
nested with two layers. Each outer slot represents a single call, and is
named with a unique identifier for that call, which is comprised of the
date (YYYYMMDD) and time (HHMM) of the original recording as well as the
Raven selection box number of the call in question. The inner slots
represent individual acoustic contours, each of which is a vector of
values. The acoustic contours are as follows: aenv = Hilbert amplitude
envelope aenvtime = time points (sec) for aenv (not used in analysis)
pmsband(1-26): energy in each of the 26 bands of a mel spectrogram,
calculated from 0-500 Hz. Each energy value within a given band
corresponds to one of the time points in aenvtime pmsD(1-26):
differences between successive values of the corresponding pmsband
pmsDD(1-26): differences between successive values of the corresponding
pmsD cc(1-12): the first 12 cepstral coefficient contours, calculated
from a mel spectrogram with 40 bands from 0-500 Hz. Each cepstral
coefficient value in a given contour corresponds to a particular time
point in aenvtime ccD(1-12): differences between successive values of
the corresponding cc contour ccDD(1-12): differences between successive
values of the corresponding ccD contour

The features referred to in the paper and the R code as “spectral
features” are derived from the following contours: aenv, pmsband(1-26),
pmsD(1-26), and pmsDD(1-26)

The features referred to in the paper and the R code as “cepstral
features” are derived from the following contours: aenv, cc(1-12),
ccD(1-12), and ccDD(1-12)

“20231225_dfall_obs_data_and_spectral_features_revision1_n469.csv”
This file is a spreadsheet containing the observational data (caller ID,
receiver ID, etc.) together with the derived SPECTRAL acoustic features.
This spreadsheet is generated in section 4 of the code using the data
files “20231225_observational_data_for_upload_revision1_n469.csv”
and “20231225_acoustic_contours_for_upload_revision1_n469.rds”.
However, instead of recreating it de-novo the reader may skip sections
2-4 of the code and instead load this file in section 5 of the code. The
variables in this spreadsheet are all described in the file
“20240123_variable_definitions_observational_data_for_upload_revision1.csv”.
Values of “NA” in this spreadsheet indicate that the information in
question was unknown for a given observation.

“20240106_dfall_obs_data_and_cepstral_features_revision1_n469.csv”
This file is a spreadsheet containing the observational data (caller ID,
receiver ID, etc.) together with the derived CEPSTRAL acoustic features.
This spreadsheet is generated in section 4 of the code using the data
files “20231225_observational_data_for_upload_revision1_n469.csv”
and “20231225_acoustic_contours_for_upload_revision1_n469.rds”.
However, instead of recreating it de-novo the reader may skip sections
2-4 of the code and instead load this file in section 5 of the code. The
variables in this spreadsheet are all described in the file
“20240123_variable_definitions_observational_data_for_upload_revision1.csv”.
Values of “NA” in this spreadsheet indicate that the information in
question was unknown for a given observation.

“20230409_playback_data_for_upload.csv” This spreadsheet contains
the results of the playback experiment. In this experiment, wild
elephants were exposed to playbacks of a call originally addressed to
them (test) and a call from the same caller originally addressed to
someone else (control). We filmed the elephants’ responses and recorded
from the video and audio recording the latency to approach the speaker
within 10 min of the playback, the latency to vocalize within 10 min of
the playback, the latency to vigilance within 10 min of the playback,
the number of vocalizations produced within 10 min of the playback, and
the amount of time spent being vigilant in 1 min following the playback
minus the amount of time spent being vigilant in 1 min before the
playback. The variables in this spreadsheet are all listed and defined
in the spreadsheet named
“20230409_variable_definitions_playback_data_for_upload.csv”.
Values of “NA” in this spreadsheet indicate that the information in
question could not be assessed for a given observation.

“20230409_variable_definitions_playback_data_for_upload.csv” This
spreadsheet defines each of the variables in
“20230409_playback_data_for_upload.csv”.

“20240104_vocal_labeling_code_for_upload_revision1.R” This R
script contains the code for all analyses in this study. Each section of
the code is labeled to describe what it does. This script replaces an
earlier version of the script called
“20230409_vocal_labeling_code_for_upload.R”.

Code/Software

“20240104_vocal_labeling_code_for_upload_revision1.R” This R
script contains the code for all the analyses in this study. Each
section of the code is labeled to describe what it does. The packages
required for the code to run are all listed in the first section. The
sections in this script are as follows:

1. Load the required packages for the whole script.

2. Code used to segment the calls from the original raw sound files and
   save them as clips for further analysis. The original sound files
   themselves are not included in this upload so this section cannot be
   run, but it is included to illustrate how we segmented the
   recordings

3. Code used to measure acoustic contours on the segmented sound clips.
   The sound clips are not included in this upload so this section can
   also not be run, but it is included to illustrate how we performed
   the acoustic contour measurements

4. Code to derive acoustic features from the acoustic contours and
   append them to the observational data. This section requires two
   files that are included in this upload:
   “20231225_observational_data_for_upload_revision1_n469.csv”
   and
   “20231225_acoustic_contours_for_upload_revision1_n469.rds”. We
   ran all of our analyses twice: once with the spectral features and
   once with the cepstral features (see above). Starting in this
   section of the script, where it is relevant we have included lines
   of code to allow the user to use either the spectral features or the
   cepstral features, with comments to indicate which is which. The
   code pertaining to cepstral features is currently commented out
   while the code pertaining to spectral features is not. To run the
   analysis with the cepstral features, the user should remove the #
   in front of the relevant lines of code for using the cepstral
   features, and comment out the corresponding lines for the spectral
   features.

5. Set the dataset that you wish to use for subsequent analyses
   (spectral features or cepstral features)

6. Code to calculate which calls belong to the same bout and the ranked
   position of each call within its bout. Bout defined as successive
   calls from the same caller in the same original raw recording with
   no more than 30 sec between successive calls.

7. Analysis to compare the mean similarity of different categories of
   call pairs. Similarity is measured using random forest proximity
   scores calculated with the algorithm devised by Rhodes, Cutler, and
   Moon 2022. We run two comparisons: pairs of calls with the same
   caller and same receiver vs. pairs of calls with the same caller and
   different receivers, and pairs of calls with different callers and
   the same receiver vs. pairs of calls with different callers and
   different receivers.

8. Code to add another variable to the dataset which indicates whether
   a given call was “imitative” of the receiver’s calls or “not
   imitative” of the receiver’s calls. This is determined by first
   calculating the mean proximity score between the focal call and each
   of the calls made by the receiver (meanElic1), and then calculating
   the mean proximity score between each of the calls made by the focal
   caller to someone other than the focal receiver and each of the
   calls made by the focal receiver (meanOther). The idea is to
   determine whether the focal call is more similar to the receiver’s
   calls than is typical for the focal caller. If meanElic1 >
   meanOther, the focal call is deemed to be more similar to the
   receiver’s calls than typical for that caller, and is therefore
   considered “imitative”. Otherwise, it is considered “non-imitative”.

9. Analysis to determine if calls are more similar to the receiver’s
   calls than to calls made by other individuals. We compare pairs of
   calls in which the caller of one call was the receiver of the other
   call to pairs of calls in which this was not the case, while
   controlling for the identity of the two callers.

10. Random forest models using the ranger package to predict receiver ID
    as a function of the derived acoustic features. The models are run
    with 7-fold cross validation, and each 7-fold model is run 2000x and
    the classification accuracy is averaged across the 2000 runs. The
    first model stratifies the cross-validation folds by caller and
    receiver ID to ensure that the same number of calls from each
    caller-receiver dyad are included in each fold to the extent
    possible. This tests the hypothesis that receiver ID can be
    predicted from call structure within a given caller. The second
    model partitions the cross-validation folds such that all calls from
    the same caller-receiver dyad are always included in the same fold.
    This tests the hypothesis that receiver ID can be predicted from
    call structure across callers.

11. Run the same random forest models from section 10, but this time
    randomly permute the acoustic feature variables so they are
    uninformative. Run each model 10,000x instead of 2000x, and save the
    classification accuracy from each run in a data table. This data
    table will be used in the following section as a null distribution
    to test the hypothesis that the corresponding model from step 10 can
    predict receiver ID better than chance.

12. Calculate a p-value for the classification accuracy of each random
    forest model by comparing it to the null distribution generated in
    section 11

13. Run logistic regression model to assess the factors influencing the
    probability of a call’s receiver ID being correctly/incorrectly
    predicted in >95% of runs of the first random forest model in
    section 10. Calls that were not correct in >95% of runs or
    incorrect in >95% of runs are excluded.

14. Calculate variable importance scores from two random forest models:
    one to predict receiver ID from acoustic features and one to predict
    caller ID from acoustic features. We use a conditional inference
    random forest in the package party rather than using ranger to
    obtain unbiased estimates of variable importance. This code
    generates a table of the variable importance scores for the derived
    acoustic features (Supplementary Table S2) and also calculates the
    relative importance of the original acoustic contours and puts them
    in a separate table (Supplementary Table S3. We only calculate
    variable importance scores for the spectral features because they
    are easier to interpret than the cepstral features.

15. Code to generate histograms and boxplots illustrating the results of
    the analyses in sections 7, 9, and 10 (Figs 1-3 in the main
    manuscript.

16. Code to generate figure illustrating how spectral acoustic contours
    were measured (Extended Data Fig. 1)

17. Code to generate multi-panel 3-D scatterplots showing similarity
    among calls addressed to same vs. different receivers within a given
    caller (Extended Data Fig. 2), and code to generate a single-panel
    scatterplot showing similarity among calls addressed to the same
    receiver across callers (Extended Data Fig. 3)

18. Analysis of playback results. Models to test the hypotheses that
    latency to approach, latency to vocalize, latency to vigilance,
    number of vocalizations, and time spent vigilant after playback
    minus before playback differ between test and control trials

19. Code to generate plots illustrating the results of the analyses in
    section 19 (Fig. 4 in main manuscript)

20. Code to calculate sample sizes for all the analyses in this script

21. Code to generate table of number of calls for each caller-receiver
    pair (Supplementary Table S1)