Modelling ocean temperatures from bio-probes under preferential sampling

TL;DR

This paper extends geostatistical models to account for preferential sampling in mobile bio-probe data, improving the reliability of ocean temperature predictions from animal-mounted sensors.

Contribution

It introduces a novel framework for modeling spatial fields with moving, stochastic sampling locations, addressing a gap in existing preferential sampling methodologies.

Findings

Preferential sampling significantly affects SST predictions.

The proposed model outperforms standard models when sampling is preferential.

Incorporating preferential sampling alters the conclusions drawn from SST data analysis.

Abstract

In the last 25 years there has been an important increase in the amount of data collected from animal-mounted sensors (bio-probes), which are often used to study the animals' behaviour or environment. We focus here on an example of the latter, where the interest is in sea surface temperature (SST), and measurements are taken from sensors mounted on Elephant Seals in the Southern Indian ocean. We show that standard geostatistical models may not be reliable for this type of data, due to the possibility that the regions visited by the animals may depend on the SST. This phenomenon is known in the literature as preferential sampling, and, if ignored, it may affect the resulting spatial predictions and parameter estimates. Research on this topic has been mostly restricted to stationary sampling locations such as monitoring sites. The main contribution of this manuscript is to extend this methodology to observations obtained by devices that move through the region of interest, as is the case with the tagged seals. More specifically, we propose a flexible framework for inference on preferentially sampled fields, where the process that generates the sampling locations is stochastic and moving over time through a 2-dimensional space. Our simulation studies confirm that predictions obtained from the preferential sampling model are more reliable when this phenomenon is present, and they compare very well to the standard ones when there is no preferential sampling. Finally, we note that the conclusions of our analysis of the SST data can change considerably when we incorporate preferential sampling in the model.

Figures

Equations

$$Y_i=\mu +S({\text{X}}_i)+Z_i,{\text{X}}_i\in \mathcal{D},i=1,\dots ,n,$$

$$C(r)={\sigma }^2\frac{2^{1-\kappa }}{\Gamma (\kappa )}{\left(\frac{r}{\varphi }\right)}^{\kappa }{K}_{\kappa }\left(\frac{r}{\varphi }\right),r>0,$$

$$L(\mathbf{\theta })=\left[\text{X},\text{Y};\mathbf{\theta }\right]=\int \left[\text{X},\text{Y},S;\mathbf{\theta }\right]\mathrm{d}S.$$

$$L(\mathbf{\theta })=\int \left[\text{X},\text{Y},S;\mathbf{\theta }\right]\mathrm{d}S=\int \left[\text{Y}|S,\text{X};{\mathbf{\theta }}_F\right]\left[\text{X}|S;{\mathbf{\theta }}_L\right]\left[S;{\mathbf{\theta }}_F\right]\mathrm{d}S.$$

$$\int \left[\text{Y}|\text{S},\text{X};{\mathbf{\theta }}_F\right]\left[\text{X}|\text{S};{\mathbf{\theta }}_L\right]\left[\text{S};{\mathbf{\theta }}_F\right]\mathrm{d}\text{S},$$

$$f(\text{S},\mathbf{\theta })=-\log \left(\left[\text{Y}|\text{S},\text{X};{\mathbf{\theta }}_F\right]\left[\text{X}|\text{S};{\mathbf{\theta }}_L\right]\left[\text{S};{\mathbf{\theta }}_F\right]\right),$$

$$\hat{\text{S}}(\mathbf{\theta })=\underset{\text{S}}{\mathrm{argmin}}f(\text{S},\mathbf{\theta }),$$

$$\text{X}(t_{k+1})=g_{\text{NP}}(\text{X}(t_{1:k}),{\mathbf{\theta }}_L)+ϵ({\text{t}}_{1:k},{\mathbf{\theta }}_L)$$

$$\text{X}(t_{k+1})=g_{\text{P}}(\text{X}(t_{1:k}),S,{\mathbf{\theta }}_L)+ϵ({\text{t}}_{1:k},{\mathbf{\theta }}_L)$$

$$\text{X}(t_{k+1})=\text{X}(t_k)+\mathbf{\mu }(\text{X}(t_{1:k}),S,{\mathbf{\theta }}_L)(t_{k+1}-t_k)+\mathbf{\Sigma }({\mathbf{\theta }}_L){\mathbf{A}}_k\sqrt{t_{k+1}-t_k},$$

$$\mathbf{\mu }(\text{X}(t_{1:k}),S,{\mathbf{\theta }}_L)=f({\beta }_{t_k})\mathbf{\varphi }(\text{X}(t_k),S,{\mathbf{\theta }}_L)+(1-f({\beta }_{t_k}))\text{v}(\text{X}(t_{1:k})),$$

$$\text{v}(\text{X}(t_{1:k}))=\frac{\text{X}(t_k)-\text{X}(t_{k-1})}{t_k-t_{k-1}},$$

$$f({\beta }_{t_k})=\frac{\exp ({\beta }_{t_k})}{1+\exp ({\beta }_{t_k})},$$

$$\begin{array}{rl}\mathbf{\mu }(\text{X}(t_k),S,{\mathbf{\theta }}_L)& =\frac{\exp ({\beta }_{t_k})}{1+\exp ({\beta }_{t_k})}\mathbf{\varphi }(\text{X}(t_k),S,{\mathbf{\theta }}_L)+\frac{1}{1+\exp ({\beta }_{t_k})}\text{v}(\text{X}(t_{1:k})),\\ {\beta }_{t_{k+1}}& ={\beta }_{t_k}+{\sigma }_{\beta }{B}_k\sqrt{t_{k+1}-t_k},\end{array}$$

$$\mathbf{\varphi }(\text{X},S,{\mathbf{\theta }}_L)=-\alpha S(\text{X})\nabla S(\text{X}),$$

$$\mathbf{\varphi }(\text{X},S,{\mathbf{\theta }}_L)=-\alpha S^{\ast }(\text{X})\nabla S(\text{X}),$$

$${\mathrm{RMSPE}}_i=\frac{1}{M}\sum _{j=1}^M\sqrt{{\left({\text{S}}_{j,i}-{\hat{\text{S}}}_{j,i}\right)}^2},$$

$${\mathrm{MIGN}}_j=\frac{1}{N}\sum _{i=1}^N\left\{\frac{({\text{S}}_{j,i}-{\hat{\text{S}}}_{j,i}{)}^2}{2{\hat{\sigma }}_{j,i}^2}+\log {\hat{\sigma }}_{j,i}\right\}j=1,\dots ,M,$$

$${\mathrm{LIGN}}_i=\frac{1}{M}\sum _{j=1}^M\left\{\frac{({\text{S}}_{j,i}-{\hat{\text{S}}}_{j,i}{)}^2}{2{\hat{\sigma }}_{j,i}^2}+\log {\hat{\sigma }}_{j,i}\right\}i=1,\dots ,N.$$

$$\begin{array}{rl}{\mathrm{RMSPE}}_i^{\text{Diff}}& ={\mathrm{RMSPE}}_i^{\text{P}}-{\mathrm{RMSPE}}_i^{\text{NP}},\\ {\mathrm{MIGN}}_j^{\text{Diff}}& ={\mathrm{MIGN}}_j^{\text{P}}-{\mathrm{MIGN}}_j^{\text{NP}},\\ {\mathrm{LIGN}}_i^{\text{Diff}}& ={\mathrm{LIGN}}_i^{\text{P}}-{\mathrm{LIGN}}_i^{\text{NP}},\end{array}$$

$${\text{D}}_i=({\hat{S}}_{1,i}^{\text{P}}-{\hat{S}}_{1,i}^{\text{NP}},\dots ,{\hat{S}}_{50,i}^{\text{P}}-{\hat{S}}_{50,i}^{\text{NP}}{)}^T,$$

$${\text{C}}_{i,j}$$

$${\text{G}}_{i,j}$$

$$({\text{K}}_{\varphi }{)}_{i,j}$$

$${\text{Q}}_{1,\varphi }$$

$${\text{Q}}_{2,\varphi }$$

$${\text{Q}}_{\alpha ,\varphi }$$

$${\text{Q}}_{2,\varphi }$$

$$=({\varphi }^{-2}\text{C}+\text{G}){\text{C}}^{-1}({\varphi }^{-2}\text{C}+\text{G}),$$

$$={\varphi }^{-4}\text{C}+2{\varphi }^{-2}\text{G}+\text{G}{\text{C}}^{-1}\text{G},$$

$$={\varphi }^{-4}{\text{M}}_0+2{\varphi }^{-2}{\text{M}}_1+{\text{M}}_2,$$

$${\text{Q}}_{3,\varphi }$$

$$={\text{K}}_{\varphi }{\text{C}}^{-1}{\text{K}}_{\varphi }{\text{C}}^{-1}{\text{K}}_{\varphi },$$

$$=({\varphi }^{-2}\text{C}+\text{G}){\text{C}}^{-1}({\varphi }^{-2}\text{C}+\text{G}){\text{C}}^{-1}({\varphi }^{-2}\text{C}+\text{G}),$$

$$={\varphi }^{-6}\text{C}+3{\varphi }^{-4}\text{G}+3{\varphi }^{-2}\text{G}{\text{C}}^{-1}\text{G}+\text{G}{\text{C}}^{-1}\text{G}{\text{C}}^{-1}\text{G},$$

Full text

\startlocaldefs\endlocaldefs

and

1 Introduction

The use of animal mounted sensors (bio-probes) to analyse population patterns has grown quickly in the last 25 years (Fedak, 2004; Ungar et al., 2005; Evans, Lea and Patterson, 2013), with tags attached to both marine and land based animals. These tags can be used to provide valuable information by collecting data on the environment where the animals live, particularly in regions that are difficult to observe otherwise. One example is given by the use of marine mammal tags to measure oceanographic data, such as water temperature, salinity and others. Fedak (2013) highlight the usefulness of such tags in profiling oceanographic data in polar regions, where data is typically difficult to obtain.

Although the methodology described in this paper is applicable to a range of different problems, we will focus here on data collected from CTD (Conductivity-Temperature-Depth) bio-probe tags attached to Elephant Seals in the Southern Indian ocean. These data were collected and made freely available as part of the MEOP (Marine Mammals Exploring the Oceans Pole to Pole) database (Roquet et al., 2013) and we utilise the South Indian ocean data subset described by Roquet et al. (2014). This data set was collected to supplement the Advanced Research and Global Observation Satellite (Argos) float and ship based measures of water masses in typically under sampled areas of the Southern oceans, an area which is drastically changing and needs to be further understood (Jacobs, 2006). While Northern oceans have been regularly sampled since the early 2000s using Argo profilers (Gould et al., 2004), utilising this method in the Southern oceans is typically complicated by the presence of sea ice.

The data consist of location coordinates (longitude and latitude which are only available when the animal is surfaced) and corresponding sea surface temperature (SST) measurements. The animal locations in our applied example are determined using Argos and are typically accurate within $\pm 5$ kilometers, whilst the temperatures are accurate within $\pm 0.03^{\circ}$C (Roquet et al., 2014). We use only the temperature data recorded at a depth of 6 meters, to represent the SSTs as closely as possible, and restrict ourselves to the region between -45 and -65 degrees latitude, between 60 and 120 degrees in longitude over the months of July to September 2012 and use tracks with 50 observations or more. The final data set consists of 9 separate tracks with 1630 observations in total, which can be seen in Figure 1.

Geostatistical models and methods (Diggle and Ribeiro, 2007) provide a natural framework to analyse these data. While standard spatial statistical methods consider sampling locations that are chosen independently from the response variable of interest, it is important to note that the movement of the tagged animals (and thus, the locations at which measurements are taken) may depend on the variable of interest (SST). Elephant Seals are likely to adjust their foraging due to the warming of ocean temperatures (McIntyre et al., 2011), because their prey tends to prefer cooler and deeper waters. It also appears that there are less successful forages when diving in warmer water (Guinet et al., 2014).

Preferential sampling (Diggle, Menezes and Su, 2010) refers to the situation where the process that determines the data locations and the spatial field of interest may depend on each other (as can be the case with data collected from animal-mounted tags). The effect of preferential sampling on subsequent inference can be two-fold. On the one hand, data collected in this way may tend to include a reduced range of the response variable. For example, if seals follow prey towards colder water, their water temperature measurements will tend to not include observations from warmer regions of their environment. On the other hand, the locations at which data are collected might carry useful information about the variable of interest. In the previous example, one may expect that regions that were not sampled had slightly warmer waters than those visited by the seals. The first issue refers to the fact that the sample of responses may not be a reliable representation over the area of interest, while the second one implies that it may be advantageous for the analysis to take into account the observed locations rather than conditioning on them, as it is done in standard spatial statistical methods.

The impact of preferential sampling on estimation and prediction has been discussed in various recent papers. Shaddick and Zidek (2014) highlighted the preferential nature by which air pollution monitoring sites were added and removed from a UK monitoring network from the 1960s until 2006. Particularly in the early years of the study, there is evidence that removed sites had a lower annual mean pollution reading than those which were added. Gelfand, Sahu and Holland (2012) showed the prediction effect of preferentially chosen ozone monitoring sites in California, whilst Diggle, Menezes and Su (2010) and Dinsdale and Salibian-Barrera (2018) illustrated the effect that preferential sampling of lead concentration in Galicia may have on the resulting predictions. Pati, Reich and Dunson (2011) studied the effect of preferential selection of monitoring sites measuring ozone levels in Eastern U.S.A. and more recently Conn, Thorson and Johnson (2017) showed that preferential sampling in animal population surveys may cause large biases in the animal density estimates, giving an example using aerial survey counts of bearded seals in the Eastern Bering sea.

Research on the issue of preferential sampling has so far mostly been restricted to stationary sampling locations such as monitoring stations. The main contribution of this manuscript is to extend this methodology to observations obtained by devices that move through the region of interest, such as those mounted on animals or people. In this paper the animal locations are modelled using a correlated random walk (Jonsen, Flemming and Myers, 2005; Johnson et al., 2008), where, to allow for the animal movement to depend on the variable of interest, we allow the drift function (which represents the direction preference at each location) to depend on the SST field. If this field can be assumed to be differentiable we include a term depending on its gradient to account for the animals’ possible preference for particular water temperatures. Finally, we use a behavioural state component as in Auger-Méthé et al. (2017) to allow for a combination of “momentum” and environmental preferences in the animal’s movement. Although we believe this movement model to be useful for the analysis of the CTD data, the focus of our work is on studying how to incorporate the possibility of preferential sampling to this type of models, rather than advocating for the superiority of a specific movement model over others.

It is interesting to note that the effect of preferential sampling in these “dynamic” spatial models (where sampling locations move through the domain) can be different from what happens when sampling locations are stationary. In the latter case, when there is no preferential sampling, observed locations are usually assumed to have a non-informative point pattern sampling distribution over the spatial domain. Thus, when preferential sampling is present it often results in distinct and potentially informative patterns in the monitoring locations (Diggle, Menezes and Su, 2010; Dinsdale and Salibian-Barrera, 2018). In contrast, even without preferential sampling, the locations visited by the animals in our study would typically not be expected to be evenly distributed over the area of interest (due to the dynamic nature of their movement). As a result, the information about SST contained in the locations the animals visited (or in the regions without observations) may be less apparent for these dynamical models than it is for stationary locations. Nevertheless, we show below that not taking into account preferential sampling may still negatively affect the quality and reliability of the estimated parameters and predictions.

To estimate the parameters of interest in the model, we utilise a likelihood approach based on a standard discretisation of the movement equations. The dependence between the sampling locations and the field of interest (SST) results in a likelihood that is computationally challenging to optimize. We follow Dinsdale and Salibian-Barrera (2018) in using the flexible Template Model Builder R package TMB (Kristensen et al., 2016) to deal with the computational complexity of the preferential sampling analysis.

In many related models discussed in the literature, the underlying field is only assumed to be continuous (but not differentiable). This choice appears to originate on mathematical and computational complexity considerations rather than on the specifics of the phenomenon being studied (Rue, Martino and Chopin, 2009; Lindgren, Rue and Lindström, 2011; Simpson, Lindgren and Rue, 2012). Since our model involves the gradient of the underlying SST field (see equation (3.7)), we will consider differentiable Gaussian random fields, which are relatively easy to use in the likelihood function when the latter is approximated using using TMB. Details can be found in Appendix A.

Our numerical experiments confirm that better predictions and parameter estimates can be obtained when the model appropriately reflects the potential presence of preferential sampling. Moreover, if the sampled locations are not preferentially chosen, the resulting predictions and parameter estimates are very close to those obtained with the standard model that conditions on the locations. Interestingly, our analysis also shows that the predicted SST in our motivating example obtained with and without a preferential sampling model are in fact different, which suggests that the former may be more reliable.

The rest of the paper is organized as follows. Section 2 briefly reviews preferential sampling spatial models and the methods used to perform inference and prediction based on them. Our model for the SST data is discussed in Section 3. The simulations results reported in Section 4 illustrate the benefits of accounting for preferential sampling when it may be present. Section 5 contains the analysis of Southern Indian ocean temperatures from tagged Elephant Seals, where we compare how accounting for possible preferential sampling might alter the conclusions reached by researchers. A final discussion can be found in Section 6.

2 The Preferential Sampling Problem

2.1 Standard Model Framework

Our Elephant Seal data consists of a response variable of interest (SST with possible measurement error, which we will denote by $Y$), a time stamp and the corresponding location in latitude and longitude coordinates $({\textbf{X}}\in{\mathbb{R}}^{2})$ which we assume contains no measurement error, see Section 2.4 for a discussion on this assumption. Since these tags measure water temperature at regular intervals, we consider the temperature measurement taken immediately before a location was obtained (which happens when the animal surfaces). We assume that the measurements taken from each trip are independent from each other.

To model these data we follow the geostatistical framework and notation of Diggle and Ribeiro (2007). More specifically, we assume that the data consist of a finite set of observations from a spatially continuous phenomenon $\left\{S({\textbf{x}}):{\textbf{x}}\in\mathcal{D}\subseteq\mathbb{R}^{2}\right\}$. In our application, $S({\textbf{x}})$ denotes the true SST at location x, and ${\cal D}$ is the region of the Southern Indian ocean. Note that $S(\cdot)$ does not vary over time, we discuss this further in Section 2.4. The model for the measurements $Y_{1}$, …, $Y_{n}$ obtained in one trip is:

[TABLE]

where $\mu\in\mathbb{R}$ is a constant mean parameter over $\mathcal{D}$ and ${\textbf{X}}_{i}$ is the measurement location. The $Z_{i}$ is included in the model above to account for measurement errors in the sea surface temperatures, and are assumed to be mutually independent random variables with mean 0 and so called “nugget variance” $\tau^{2}$.

We assume that the SST field $S$ is a Gaussian Process with mean 0 and Matérn covariance function given by

[TABLE]

where $r$ is the distance between two points and $K_{\kappa}$ is the modified Bessel function of the second kind. The scale (range) $\phi$ and smoothness $\kappa$ parameters control the rate of correlation decay over distance and $\sigma^{2}$ the marginal variance of the process (partial sill). Furthermore, $\kappa$ controls the smoothness of the process realizations. Since our model involves the gradient $S$ (see (3.7)), in what follows we will assume that the surface of SST is differentiable (in mean-square sense), which corresponds to $\kappa>1$ (Diggle and Ribeiro, 2007). More specifically, we take $\kappa=2$ which results in a spatial process $S$ that is mean-square differentiable (Banerjee, Gelfand and Sirmans, 2003; Banerjee and Gelfand, 2003).

The main goal of our analysis is to obtain predictions for the SST field $S$ over a grid of unobserved locations. Given estimates for the unknown field parameters in the model one can use standard spatial prediction techniques (e.g. kriging). However, as discussed in Diggle, Menezes and Su (2010) and Dinsdale and Salibian-Barrera (2018), when the locations X may depend on the field of interest $S$, one can obtain better predictions by including the information on $S$ contained in X. In the rest of this section we show how to construct an appropriate likelihood function that incorporates the possibility of preferential sampling, and how it can be maximized to obtain more accurate parameter estimates and predictions for the field $S$.

2.2 Maximum Likelihood Estimation

We use $[{\textbf{A}};\boldsymbol{\theta}]$ to denote the density or probability mass function of the random object A, which depends on a vector of parameters $\boldsymbol{\theta}$. In our case, $\boldsymbol{\theta}$ is the vector of all parameters in the model, which can be partitioned as $\boldsymbol{\theta}=(\boldsymbol{\theta}_{F}^{\top},\boldsymbol{\theta}_{L}^{\top})^{\top}$, where $\boldsymbol{\theta}_{F}$ are the parameters of the latent field process (e.g. $\boldsymbol{\theta}_{F}=(\mu,\tau,\kappa,\phi,\sigma^{2})^{\top}$ when using the model from (2.1) and (2.2)), and $\boldsymbol{\theta}_{L}$ are the parameters of the sampling location distribution. This last vector is divided into $\boldsymbol{\theta}_{L1}$, which are parameters that account for dependence between sampling locations and the latent field, and $\boldsymbol{\theta}_{L2}$, which do not relate to this dependence: $\boldsymbol{\theta}_{L}=(\boldsymbol{\theta}_{L1}^{\top},\boldsymbol{\theta}_{L2}^{\top})^{\top}$.

We then consider the likelihood function $L(\boldsymbol{\theta})$ based on the observed data Y and X:

[TABLE]

Typically, one has $\left[{\textbf{X}},{\textbf{Y}},S;\boldsymbol{\theta}\right]=\left[\left.{\textbf{Y}}\right|S,{\textbf{X}};\boldsymbol{\theta}_{F}\right]\left[\left.{\textbf{X}}\right|S;\boldsymbol{\theta}_{L}\right]\left[S;\boldsymbol{\theta}_{F}\right]$. Standard geostastical models assume that the process that selects the measurement locations X is independent from the response process $S$ (in symbols: $\left[\left.{\textbf{X}}\right|S;\boldsymbol{\theta}_{L}\right]=\left[{\textbf{X}};\boldsymbol{\theta}_{L2}\right]$), and hence $\left[{\textbf{X}},{\textbf{Y}},S;\boldsymbol{\theta}\right]=\left[\left.{\textbf{Y}}\right|S,{\textbf{X}};\boldsymbol{\theta}_{F}\right]\,\left[{\textbf{X}};\boldsymbol{\theta}_{L2}\right]\,\left[S;\boldsymbol{\theta}_{F}\right]$. In this case it follows that $L(\boldsymbol{\theta})=\left[{\textbf{X}};\boldsymbol{\theta}_{L2}\right]\left[\left.{\textbf{Y}}\right|{\textbf{X}};\boldsymbol{\theta}_{F}\right]$, and inference about $\boldsymbol{\theta}_{F}$ can be carried out conditionally on the observed locations X.

Preferential sampling refers to the situation where the observed locations X may depend on the unobserved process $S$. When $\left[\left.{\textbf{X}}\right|S;\boldsymbol{\theta}_{L}\right]\neq\left[{\textbf{X}};\boldsymbol{\theta}_{L2}\right]$ care must be taken when constructing the likelihood function in (2.3). In this case we cannot simply condition on the sampling locations, but should rather use the full likelihood function:

[TABLE]

2.3 Preferential Sampling Using Template Model Builder

Evaluation of the integral (2.4) is computationally challenging and hence optimising the function is a difficult problem.  Diggle, Menezes and Su (2010) proposed a Monte Carlo (MC) approximation to a discrete version of the integral, namely

[TABLE]

where S is a set of values of $S$. The exact locations used in the discretisation depends on the model used for the sampling locations and is discussed in more detail in Section 3.2. A direct MC approximation via simulated instances of the vector S is particularly inefficient, since many of the realisations of S may not be compatible with the observed measurements Y. Alternative representations of the likelihood function require sampling from the distribution of the discretized field S conditional on the observed locations and measurements (i.e.​ ${\textbf{S}}|{\textbf{Y}},{\textbf{X}}$), which is generally intractable (Dinsdale and Salibian-Barrera, 2018). Other alternatives to inference under preferential sampling have been proposed. For example, Pati, Reich and Dunson (2011) considered a Bayesian alternative and one could also use the R-INLA package (Rue, Martino and Chopin, 2009), which utilises integrated nested Laplace approximation.

A common assumption made in the literature about preferential sampling is that, conditional on the random field $S$, the sampling locations are static, often modelled via an inhomogeneous Poisson or similar process where the intensity function of $[{\textbf{X}}|S;\boldsymbol{\theta}_{L}]$ depends on $S$. However, in our application the sampling locations are obtained from a process continuously moving through the 2-dimensional domain. For this reason, we wish to use a flexible modelling framework in which we can evaluate the likelihood (2.4) efficiently for more complex forms of $[{\textbf{X}}|S;\boldsymbol{\theta}_{L}]$. Although R-INLA provides a highly efficient computational framework, it did not accommodate our relatively complex models for $[{\textbf{X}}|S;\boldsymbol{\theta}_{L}]$ with an underlying smooth process $S$ (e.g. a mean square differentiable SST surface).

We use the R package Template Model Builder (TMB) (Kristensen et al., 2016) to maximise (2.4) for the dynamic movements models discussed in the next Section. This package uses Automatic Differentiation (AD) (Griewank and Walther, 2008) of a Laplace Approximation to the likelihood to efficiently maximise it with respect to the full parameter vector $\boldsymbol{\theta}$. We define the joint negative log-likelihood function

[TABLE]

and TMB computes an approximation to $\int\exp[-f({\textbf{S}},\boldsymbol{\theta})]\,\mathrm{d}{\textbf{S}}$, which can be optimized numerically with respect to $\boldsymbol{\theta}$.

The dimension of the integral in (2.5) grows rapidly with the size of the grid that is used to discretize the field $S$. Important efficiencies can be obtained by using the stochastic partial differential equation (SPDE) approximations for Gaussian fields (Lindgren, Rue and Lindström, 2011), which are also exploited by R-INLA. More specifically, these SPDE approximations allow the use of sparse precision matrices to more efficiently evaluate the high dimensional integral in (2.5). Although the R-INLA package currently only allows continuous but not differentiable fields $S$ in (2.4), it is not difficult to extend the same approach for smoother random fields when using TMB to approximate (2.5). In particular, we work with mean-square differentiable random fields (Banerjee, Gelfand and Sirmans, 2003; Banerjee and Gelfand, 2003). Details can be found in Appendix A.

An important goal of this type of analyses is the prediction of SST on nearby locations that were not sampled. To construct predictions that take into account the preferential nature of the data, it is generally not sufficient to use kriging, even with parameter estimates obtained through a corrected likelihood function as in (2.4). Such an approach would effectively ignore the dependency between X and S (da Silva Ferreira et al., 2015; Dinsdale and Salibian-Barrera, 2018). Although the true predictive distribution of S is intractable in most cases when X and S are dependent, TMB provides point predictions and prediction variances from the estimated mode of $[{\textbf{S}}|{\textbf{Y}},{\textbf{X}};\boldsymbol{\theta}]$ at $\boldsymbol{\theta}=\boldsymbol{\theta}_{\text{opt}}$, where $\boldsymbol{\theta}_{\text{opt}}$ is the vector of parameters that maximises (2.5). Specifically, let

[TABLE]

then $\hat{{\textbf{S}}}(\boldsymbol{\theta}_{\text{opt}})$ is a predictor for the discretised version S of the process $S$ based on the preferential sampling model (2.4).

2.4 Assumptions

There are two key assumptions made so far in this section that should be noted. First is the assumption that the observed sampling locations X are the true positions of the animals. In reality there will be some degree of measurement error attached to the sampling locations, which depends on the tag type. Argos locations tend to include significant measurement error in comparison to Global Positioning System (GPS) locations which is common in land based tracking such as polar bears (Auger-Méthé et al., 2016) and various birds such as albatrosses (Weimerskirch et al., 2002) and gannets (Votier et al., 2010). However such systems are less appropriate for marine systems since GPS requires several seconds of exposure to obtain a location estimate (Dujon, Lindstrom and Hays, 2014). Recently Fastloc-GPS tags (http://www.wildtracker.com) have become more popular due to their improved accuracy compared to Argos and that these systems only require a fraction of a second to obtain a location estimate. When using such data Auger-Méthé et al. (2017) consider the location measurement error negligible enough to be ignored.

For the purpose of this paper, which is to emphasize the preferential sampling problem, we also decided to ignore measurement error in order to provide realistic but not overly complex models. However, this particular compromise between model complexity and computational efficiency by neglecting measurement error may have a serious impact on the analysis in certain situations. In particular, large unaccounted sampling location errors may lead to erroneous conclusions on the animal movement and consequently the preferential sampling effect. A possible strategy to incorporate measurement errors in the locations into our model is to include a latent state of true but unobservable locations (see for example Albertsen et al. (2015); Johnson et al. (2008)), which would add further latent states to the integral in (2.3). Note that in this case, care will be needed when considering the interplay between $S$ and X, since both objects will be unobservable.

The second assumption is that although the samples are taken at various time points, that SST depends only on location and not time. Therefore we can view the continuous SST field $S$ as constant over time. Further research in this area relaxing this assumption would be valuable. Enabling $S$ to vary over time, as it would in real life, would allow for analysis of data over longer periods of time with the model adapting to SST over various seasons and years. In this case one would consider the data to be of the form $Y_{t}=S({\textbf{X}}_{t},t)+Z_{t}$ so that $S(\cdot)$ is a function of both location and time. For the real data analysis in Section 5 we consider data across only 3 months to reduce the impact of a changing temperature field.

3 A Preferential Movement Model

To account for preferential sampling of ocean temperatures, we need to define a model for the location of marine mammals that takes into account possible relationships between movement velocity and ocean temperature. We will use discretised models as functions of the observed locations ${\textbf{X}}=({\textbf{X}}(t_{1}),\ldots,{\textbf{X}}(t_{n}))$, in which movement may also depend on previous locations. In these cases, we can write the location at time $t_{k+1}$ as

[TABLE]

where ${\textbf{X}}(t_{1:k})=({\textbf{X}}(t_{1}),\ldots,{\textbf{X}}(t_{k}))$ and $\boldsymbol{\theta}_{L}$ is now the vector of all movement parameters. The function $g_{\text{NP}}(\cdot)$ is some deterministic movement function where NP stands for non-preferential and $\epsilon(\cdot)$ is an error term.

Using the representation in (3.1), under preferential sampling we need to define a movement model of the form

[TABLE]

where the function $g_{\text{P}}$ is now a function of $S$, therefore enabling the movement model to depend on the temperature field.

3.1 A “Preferential” CRW Model for Marine Mammal Movement

We consider a model similar to the first-difference correlated random walk (DCRW) model (Jonsen, Flemming and Myers, 2005). We wish to include non-regular time intervals to account for possible irregularity in the surfacing of the marine mammals. Maintaining a constant time step through interpolation of the data, as discussed in Jonsen, Flemming and Myers (2005); McClintock et al. (2012); Hooten et al. (2017) among others, is not possible in our preferential sampling framework. This is due to the necessity of maintaining the link between the sampling locations (X) and the corresponding latent field measurements (Y). If we interpolate the trajectory it is not clear how we would obtain the corresponding $Y_{i}$ measurements at these interpolated locations, other than using a method such as kriging, which may dilute any preferential effect that was present in the original data. Another option may be to use thinning (Gurarie et al., 2017), however we wanted to avoid this in our application in this paper, due to the limited temporal resolution of the data to which we have access.

We term this model the “preferential correlated random walk” (PCRW) model and assume that the sampling locations ${\textbf{X}}(t_{1}),\ldots,{\textbf{X}}(t_{n})$ follow

[TABLE]

where ${\textbf{A}}_{k}$ denotes a standard bivariate normal random vector, $\boldsymbol{\Sigma}$ is a $2\times 2$ matrix that corresponds to the variance of the diffusion terms, and $t_{k}$ are the observation times.

To capture various movement patterns such as foraging and directed movement, rather than using discrete behavioural states (Morales et al., 2004; Breed et al., 2009; McClintock et al., 2012), we propose a continuous behavioural state system similar to Auger-Méthé et al. (2017); Breed et al. (2012). This method was chosen to obtain a differentiable likelihood function through the Laplace approximation outlined in Section 2.3, which would be invalidated with the more commonly used discrete states (Bolker et al., 2013). An alternative approach could be to estimate the movement parameters in TMB, then follow this with behavioural state estimation using the Viterbi algorithm (Whoriskey et al., 2017).

The drift function and behavioural states at the measured locations and times satisfy:

[TABLE]

where $\beta_{t_{k}}\in\mathbb{R}$ for all $t_{k}$, ${\textbf{v}}:\mathbb{R}^{2}\to\mathbb{R}^{2}$ represents the “velocity” of the animal, and $\boldsymbol{\phi}:\mathbb{R}^{2}\times\mathbb{R}^{2}\to\mathbb{R}^{2}$ can be thought of as the foraging movement function that depends on the location and the latent temperature field $S$, which models the possible preference of the animals for different water temperatures. The behavioural state function $f:\mathbb{R}\to[0,1]$ depends on $\beta_{t_{k}}$ and controls the auto-correlation of the movement at each time point. This ensures that when $f(\beta_{t_{k}})\approx 1$ then $\boldsymbol{\phi}({\textbf{X}}(t_{k}),S,\boldsymbol{\theta}_{L})$ becomes the expected drift direction, whilst when $f(\beta_{t_{k}})\approx 0$ movement tends in the direction of the current velocity ${\textbf{v}}({\textbf{X}}(t_{1:k}))$.

It is important to note that the velocity function may depend on more than just the previous sampling location. Consequently, like the DCRW, our PCRW model is not Markovian. Although one may consider including a latent velocity state (v) similar to the continuous time correlated random walk model (CTCRW) of Johnson et al. (2008), if this velocity was to depend on the locations X and also $S$, such a model becomes drastically more complicated. We discuss this in more detail in Appendix B. Alternative continuous-time models may also be adaptable to model preferential movement, however. For example, the correlated velocity model (CVM) (Gurarie and Ovaskainen, 2011) and functional movement models (FMMs) (Buderman et al., 2016; Hooten and Johnson, 2017).

We consider a non-latent velocity state approximation, taken to be

[TABLE]

and specify a behavioural function

[TABLE]

so that as $\beta_{t_{k}}$ increases, so does the influence of $\boldsymbol{\phi}$, whereas when $\beta_{t_{k}}$ decreases the current velocity ${\textbf{v}}({\textbf{X}}(t_{k}))$ becomes more of a factor in the movement. Therefore, our PCRW model can be written

[TABLE]

where $B_{k}$ are univariate standard normal random variables and $\sigma_{\beta}>0$ determines the evolution of the random states $\beta_{t_{k}}$. Note that the inclusion of the random $\beta$ states means we need to re-specify the preferential likelihood for our Laplace approximation, which we show in Appendix C. Initial values for ${\textbf{X}}(t_{0})$ and $\beta_{t_{0}}$ used in our simulation studies are described in Appendix D.1.

Depending on the application of the PCRW model, a variety of forms for the foraging function $\boldsymbol{\phi}$ may be used. In our case, to model the possible tendency of an animal to move towards particular water temperatures when searching for prey we propose the following:

[TABLE]

where $\alpha\in\mathbb{R}$, $S({\textbf{X}})$ is the value of the random field at location X and $\nabla S({\textbf{X}})$ is the gradient of $S$ at X. Although the parameter $\alpha$ above may appear to be unidentifiable, this is in fact not the case when you consider the full likelihood function which also includes the density functions $[{\textbf{Y}}|S,{\textbf{X}};\boldsymbol{\theta}_{F}]$ and $[S;\boldsymbol{\theta}_{F}]$.

The form of $\boldsymbol{\phi}$ in (3.7) defines the expected drift as descending (or ascending if $\alpha<0$) along the gradient of the SST field, with a velocity that depends both on the temperature at the present location and a scalar $\alpha$. This is somewhat similar to the varying motility surface used by Russell et al. (2018) in an stochastic differential equation (SDE) model to allow the magnitude of the velocity vector to depend on the location. More specifically, we can view the latent field $S$ as a scaled potential surface in which the gradient of $S$ directs the expected movement with a velocity that also depends on the value of $S$ at that location. Potential surfaces have previously been used to model the movement of animals including monk seals (Brillinger, Stewart and Littnan, 2008), elk (Brillinger et al., 2002; Preisler, Ager and Wisdom, 2013), and ants (Russell et al., 2018), with various estimation methods for potential surface SDE models compared by Gloaguen, Etienne and Le Corff (2018).

It should be noted that appropriate forms of $\boldsymbol{\phi}$ for particular species may require specialist input and we do not claim that the one shown in (3.7) is necessarily the best model for seals. The preferential sampling effect may vary over species, locations and possibly even individuals. However, the form in (3.7) may identify preferential movement and we use it as an example for integrating a foraging function into the preferential sampling framework. Furthermore, in application we may wish to adjust (3.7) to

[TABLE]

where $S^{*}=S+c$ for some constant $c\in\mathbb{R}$ specified by the user. We would do this to ensure that $S^{*}$ is the same sign across the domain, to prevent the switching of movement patterns when $S$ goes from negative to positive or vice versa.

Finally, using (2.1), (3.3)-(3.7) and a finite-differences approximation to the gradient of the field $S$ we construct a likelihood function that can be maximized numerically, as described in Section 2.3 above. Predicted values for the SST field can be obtained by using (2.7) at the vector of optimal model parameters.

3.2 Relationship Between the Movement Model and Likelihood Integral

In this section we discuss the form of the discrete grid S used to approximate the integral in (2.5). The locations at which we need to integrate over depends entirely on the distribution of ${\textbf{X}}|S$. Previous literature, in which the sampling locations were point patterns, suited a finely spaced lattice covering the entire domain (Dinsdale and Salibian-Barrera, 2018). This was because every point in the domain was a factor in determining the distribution of sampling locations.

In the case of a moving animal measuring a temperature field, this may or may not be required. Take for example, a model that assumes the animal might be knowingly moving towards distant points of attraction that are related to the latent field, for example a high prey region with low water temperature. In this case, locations far away from the animal might impact on the movement, hence requiring an approach similar to the point pattern integrals in which we require using a finely spaced lattice covering the entire, or majority, of the domain.

On the other hand, the preferential CRW model we are proposing in this paper assumes movement only depends on the animal’s immediate vicinity. This can be seen by observing that the only influence of $S$ on the movement is in (3.8), in which the current temperature and gradient of temperature field impacts movement. Hence, it would be more efficient to use a smaller grid for S, which contains only the sampling locations and those areas nearby which can be used to calculate the temperature gradient at sampled locations using a finite differences approach.

4 Simulation Experiments

In this section we discuss the results of a simulation study conducted to illustrate the effect of preferential sampling on the analysis of spatial data where sampling locations are moving through space, such as is the case in our SST data. The goal is to show to what extent incorporating preferential sampling in the model may improve the resulting SST predictions and parameter estimates. Furthermore, our results indicate that when preferential sampling is not present, there is almost no difference between using a model that incorporates preferential sampling and the usual geostatistical model that conditions on the locations. The non-preferential analysis can be found in Appendix D.3.

We generated 100 data sets, each of them consisting of up to 300 observations on 3 animal tracks following the Preferential-CRW movement model described in Section 3.1. Data were first generated on a fine time grid, and a subsample selected to form each of the 100 data sets. Details on the data generation process be found in Appendix D.1 and the corresponding code found in LABEL:suppA. With our simulation parameters we expect the tracks to oversample cooler regions. Figure 2 shows one simulated data set where the preferential sampling effect is apparent.

4.1 Parameter Estimates

Note that because of the way the data was generated (by subsampling trajectories created on a relatively fine time scale), the estimated movement parameters may not correspond to their nominal values used to create the data (Gurarie et al., 2017). Hence, we report here results for the estimates of the parameters of the spatial process.

Figure 3 shows the boxplots of the 100 estimated parameters for the spatial field process ($\boldsymbol{\theta}_{F}$) using each of the two likelihoods (standard and accounting for preferential sampling). The grey horizontal lines represent the true values. As expected, when the model does not account for preferential movement, the prevalence of lower temperatures in the sample introduces a negative bias to the estimates for the mean parameter ($\mu$). A similar pattern is observed also for the scale ($\phi$) and marginal variance ($\sigma^{2}$). The bias in the variance estimates is likely due to the tracks avoiding higher temperature regions and recording temperatures with a reduced range than they would otherwise. This may also explain the negative bias in the scale estimates. In contrast, using the Preferential-CRW model with TMB results in better parameter estimates. This is particularly noticeable for the estimates of $\mu$. Movement parameter estimates are discussed in Appendix D.2.

4.2 Prediction

We now turn our attention to the predictions for the underlying spatial field $S$. The first set of predictions are computed via kriging with parameter estimates obtained from the standard model that conditions on the observed locations, while the preferential sampling ones correspond to the estimated mode of $[{\textbf{S}}|{\textbf{X}},{\textbf{Y}};\boldsymbol{\theta}_{\text{opt}}]$, where $\boldsymbol{\theta}_{\text{opt}}$ are the parameter estimates from our preferential likelihood (see equation (2.7)).

Predictions were computed on a $26\times 26$ lattice for each of the $M=100$ data sets and we used two different measures of their quality. The root mean square prediction error (RMSPE) over the discrete domain is given by

[TABLE]

for each location $i=1,\ldots,N=26^{2}$ on the prediction grid. Here ${\textbf{S}}_{j,i}$ is the true value of field for the $j$-th simulated data set ${\textbf{S}}_{j}$ at the $i$th prediction location (${\textbf{S}}_{j}({\textbf{x}}_{i})$) and $\hat{{\textbf{S}}}_{j,i}$ is the corresponding predicted value.

To compare the resulting predictions whilst accounting for their variances, we used Ignorance Scores (Roulston and Smith, 2002), which are given by $\operatorname{IGN}(x)=-\log(p(x))$, where $p$ is the predictive density and $x$ is the target forecast which would be the true SST at that location (Siegert, Ferro and Stephenson, 2014; Gneiting and Raftery, 2007). For each of the $j=1,\ldots,M=100$ preferential simulations we calculated the Mean Ignorance Score (MIGN) of our latent field predictions as

[TABLE]

where $\hat{\sigma}_{j,i}^{2}$ is the prediction variance of $\hat{{\textbf{S}}}_{j,i}$. This measure gives an indication of model performance for each simulation averaged over the entire domain. We also calculated location-specific Ignorance Scores, averaging the IGN of each location over the 100 samples, we call them Location Ignorance Scores (LIGN):

[TABLE]

The LIGN gives an assessment of model prediction across each region of the domain.

To compare the “standard” predictions with the “preferential sampling” ones we computed the corresponding differences of the above 3 measures:

[TABLE]

where NP and P indicate the values of the scoring functions for the non-preferential (standard) and preferential models, respectively.

The first panel in Figure 4 shows the values of $\operatorname{RMSPE}^{\text{Diff}}_{i}$ at each location, colouring the areas in blue for which this measure was negative, which correspond to regions where the RMSPE’s for preferential sampling predictions were better. Nearly all regions were predicted more accurately on average using the preferential sampling model, with only a small number of positive (red) locations.

The boxplot of the $M=100$ differences in MIGN is displayed in the second panel of Figure 4. Negative values of this difference imply that the P model had a lower average Ignorance Score (over the region) for that particular simulation run. Again we see that the preferential sampling model performs better than standard methods across the majority of simulated data sets. This conclusion is further supported by the third panel in Figure 4 which shows $\operatorname{LIGN}^{\text{Diff}}_{i}$ for each lattice prediction location. All locations had a smaller Location Ignorance Score on average when using the preferential sampling model, which interestingly includes the locations where point predictions were slightly inferior. This suggests that, even in this case, the prediction variances were more reasonable.

5 Real Data Example

In this section we analyse Southern Indian ocean temperatures at a depth of 6 meters, which we will be calling sea surface temperatures (SSTs), using CTD (Conductivity-Temperature-Depth) sensor data from tags attached to Elephant Seals. These data were collected and made freely available as part of the MEOP (Marine Mammals Exploring the Oceans Pole to Pole) database (Roquet et al., 2013, 2014) and is described in Section 1 and shown in Figure 1.

Since these data are relatively near to the South Pole, there is a large change in true distance between latitudinal and longitudinal degrees over the domain. For example, we can see in Figure 1 that the change in latitude is not constant in distance between $-45^{\circ}$ and $-65^{\circ}$. This is due to the difficulties of representing locations from a 3-dimensional sphere on a 2-dimensional surface, which is of particular importance near the poles. One option is to use the Haversine formula to calculate the great-circle distance between two points in the latitude/longitude space (Robusto, 1957). However, as is discussed by Jeong and Jun (2015); Gneiting (2013) among others, using the Matérn covariance along with great circle distances requires the use of a non mean-square differentiable processes ($0<\kappa\leq 0.5$), since the Matérn class is not isotropic otherwise.

To more accurately represent distance over our domain but remaining in the Euclidean space, we transformed the sampling locations. Our transformation used a scaled version of the Universal Transverse Mercator (UTM) projection (zone 43) which can be seen in Figure 5. This transformation provides us with a more accurate representation of distance than using raw latitudes and longitudes, whilst remaining in the 2-dimensional Euclidean space and retaining the isotropic properties of our correlation model. The values of our scale have no real-world interpretation, other than to provide more realistic scaled distances between sampled locations.

We compare both the parameter estimates and corresponding field predictions using the PCRW model described in Section 3.1, with those obtained from a standard geostatistical model in which the sampling locations were considered independent of the temperature field. In order to control the computational complexity of the analysis and also to explore the sampling distribution of the field parameter estimators, we use 50 subsamples from the data. For each of these 50 replications we randomly sampled 40 observations from each of the 9 tracks and estimated the parameters on the sub-sampled data. The resulting parameter estimates are displayed in Figure 6.

We note that the difference between the two sets of estimates is very small in the scale and variance parameter estimates. However, there was an increase in the mean parameter estimate using the PCRW model. The consistently positive $\hat{\alpha}$ estimates suggest that there was a tendency of the animals to move towards cooler regions according to our model, which explains the increased mean parameter estimates.

For predictive assessment of the ocean temperatures obtained from marine mammal tags we compare our model predictions with monthly average temperature fields obtained via the Simple Ocean Data Assimilation ocean/sea ice reanalysis (SODA) (Carton and Giese, 2008), specifically the SODA version 3 (SODA3) reanalysis (Carton et al., 2018), which uses all temperature and salinity profiles from the World Ocean Database. This data is available for depths between 5 and 5000 meters below the surface at a spatial resolution of 0.5 degrees latitude and longitude. Data is available from January 1980 to present.

The first panel of Figure 7 shows the SODA3 monthly average ocean temperature field for August for a depth of 5 meters. We compare our predictions with these monthly averages at $N=461$ locations which were chosen as those on the original $26\times 26$ lattice points which were close enough to sampling locations to obtain kriging predictions that did not simply tend to the mean trend parameter $\hat{\mu}$. The locations on the lattice which were not used are shown in grey. To compare point predictions we consider the quantiles of the difference between the two prediction methods on each of the data sets at each prediction location. In other words we consider the quantiles of each coordinate of the vector ${\textbf{D}}=(\mathbb{\textit{D}}_{1},\ldots,\mathbb{\textit{D}}_{N})$ such that

[TABLE]

for $i=1,\ldots,461$, where $\hat{S}^{\textit{P}}_{j,i}$ is the prediction at location $i$ in simulation $j$ using the preferential model and $\hat{S}^{\textit{NP}}_{j,i}$ the equivalent using the non-preferential model. We plot the $50\%$ quantile of D in the second panel of Figure 7. As we would expect with the positive $\alpha$ estimates, most areas further from sampling locations had an increased SST prediction using the PCRW model. Interestingly however, the PCRW model actually tended to decrease SST prediction in certain areas of the West and North-East predictive region. We discuss the magnitude of these prediction differences further in Appendix E.

To assess the prediction accuracy compared to the SODA3 data set we consider ignorance scores and RMSPE’s in the same manner as in Section 4.2. We plot the difference in RMSPE, MIGN and LIGN defined in (4.4) in Figure 8. The first panel shows that the PCRW model tended to reduce RMSPE in comparison to non-preferential prediction in the Northern regions, whilst underperforming in areas of the South. However, the LIGN in the third panel shows improved ignorance scores in general using the preferential model. An indication of superior prediction performance is through the second panel of Figure 8 which shows the majority of MIGN’s across the 50 subsamples were negative, and indication of better predictions using the PCRW model.

To summarise, the PCRW model applied using TMB showed differences in mean parameter estimation when compared to the naive method and identified possible tendency of sampler movement towards the cooler regions. This translated to many areas of increased SST prediction when using the preferential model in spatial prediction, but not across the entire spatial region. Small changes in SST prediction as shown by our model may be of increased importance when dealing with complex systems as ocean temperatures. Although this simulation used a simplistic relationship between the sampler and the process being modelled, our results suggest that developing more realistic models for this type of data can improve the resulting statistical inference about our ocean climates.

6 Discussion

We have shown the damaging effect that preferential sampling may have on statistical inferences based on spatial models where monitoring locations are not stationary. The simulation experiments reported in Section 4 illustrate how predictions may be improved when accounting for the preferential nature of movement in the sampling model. This is evident in the parameter estimates and also in the predictive performance. Here the combination of corrected parameters and a predictive distribution that accounts for the relationship between the sampling locations and the spatial field of interest improves notably upon the simple extrapolation of kriging. When we compare the results of the standard method (that conditions on the observed locations) and the Preferential-CRW model in our motivating example we observe increased mean parameter estimates for the SST field, whilst the predicted SST fields show consistent differences.

These analyses highlight the importance of expanding preferential sampling methods to the case of non-stationary sampling locations, which is becoming a more prevalent situation. We show how implementing a Laplace approximation to the likelihood function via the R package TMB allows for flexible movement model specification. This method can expand beyond the typical point process models to sampling processes derived from movement models which may depend on our field of interest $S$. This is just the first step in incorporating preferential sampling into the statistical analyses of tagged marine mammals. Considering the observed measurement locations to contain some measurement error by assuming there is an latent true location state would help to account for the uncertainty in the Argos location estimates, whilst there are many other applications such as land animals and other non-stationary sampling processes as the next step in expanding our applications and methodology. Furthermore, a natural extension of this research is to consider data at any depth. Although SST analysis is a useful first step, an important application of these methods would be to aid in mapping the water masses at depths unobservable by satellites. Theoretically little would change in comparison to the SST data we consider here in analyses at specific depths, however the location accuracy underwater is reduced. In this case, location estimates could inferred based on accelerometer data using methods such as the dead reckoning algorithm (Wilson and Wilson, 1988; Liu et al., 2015) and ocean temperature at specific depths analysed to build a water temperature profile.

Appendix A Using Smooth Random Fields with TMB

Many geostatistical models proposed to analyse spatial data assume that the underlying random field is continuous (but not necessarily differentiable). In what follows we show how TMB can be used to efficiently employ mean-square differentiable Gaussian fields (GFs) in our models, using the approach by Lindgren, Rue and Lindström (2011) who show how Gaussian Markov random field (GMRF) representations of GFs with Matérn covariance (as defined in (2.2)), can be constructed through the solution to a stochastic partial differential equation (SPDE) when driven by white noise.

Methods that utilise GMRF representations of GFs require constructing a sparse precision matrix Q that closely represents the covariance of the GF. We follow directly after (9) in Lindgren, Rue and Lindström (2011) but note the change in notation from $\kappa$ in the original work to $\phi$ here and that in our 2-dimensional application $\alpha=\kappa+1$ where $\kappa$ is the smoothness parameter in (2.2).

Lindgren, Rue and Lindström (2011) write the precision matrix ${\textbf{Q}}_{\alpha,\phi}$ as a combination of $m\times m$ matrices ${\textbf{C}},{\textbf{G}}$ and ${\textbf{K}}_{\phi}$,

[TABLE]

where $m$ is the number of vertices in the triangulation of the domain.

This combination depends on $\alpha$ but can be calculated recursively as

[TABLE]

Unfortunately ${\textbf{C}}^{-1}$ is dense, but Lindgren, Rue and Lindström (2011) show that C can be replaced by the diagonal matrix $\tilde{{\textbf{C}}}$, where $\tilde{C}_{i,i}=\langle\psi_{i},1\rangle$. Hence $\tilde{{\textbf{C}}}$ is sparse and the resulting precision matrix ${\textbf{Q}}_{\alpha,\phi}$ also sparse. Further details can be found in Lindgren, Rue and Lindström (2011); Simpson, Lindgren and Rue (2012); Lindgren and Rue (2015).

A.1 Implementation of $\kappa=2$ in TMB

R-INLA only implements the SPDE approximation for $\alpha\in(0,2]$. At the largest smoothness of $\alpha=2$ ($\kappa=1$), the corresponding Matérn field is mean square continuous but not mean square differentiable. In some applications this might not suit the problem at hand (for example, when derivatives of the spatial field may be part of the model). Below we extend the current approximations to the case $\alpha=3$ ($\kappa=2$) and show how this can be implemented in TMB via the built-in compatibility with R-INLA.

To begin, note that for $\alpha=2$ we can expand the formula for ${\textbf{Q}}_{2,\phi}$ above as follows:

[TABLE]

where ${\textbf{M}}_{0}={\textbf{C}},{\textbf{M}}_{1}={\textbf{G}}$ and ${\textbf{M}}_{2}={\textbf{G}}{\textbf{C}}^{-1}{\textbf{G}}$. Note that the matrices ${\textbf{M}}_{1},{\textbf{M}}_{2}$ and ${\textbf{M}}_{3}$ do not depend on $\phi$ or $\kappa$ and can be computed with the R-INLA function inla.spde2.matern. Now we expand ${\textbf{Q}}_{3,\phi}$:

[TABLE]

Hence we can use ${\textbf{M}}_{0},{\textbf{M}}_{1}$ and ${\textbf{M}}_{2}$ provided by inla.spde2.matern to construct ${\textbf{Q}}_{3,\phi}$ for a solution to the SPDE approximation for $\alpha=3$ ($\kappa=2$). We simply need to take a new combination of these matrices which in TMB can be done easily within the function that computes the likelihood.

Appendix B Preferential-CRW Models with a Field-Dependent Velocity Term

If one considers a velocity term v defined in (3.4) that may depend on both the locations X and the underlying field of interest $S$, the likelihood quickly becomes intractable. Considering a latent velocity state in the model, the full likelihood is

[TABLE]

Notice that this factorisation includes the term $[{\textbf{v}}|{\textbf{S}};\boldsymbol{\theta}_{L}]$. If the velocity was to only depend on the latent field S and not on the location X, then we could evaluate the likelihood. However, this makes little biological sense in practice. Preferential movement induced by a latent velocity would require this velocity to depend on the field S at the current location X of the animal. In other words the velocity would be dependent not just on the latent field but the current location within this field. However, this means that we would only have a tractable form for $[{\textbf{v}}|{\textbf{S}},{\textbf{X}};\boldsymbol{\theta}_{L}]$ and would therefore need to evaluate

[TABLE]

which is difficult to compute due to the complexity of $[{\textbf{X}}|{\textbf{S}};\boldsymbol{\theta}_{L}]$. The factorisation of the form in (B.1) appears intractable and therefore the Preferential-CRW model in which velocities depend on location and the latent field become difficult to implement using likelihood methods when preferential sampling may be present.

Appendix C Preferential CRW Likelihood

In Section 2.3 we outlined how to evaluate the preferential model using a Laplace approximation to the likelihood function. In previous cases this factorisation was of the form $[{\textbf{X}},{\textbf{Y}};\boldsymbol{\theta}]=\int[{\textbf{X}},{\textbf{Y}},{\textbf{S}};\boldsymbol{\theta}]\mathrm{d}{\textbf{S}}$, however in the case of the Preferential-CRW model we have a second latent vector in the behavioural states $\boldsymbol{\beta}=(\beta_{t_{1}},\ldots,\beta_{t_{n}})$. Therefore we need to re-specify the full likelihood for our Laplace approximation routine.

[TABLE]

Notice that in the Preferential-CRW model we have $[\boldsymbol{\beta}|{\textbf{S}};\boldsymbol{\theta}_{L}]=[\boldsymbol{\beta};\boldsymbol{\theta}_{L2}]$ so therefore in the Laplace approximation implemented using TMB we will need to redefine the joint negative log-likelihood as

[TABLE]

Appendix D Simulation Details

D.1 Data Generation

To generate each track, we initialised a starting location chosen uniformly at random over the 2-dimensional domain $[-150,150]\times[-150,150]$, then simulated 360 observations with the time between consecutive observations following an exponential distribution with rate parameter $\lambda=10$. The first 60 positions of each animal were considered a burn-in period and discarded, resulting in 300 remaining observations. Finally the track was thinned by taking every 3rd observation, to retain a final track of 100 observations.

For each data set we simulated a random field $S$ over the domain and discretised it on a $51\times 51$ grid. We used a Matérn covariance function as in (2.2) with smoothness parameter $\kappa=2$, scale $\phi=25$, marginal variance $\sigma^{2}=1.5$ and a constant mean $\mu=5$. Since we know the true field in generating the tracks, the gradient of the field used to direct movement was approximated using finite differences from points not necessarily on the grid. For the movement model we set $\sigma_{\beta}=0.1,\boldsymbol{\Sigma}=3\,{\textbf{I}}_{2}$, where ${\textbf{I}}_{2}$ denotes the $2\times 2$ identity matrix, $\alpha=100$ and initiated the behavioural states at $\beta_{0}=-1.5$. From (3.6) we see that this choice initialises primarily directed movement ($f(\beta_{0})=0.18$), but with a slight influence of the foraging (preferential) function $\phi(\cdot)$. We assume that $\tau^{2}$ is a known parameter due to the assumption that the measuring device will have a known sampling error and as it is commonly done in the literature, we also assume that the smoothness parameter $\kappa$ is known.

D.2 Movement Parameter Estimates

For the preferentially generated data, Figure 9 shows that the movement parameters estimates, with the positive $\hat{\alpha}$ estimates accounting for the tendency of the sampler to avoid warmer warmers, which explains the correction to the mean parameter estimates observed in Figure 3. Although these estimates cannot in general be compared with the values used to generate the fine-time-scale data (Gurarie et al., 2017), the boxplots in Figure 9 show that the estimators have a reasonable sampling distribution which seem to be unimodal and symmetric around their means.

D.3 Results for Non-Preferential Data

We also generated 100 data sets with non-preferentially sampled data (i.e. setting $\alpha=0$ in (3.7)), and estimated the parameters of the spatial process using the Preferential-CRW likelihood and the standard one (that conditions on the observed locations). The results in Figure 10 show that, as expected, there is no practical difference between the parameter estimates obtained using either likelihood when no preferential sampling is present.

We also performed the same comparisons as above over the 100 non-preferentially sampled data sets to verify that in this case, as expected, there was little qualitative difference between the predictions obtained using either method. The results are displayed in Figure 11 and show that, although RMSPE’s are often larger with the preferential model, the difference between the two models is minor compared to the preferentially sampled data. This can be observed by comparing the scale of the differences in the plots, which are considerably smaller in the non-preferentially sampled analysis.

Appendix E Magnitude of Prediction Differences

We have tried different ways to assess the significance of the prediction differences in Figure 7. Panel (b) shows the median difference, while the first (25%) and 3rd (75%) quartiles of the prediction differences (over the 50 subsamples) are shown in panels (a) and (b) of Figure 12. Note that in this plot, the white areas are locations of no-significance (negative valued regions in the 25% quantile and positive regions in the 75% quantile). In other words, the coloured regions of the 25% quantile plot are regions in which prediction difference were consistently positive, and coloured regions in the 75% quantile plot consistently negative.

As it may be expected, in large parts of the domain the differences in predictions do not show a clear trend either way (white locations), but some conclusions could be drawn from these figures. The coloured regions in panel (a) are zones were most of the predictions for the preferential sampling model were higher than from the standard one. We see that they are in regions outside the area covered by the data, which is consistent with the positive estimates for $\alpha$ (see Figure 8 on page 24). Similarly, the coloured regions in panel (b) are zones where the predictions were mostly lower. Note that these areas are almost all among the observed locations, which is also consistent with our conclusions above.

{supplement}

[id=suppA] \snameSupplement A \stitleSimulation Code \slink[url]https://github.com/msalibian/PreferentialMovement.git \sdescriptionCode used for the simulations in this paper, with an example shown in the README.