GoMx Cuvier’s beaked whale habitat models: Density with GAMs
Kaitlin Frasier
Scripps Institution of Oceanography, UC San Diego19 June 2020
1. Exploratory analysis
1.1 Data Inputs
NOAA SEFSC visual data goes back to 1992, but as shown in the figure below, many predictor variables are only available starting in 2003, therefore earlier visual data is currently excluded from further analyses.
Note: Future work could use monthly climatologies (averages) so that older sightings data could be used. Some dynamic drivers like eddy and front locations would not be able to be considered using that approach.
Visual data predictor variable availability:
1.1.1 Splitting into testing and training sets
The data are split into training and testing sets. In this case, visual data from 2009 and acoustic data from 2013 were used only for testing. Only observations from 2003 or later were used for modeling due to covariate limitations.
1.1.2 Map of visual sightings data
The visual data selected for modeling are displayed on the map below. Data from 2009 were held back for testing. Blue markers indicate HARP locations.
1.1.3 Time series of acoustic data
The time series below show timeseries of estimated densities from passive acoustic data used for modeling (Densities were calculated following methods detailed in [1]). Data from 2011 and 2012 were used for training, and 2013 data was held back for testing.
Acoustic Timeseries:
1.2 Examination of covariates
1.2.1 Covariate distribution check
Distributions of covariates from acoustic observations (training data only):
Distributions of covariates from the visual observations (training data only):
Some of these covariates are more or less interrelated. Correlations are examined in the figure below. Numbers closer to 1 above the diagonal in the figure below represent correlation coefficients. If a pair of covariates is highly-correlated only one should typically be used in the model.
Covariate Correlations:
1.2.2 Transformation of predictor variables
Some variables, including chlorophyll, mixed layer depth and distance to fronts are highly skewed and were log-transformed for input to GAMs.
Below, the two sets of covariates have been combined and transformed:
1.2.3 Preliminary check of predictive power
To get an idea of the basic predictive power of these covariates, we can look at presence/absence relative to each variable. This also provides an opportunity to look at the range of values observed for each covariate in the visual and acoustic datasets. In the plots below dotted lines indicate the distribution of each covariate when Cuvier’s beaked whale were present, and solid lines indicate the distribution when Cuvier’s beaked whale were absent. Note that these plots do not account for effort.
Acoustic kernel densities:
Visual kernel densities:
1.2.4 Estimation of relative weights
To train the model, we need to know how much power the various data points have relative to one another. This is important because the duration, spatial coverage, and detection probabilities are quite different between the visual and acoustic data sets. If an animal is seen or heard, we know for certain that the species was present. However, if it was not heard, then it either wasn’t present, or it was present but missed.
“Zero inflation” or an excess of false zeros more common in the visual survey data because each data point represents a 10km transect section, traversed at survey speed (>10 knots) or approximately 30 minutes of observation effort. In contrast, the acoustic data are binned by day with a stationary instrument, therefore the probability of missing a group over the course of a day is lower.
For each data type, we estimated the probability of a missed detection to account for differences in zero-inflation, downweighting zeros according to the probability of a recording false negative.
The visual data represent wether or not Cuvier’s beaked whale were seen during each transect segment. The probability of missing a sighting of Cuvier’s beaked whale was estimated as
\[P_{V}(detect|present) = \mu_{det} * g0\ = 0.46 * 0.29 = 0.13\]
where _{det} is the mean detection probability as estimated by a model fit using the mrds package, and g0 is the probability of observing an animal on the transect line [2]. We assume that reported absences are likely to be true absences 13 % of the time, therefore zeros are given a weight of 0.13 on a scale of [0,1].
The acoustic data represent presence or absence of Cuvier’s beaked whale in one day bins. Given that a group of animals is present near the sensor, the probability of detecting them in a 5 minute period within a 4 km range is estimated at 0.36, therefore the probability of missing an encounter is 1 - 0.36 = 0.64 [1,3]. Given that animals were present, the probability of missing a group for a full day (288 5-minute periods) is estimated as \[P_{A}(detect|present) = (1-0.36)^{288} \approx 1\]
Therefore we assume that there are no false negative days in the passive acoustic timeseries, and all acoustic observations are given weight = 1.
Best visual detection probability model for Cuvier’s beaked whale:
2. Model Fitting
Models were fit using gam from the mgcv package in R [4].
2.1 Set up weighting
2.2 Run GAMs
Run GAMs on acoustic-only, visual-only, and joint acoustic/visual datasets.
varNames <- names(transformedCovars_AcOnly.train)
varNamesFormula <- varNames[2:10]
formulaAcOnly_allVars <-paste('s(', varNamesFormula, ", bs = 'ts', k = kVal)", sep = "", collapse = ' + ')
formulaAcOnly_allVars <- as.formula(paste('yAcOnly ~',formulaAcOnly_allVars))
gam_full_AcOnly <- gam(formulaAcOnly_allVars,
data = transformedCovars_AcOnly.train,
na.action = na.omit,family = tw())formulaVisOnly_allVars <-paste('s(', varNamesFormula, ", bs = 'ts', k = kVal)", sep = "", collapse = ' + ')
formulaVisOnly_allVars <- as.formula(paste('yVisOnly ~',formulaVisOnly_allVars))
gam_full_VisOnly <- gam(formulaVisOnly_allVars,
data = transformedCovars_VisOnly.train,
na.action = na.omit,
family = Tweedie(p = 1.5, link = log),
weights = VisOnly.train_weightsG0)formulaJoint_allVars <-paste('s(', varNamesFormula, ", bs = 'ts', k = kVal)", sep = "", collapse = ' + ')
formulaJoint_allVars <- as.formula(paste('y ~',formulaJoint_allVars))
gam_full_Joint <- gam(formulaJoint_allVars,
data = transformedCovars.train,
weights = joint_train_weightsG0,
na.action = na.omit,
family = Tweedie(p = 1.5, link = log))3. Model Summaries
Model summaries show that some terms were not significant.
3.1 Acoustic-only model
Model summary:
##
## Family: Tweedie(p=1.298)
## Link function: log
##
## Formula:
## yAcOnly ~ s(SST, bs = "ts", k = kVal) + s(SSH, bs = "ts", k = kVal) +
## s(log10_HYCOM_MLD, bs = "ts", k = kVal) + s(HYCOM_SALIN_0,
## bs = "ts", k = kVal) + s(HYCOM_UPVEL_50, bs = "ts", k = kVal)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.13197 0.05806 2.273 0.0231 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(SST) 1.8664 2 198.259 < 2e-16 ***
## s(SSH) 1.8550 2 42.665 < 2e-16 ***
## s(log10_HYCOM_MLD) 1.7893 2 23.756 6.32e-12 ***
## s(HYCOM_SALIN_0) 1.7274 2 232.812 < 2e-16 ***
## s(HYCOM_UPVEL_50) 0.8236 2 1.944 0.0277 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.294 Deviance explained = 45%
## -REML = 2770.6 Scale est. = 4.2359 n = 1789Smooths of predictor variables:
## png
## 2
Plot of residuals:
3.2 Visual-only model
Model summary:
##
## Family: Tweedie(1.5)
## Link function: log
##
## Formula:
## yVisOnly ~ s(SST, bs = "ts", k = kVal) + s(log10_CHL, bs = "ts",
## k = kVal) + s(log10_HYCOM_MLD, bs = "ts", k = kVal) + s(HYCOM_SALIN_0,
## bs = "ts", k = kVal) + s(log10_HYCOM_MAG_0, bs = "ts", k = kVal) +
## s(HYCOM_UPVEL_50, bs = "ts", k = kVal) + s(Neg_EddyDist,
## bs = "ts", k = kVal) + s(Pos_EddyDist, bs = "ts", k = kVal)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.308 0.190 -6.885 7.8e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(SST) 1.975 2 48.875 < 2e-16 ***
## s(log10_CHL) 1.984 2 28.092 7.23e-13 ***
## s(log10_HYCOM_MLD) 1.992 2 42.951 < 2e-16 ***
## s(HYCOM_SALIN_0) 1.998 2 41.204 < 2e-16 ***
## s(log10_HYCOM_MAG_0) 1.951 2 3.015 0.0456 *
## s(HYCOM_UPVEL_50) 1.002 2 10.506 3.53e-06 ***
## s(Neg_EddyDist) 2.000 2 22.820 1.29e-10 ***
## s(Pos_EddyDist) 1.997 2 35.125 8.52e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = -3.71 Deviance explained = 36.5%
## GCV = 0.83949 Scale est. = 2.1197 n = 1956
Smooths of predictor variables:
## png
## 2
Plot of weighted residuals:
3.3 Joint model
Model summary:
##
## Family: Tweedie(1.5)
## Link function: log
##
## Formula:
## y ~ s(SST, bs = "ts", k = kVal) + s(SSH, bs = "ts", k = kVal) +
## s(log10_CHL, bs = "ts", k = kVal) + s(log10_HYCOM_MLD, bs = "ts",
## k = kVal) + s(HYCOM_SALIN_0, bs = "ts", k = kVal) + s(log10_HYCOM_MAG_0,
## bs = "ts", k = kVal)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.62601 0.05404 11.58 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(SST) 1.8507 2 69.536 < 2e-16 ***
## s(SSH) 0.9828 2 2.573 0.01841 *
## s(log10_CHL) 1.0205 2 3.814 0.00431 **
## s(log10_HYCOM_MLD) 1.9935 2 15.106 2.42e-07 ***
## s(HYCOM_SALIN_0) 1.9851 2 133.529 < 2e-16 ***
## s(log10_HYCOM_MAG_0) 1.8213 2 7.763 0.00021 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0726 Deviance explained = 22.6%
## GCV = 1.5516 Scale est. = 3.6868 n = 3746Smooths of predictor variables:
## png
## 2Plot of weighted residuals:
3.4 Model prediction on joint dataset
## [1] "Mean absolute error scores (lower is better)"## MAE
## [1,] 2.826
## [2,] 3.510
## [3,] 2.637## [1] "Weighted mean absolute error scores (lower is better)"## wMAE
## Acoustic - Test Joint Data 0.990
## Visual - Test Joint Data 1.074
## Joint - Test Joint Data 0.909## [1] "RMSE scores (lower is better)"## RMSE
## Acoustic - Test Joint Data 4.966
## Visual - Test Joint Data 14.375
## Joint - Test Joint Data 4.373## [1] "Weighted RMSE scores (lower is better)"## wRMSE
## Acoustic - Test Joint Data 2.889
## Visual - Test Joint Data 5.831
## Joint - Test Joint Data 2.6194. Model Predictions
4.1 Temporal predictions
Predictions were made on the acoustic test dataset, and compared with actual observations for 2013. The predicted density of animals at each site was compared with the estimated daily density from the pasive acoustic record.
4.1.1 Acoustic-only prediction
Predicted and observed densities at passive acoustic monitoring sites using the acoustic-only model:
## geom_custom_ann: grob = list(name = "GRID.gTree.11", gp = NULL, vp = NULL, children = list(GRID.text.12 = list(label = "MC", x = 0.01, y = 0.9, just = "centre", hjust = 0, vjust = NULL, rot = 0, check.overlap = FALSE, name = "GRID.text.12", gp = list(col = "black", fontsize = 13, fontface = "bold", font = c(bold = 2)), vp = NULL)), childrenOrder = "GRID.text.12"), xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf
## stat_identity: na.rm = FALSE
## position_identity## png
## 24.1.2 Visual-only prediction
Predicted and observed densities at passive acoustic monitoring sites using the visual-only model:
## png
## 24.1.3 Joint prediction
Predicted and observed encounter probabilities at passive acoustic monitoring sites using the joint model:
## png
## 24.2 Spatial predictions
Models were evaluated for summer (July 2009) and winter(January 2009) across the entire Gulf of Mexico (US EEZ beyond the 200m contour).
Note: As an alternative to training joint models, an average surface was computed across the acoustic and visual-only predictions.
Each grid cell in the following maps represents a 10x10 km square. Densities are therefore shown as estimated number of animals per 100 km2.
4.2.1 Summer 2009 predictions
Summer 2009 predicted distribution and test sightings:
4.2.2 Winter 2009 predictions
Winter 2009 predicted distribution:
5. Monthly model predictions
Spatial model predictions were generated using climatological means of oceanographic variables, averaged by month between 2003 and 2015.
#References
1. Frasier KE, Wiggins SM, Harris D, Marques TA, Thomas L, Hildebrand JA. Delphinid echolocation click detection probability on near-seafloor sensors. The Journal of the Acoustical Society of America. 2016;140: 1918–1930. doi:http://dx.doi.org/10.1121/1.4962279
2. Palka DL. Summer abundance estimates of cetaceans in us north atlantic navy operating areas. Northeast Fish Sci Cent Ref Doc. 2006; 06–03.
3. Hildebrand J, Baumann-Pickering S, Frasier K, Tricky J, Merkens K, Wiggins S, et al. Passive acoustic monitoring of beaked whale densities in the gulf of mexico during and after the deepwater horizon oil spill. Nature Scientific Reports. 2015;5: 16343.
4. Wood S. Generalized additive models: An introduction with r. CRC press; 2006.