Analyzing whale calling through Hawkes process modeling

The single-channel dataset was collected and processed by the National Oceanic and Atmospheric Administration (NOAA) (DCLDE 2013,, 2023). This comprised recordings from a week-long deployment of a bottom-mounted hydrophone system, referred to as a marine autonomous recording unit or MARU (Calupca et al.,, 2000); the MARU was deployed to the west of Race Point on the northern edge of Cape Cod Bay (Figure 1). MARUs were programmed to record sound continuously in the 10-585 Hz low-frequency band (Parks et al.,, 2009) including biotic, abiotic, and anthropogenic sources. Post-processing was done with a software routine designed to automatically detect and classify low-frequency calls of baleen whales, including NARW (Baumgartner and Mussoline,, 2011). Following automated detection, a trained analyst manually reviewed and verified possible NARW up-calls, and the processed dataset was made available to the acoustic processing community as a benchmark dataset to be used in subsequent algorithm development. There were 5,302 unique and verified up-calls detected across the seven days (2009-03-28 to 2009-04-03), and we extracted ambient noise levels (dB re 1 µPa) across the same frequency band (60-400 Hz) as the NARW up-call using the PAMGuard software program (https://www.pamguard.org). Ambient sound was in dB re 1 µ-Pa, i.e., decibel units of Sound Pressure Level, root mean square, with reference to 1 µ-Pascal—a standard metric in the bioacoustics community (Au and Hastings,, 2008). This noise metric was used as a covariate in the modeling.

The multi-channel dataset consists of recordings from the distributed network of 10 MARUs located within Cape Cod Bay, MA. The distances between MARUs are found in the supplement Figure S1. This bay is seasonally inhabited by NARWs (Mayo et al.,, 2018), with upwards of 70% of the entire species sometimes being observed during aerial surveys and vessel-based surveys (C. Hudak, pers comm). Starting in the early 2000s, researchers deployed passive acoustic arrays of different size and configuration in CCB coincident with the season when NARWs occur in this habitat (Clark et al.,, 2010; Urazghildiiev,, 2014). These arrays use the same type of MARU used in the single-channel dataset described earlier, however the data processing is slightly different. For our multi-channel analysis, we used acoustic data collected from 10 MARUs in a distributed network across 9 days (2010-04-02 to 2010-04-10). The data are archival, so the buoys must be recovered to obtain the sound files (Calupca et al.,, 2000).

We used the software PAMGuard to process the raw sound files, detect and classify NARW up-calls, extract ambient noise metrics, and localize individual up-calling NARWs (Gillespie,, 2004; Gillespie et al.,, 2020). Localization algorithms use the time-difference-of-arrival of three or more matched up-calls (Urazghildiiev,, 2014). To determine the times of unique up-calls, i.e., those arriving from a single whale, we found all up-calls associated with each unique localization; from this associated set we extracted the first arrival time, as well as its associated MARU. Across the 9 days, a total of 2,750 unique up-calls were extracted from the data. The arrival times on the MARUs are the data that comprise the analysis. Figure 2 (a) shows the ambient noise for two illustrative MARUs (MARU 5 and MARU 10 in Figure 1)) highlighting variation across MARUs as well as within MARUs across days. Figure 2 (b) depicts the event sequence of unique up-calls received on each of the MARUs within CCB for the 8-day period. The up-call sequences show diel patterns as well as possible excitement in up-call rates, motivating the inclusion of diel harmonics in the intensity of the process as well as Hawkes process modeling.

Finally, we investigate possible spatial dependence in the up-calling behavior captured by the array of MARUs. We might expect nearby hydrophones to capture similar up-calling behavior as a response to noise and temporal harmonics. In addition, we might expect the number of counter-calls received at MARU $k$ that were excited by up-calls at MARU $l$ to decrease with the distance between MARUs $l$ and $k$. For each MARU, we compute the kernel intensity estimate for the total event sequence of up-calls. Then, we compute pairwise correlations between these estimates. In Figure 3, we show these pairwise correlations spatially between MARU 10 and each of the other MARUs. Positive correlation is detected between MARU 10 and each of the nearby MARUs, suggesting evidence of spatial dependence across the array.

We propose four different multivariate Hawkes process models for studying the behavior of NARW up-calls. Let $\mathcal{T}=\{(t_{i},m_{i}),i=1,\dots,n\}$ be the sequence of paired observations where $t_{i}\in(0,T]$ is an up-call time with $t_{i}<t_{i+1}$ and $m_{i}\in\{1,\dots,K\}$ indicates the MARU which received up-call $i$ at time $t_{i}$. Let $\mathcal{H}_{t}$ be the history of all up-calls up to time $t$, i.e., $\mathcal{H}_{t}=\{(t_{i},m_{i}):t_{i}<t\}$. Such a point process in time can be characterized by a conditional intensity defined as

$\displaystyle\lambda(t\mid\mathcal{H}_{t})=\lim_{\Delta t\to 0}\frac{\textup{E}\left\{N\left((t,t+\Delta t]\right)\mid\mathcal{H}_{t}\right\}}{\Delta t},$

where $N(A)$ is the number of up-calls recorded over $A\subseteq(0,T]$.

We label our most general model as NHPP+GP+CC (nonhomogeneous Poisson process + Gaussian process + counter-call process). For this model we define the conditional intensity of up-calls at MARU $k$ at time $t$ as

$\displaystyle\lambda_{k}(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})$

$\displaystyle=\mu_{k}(t)+\sum_{i:t_{i}<t,t_{i}\in\mathcal{T}}\alpha_{m_{i}}e^{-\eta(t-t_{i})}e^{-\phi d_{m_{i},k}},$

(1)

where $\theta$ is the collection of model parameters. The conditional intensity of up-calls across $K$ MARUs at time $t$ is $\lambda(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})=\sum_{k=1}^{K}\lambda_{k}(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})$.

The first term $\mu_{k}(t)$ of (1) explains the rate of contact calls at MARU $k$ at time $t$. We use minutes and kilometers (km) for the units of time and space. To capture the impact of noise on contact calls as well as the potential daily pattern of contact calls, we assume

	$\displaystyle\mu_{k}(t)=\exp\bigg{\{}$	$\displaystyle\beta_{0k}+\beta_{1k}\text{Noise}_{k}(t)+\beta_{2k}\sin\left(\frac{2\pi t}{8\times 60}\right)+\beta_{3k}\cos\left(\frac{2\pi t}{8\times 60}\right)+\beta_{4k}\sin\left(\frac{2\pi t}{12\times 60}\right)$
		$\displaystyle+\beta_{5k}\cos\left(\frac{2\pi t}{12\times 60}\right)+\beta_{6k}\sin\left(\frac{2\pi t}{24\times 60}\right)+\beta_{7k}\cos\left(\frac{2\pi t}{24\times 60}\right)+\delta_{k}w(t)\bigg{\}},$

where Noise${}_{k}(t)$ is the ambient noise level at MARU $k$ at time $t$, and the three sine and cosine pairs are the 8-hour, 12-hour, and 24-hour harmonics, respectively. The GP, $w(\cdot)$, is mean 0 and variance 1, and is used to account for variability in the contact call process not captured by the other covariates. We assume an exponential correlation function with an effective range of 3 hours to capture local temporal dependence. The $\delta_{k}$’s are MARU-specific coefficients for the GP. Let $\mbox{\boldmath$\beta$}_{j}=\left(\beta_{j1},\dots,\beta_{jK}\right)^{\top}$ denote the length $K$ variable-specific vector of coefficients, where $j=0$ denotes the intercept. Let $\mbox{\boldmath$\delta$}=\left(\delta_{1},\dots,\delta_{K}\right)^{\top}$ denote a vector of the GP coefficients. To account for spatial dependence, we assume that $\beta_{jk}$’s are dependent across $k$. Let $d_{k,\ell}$ denote the distance between MARUs $k$ and $\ell$. To improve Markov chain mixing, we incorporate customary hierarchical centering (Gelfand et al.,, 1996). For $j=0,1,\dots,p$, we specify $\mbox{\boldmath$\beta$}_{j}\sim\text{MVN}(\tilde{\beta}_{j}\textbf{1},\tau_{j}V)$ where $V_{k,\ell}$ is $e^{-3d_{k,\ell}/\max_{k^{\prime},\ell^{\prime}}\{d_{k^{\prime},\ell^{\prime}}\}}$ so that its effective range entirely covers the array of MARUs. The GP coefficients $\delta_{k}$’s are also assumed to be spatially dependent and given by $\log(\mbox{\boldmath$\delta$})\sim\text{MVN}(\tilde{\delta}\textbf{1},\tau_{\delta}V)$.

The second term of (1) captures the conditional intensity component for counter-calls at MARU $k$ at time $t$. It is natural to assume that the intensity of counter-calls caused by an up-call at $t_{i}$ is higher in time near $t_{i}$ and higher in space near $m_{i}$ and decreases as we move forward in time and away in space. Some motivation for this specification is provided by Figure 3. We assume the $\alpha_{k}$’s are independent with a Gamma prior.

In this manuscript we consider the following four models.

1. (i)

   NHPP: $\lambda_{k}(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})=e^{\beta_{0k}+\textbf{x}_{k}(t)^{\top}\boldsymbol{\beta}_{k}}$

2. (ii)

   NHPP+GP: $\lambda_{k}(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})=e^{\beta_{0k}+\textbf{x}_{k}(t)^{\top}\boldsymbol{\beta}_{k}+\delta_{k}w(t)}$

3. (iii)

   NHPP+CC: $\lambda_{k}(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})=e^{\beta_{0k}+\textbf{x}_{k}(t)^{\top}\boldsymbol{\beta}_{k}}+\sum_{i:t_{i}<t}\alpha_{m_{i}}e^{-\eta(t-t_{i})}e^{-\phi d_{m_{i},k}}$

4. (iv)

   NHPP+GP+CC: $\lambda_{k}(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})=e^{\beta_{0k}+\textbf{x}_{k}(t)^{\top}\boldsymbol{\beta}_{k}+\delta_{k}w(t)}+\sum_{i:t_{i}<t}\alpha_{m_{i}}e^{-\eta(t-t_{i})}e^{-\phi d_{m_{i},k}}$

Models (i)-(iii) are submodels of Model (iv). Models (i) and (ii) assume all up-calls are contact calls and that there are no counter-calls. Models (i) and (iii) are suitable when the covariate effects are sufficient to explain the rate of contact calls. The likelihood function (White and Gelfand,, 2021) associated with Model (iv) is given by

$\displaystyle L(\mbox{\boldmath$\theta$}\mid\mathcal{T})=\exp\left\{-\sum_{k=1}^{K}\int_{0}^{T}\lambda_{k}(t\mid\mathcal{H}_{t})dt\right\}\prod_{i=1}^{n}\lambda_{m_{i}}(t_{i}\mid\mathcal{H}_{t_{i}}).$

(2)

Model fitting details are found in the supplementary material.

The compensator for a point process provides the expected number of events under the intensity in the time window $(0,t]$. The compensator for a process on $(0,T]$ with intensity $\lambda(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})$ is defined as $\Lambda(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})=\int_{0}^{t}\lambda(u\mid\mathcal{H}_{u};\mbox{\boldmath$\theta$})du$. For any $t$, the compensator is a function of the data, i.e., it is a random variable. Here we focus on inference using the compensator. We will employ the compensator for assessing model adequacy in Section 4.1. If we evaluate the compensator at $t=T$, i.e., $\Lambda(T\mid\mathcal{H}_{T};\mbox{\boldmath$\theta$})$, we obtain the expected total number of events given the data. If we separate the integral of the conditional intensity, we obtain the expected number of contact calls and counter-calls received, respectively, in the time window, $(0,T]$, given the data. Using an earthquake analogy, this separates the expected number of earthquakes and aftershocks.

First, suppose that there is a single MARU. If we integrate the counter-call component, i.e., $\sum_{i:t_{i}\in\mathcal{T}}\int_{t_{i}}^{T}\alpha e^{-\eta(t-t_{i})}dt$, we obtain $\alpha/\eta\sum_{i:t_{i}\in\mathcal{T}}[1-e^{-\eta(T-t_{i})}]$. We can interpret the $i$th term in the sum as the contribution of the event at time $t_{i}$ to the expectation. From the Bayesian perspective, given the data, using posterior samples of the model parameters, we can obtain the posterior distributions of these functions, i.e., the joint distribution of the expected number of contact calls and the expected number of counter-calls, their marginal distributions, and the distribution of the sum.

Consider an array of $K$ MARUs. The intensity of up-calls received by each MARU $k$ is specified by $\lambda_{k}(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})=\mu_{k}(t)+\sum_{i:t_{i}<t}\alpha_{m_{i}}e^{-\eta(t-t_{i})}e^{-\phi d_{m_{i}k}}$, which splits the conditional intensity into an intensity for contact calls received by MARU $k$ and an intensity for counter-calls received by MARU $k$. An up-call received at any MARU at any time prior to $t$ provides a contribution to the counter-call intensity for MARU $k$ at time $t$. This is the interpretation of $\sum_{i:t_{i}<t}\alpha_{m_{i}}e^{-\eta(t-t_{i})}e^{-\phi d_{m_{i}k}}$. The expected number of contact calls received by MARU $k$ over $(0,T]$ is $\int_{0}^{T}\mu_{k}(t)dt$. The expected number of counter-calls received at MARU $k$ over $(0,T]$ is $\int_{0}^{T}$ $\sum_{t_{i}<t}$ $\alpha_{m_{i}}$ $e^{-\eta(t-t_{i})}$ $e^{-\phi d_{m_{i},k}}dt=\sum_{i:t_{i}\in\mathcal{T}}\int_{t_{i}}^{T}\alpha_{m_{i}}e^{-\eta(t-t_{i})}e^{-\phi d_{m_{i},k}}dt=\sum_{i:t_{i}\in\mathcal{T}}(\alpha_{m_{i}}/\eta)[1-e^{-\eta(T-t_{i})}]e^{-\phi d_{m_{i},k}}$. With multiple MARUs, if we sum the contact call component over $k$ we obtain the overall expected number of contact calls. If we sum the counter-call component over $k$ we obtain the overall expected number of counter-calls. As above, we can obtain the joint posterior distribution of the expected number of contact calls and the expected number of counter-calls, their marginal distributions, and the posterior distribution of total expected number of up-calls. We can also obtain the expected number of counter-calls at MARU $k$ that were excited by up-calls at a specific MARU $\ell$ in the time window $(0,T]$. This is given by $\sum_{i:t_{i}\in\mathcal{T}}\mathbb{I}(m_{i}=\ell)(\alpha_{m_{i}}/\eta)[1-e^{-\eta(T-t_{i})}]e^{-\phi d_{m_{i},k}}$ where $\mathbb{I}(\cdot)$ is the indicator function.

We define the compensator associated with our model as $\Lambda(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})=\int_{0}^{t}\lambda(u\mid\mathcal{H}_{u};\mbox{\boldmath$\theta$})du$ as in Section 3.2. By the random time change theorem (RTCT) (Brown and Nair,, 1988), if $\theta$ is “true” and $t^{*}_{i}=\Lambda(t_{i}\mid\mathcal{H}_{t_{i}};\mbox{\boldmath$\theta$})$ for $i=1,2,...,n$, then the set $\{t^{*}_{i}\}$ are a realization from a homogeneous Poisson process with rate 1 (HPP(1)) and the associated $d^{*}_{i}=t^{*}_{i}-t^{*}_{i-1}$ with $t^{*}_{0}=0$ are an i.i.d. sample from an exponential distribution with rate 1 (Exp(1)). For example, consider data $\mathcal{T}$ generated from an HPP$(\lambda)$. The associated waiting times $d_{i}=t_{i}-t_{i-1}$ with $t_{0}=0$ are independent and have an Exp$(\lambda)$ distribution. The compensator associated with HPP($\lambda$) is $t^{*}_{i}=\int_{0}^{t_{i}}\lambda du=\lambda t_{i}$. Then the associated $d^{*}_{i}=t^{*}_{i}-t^{*}_{i-1}=\lambda(t_{i}-t_{i+1})=\lambda d_{i}$ has an Exp(1) distribution.

The compensator $\Lambda(t\mid\mathcal{H}_{t};\mbox{\boldmath$\theta$})$, as a function of the data, is a random function. In the frequentist setting, there is a “true” $\theta$ and only $\mathcal{T}$ is random. We have to estimate $\theta$ and then use the estimated compensator $\Lambda(t\mid\mathcal{H}_{t};\hat{\mbox{\boldmath$\theta$}})$. This yields $\{d^{*}_{i}(\mathcal{T};\hat{\mbox{\boldmath$\theta$}})\}$. The RTCT does not formally apply to $\{d^{*}_{i}(\mathcal{T};\hat{\mbox{\boldmath$\theta$}})\}$ but the argument is that, if $\hat{\mbox{\boldmath$\theta$}}$ is close to the true $\theta$, the result holds approximately.

We consider the RTCT in a fully Bayesian modeling setting. From a frequentist perspective the compensator is a random variable as a function of the data. From the Bayesian perspective it is a random variable as a function of the parameters given the data and we examine it through posterior inference. Consider the simple HPP($\lambda$) case. If we draw $\lambda_{b}\sim p(\cdot)$ a priori and $\mathcal{T}_{b}\sim$ HPP$(\lambda_{b})$, then the set $\{d_{b,i}^{*}\}$ are an i.i.d. sample from Exp(1). This is also true if we generate $\mathcal{T}\sim$ HPP$(\lambda)$ and define $\lambda_{b}\sim\pi(\cdot\mid\mathcal{T})$ as the posterior distribution. If the model is correctly specified, the posterior sample $\{d_{b,i}^{*}\}$ are a realization from Exp(1). In other words, if the posterior sample $\{d_{b,i}^{*}\}$ do not look like Exp(1) samples, then the model is not well specified. The same argument applies to a general point process model.

For the multivariate Hawkes process model, model assessment using RTCT only requires evaluating $d_{b,i}^{*}=t_{b,i}^{*}-t_{b,i-1}^{*}=\int_{t_{i-1}}^{t_{i}}\lambda(u\mid\mathcal{H}_{u};\mbox{\boldmath$\theta$}_{b})du$ for $i=1,\dots,n$ and for each posterior sample $\mbox{\boldmath$\theta$}_{b}$. We can obtain the estimated posterior mean $\hat{d}_{i}^{*}$ of $d_{i}^{*}$ with associated uncertainty and examine if it resembles a sample from Exp(1).

To examine model adequacy, we use the Q-Q plot for $\{d_{i}^{*}\}$ against an Exp(1) distribution. Specifically, we obtain the order statistics $\{\hat{d}_{(i)}^{*}\}$ of the estimated posterior means $\{\hat{d}_{i}^{*}\}$ and plot them against the theoretical quantiles for an Exp(1), i.e., $\log\frac{n}{n-i-0.5}$. The Q-Q plot can be overlaid with the associated 95% credible band. We look for fidelity to a $45^{o}$ reference line as an assessment of model adequacy.

In this manuscript we consider four different models: (i) NHPP, (ii) NHPP+GP, (iii) NHPP+CC, (iv) NHPP+GP+CC. Visual comparison of the Q-Q plots, with associated uncertainties provides informal model comparison. Using the metric $\frac{1}{n}\sum_{i}\big{(}\hat{d}_{(i)}^{*}-\log\frac{n}{n-i-0.5}\big{)}^{2}$, which quantifies the mean squared difference between sample and theoretical quantiles, we can assign a measure of model performance associated with the Q-Q plot.

For model comparison we consider DIC = $\widehat{-2\log L}+p_{D}$. The first term is the estimated posterior mean of $-2\log L(\mbox{\boldmath$\theta$}\mid\mathcal{T})$ and rewards model fit. The second term $p_{D}=\widehat{-2\log L}+\log L(\hat{\mbox{\boldmath$\theta$}}\mid\mathcal{T})$, where $\hat{\mbox{\boldmath$\theta$}}$ is the estimated posterior mean of $\theta$, is the estimated effective number of parameters and penalizes model complexity. The DIC is easily evaluated using posterior samples in the Bayesian framework.

First we analyze NARW up-calls observed by a single recorder (Section 2.1). We set the number of MARUs $K=1$ and fit Models (i)-(iv) described in Section 3.

We assess model adequacy using the RTCT. Figure 6 displays the posterior mean empirical Q-Q plots for $\{d^{\ast}_{i},i=1,\dots,n\}$ against an Exp(1) distribution. For Model (i) the points significantly deviate from the 45^∘ reference line, which is not captured by the 95% credible band (CB). This implies that the NHPP is not adequate for this dataset. The 95% CBs of Models (ii) and (iv) are wider than the other models, attributable to the extra uncertainty from the Gaussian process, and do include the reference line. The points of Models (iii) and (iv) follow the line best although they slightly curve away near the upper tail. As for the mean squared difference (MSD) between sample and theoretical quantiles, Model (iv) provides the smallest MSD = 0.006, followed by Model (iii) with the second smallest MSD = 0.012, Model (ii) with the third smallest MSD = 0.08, and Model (i) with the largest MSD = 24.02.

We use DIC for model comparison. Table 2 shows that Model (iv) has a significantly small value of $-2\log L$ compared to Models (i) and (iii). This implies that Model (iv) fits the data significantly better than Models (i) and (iii). Model (ii) fits the data similarly well to Model (iv) but has more than twice as large of a model complexity penalty as Model (iv). This may indicate that the GP component of Model (ii) inappropriately accounts for variability from excitement. DIC provides consistent results with MSD. Model (iv) has the smallest DIC, followed by Model (iii) with the second smallest DIC, Model (ii) with the third smallest DIC, and Model (i) with the largest DIC. Collectively, these metrics support our choice of Model (iv) for this dataset and for which we present further model inference.

The estimated posterior mean of the coefficient $\beta_{1}$ of the ambient noise variable is -0.49 with 95% HPD interval of (-0.71, -0.28). This significant effect implies that NARWs are less likely to make up-calls during periods with high levels of ambient noise. The observed responses of NARWs to increased levels of noise include changes in sound level (Parks et al.,, 2010; Palmer et al.,, 2022), frequency content (Urazghildiiev,, 2014) and up-calling rate (CWC pers obs.). Here we quantify the reduction in up-calling behavior as background noise levels increases, a response that is consistent with previous work on NARWs in CCB (Urazghildiiev,, 2014) and bowhead whales (Balaena mysticetus) in the Arctic (Blackwell et al.,, 2013). Figure 7 shows the resulting estimated overall daily effect as captured by the harmonics over time. Shown is the posterior estimate of $\hat{\beta}_{2}\sin(\frac{2\pi t}{8\times 60})+\hat{\beta}_{3}\cos(\frac{2\pi t}{8\times 60})+\hat{\beta}_{4}\sin(\frac{2\pi t}{12\times 60})+\hat{\beta}_{5}\cos(\frac{2\pi t}{12\times 60})+\hat{\beta}_{6}\sin(\frac{2\pi t}{24\times 60})+\hat{\beta}_{7}\cos(\frac{2\pi t}{24\times 60})$ over the 24-hour window, which better depicts the diel behaviors than any of the individual coefficient estimates separately. Contact calls tend to increase during the twilight hours between dusk and sunset, as well as around midnight, which has been well documented in NARW across a broad range of habitats (Matthews et al.,, 2001; Matthews and Parks,, 2021).

From this analysis of NARW up-call times of occurrence, we find strong evidence of excitement. The estimated posterior mean of the excitement parameter $\alpha$ is 0.34 with 95% HPD interval of (0.31, 0.38), which is much greater than 0. The estimated posterior mean of the decay parameter $\eta$ is 0.51 with 95% HPD interval of (0.44, 0.58). This implies that the expected median response time for NARW up-calls is approximately $\frac{\log(2)}{0.51}=1.36$ minutes, consistent with direct measurements of up-calls from animals tagged with DTAGs in other NARW habitats (Parks et al.,, 2010). During the study period, the estimated posterior mean for the expected total number of up-calls is 5301 with 95% HPD interval of (5162, 5448), which contains the observed total number of up-calls of 5302. The estimated posterior means of the expected number of contact calls and counter-calls are 1735 with 95% HPD interval of (1364, 2066) and 3567 with 95% HPD interval of (3223, 3947), respectively. We observe approximately twice as many counter-calls as contact calls on average.

To study spatio-temporal patterns of NARW up-calls across CCB, we apply Models (i)-(iv) described in Section 3 to the sequence of up-call times received by the array of ten MARUs within the region (Section 2.2).

Model adequacy for the multi-channel data is assessed using the RTCT. Figure 8 shows that Model (i) is again inadequate for these NARW up-call data, as evident by the points and 95% credible band (CB) that significantly deviate from the reference line. The points of Models (ii) to (iv) follow the line better, but deviate slightly near the upper tail. As with the results for the single-channel data, the 95% CBs of Models (ii) and (iv) are wide due to the extra uncertainty from the Gaussian process; only these two models completely include the reference line. As for the mean squared difference (MSD) between sample and theoretical quantiles, Models (ii) to (iv) provide the smallest MSD = 0.02 and Model (i) provides the largest MSD = 1.87.

Table 3 reports the DIC values for model comparison. Notably, Model (iv) has a significantly smaller value of $-2\log L$ than the other models. Model (iv) also achieves the smallest DIC, indicating this is the preferred model for these data. We again select Model (iv), that with the GP and CC, as our final model for this dataset and present the associated results.

Figure 9 (a) displays the estimated posterior mean and HPD intervals of the noise coefficient $\beta_{1k}$ for each MARU $k$. The posterior mean estimate is negative for each MARU, and the coefficient is significantly negative as evident by the 95% HPD for MARUs located in the bottom edge of the array, i.e., MARUs 3, 4, 5, 9, and 10. The 95% HPD intervals of MARU 1 and 7 include zero but have upper bounds very close to zero. MARU 6 has a very wide HPD interval, which we attribute to the small sample size for this MARU; the expected number of unique up-calls at MARU-6 is estimated to be only 26. Figure 9 (b) shows the estimated overall daily effect as captured by the harmonics over the 24-hour window for each MARU. In general, contact calls tend to increase during the twilight hours between dusk and sunset. We observe a significant increase in up-calling during the night for MARUs 1, 7, 8, and 10. Since MARU 10 is close to a shipping lane, it’s possible that up-calling increases at nighttime with lower shipping. NARW have been known to up-call more at night (Matthews et al.,, 2001; Matthews and Parks,, 2021); however because we lack visual observations during this period, we cannot say for sure if this is in response to non-foraging aggregations or distribution changes within and across the bay.

Except for MARU 10, we find strong evidence of excitement in NARW up-calls; the 95% HPD interval of $\alpha_{k}$ is greater than zero for $k=1,\dots,9$. MARU 10 is located in a region with the highest shipping traffic. This might prevent whales from hearing and responding to up-calls and disturb communication between whales. The estimated posterior mean of the temporal decay parameter $\eta$ is 0.151 with 95% HPD interval of (0.150, 0.154). This implies that the expected median response time for NARW up-calls is approximately $\frac{\log(2)}{0.151}=4.59$ minutes. The estimated posterior mean of the spatial decay parameter $\phi$ is 0.32 with 95% HPD interval of (0.29, 0.35). The shortest distance between MARUs is 7.1 km. This implies that the probability of a counter-call being received at a MARU at least 7.1km away is approximately $P(\text{Exp}(0.32)>7.1)=0.10$.

During the study period, the estimated posterior mean for the expected total number of up-calls is 2716 with 95% HPD interval of (2613, 2820) which contains the observed total number of up-calls of 2750. The model also recovers the observed number of up-calls recorded by each MARU within the 95% HPD intervals. Figure 10 (a) displays the estimated posterior mean for the expected number of within-recorder counter-calls, i.e., counter-calls at MARU $k$ that were excited by calls at the MARU. Figure 10 (b) shows the estimated posterior mean for the expected number of cross-recorder counter-calls, i.e., counter-calls at MARU $k$ that were excited by up-calls at MARU $\ell$ for $\ell\neq k$. We see that MARU 10 received a relatively small number of within- and cross-recorder counter-calls. MARU 10 is closest to the path of vessels that transit from the Cape Cod Canal to Boston Harbour. It is the loudest MARU (Figure 2), which may indicate that underwater noise caused by shipping can inhibit whales from being able to communicate with each other. Elevated ambient noise can also impact our ability to detect up-calls (Palmer et al.,, 2022). MARUs 1, 4, 5, and 8 received a substantial number of within-recorder counter-calls compared to the others. In contrast to MARU 10, MARU 5 was one of the quietest MARUs (Figure 2). MARUs 1 and 8 are at the top of the array (Figure 1), and thus are likely hearing up-calls that are made outside of CCB as well. CCB is a relatively small area, and sound can be heard across much of the Bay; the average swim speed of NARW is approximately 3.7 km per hour, and the shortest distance between MARUs is 7.1 km. Thus, the excitement across recorders likely indicates communication between different whales, which is consistent with the role of the contact call (Clark and Clark,, 1980; Clark,, 1982).

Figure 11 shows the joint posterior distribution of the expected number contact calls and counter-calls received at each MARU. The majority of MARUs received a similar or greater number of counter-calls compared to contact calls. Uniquely, MARU 3 received twice as many contact calls as counter-calls and MARU 10 received mostly contact calls and very few counter-calls. More detailed evidence of this behavior is presented in the supplementary material. Figure S5 displays the posterior mean estimates of the background intensity and the counter-call intensity. The panel for MARU 10 further depicts the sparsity of expected counter-calls there, reflecting the estimate of $\alpha_{10}$ which is essentially $0$.

Figure S6 of the Supplementary Material aggregates the two intensity components for each MARU and overlays, with rug tassels, the observed event sequence (Figure 2 (b)). It offers informal/visual evidence of good model performance for the observed data.