Skip to contents

Non-uniform study areas: distance truncation and coastlines

Source: vignettes/callDensity_coast.Rmd
callDensity_coast.Rmd
library(callDensity)
library(ggplot2)
library(scales)
library(kableExtra)
library(VGAM)
#> Loading required package: stats4
#> Loading required package: splines

Further exploration and topics in callDensity package (TODO)

  1. Spatial distribution of calls
    1. Effect of land, ice, and other zero-probability areas
    2. Effect of non-uniform distribution of calls in area
  2. Truncation in SNR
  3. Detectors with different slopes
  4. Minimum number of recaptures/effort required to get a reliable estimate
  5. More realistic Transmission Loss models and effects of mismatch

As usual, we start with the density equation:

Dc=Nc(1−c)kTpaA\begin{equation} D_c = \frac{N_c(1-c)}{kTp_aA} \end{equation}

where: DcD_c is call density, NcN_c is number of calls, cc is false discovery rate, kk is number of sensors (here always 1), TT is the duration of data analysed (in hours), pap_a is the probability of detection in the study area, and AA is the study area (in km^2).

1A) Effect of land and other zero-probability areas

Scenario: Antarctic blue whales in a quiet ocean with spherical spreading of sounds. Instead of open-ocean, the recorder is located 144 km north of a coastline.

In this example, we provide truncation distances to the Monte-Carlo simulation to exclude the coastal regions from the simulation. We also amend our estimate of the study area to account for the area excluded by coasts.


n = 1e6; # number of simulated calls (regardless of whether detected or not); 
R = 1e6; # radius of 1000 km (but in m)
k <- 1   # number of sensors
minDate <- as.POSIXct("2025-01-01")
maxDate <- as.POSIXct("2025-12-31")
Time <- as.numeric(difftime(maxDate,minDate,unit="days"))/365 # in years

# 1) Generate n uniformly distributed calls within radius R and time period Time
sim <- simCallLocation(n=n, R=R, minDate=minDate, maxDate=maxDate)

A_circ = studyArea(R/1e3)  # Circular study area in km^2
coast= -144e3;             # Simulate coastline at y = -144 km;
sim <- subset(sim, sim$y>=coast);
n = dim(sim)[1]

# Area of a circle segment given radius r and segment height, h
# A =    r^2 *   acos( (r-h)/r ) - (r-h)*sqrt(2*r*h-h^2)
h = R-abs(coast) # Height of segment

# Area of excluded segment
A_exclude = (R^2 * acos( (R-h)/R) - (R-h)*sqrt(2*R*h-h^2) )/ 1e3^2 # in km^2
A = A_circ - A_exclude # Subtract the excluded area from the circle (in km^2)

TrueCallDensity = n/(A*Time)

# Map showing true location of all calls
ggplot(data=sim, aes(x=x/1e3,y=y/1e3))+
  # geom_point(colour = "black",alpha = 0.1, size=0.1)+
  geom_bin_2d(alpha=1,binwidth=c(10,10))+
  coord_equal()+
  xlab("X location (km)")+
  ylab("Y location (km)")

Sonar equation parameters for this scenario:

TL=20log10d\begin{equation} TL = 20 log_{10} d \end{equation}SNR=SL−TL−NL\begin{equation} SNR = SL-TL-NL \end{equation}

In previous simulations we had four radial transects. Here we create four additional transect radials (eight in total) in order to approximate the study area with more fidelity.


SL <- data.frame(mean=190, sd=4, sampleSize=350) # True SL distribution for sim
NL = data.frame(mean=84, sd = 4, sampleSize = n)  # True NL distribution for sim
tlFunc <- function(r)20*log10(r)

# Simulate the true acoustic properties of each call and false positive 
sim <- simCallAcoustics(sim, SL, NL, TL = tlFunc)

# Per-transect estimate of TL (used by callDensity package to estimate pDet)
numTransects=8
TL <- simTLradials_20logR(maxRange=R, rangeStep=100, numTransects=numTransects)

Estimate Truncation distances

The TL models in these transects assume spherical spreading (without regard for the coast).

Regardless of whether the TL models are a good match for the physical environment or not, we need to exclude from the simulation (and study area) portions of the transects where the animals cannot be located (e.g. on land).

# Work out the range for each radial over which to truncate. 
radials = seq(from=0,to=315, by=360/numTransects)
range_m = seq(from=5,to=R,by=100)

# Find first r where the r*sind(radials) <= coastline
truncDist <- matrix(R,nrow=1,ncol = length(radials))

# Loop over each radial. find the range at which it encounters the coast, and
# truncate (i.e. Distances beyond this range will be excluded from the analysis)
for (i in 1:length(radials)){ 
  rangeIx <- which(cos(radials[i]*pi/180)*range_m <= coast)[1]
  truncDist[i] <- ifelse(is.na(rangeIx),R,range_m[rangeIx])
}

# callDensity package provides a function to calculate area from radials with
# different lengths
Atrunc <- callDensity::studyArea(R,truncDist)/1e6

Simulate detectors:

Detector Description Location/SNR ‘threshold’ Scale/slope False positive rate
Detector 1 Reliable and good detector, meant to emulate a human analyst. Lower Same Lower
Detector 2 Automated detector, meant to emulate a signal processing algorithm. Higher Same Higher 

#Specify parameters for detector 1
# Good detector with low threshold (location) and low false positive rate.
det1params = data.frame( 
  location=1,    # AKA intercept?
  scale=1,       # AKA slope?
  func='plogis', # logistic function
  c=0.1,         # False discovery rate
  fpMean=0,      # Mean of distribution of false positives (in dB SNR)
  fpSD = 1      # Standard deviation false positive distribution (in dB SNR)
)

# Detector 2 has a higher threshold (location), scale, and false positive rate
# than detector 1.
det2params = data.frame( 
  location=3,    # AKA intercept?
  scale=2,       # AKA slope?
  func='plogis', # logistic function
  c=0.3,         # False discovery rate
  fpMean=0,      # Mean of distribution of false positives (in dB SNR)
  fpSD = 1      # Standard deviation false positive distribution (in dB SNR)
)

# Simulate two detectors each with the different parameters
simDet1<- simulateDetector(det1params,sim)
simDet2<- simulateDetector(det2params,sim)

View the spatial detection density (which is not the same as spatial call density since it accounts for neither detection probability nor time).


p1<-plotSpatialDetections(simDet1)
p2<-plotSpatialDetections(simDet2)
gridExtra::grid.arrange(p1,p2,nrow=1)

Subsample data for detector characterisation

Castro et al. (2024) subsample approximately 200 hours evenly spaced throughout the year to create detector characterisation curves (detection vs SNR). Here we subsample 1 hour in every 41 hours yielding 214 subsamples that are 1 hour in duration.


subsampleDet1 <- subsampleSimInTime(simDet1)
subsampleDet2 <- subsampleSimInTime(simDet2)

Detection matching between detector 1 and 2

Both detectors (simulated subsets) operate on the same set of calls (i.e. they are different subsets of the same underlying simulation). So we can find detections from one detector that match those of the other detector by merging the detections into a capture history table.

Presently, false positives have randomly generated date-times, so have a negligible chance of matching between the two detectors.

TODO: Create false positives in a way that better approximates their real-world occurrence (i.e. in a way that allows them to be matched between detectors).


capHistTab<- simsTocaptureHistoryTable(subsampleDet1,subsampleDet2)

# simsTocaptureHistoryTable() keeps every observer's own signal/noise/SNR
# columns separately, suffixed (matchbox-native), rather than
# auto-consolidating them the way it used to. cde()'s own
# falseDiscoveryRate() call needs a single SNR column for its truncation
# handling (a separate mechanism from chtToSNRinfo()'s signalCol/noiseCol
# averaging below), so that's built explicitly here.
capHistTab$SNR <- rowMeans(capHistTab[, c("snr_observer1", "snr_observer2")], na.rm = TRUE)

# Total number of true positive detections for detector1 OR detector2
nDetectedSubset <- with(capHistTab,sum(detect_observer1 & groundTruth |   
                                        detect_observer2 & groundTruth) )

# Total number of positive detections for detector1 OR detector2
nPositiveSubset <- with(capHistTab,sum(detect_observer1 | detect_observer2))

# Number positive detections in FULL set
Nc <- sum(simDet2$detect_table)       

# SNRinfo is used to calculate snrDetFun for detector2 assuming detector1 is the
# ground truth. This is an observer ground (OG) so includes false positives from
# detector1, but does not include detections from detector2 not detected on
# detector1. signalCol/noiseCol point at each observer's own suffixed
# signal/noise columns, which chtToSNRinfo() averages per event -- there is
# no single consolidated signalRMSdB/noiseRMSdB column to fall back on.
SNRinfo <- chtToSNRinfo(capHistTab, groundTruth = "observer1", observers = "observer2",
                        signalCol = c("signalRMSdB_observer1", "signalRMSdB_observer2"),
                        noiseCol  = c("noiseRMSdB_observer1",  "noiseRMSdB_observer2"))

Calculate call densities using callDensity package

We fit a capture-recapture model to these adjudicated data using the VGAM package with the function ‘vglm’ (Yee et al. 2015).


snrDetFun.glm <- callDensity::fitDetFun(SNRinfo, modelType = "glm")

results.glm <- cde(Nc=Nc, capHistTab = capHistTab, snrDetFun = snrDetFun.glm,
                   SL = SL, TL = TL, A = A, k = k, T = Time, modelType = 'glm', 
                   truncationDistance = truncDist,
                   groundTruthCol = "detect_observer1", observerCol = "detect_observer2",
                   signalCol = c("signalRMSdB_observer1", "signalRMSdB_observer2"),
                   noiseCol  = c("noiseRMSdB_observer1",  "noiseRMSdB_observer2"))

### Adjudicated capture recapture model
adjudicated <- subset(capHistTab,capHistTab$groundTruth &
                (capHistTab$detect_observer1 | capHistTab$detect_observer2))

# Combine SNR of both detectors (shouldn't be necessary unless we've somehow
# made them different)
adjudicated$SNR <-   rowMeans(
  subset(adjudicated,select=c('snr_observer1','snr_observer2')) ,na.rm=T) #dB

observerNames = c("detect_observer1", "detect_observer2")
snrDetFun.vglm <- fitDetFun(adjudicated, modelType = "vglm", yColNames = observerNames,
                             whichObserver = "detect_observer2")

# Best estimate of NL 
NLadj <- nlFromSnrInfo(chtToSNRinfo(capHistTab, groundTruth = "observer1", observers = "observer2",
                                    signalCol = c("signalRMSdB_observer1", "signalRMSdB_observer2"),
                                    noiseCol  = c("noiseRMSdB_observer1",  "noiseRMSdB_observer2")),
                       snrDetFun.vglm)

# groundTruthCol/observerCol point cde() at the simulated ground truth
# (groundTruth) directly, matching the manual fp/precision calculation
# this replaces (fp/(fp+tp), the same definition falseDiscoveryRate()
# already computes internally) -- no separate pDetInArea() call or
# hand-derived Dc formula needed. returnPdetDetail=TRUE keeps the
# per-transect detail available below for the same diagnostic plots
# pDetInArea() itself would have given.
results.vglm <- cde(Nc=Nc, capHistTab = capHistTab, snrDetFun = snrDetFun.vglm,
                   SL = SL, TL = TL, NL = NLadj, A = Atrunc, k = k, T = Time,
                   modelType = 'vglm', truncationDistance = truncDist,
                   groundTruthCol = "groundTruth", observerCol = "detect_observer2",
                   signalCol = c("signalRMSdB_observer1", "signalRMSdB_observer2"),
                   noiseCol  = c("noiseRMSdB_observer1",  "noiseRMSdB_observer2"),
                   returnPdetDetail = TRUE)

resultsTrue <- data.frame(season='year', siteCode='', Nc=n, c=det2params$c,
                          k=k, T=Time, A=A, pa=mean(simDet2$p_det,na.rm=TRUE),
                          SLmean=SL$mean, SLsd=SL$sd, NLmean=NL$mean, NLsd=NL$sd,
                          modelType=det2params$func, CV.Nc=0, CV.c=0, CV.pa=0,
                          Dc=TrueCallDensity, CV.Dc=0)

results <- rbind(resultsTrue, results.glm, results.vglm)

kableExtra::kbl(results[, c('Dc','Nc','A','NLmean','NLsd','modelType','c','pa','CV.Dc')],
                digits = c(4,0,0,1,2,NA,3,4,4),
                col.names= c('Dc','Sample Size','A','NL_mean','NL_sd',
                             'Model','c','P_a','CV.Dc')) %>% 
  kableExtra::kable_classic(full_width=FALSE)
Dc Sample Size A NL_mean NL_sd Model c P_a CV.Dc
0.3195 591874 1857798 84.0 4.00 plogis 0.300 0.0770 0.0000
0.2699 65240 1857798 83.6 4.03 glm 0.456 0.0710 0.4137
0.3050 65240 2004230 84.0 4.12 vglm 0.297 0.0752 0.3899

The estimate of Dc derived from the vglm model is very close to the true value for this simulation.


plotPDetRadials(attr(results.vglm, "pDetResults"))

The spatial exclusions and differing slopes of the detectors appear to have been modelled with acceptable accuracy.

References

Castro, Franciele R., Danielle V. Harris, Susannah J. Buchan, Naysa Balcazar, and Brian S. Miller. 2024. “Beyond Counting Calls: Estimating Detection Probability for Antarctic Blue Whales Reveals Biological Trends in Seasonal Calling.” Frontiers in Marine Science 11 (July). https://doi.org/10.3389/fmars.2024.1406678.
Yee, Thomas W., Jakub Stoklosa, and Richard M. Huggins. 2015. “The VGAM Package for Capture-Recapture Data Using the Conditional Likelihood.” Journal of Statistical Software 65 (June): 1–33. https://doi.org/10.18637/jss.v065.i05.