11institutetext: Federal Office for Information Security, Bonn, Germany
11email: ¡first name¿.¡last<EMAIL_ADDRESS>22institutetext: LatticeFlow,
Zurich, Switzerland
22email: ¡first<EMAIL_ADDRESS>
# Robustness testing of AI systems: A case study for traffic sign recognition

Christian Berghoff Authors are listed in alphabetical order.11 Pavol Bielik 22
Matthias Neu 11 Petar Tsankov 22 Arndt von Twickel 11

###### Abstract

In the last years, AI systems, in particular neural networks, have seen a
tremendous increase in performance, and they are now used in a broad range of
applications. Unlike classical symbolic AI systems, neural networks are
trained using large data sets and their inner structure containing possibly
billions of parameters does not lend itself to human interpretation. As a
consequence, it is so far not feasible to provide broad guarantees for the
correct behaviour of neural networks during operation if they process input
data that significantly differ from those seen during training. However, many
applications of AI systems are security- or safety-critical, and hence require
obtaining statements on the robustness of the systems when facing unexpected
events, whether they occur naturally or are induced by an attacker in a
targeted way. As a step towards developing robust AI systems for such
applications, this paper presents how the robustness of AI systems can be
practically examined and which methods and metrics can be used to do so. The
robustness testing methodology is described and analysed for the example use
case of traffic sign recognition in autonomous driving.

###### Keywords:

Neural networks Robustness Autonomous driving

##  1 Introduction

AI systems based on neural networks have tremendously increased their
performance over the course of the last years and are now used in a plethora
of very diverse areas. The improved performance is in particular based on a
strong increase in available computing power and data, but also on theoretic
advances in the area of AI in general and in the design of the networks in
particular. Current applications of neural networks range from predictions
based on structured data to natural language processing and computer vision
tasks. In the latter domain, AI systems are for instance used in biometrics
for identifying people, and as an important building block of (partially)
autonomous driving for processing and analysing sensor data.

However, unlike the classical systems of symbolic AI, neural networks can, in
most cases, not be manually developed by experts, but are instead trained on
the basis of large data sets so as to provide the desired functionality. As a
result of this training procedure, several problems arise when evaluating the
behaviour of AI systems after they have been deployed to live operation. On
the one hand, the systems are strongly dependent on the quality and
representativeness of training data, whose amount is necessarily limited and
which, in most realistic application scenarios, cannot cover all possible
inputs [2]. Furthermore, it is becoming increasingly evident that even if we
could easily obtain much larger data sets, this alone will not be sufficient
to ensure the safe operation of these models [12, 7].

As a result, significant effort has been recently devoted to developing
techniques capable of providing formal robustness guarantees [14, 15, 11, 19].
Such guarantees are of paramount importance, especially when considering use
cases where severe damages may potentially arise. As an example, malfunction
of an AI system used in autonomous driving can lead to significant material
loss and in extreme cases to fatalities. Unfortunately, such techniques are
currently not applicable to state-of-the-art computer vision models due to
their limited scalability and the limited type of robustness properties they
support.

This article focuses on the complementary approach of developing a methodology
for testing the robustness of AI systems. This allows the developers of the AI
systems to assess the robustness of any state-of-the-art system in a
principled way and to identify failure modes that need to be explored in more
detail and remedied before the system can be deployed. While this empirical
approach cannot provide absolute guarantees, the amount of robustness testing
carried out can be adjusted to achieve the required confidence level about the
model correctness. It allows systematically taking into account prior human
knowledge by defining task-specific robustness properties, and thus aids in
making the development of robust neural networks more efficient. In what
follows, we will briefly describe the main components of our work and
instantiate the approach to the concrete use case of traffic sign recognition.

##  2 Related Work

The robustness of AI systems has been extensively studied in the literature,
but the large majority of research has focused on their robustness to attacks
specially designed to break the system, an area commonly referred to as
adversarial machine learning (AML). Starting with the first publication [21]
applying adversarial attacks (more generally known as evasion attacks years in
advance [6]) to neural networks, many approaches to mitigating the
vulnerability of neural networks to these attacks have been developed and many
of them have been broken subsequently [4]. Another similar, albeit somewhat
less popular subject of research are poisoning attacks [4]. One common feature
of research in AML is the lack of a standardised and agreed-upon framework for
evaluating robustness, although some progress has been made towards this end
[5].

When considering research on the robustness of AI systems to perturbations
that may occur naturally during operation and should thus be much easier to
cope with than malicious attacks, results are much more scarce. Some
publications deal with the robustness to natural perturbations and propose
measures for increasing it [17, 18, 16], but they do not set out to perform a
fine-grained and systematic assessment of the phenomenon. Probably the work
most closely related to ours is [23, 22]. Similar to our work, the authors
define a set of robustness properties and assess the model robustness against
them. However, while the main goal of their work is to introduce a new
benchmark, the goal of our work is a thorough assessment of state-of-the-art
models. As a result, (i) we assess the robustness of significantly better
models that achieve 99% standard accuracy (compared to models with at most 90%
accuracy used in [23]), (ii) we show the need of considering robustness
properties jointly, rather than in isolation, and (iii) we focus on an
iterative methodology to assess and improve model robustness (which includes
eliciting new properties and discovering model failure modes), rather than
proposing a fixed data set.

One line of work that can yield results for this research question is that of
formal verification [14, 15, 11, 19]. The main limitations of formal
verification are that the methods developed so far do not scale to large
neural networks used in practice and that the type of robustness properties
that can be efficiently encoded is limited. For example, the majority of the
current works consider only norm-constrained pixel perturbations, and
verifying even simple geometric transformations such as rotations is highly
non-trivial [1]. Further, formal verification approaches suffer from the issue
of incomplete specification, since the number of robustness properties they
can encode is limited (mostly pixel perturbations). As a result, one has to be
careful when interpreting their results and make sure not to incur a false
sense of security by performing a thorough assessment of a wide range of the
robustness properties, not only a subset of those for which formal
verification can provide strong guarantees.

Figure 1: Robustness assessment of neural networks. The testing approach takes
as input a neural network and its test data set, along with robustness
properties that capture important environmental conditions (e.g., camera
position, sign stickers), and assesses the neural network against these
together with its common failure modes.

##  3 Approach

This section presents details about our case study on the robustness testing
of traffic sign classifiers. Figure 1 provides an overview of the approach.
Robustness testing is performed on a neural network based on a test data set
and a specification of certain properties. A metric is used to compute the
robustness scores. Apart from identifying common failure modes of the model,
the results can be used to understand how to enhance the data set and to
define additional custom properties, which improve both the neural network and
the testing process. More detailed information on the individual components
for the case study on traffic sign classifiers is described in sections
3.2–3.4. The general methodology depicted in Figure 1 can be applied to other
use cases by exchanging and adjusting the components in a suitable way.

###  3.1 Models

The following models, which were trained by the authors, were used for the
robustness testing:

  * •

The first neural network, called pre-trained, has 99.0% accuracy. It uses the
pre-trained Inception-v3 model [20], which is trained on the ImageNet data set
[8] and fine-tuned for the GTSRB data set [13].

  * •

The second one, called self-trained, has 97.4% accuracy. It uses an
architecture based on Inception-v3 but with reduced size and without pre-
training.

Both networks were trained on the GTSRB data set with a data augmentation
policy that applies the following random transformations: random cropping of a
portion of the original image with a side length between 0.6 and 1.0 of its
original length, rotation by an angle within $-15$ and $+15$ degrees, colour
changes – brightness, contrast, saturation, and hue – by a factor of up to
0.1, and random change to a gray-scale format with a probability of 0.1. The
transformed images are then scaled to the neural networks’ expected resolution
– $32\times 32$ pixels for the self-trained network and $299\times 299$ pixels
for the pre-trained network.

###  3.2 Data Set

We use the German Traffic Sign Recognition Benchmark data set (GTSRB) [13], as
illustrated in Figure 2. This data set consists of 39,209 colour images used
for training and 12,630 images for testing, each assigned to one of 43
classes. The data set images have different resolutions, with height and width
ranging between 25–266 and 25–232 pixels, respectively.

Figure 2: Example traffic signs drawn from the GTSRB data set [13], one for
each class.

###  3.3 Robustness properties

The robustness properties used were chosen according to the ambient conditions
of the use case. Traffic sign recognition, and autonomous driving in
particular, is a very challenging but also relevant use case in this respect.
Since autonomous driving technologies are deployed in the real world and on
the roads, ambient conditions can change dramatically depending on the time of
day, season or weather. Traffic signs may also to some extent be rotated,
occluded and damaged. The sensors used to capture images of the traffic signs
may themselves be degraded from long use, damaged or soiled, and the quality
of the images captured can strongly vary depending on the specific lighting
conditions. This situation is in (stark) contrast to other sensitive use cases
such as medical image classification or biometrics used in border control,
where ambient conditions can be standardised to a (much) larger extent, and
the conditions after deployment can be made reasonably close to the ones
considered during training.

The following list provides a brief overview of the robustness properties
used. For the complete list of properties along with the used robustness
bounds, we refer the reader to our detailed technical report [3].

  * •

Image noise includes Gaussian noise, uniform noise and impulse noise, each of
them modelling different deviations from optimal conditions in image capture
or pre-processing.

  * •

Pixel perturbations define the robustness over individual pixels, each of
which is allowed to change independently from other pixels subject to a total
perturbation budget measured in the $L_{0}$ or $L_{\infty}$ norm. The $L_{0}$
norm allows modifying a limited amount of pixels to any extent, while the
$L_{\infty}$ norm allows changing all pixels by less than a certain threshold.

  * •

Geometric transformations model different orientations or positions of the
traffic signs as well as faults in image acquisition and pre-processing. The
perturbations considered are rotation, translation, scaling, shearing,
blurring, sharpening and flipping (flipping was only considered on images
whose class was not changed by applying it).

  * •

Colour transformations model changes in lighting conditions or colour post-
processing performed on the images. The colour transformations considered are
brightness, contrast, saturation, hue, gray-scale and colour-depth.

###  3.4 Metric

The metric used to assess the robustness of a model to a specified property is
the robustness score, defined as the fraction of robust samples in the data
set. A sample is called robust to a given property if the network produces the
correct label for all transformations captured by the property. For example,
suppose the model achieves 86.5% robustness to rotations up to 15 degrees.
This indicates that for 13.5% of the input samples, a rotation of the original
image within 15 degrees was found which causes the model to predict the wrong
label. Robustness is defined only over input images for which the model
already predicts the correct label as all incorrect inputs are trivially non-
robust.

SELF-TRAINEDPRE-TRAINED95.5%98.8%Gaussian Noise RobustnessINPUTMORE NOISE
$\rightarrow$

SELF-TRAINEDPRE-TRAINED61.3%96.8%Pixel $L_{0}$ RobustnessINPUTHIGHER NUMBER OF
PERTURBED PIXELS $\rightarrow$

SELF-TRAINEDPRE-TRAINED53.5%2.1%Pixel $L_{\infty}$ RobustnessINPUTHIGHER
MAXIMUM CHANGE TO EACH PIXEL $\rightarrow$

SELF-TRAINEDPRE-TRAINED97.6%99.5%Rotation Robustness$\leftarrow$ NEGATIVE
ROTATIONPOSITIVE ROTATION $\rightarrow$

SELF-TRAINEDPRE-TRAINED89.0%96.9%Scaling Robustness$\leftarrow$ SMALLER
SCALELARGER SCALE $\rightarrow$

SELF-TRAINEDPRE-TRAINED99.6%99.9%Hue Robustness$\leftarrow$ LOWER HUEHIGHER
HUE $\rightarrow$

Figure 3: Overview of the models’ robustness to different transformations.

##  4 Results

This section provides a summary of the results obtained using the data set,
models, robustness properties and metric from section 3. For comprehensive
results, as well as exact definitions of the tested properties, the reader is
referred to [3].

###  4.1 Basic Robustness Tests

Figure 3 shows the robustness of the models with respect to the robustness
properties as outlined in section 3.3. Each part of Figure 3 is organised as
follows: The robustness score of both the self-trained and pre-trained model
to the respective property is given, and the result of applying different
intensities of the property on the original image is provided to show the
visual impact for easy human inspection. Figure 3 shows that the robustness
scores significantly differ, both between models and properties. Whereas both
models are relatively robust (score $\geq 93.9\%$) to Gaussian noise,
rotations and colour transformations such as changes in brightness (Figure 4)
and hue, the robustness of the self-trained model to rescaling the input is
lower and its robustness to pixel $L_{0}$ perturbations even much lower, with
more than one-third of the inputs being non-robust. Of the properties selected
for Figure 3, the robustness to pixel $L_{\infty}$ perturbations exhibits by
far the lowest scores, with only about half the inputs being robust for the
self-trained model and virtually none for the pre-trained one.

SELF-TRAINEDPRE-TRAINED93.9%96.6%Brightness+=

SELF-TRAINEDPRE-TRAINED37.4%56.1%Adaptive Brightness+=

Figure 4: Robustness comparison of the basic (top) and adaptive brightness
(bottom). Basic brightness applies the same transformation to the whole image,
whereas adaptive brightness allows different transformations for different
parts of the image.

SELF-TRAINEDPRE-TRAINED33.8%27.2%Traffic Sign Stickers+=

Figure 5: Task-specific property that inserts a single sticker of varying
position, size and orientation on the traffic sign.

###  4.2 Stronger and Task-specific Properties

In addition to the generic and basic computer vision properties considered in
section 4.1, many of which naturally transfer to other use cases and
application domains, the robustness of the models for traffic sign recognition
was assessed with respect to more-task specific and stronger properties.
Generally speaking, these properties better reflect specific ambient
conditions that may occur in this use case. Figures 4 and 5 show the results
for two task-specific properties. On the one hand, the adaptive brightness
property is defined as a generalisation of the basic brightness
transformation. While the basic brightness transformation applies brightness
changes uniformly to the whole image, the adaptive brightness transformation
makes fine-grained changes, allowing some parts of the image to become
brighter, and others darker. This better reproduces situations that can occur
in practice, as a result of the material properties of object surfaces with
respect to absorption and reflection of light. On the other hand, the
robustness to generic stickers placed on the traffic signs is assessed.
Especially in metropolitan areas, it is not uncommon for such stickers to be
placed on traffic signs. As can be seen from Figures 4 and 5, the models
become much more brittle when generalising the basic brightness transformation
to its adaptive version, and their robustness to stickers is very low indeed.

###  4.3 Testing Combinations

So far, the robustness properties were analysed independently from each other.
However, a combination of all of the properties with different strengths
reflects real-world conditions more realistically. As a first step towards
evaluating these situations, the properties can be composed together (e.g.,
combining brightness with rotation), effectively creating a new set of
properties with additional degrees of freedom that should be assessed for
robustness. Figure 6 provides an overview of robustness scores of the models
when tested against combinations of two properties. In most cases, the
robustness scores for combinations decrease. For instance, combining
translation and scaling yields a robustness score of 85.5%, whereas the
robustness to translation and scaling only is 94.7% and 93.0%, respectively.
The combination thus decreases robustness by an additional 7.5% as compared to
the lower of the two individual robustness scores.

$\leftarrow$ WORSE $\mid$ BETTER $\rightarrow$Robustness decrease of combined
propertiesSELF-TRAINED$\leftarrow$ WORSE $\mid$ BETTER $\rightarrow$Robustness
decrease of combined propertiesPRE-TRAINED Translate: 94.7%
Scale: 93.0%
Translate + Scale: 85.5% (-7.5%) Rotation: 96.1%
Blur: 96.9%
Rotation + Blur: 93.5% (-3.4%) Figure 6: Effect of assessing model against
combinations of robustness properties. Each point corresponds to the relative
change in robustness of a selected properties combination compared to
considering the properties independently. For most properties, their
combination decreases the robustness (shown as negative change).

##  5 Discussion

The results in section 4 show that significant differences in robustness exist
both between models and properties. Variations in robustness scores can be
partly ascribed to differences in the structure of the models and in the
resolution of their input images (which were taken into account in testing in
order not to distort the results). The different levels of brittleness the
models exhibit to different properties at least partly correlate with the
degree of representation of these properties in the training data set or the
data augmentation process. It is assumed that, hereby, the models are actively
endowed with a certain baseline robustness against such properties. It is
concerning to see that the models’ robustness to some properties is very low
indeed, and using such brittle models in safety- and security-relevant use
cases seems audacious.

The more detailed analyses of stronger and task-specific properties as well as
of combinations of properties outlined in sections 4.2 and 4.3 show that
performing a standard set of basic robustness tests provides a good baseline,
but more comprehensive tests can yield higher levels of confidence and better
insights as to under what conditions the models are robust.

As an additional refinement, the aggregate robustness scores can be replaced
with more precise plots of model robustness across the perturbation parameter
space (e.g., [9, 10]). This can provide detailed insights into the specific
regions where models may fail, which is especially beneficial in safety- and
security-critical applications. Beyond that, a careful, fine-grained
inspection of the robustness results can help pinpoint failure modes, i.e.
specific conditions where model robustness is lower than in others (e.g.,
direct sunlight or very low brightness, as shown in Figure 7). Improving the
training data sets based on the results obtained in the robustness tests, e.g.
by using specially selected, augmented or synthesised data, seems a promising
direction for alleviating their shortcomings.

Priority RoadNon-robust samples per property$\dots$ similar background direct
sunlight object occlusion very low brightness Figure 7: Example of failure
modes for the self-trained model. For a given property, each square denotes a
non-robust sample (or a set of samples). Analysis of the robust (not shown)
and non-robust samples identifies failure modes such as direct sunlight.

The results presented in section 4 are those obtained using the specific use
case of traffic sign recognition, a specific data set and two specific models
trained with very limited resources. It is to be expected that models deployed
in the real-world by car vendors will have been trained using much larger
resources and will, therefore, exhibit higher robustness on average.
Nevertheless, a principled procedure for testing their robustness is necessary
and the approach discussed in this article yields one possible solution.

##  6 Conclusion and Outlook

The tests that were performed in the scope of this work mainly focus on the
robustness of models against naturally occurring, random properties. They only
address to a limited extent the stability of the models with respect to
perturbations crafted by an attacker in a targeted way. On the one hand, it is
to be expected that the models’ robustness can be further decreased by using
specialised attacks and that remedying this problem is much harder than
improving robustness to certain naturally occurring failure modes as
identified by the tests. On the other hand, the effort required for performing
different kinds of attacks strongly varies and, therefore, not all attacks are
realistic depending on the use case, the motivation and the capabilities of
possible attackers.

The testing methodology was applied to the special use case of traffic sign
recognition. The results show that performing such tests in a principled way
is feasible, and they yield valuable insights into the existing limitations of
the models. As a next step, such tests should also be applied to other use
cases, in order to obtain the respective results. In doing so, the questions
to what extent the testing methodology from traffic sign recognition can be
transferred to other use cases and where it needs to be adapted should be
addressed. Whereas the basic tests do not have any specific connection to the
example use case considered and can probably be transferred easily to check
models solving different computer vision problems, the task-specific
properties may well need to be adjusted to reflect the specific ambient
conditions that may arise in such new use cases. More generally, the abstract
testing methodology might also be generalised to AI models solving problems
not based on image data, but on data from other input domains, e.g. acoustic
signals or even structured tabular data, and to other model architectures
beyond neural networks.

For AI models to be used in safety- and security-critical areas such as
(partly) autonomous driving in the years to come, a standardised methodology
and concrete test criteria will be required in order to assess and evaluate
the robustness of these models with respect to random as well as targeted
perturbations. For ensuring an adequate level of safety and security, such
criteria must be developed, and checking compliance with them must be made
compulsory.

Another question that should be further studied is under what conditions and
to what extent it is feasible not only to test the robustness of the models in
an empirical, though principled, way but also to obtain formally verified
statements, which could guarantee the required level of safety and security.
Intense research should be devoted to this question, while aiming at practical
applicability of solutions at least as a mid-term goal. Combining robustified
models with a reject option for specific, non-verifiable regions of the input
space might be a first step in this direction. In addition to that, further
research effort should be devoted to investigating defensive measures that are
able to ensure a high security level even when facing strong attackers, i.e.
adaptive attackers with knowledge of the security measures in place (e.g.,
[24]), and to developing fundamental approaches for increasing the robustness
of AI models.

##  7 Acknowledgements

The authors would like to thank the reviewers for their helpful suggestions.

## References

  * [1] Balunovic, M., Baader, M., Singh, G., Gehr, T., Vechev, M.: Certifying geometric robustness of neural networks. In: Wallach, H., Larochelle, H., Beygelzimer, A., d’Alché Buc, F., Fox, E., Garnett, R. (eds.) Advances in Neural Information Processing Systems. vol. 32. Curran Associates, Inc. (2019)
  * [2] Berghoff, C., Neu, M., von Twickel, A.: Vulnerabilities of Connectionist AI Applications: Evaluation and Defense. Frontiers Big Data 3, 23 (2020). https://doi.org/10.3389/fdata.2020.00023
  * [3] Bielik, P., Tsankov, P., Krause, A., Vechev, M.: Reliability Assessment of Traffic Sign Classifiers. Tech. rep., Bundesamt für Sicherheit in der Informationstechnik (2020), https://www.bsi.bund.de/ki
  * [4] Biggio, B., Roli, F.: Wild patterns: Ten years after the rise of adversarial machine learning. Pattern Recognit. 84, 317–331 (2018). https://doi.org/10.1016/j.patcog.2018.07.023
  * [5] Carlini, N., Athalye, A., Papernot, N., Brendel, W., Rauber, J., Tsipras, D., Goodfellow, I.J., Madry, A., Kurakin, A.: On evaluating adversarial robustness. CoRR abs/1902.06705 (2019)
  * [6] Dalvi, N.N., Domingos, P.M., Mausam, Sanghai, S.K., Verma, D.: Adversarial classification. In: Kim, W., Kohavi, R., Gehrke, J., DuMouchel, W. (eds.) Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 99–108. ACM (2004). https://doi.org/10.1145/1014052.1014066
  * [7] D’Amour, A., et al.: Underspecification presents challenges for credibility in modern machine learning. CoRR abs/2011.03395 (2020)
  * [8] Deng, J., Dong, W., Socher, R., Li, L., Li, K., Li, F.: ImageNet: A large-scale hierarchical image database. In: 2009 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. pp. 248–255. IEEE Computer Society (2009). https://doi.org/10.1109/CVPR.2009.5206848
  * [9] Engstrom, L., Tran, B., Tsipras, D., Schmidt, L., Madry, A.: Exploring the landscape of spatial robustness. In: Proceedings of the 36th International Conference on Machine Learning. vol. 97, pp. 1802–1811. PMLR (2019)
  * [10] Fawzi, A., Moosavi-Dezfooli, S.M., Frossard, P., Soatto, S.: Empirical study of the topology and geometry of deep networks. In: 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 3762–3770 (2018). https://doi.org/10.1109/CVPR.2018.00396
  * [11] Gehr, T., Mirman, M., Drachsler-Cohen, D., Tsankov, P., Chaudhuri, S., Vechev, M.T.: AI2: Safety and Robustness Certification of Neural Networks with Abstract Interpretation. In: 2018 IEEE Symposium on Security and Privacy. pp. 3–18. IEEE Computer Society (2018). https://doi.org/10.1109/SP.2018.00058
  * [12] Geirhos, R., Jacobsen, J.H., Michaelis, C., Zemel, R., Brendel, W., Bethge, M., Wichmann, F.A.: Shortcut learning in deep neural networks. Nature Machine Intelligence 2, 665–673 (Nov 2020)
  * [13] Houben, S., Stallkamp, J., Salmen, J., Schlipsing, M., Igel, C.: Detection of traffic signs in real-world images: The German traffic sign detection benchmark. In: The 2013 International Joint Conference on Neural Networks. pp. 1–8. IEEE (2013). https://doi.org/10.1109/IJCNN.2013.6706807
  * [14] Huang, X., Kwiatkowska, M., Wang, S., Wu, M.: Safety Verification of Deep Neural Networks. In: Majumdar, R., Kuncak, V. (eds.) Computer Aided Verification - 29th International Conference. Lecture Notes in Computer Science, vol. 10426, pp. 3–29. Springer (2017). https://doi.org/10.1007/978-3-319-63387-9_1
  * [15] Katz, G., Barrett, C.W., Dill, D.L., Julian, K., Kochenderfer, M.J.: Reluplex: An Efficient SMT Solver for Verifying Deep Neural Networks. In: Majumdar, R., Kuncak, V. (eds.) Computer Aided Verification - 29th International Conference. Lecture Notes in Computer Science, vol. 10426, pp. 97–117. Springer (2017). https://doi.org/10.1007/978-3-319-63387-9_5
  * [16] Kim, Y., Hwang, H., Shin, J.: Robust object detection under harsh autonomous-driving environments. IET Image Processing (2021). https://doi.org/10.1049/ipr2.12159
  * [17] Michaelis, C., Mitzkus, B., Geirhos, R., Rusak, E., Bringmann, O., Ecker, A.S., Bethge, M., Brendel, W.: Benchmarking robustness in object detection: Autonomous driving when winter is coming. CoRR abs/1907.07484 (2019)
  * [18] Ponn, T., Kröger, T., Diermeyer, F.: Identification and explanation of challenging conditions for camera-based object detection of automated vehicles. Sensors 20(13), 3699 (2020). https://doi.org/10.3390/s20133699
  * [19] Singh, G., Gehr, T., Püschel, M., Vechev, M.T.: An abstract domain for certifying neural networks. Proc. ACM Program. Lang. 3(POPL), 41:1–41:30 (2019). https://doi.org/10.1145/3290354
  * [20] Szegedy, C., Vanhoucke, V., Ioffe, S., Shlens, J., Wojna, Z.: Rethinking the Inception Architecture for Computer Vision. In: 2016 IEEE Conference on Computer Vision and Pattern Recognition. pp. 2818–2826. IEEE Computer Society (2016). https://doi.org/10.1109/CVPR.2016.308
  * [21] Szegedy, C., Zaremba, W., Sutskever, I., Bruna, J., Erhan, D., Goodfellow, I.J., Fergus, R.: Intriguing properties of neural networks. In: Bengio, Y., LeCun, Y. (eds.) 2nd International Conference on Learning Representations (2014), http://arxiv.org/abs/1312.6199
  * [22] Temel, D., Chen, M., AlRegib, G.: Traffic sign detection under challenging conditions: A deeper look into performance variations and spectral characteristics. IEEE Transactions on Intelligent Transportation Systems pp. 1–11 (2019). https://doi.org/10.1109/TITS.2019.2931429
  * [23] Temel, D., Kwon, G., Prabhushankar, M., AlRegib, G.: CURE-TSR: Challenging unreal and real environments for traffic sign recognition. In: Neural Information Processing Systems (NeurIPS) Workshop on Machine Learning for Intelligent Transportation Systems (2017)
  * [24] Tramer, F., Carlini, N., Brendel, W., Madry, A.: On adaptive attacks to adversarial example defenses. In: Advances in Neural Information Processing Systems. vol. 33, pp. 1633–1645. Curran Associates, Inc. (2020)


# Understanding step selection analysis
through numerical integration

Théo Michelot1111Email<EMAIL_ADDRESS>1Dalhousie University,
2University of Sheffield, 3University of Minnesota  Natasha J. Klappstein1
1Dalhousie University, 2University of Sheffield, 3University of Minnesota
Jonathan R. Potts2 1Dalhousie University, 2University of Sheffield,
3University of Minnesota  John Fieberg3 1Dalhousie University, 2University of
Sheffield, 3University of Minnesota

###### Abstract

Step selection functions (SSFs) are flexible statistical models used to
jointly describe animals’ movement and habitat preferences. The popularity of
SSFs has grown rapidly, and various extensions have been developed to increase
their utility, including the ability to use multiple statistical distributions
to describe movement constraints, interactions to allow movements to depend on
local environmental features, and random effects and latent states to account
for within- and among-individual variability. Although the SSF is a relatively
simple statistical model, its presentation has not been consistent in the
literature, leading to confusion about model flexibility and interpretation.
We believe that part of the confusion has arisen from the conflation of the
SSF model with the methods used for statistical inference, and in particular,
parameter estimation. Notably, conditional logistic regression can be used to
fit SSFs in exponential form, and this model fitting approach is often
presented interchangeably with the actual model (the SSF itself). However,
reliance on conditional logistic regression reduces model flexibility, and
suggests a misleading interpretation of step selection analysis as being
equivalent to a case-control study. In this review, we explicitly distinguish
between model formulation and inference technique, presenting a coherent
framework to fit SSFs based on numerical integration and maximum likelihood
estimation. We provide an overview of common numerical integration techniques
(including Monte Carlo integration, importance sampling, and quadrature), and
explain how they relate to popular methods used in step selection analyses.
This general framework unifies different model fitting techniques for SSFs,
and opens the way for improved inferential methods. In this approach, it is
straightforward to model movement with distributions outside the exponential
family, and to apply different SSF model formulations to the same data set and
compare them with AIC. By separating the model formulation from the inference
technique, we hope to clarify many important concepts in step selection
analysis.

##  1 Introduction

The increased availability of animal tracking data has led to the widespread
use of statistical methods to estimate habitat selection at the scale of the
observed movement step. Perhaps the most common model is the step selection
function (SSF; Rhodes et al.,, 2005; Fortin et al.,, 2005), whereby the
likelihood of moving to the spatial location $\bm{x}_{t+1}$ given previous
locations $\bm{x}_{1:t}=\\{\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{t}\\}$ is in
the following form,

$p(\bm{x}_{t+1}\mid\bm{x}_{1:t})=\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\int_{\Omega}w(\bm{x}_{t},\bm{z})\phi(\bm{z}\mid\bm{x}_{1:t})d\bm{z}},$
(1)

where $\Omega$ is the study region. The function $w$ describes the effects of
environmental variables (e.g., resources, risks, and environmental conditions;
Matthiopolous et al.,, 2023), and $\phi$ accounts for the effects of movement
constraints (e.g., on the range of observed step lengths). The habitat
selection function is often assumed to take an exponential (or “log-linear”)
form, i.e.,
$w(\bm{x}_{t},\bm{x}_{t+1})=\exp\\{h(\bm{x}_{t},\bm{x}_{t+1})\bm{\beta}^{\intercal}_{h}\\}$,
where $h(\bm{x}_{t},\bm{x}_{t+1})$ is a vector of habitat variables for the
step, and $\bm{\beta}_{h}$ is the vector of associated selection parameters.
The form of $\phi$ reflects assumptions about movement patterns of the animal,
and it is often written as a function of the step length and turning angle to
capture movement speed and tortuosity of an animal’s movement. We call “step
selection function” the numerator of Equation 1, but the terminology is not
consistent across the literature, and the term has been used variously to
refer to $w$, to $w\times\phi$, or to the whole right-hand side of Equation 1.
We choose to define $w\times\phi$ as the SSF to reflect the fact that an
animal’s selection of a step is based on both habitat preferences and movement
constraints. Figure 1 shows an example step selection function with the two
model components.

Figure 1: Example model components at some time $t$, where the last point
$\bm{x}_{t}$ is at the centre of each panel and the previous steps are shown
as black segments: (a) habitat selection function $w(\bm{x}_{t},\bm{z})$, (b)
movement kernel $\phi(\bm{z}|\bm{x}_{t})$ based on distance from $\bm{x}_{t}$,
and (c) resulting step selection function
$w(\bm{x}_{t},\bm{z})\phi(\bm{z}|\bm{x}_{t})$. The integral on the denominator
of Equation 1 is the volume under the step selection function, which is
required to transform the step selection function into a probability
distribution (sometimes called the step density).

Step selection analysis refers to a wide range of methods for applying SSFs to
animal tracking data, with the aim to estimate the parameters of the habitat
selection function $w$ and the movement kernel $\phi$. Although the data-
generating mechanism for this model is described by Equation 1, this is a
difficult statistical problem due to the presence of the integral in the
denominator. This integral is required in Equation 1 to ensure that
$p(\bm{x}_{t+1}\mid\bm{x}_{1:t})$ integrates to 1, i.e., that it is a valid
probability density function with respect to $\bm{x}_{t+1}$ (Rhodes et al.,,
2005; Forester et al.,, 2009). It also has a more intuitive interpretation: to
evaluate the likelihood of a given step, we weigh its suitability against the
suitability of all other possible steps in the study region. Here,
“suitability” refers both to the habitat quality of a location (as captured by
$w$) and to its accessibility (as captured by $\phi$). This integral cannot
generally be calculated analytically, because the integrand (i.e., the
expression that is integrated) depends on $w$, which is usually a function of
environmental covariates with no mathematically-convenient functional forms.
That is, the integral cannot be rewritten in terms of simple functions that
could be directly implemented with a computer, and so Equation 1 cannot
generally be evaluated for a given movement track and set of parameters.

Although the integral in Equation 1 cannot be computed directly, methods have
been developed to replace the expression by a tractable approximation. In some
cases, such approximations are equivalent to applying conditional logistic
regression (CLR) to a case-control data set, where each observed location
(“case”) is associated with a set of locations from the landscape (“controls”;
Forester et al.,, 2009). This has been a popular framework for step selection
analysis, because CLR can be fitted quickly and conveniently using statistical
software (e.g., using the survival package in R; Therneau,, 2023).
Consequently, SSFs are often conflated with CLR, even though the latter is
merely a convenient tool to fit the former by approximating the likelihood in
some special cases. In our view, this presentation can lead to confusion about
model interpretation, and reduces the flexibility of step selection analyses.
In particular, we can avoid the need to make strong assumptions about the
functional forms of $w$ and $\phi$ if we are willing to use numerical methods
other than CLR for parameter estimation.

In this review, we show that most methods used to estimate parameters in step
selection analyses can be viewed as applications of numerical integration. A
variety of numerical integration techniques are commonly applied in statistics
and mathematics to approximate integrals in cases where there is no known
formula to compute them exactly. We will present several numerical tools
developed for this purpose, and describe their utility in step selection
analyses. This perspective suggests we can contrast existing methods (e.g., to
identify those with lower approximation errors), and opens the way for
improved inferential techniques in step selection analysis. Lastly, we hope
that our review will motivate further exploration of numerical integration and
estimation methods, which have broad utility across a wide range of ecological
applications.

##  2 Likelihood approximation in step selection analysis

###  2.1 Maximum likelihood estimation

We consider $n$ observed locations $\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{n}$,
recorded without error at regular time intervals. The goal of step selection
analysis is to estimate the parameters $\bm{\beta}_{h}$ of the habitat
selection function $w$, which quantify the strength of selection or avoidance
of spatial covariates, and the parameters $\bm{\beta}_{m}$ of the movement
kernel $\phi$, which quantify movement tendencies (e.g. speed and tortuosity).
We denote as $\bm{\beta}=(\bm{\beta}_{h},\bm{\beta}_{m})$ the vector of all
parameters. Based on the standard assumptions that each step follows the model
of Equation 1, and that successive steps are conditionally independent (i.e.,
given past locations), the likelihood of the parameters under the step
selection model can be written as

$L(\bm{\beta};\bm{x}_{1},\dots,\bm{x}_{n})=\prod_{t=1}^{n-1}p(\bm{x}_{t+1}\mid\bm{x}_{1:t})=\prod_{t=1}^{n-1}\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\int_{\Omega}w(\bm{x}_{t},\bm{z})\phi(\bm{z}\mid\bm{x}_{1:t})d\bm{z}}.$

Note that both $w$ and $\phi$ depend on some of the parameters in $\bm{\beta}$
(specifically, $w$ depends on $\bm{\beta}_{h}$ and $\phi$ on
$\bm{\beta}_{m}$), but we do not make this explicit for notational
convenience. The likelihood (or log-likelihood) can be optimised numerically
with respect to $\bm{\beta}$, for example using the function optim or nlm in
R, to obtain maximum likelihood estimates of $\bm{\beta}$. For an overview of
maximum likelihood estimation in an ecological context, see for example
Chapters 6–7 of Bolker, (2008) or Chapter 10 of Fieberg, (2023).

Model fitting requires computing the likelihood function for the observed
data, and therefore evaluating the integral

$I=\int_{\Omega}w(\bm{x}_{t},\bm{z})\phi(\bm{z}\mid\bm{x}_{1:t})d\bm{z}$ (2)

for each time step $t\in\\{1,2,\dots,n-1\\}$. In this section, we describe
several methods to approximate $I$ by some quantity $\widehat{I}$, which can
be substituted in the likelihood formula to carry out approximate inference.
The main approaches are summarised in Table 1.

| Approximate likelihood |  $\bm{z}_{k}\sim$ ? | References
---|---|---|---
MC | $\frac{w(\bm{x}_{t},\bm{x}_{t+1})}{\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})}$ | $\phi(\cdot\mid\bm{x}_{1:t})$ | Fortin et al., (2005)
UMC | $\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t})}$ | uniform | Forester et al., (2009)
IS | $\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t})/h(\bm{z}_{k}\mid\bm{x}_{1:t})}$ | user-defined $h$ | Forester et al., (2009)
UQ | $\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t})}$ | regular grid | Rhodes et al., (2005)
Table 1: Summary of most common numerical integration approaches used in step
selection analyses: Monte Carlo with known movement kernel (MC), uniform Monte
Carlo (UMC), importance sampling (IS), and uniform quadrature (UQ). The
columns include the approximate likelihood of a step from $\bm{x}_{t}$ to
$\bm{x}_{t+1}$, method for determining the distribution of integration points
$\bm{z}_{k}$, and references to key papers presenting each approach.

###  2.2 Monte Carlo integration

We use the term “Monte Carlo integration” to refer to all forms of numerical
integration that rely on random sampling (i.e., all methods presented in this
review, except for quadrature). Monte Carlo integration is a method for
evaluating an integral of the form $\int_{\Omega}f(\bm{z})g(\bm{z})d\bm{z}$,
where $f$ is a probability density function (Section 3.2 of Robert and
Casella,, 2010). The general idea is to generate a sample from the
distribution $f$, and use it to find an unbiased estimate of the integral. It
can be shown that

$\int_{\Omega}f(\bm{z})g(\bm{z})d\bm{z}\approx\frac{1}{K}\sum_{k=1}^{K}g(\bm{z}_{k}),\quad\text{where
}\bm{z}_{k}\sim f,$ (3)

with the error of the approximation decreasing as $K$ increases. Throughout
this review, we use the notation “$\bm{z}_{k}\sim f$” to indicate that
$\bm{z}_{k}$ follows the distribution with probability density function $f$.

The most common approaches to fitting SSFs can be viewed as different forms of
Monte Carlo integration applied to the integral in Equation 2, which result
from different choices of the functions $f$ and $g$ (always with the
constraint that $f\times g=w\times\phi$). Generally, this choice may impact
the accuracy and precision of the approximation, so choosing it requires
thought (see Section 2.2.3 and Rizzo,, 2019).

Maximum likelihood estimation based on the numerical approximation of the
integral in Equation 3 defines a general framework of approximate inference
for SSFs. Methods of inference based on Monte Carlo likelihood approximations
are common in econometrics, where they are usually called _simulated_ maximum
likelihood estimation, to highlight their inherent stochasticity (Section
3.1.2 of Gourieroux and Monfort,, 1996). Even though the Monte Carlo estimator
of the integral given in Equation 3 is not biased, simulated maximum
likelihood estimators generally are, due to the log-transformation of the
likelihood before it is maximised (Gourieroux and Monfort,, 1996). As a
result, Monte Carlo-based step selection parameter estimators are biased, but
this bias decreases with the number of integration points.

Using Monte Carlo integration, Equation 2 can be approximated as
$I\approx\frac{1}{K}\sum_{k=1}^{K}g(\bm{z}_{k})$, for functions $f$ and $g$
chosen such that $\bm{z}_{k}\sim f$ and $f\times g=w\times\phi$. In this
review, we include the observed location $\bm{x}_{t+1}$ as an additional
integration point in the random sample $\\{\bm{z}_{1},\dots,\bm{z}_{K}\\}$,
and we denote $\bm{z}_{0}=\bm{x}_{t+1}$ for convenience. This slight deviation
from the formal definition of Monte Carlo integration is justified in this
context for three reasons: (1) the resulting formulas have clear links to
those presented in the literature, (2) it improves numerical stability and
decreases bias for small $K$ (as illustrated with simulations in Appendix A),
and (3) the effect of this change vanishes for large values of $K$. The
approximate likelihood of a step under Monte Carlo integration is then

$p(\bm{x}_{t+1}\mid\bm{x}_{1:t})\approx\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\frac{1}{K+1}\sum_{k=0}^{K}g(\bm{z}_{k})},\quad\text{where
}\bm{z}_{k}\sim f,$

for some general $f$ and $g$. We present several important special cases
below, which encompass most existing methods for step selection analysis
(e.g., Fortin et al.,, 2005), including extensions proposed by Forester et
al., (2009), Duchesne et al., (2015), and Avgar et al., (2016). A key
characteristic of the first approach, described in Section 2.2.1, is that it
assumes that the movement kernel $\phi$ is known prior to the step selection
analysis. It has become increasingly common to jointly estimate the movement
kernel and habitat selection, and the approaches in Sections 2.2.2–2.2.3 focus
on that situation.

####  2.2.1 Assuming that the movement kernel is known

A popular approach to step selection analysis is to define the movement kernel
prior to fitting the SSF, typically from empirical or parametric distributions
of step lengths and turning angles (Fortin et al.,, 2005). Then, since $\phi$
is assumed to be known, we can apply Monte Carlo integration to the integral
of Equation 2 by choosing $f=\phi$ and $g=w$. Equation 2 becomes

$I\approx\frac{1}{K+1}\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k}),\quad\text{where
}\bm{z}_{k}\sim\phi(\cdot\mid\bm{x}_{1:t}).$

That is, we generate the random Monte Carlo sample from the movement kernel
$\phi$, and we take the mean of the habitat selection function at those random
points to evaluate the integral. This characteristic of the typical step
selection analysis workflow has led practitioners to view the $\bm{z}_{k}$ as
a sample from the “available” landscape. The key limitation of this approach
is that the movement kernel cannot be estimated jointly with habitat selection
parameters, because the points $\bm{z}_{k}$ are sampled from it prior to model
fitting. In addition, this approach to estimating $\phi$ without consideration
of $w$ has been shown to result in biased parameter estimators since the
observed movements are a function of both processes (Forester et al.,, 2009).

Using this method, the likelihood of a step (Equation 1) is approximated by

$p(\bm{x}_{t+1}\mid\bm{x}_{1:t})\approx\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\frac{1}{K+1}\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})}\propto\frac{w(\bm{x}_{t},\bm{x}_{t+1})}{\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})},$
(4)

where $\bm{z}_{k}\sim\phi(\cdot\mid\bm{x}_{1:t})$. Note that multiplicative
constants, i.e., terms that do not depend on the parameters $\bm{\beta}$, can
be omitted from the likelihood function with no effect on inference. Here,
$\phi$ can be omitted from the numerator because it is assumed known, and
therefore does not depend on any estimated parameter, and $1/(K+1)$ is omitted
from the denominator. When $w$ is written in the usual exponential form, this
is the likelihood of a conditional logistic regression (CLR) model, and
parameters can be estimated using standard software such as the clogit
function in the survival R package (Therneau,, 2023). However, the method of
Equation 4 has also been used for non-exponential functional forms for $w$, by
using likelihood maximisation procedures other than CLR (Potts et al.,, 2014).

This approach was initially proposed by Fortin et al., (2005), and has been
widely used due to the convenience of implementation using CLR, and the
intuitive appeal of interpreting the random locations as a sample of
availability. Estimating the movement kernel separately has drawbacks,
however: as previously mentioned, this approach leads to biased parameter
estimators and does not propagate statistical uncertainty about the movement
model to the habitat selection parameters. In addition, it does not allow for
interactions between the movement parameters and local environmental features.
As a consequence, recent research in step selection analysis has focused on
formulations where $\phi$ and $w$ are estimated simultaneously (Rhodes et
al.,, 2005; Forester et al.,, 2009; Avgar et al.,, 2016), and we present those
in the next subsections. The interpretation of the random points $\bm{z}_{k}$
is different when $\phi$ is estimated as part of the SSF, and we discuss this
in Section 4.

####  2.2.2 Uniform Monte Carlo sampling

If the movement kernel $\phi$ is not known a priori, we must choose a
different probability density function $f$ from which to generate random
points. One natural choice is to use a uniform distribution over the domain
$\Omega$, because it is easy to sample from. In this case, the two functions
in Equation 3 are defined as $f(\bm{z})=1/A(\Omega)$ (where $A(\Omega)$ is the
area of $\Omega$), and
$g(\bm{z})=A(\Omega)w(\bm{x}_{t},\bm{z})\phi(\bm{z}\mid\bm{x}_{1:t})$. Then,
Equation 2 can be approximated as

$I\approx\frac{A(\Omega)}{K+1}\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t}),\quad\text{where
}\bm{z}_{k}\sim\text{Unif}(\Omega).$

Intuitively, the random points are used to estimate the mean value of the SSF
over $\Omega$, and the integral is approximated by the product of that mean
value and $A(\Omega)$ (see Figure 2a for a one-dimensional example).

Figure 2: Illustration of numerical integration in one dimension, for the
function $f\times g$ shown as a black line, over the interval $\Omega=[-5,5]$.
The orange dots are the integration points, and the blue dots are the
corresponding function evaluations. (a) Monte Carlo integration with a uniform
sample over $\Omega$; the height of the grey rectangle is the mean of function
evaluations. (b) Monte Carlo integration with a uniform sample over $[-3,3]$.
(c) Importance sampling with random points generated from a normal
distribution that roughly approximates the integrand $f\times g$. (d)
Quadrature over a regular grid using a Riemann sum. In (a), (b) and (d), the
shaded area approximates the area under the curve; there is no such simple
visualisation method for importance sampling. In this small simulated example,
the true integral is 2.53, and the approximations are (a) 2.83, (b) 2.28, (c)
2.41, and (d) 2.80.

This method, called “uniform sampling” by Forester et al., (2009), allows for
joint estimation of habitat selection and movement parameters (e.g., Schlägel
and Lewis,, 2014). However, the uniform sampling approach can be
computationally demanding, because good performance requires that the
integration points provide adequate coverage of the study region $\Omega$, and
this can often only be achieved for large values of $K$. In most step
selection analyses, the number of points can be greatly reduced based on the
observation that, at each step, the SSF decreases sharply with distance from
the start point $\bm{x}_{t}$ (due to movement constraints of the animal).
Points far from $\bm{x}_{t}$ therefore contribute a negligible amount to the
integral, and the approximation is virtually unchanged if the domain of
integration is truncated to a disc around $\bm{x}_{t}$, with radius large
enough to encompass any possible step (Boyce et al.,, 2003; Craiu et al.,,
2008). Fewer points are needed to ensure good coverage of this disc, which
reduces the computational cost of evaluating the integral (Klappstein et al.,,
2022). Uniform sampling on a truncated interval is illustrated in Figure 2b.

Using this approximation, Equation 1 becomes

$p(\bm{x}_{t+1}\mid\bm{x}_{1:t})\approx\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\frac{A(\Omega)}{K+1}\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t})}\propto\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t})},$
(5)

where $z_{k}\sim\text{Unif}(\Omega)$. Like before, we remove $A(\Omega)/(K+1)$
from the denominator because it is a constant and thus will not affect maximum
likelihood estimation. Here, we cannot omit $\phi$ from the numerator because
it is not assumed to be known, and therefore is a function of the parameters
$\bm{\beta}$ of interest. If both $\phi$ and $w$ have an exponential form,
then Equation 5 is equal to the CLR likelihood, as long as we include the
observed location as an integration point.

Note that uniform Monte Carlo sampling refers to spatially uniform points, and
this should not be confused with sampling points with uniform distances from
$\bm{x}_{t}$ (Avgar et al.,, 2016). Generating uniform distances will not
result in a spatially uniform distribution of end points $\bm{z}_{k}$, and so
the above formulas do not hold in that case. This is due to the fact that the
set of possible long steps is spread over a larger area than the set of
possible short steps; therefore if distances are uniform, points will be
relatively more concentrated around the origin than far from it (Rhodes et
al.,, 2005). The case where distances are sampled from uniform distributions
can in fact be viewed as a special case of importance sampling (Section
2.2.3).

####  2.2.3 Importance sampling

The precision of numerical integration depends on the choice of integration
points; generally, the variability in the approximation is lower if points are
concentrated in areas where the function takes large values. Importance
sampling is a method to increase the precision of Monte Carlo integration by
generating random points from a user-defined distribution $h$, called the
importance function (Section 3.3 of Robert and Casella,, 2010, and see an
illustration in Figure 2c). To apply importance sampling to an SSF, we choose
$f=h$ and $g=(w\times\phi)/h$ (so that $f\times g=w\times\phi$ as required),
and Equation 3 gives us the following approximation for the SSF integral,

$I\approx\frac{1}{K+1}\sum_{k=0}^{K}\frac{w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t})}{h(\bm{z}_{k}\mid\bm{x}_{1:t})},\quad\text{where
}\bm{z}_{k}\sim h(\cdot\mid\bm{x}_{1:t}).$ (6)

The only constraint on $h$ is that it should be strictly positive over
$\Omega$; when this is not the case, the approximation of the integral is
truncated to the support of $h$ (i.e., the geographical area over which
$h>0$). Note that, in Equation 6, the importance function $h$ is used to weigh
the contribution of each sampled point to the approximation; this is required
to correct for the preferential sampling of some points over others when
generating $\bm{z}_{k}$.

Importance sampling is useful because the function $h$ can be chosen in such a
way that the variance of the integral estimator decreases (i.e., its precision
increases). The aim is to choose a function $h$ with a shape that is as
similar as possible to $g$ (Section 6.6 of Rizzo,, 2019). This is a convenient
framework in the context of SSFs, because it is often possible to determine
where the SSF will take large values, based on the movement constraints of the
animal. Animals are likely to avoid long steps, and so the SSF often decays
rapidly as distance from the start point increases (Figure 1). The speed of
this decay can be determined approximately from the data, for example by
fitting a distribution to the observed step lengths, and this information can
be used to define the importance function $h$. For example, $h$ could be
chosen as a bivariate normal distribution centred on the last location
$\bm{x}_{t}$, or the two-dimensional spatial distribution implied by step
length and turn angle distributions estimated from the data.

Using importance sampling with function $h$, the approximate likelihood of a
step under the SSF model is

$\displaystyle p(\bm{x}_{t+1}\mid\bm{x}_{1:t})$
$\displaystyle\approx\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\frac{1}{K+1}\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t})/h(\bm{z}_{k}\mid\bm{x}_{1:t})}$
$\displaystyle\propto\frac{w(\bm{x}_{t},\bm{x}_{t+1})\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})}{\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t})/h(\bm{z}_{k}\mid\bm{x}_{1:t})},$
(7)

where $z_{k}\sim h$. The two other Monte Carlo approaches (Sections
2.2.1–2.2.2) can be viewed as special cases of importance sampling. If $h$ is
the probability density function of a uniform distribution over geographical
space, it is constant and can be omitted in Equation 2.2.3, and we obtain
Equation 5. Alternatively, if $\phi$ is assumed to be known, and we choose
$h=\phi$, then $h$ and $\phi$ cancel out in the denominator of Equation 2.2.3;
omitting $\phi$ from the numerator because is is known, this simplifies to
Equation 4.

Although the term “importance sampling” has rarely been used in the SSF
literature, this and equivalent methods have been widely advocated, starting
with the recommendation of Forester et al., (2009) to distinguish between the
sampling function $h$ and the movement model $\phi$. Forester et al., (2009)
derived a formula very similar to Equation 2.2.3, but with one small
difference: their numerator is divided by $h(\bm{x}_{t+1}\mid\bm{x}_{1:t})$.
This difference is inconsequential because $h(\bm{x}_{t+1}\mid\bm{x}_{1:t})$
does not depend on the estimated parameters, so excluding it does not affect
inference. The widely-used methods of Avgar et al., (2016), often called
integrated step selection analysis, are based on Forester et al., (2009) and
can also be viewed as importance sampling. They focus on cases where $h$ and
$\phi$ are chosen to be from the same exponential family of distributions, so
that the calculations simplify and CLR can be used, but there is no such
restriction when the approach is implemented with maximum likelihood
estimation. In the approach of Avgar et al., (2016), the parameter estimates
are adjusted after model fitting to correct for the sampling design. Here, the
bias is corrected directly by including $h$ in Equation 2.2.3 when
implementing the likelihood function. Another notable example is Johnson et
al., (2008), who explicitly suggested importance sampling for a weighted
distribution model analogous to the SSF in Equation 1, and used maximum
likelihood estimation to fit the model without the need for CLR. More
recently, Klappstein et al., (2023) and Pohle et al., (2023) also recognised
that the approach proposed by Forester et al., (2009) was a form of importance
sampling.

Figure 3 contrasts importance sampling and uniform Monte Carlo sampling
(Section 2.2.2) for an example SSF. Similar to the one-dimensional example
shown in Figure 2, the intuition is that the precision of the approximation of
the integral depends on the coverage of areas where the SSF takes high values.
More specifically, the variance of the approximation is minimised when the
distribution of integration points (i.e., the importance function $h$) is
proportional to the SSF $w\times\phi$. In practice, $w$ and $\phi$ are not
known, but a good heuristic is to approximate the movement kernel with
parametric distributions (which can then be sampled from), and use these to
define $h$.

Figure 3: Illustration of three integration designs for an example SSF. The
triangle in the centre shows the last location $\bm{x}_{t}$, and the heatmap
shows the function that needs to be integrated, i.e.,
$w(\bm{x}_{t},\bm{z})\phi(\bm{z}\mid\bm{x}_{1:t})$. This function decreases
with distance to $\bm{x}_{t}$ due to the animal’s movement constraints. The
black dots in each panel represent 100 points used for numerical integration:
(a) uniform points over a disc, (b) points generated from an importance
function based on distance to $\bm{x}_{t}$, and (c) regular quadrature grid
over a disc. Importance sampling generates more points in areas where the
function is high, and will typically lead to a better approximation of the
integral.

###  2.3 Quadrature

Quadrature is a deterministic (non-random) alternative to Monte Carlo
sampling, where the function is evaluated on a user-defined grid of points
over the domain of integration. The simplest example of quadrature is the
Riemann sum (see Figure 2d); in one dimension, the integral is approximated by
the sum of the areas of rectangles with heights determined by function
evaluations along a regular grid of points. The two-dimension extension of
this approach consists of evaluating the function at regularly-spaced points
along a two-dimensional grid, and approximating the integral by the sum of the
volumes of cuboids.

In the context of an SSF, quadrature can be applied by evaluating
$w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t})$ at the centres of
$K+1$ grid cells, $\\{\bm{z}_{0},\bm{z}_{1},\bm{z}_{2},\dots,\bm{z}_{K}\\}$,
and calculating

$I\approx
A\sum_{k=0}^{K}w(\bm{x}_{t},\bm{z}_{k})\phi(\bm{z}_{k}\mid\bm{x}_{1:t}),$

where $A$ is the area of each grid cell. Unlike Monte Carlo approaches, there
is no variance in the estimate of the integral obtained from quadrature,
because the grid is fixed (for a given $K$) rather than random. However, there
is some error in the estimate (akin to statistical bias), and this error
decreases as the spatial resolution of the quadrature grid increases (i.e., as
$K$ increases).

Rhodes et al., (2005) proposed this approach, and pointed out that the grid
cells need not cover the entire study area, and the approximation can be
truncated to the region where the function is non-negligible (within some
distance of $\bm{x}_{t}$). This can speed up computations, similar to the
observation in Section 2.2.2 that Monte Carlo integration with uniform
sampling can be performed over a disc centred on the last location, as long as
the radius of the disc is big enough. We call the approach where integration
points are on a regular grid “uniform quadrature”, and it is illustrated in
Figure 3c for an example SSF. The equation for approximating the integral
using uniform quadrature is identical to that obtained for uniform Monte Carlo
sampling, and so the expression for the SSF likelihood is given by Equation 5.

In step selection analysis, environmental covariates are often available only
over the discrete cells of a raster, making the centroids of the raster cells
a natural choice for the grid of quadrature (Rhodes et al.,, 2005), but this
is not the only possible choice. Arce Guillen et al., (2023) described a new
method of inference for SSFs, where the function to integrate is evaluated at
the nodes of a (deterministic) triangular mesh. This is another form of
quadrature, which works well for the general class of spatial point processes
implemented in the inlabru spatial modelling software that they use (Simpson
et al.,, 2016; Bachl et al.,, 2019).

Some have argued that deterministic numerical integration is preferable in
habitat selection modelling because it has better properties than Monte Carlo
integration for low-dimensional integrals like the one in Equation 2 (Warton
and Shepherd,, 2010; Arce Guillen et al.,, 2023). However, this comparison
assumed uniform Monte Carlo sampling, for which the performance can be
improved substantially using importance sampling with a well-chosen function
$h$. An interesting alternative might therefore be adaptive quadrature, where
the spatial arrangement of points in a (deterministic) spatial grid is chosen
based on the shape of the function. The general idea is to iteratively
subdivide the domain of integration, in such a way that regions where the
function is more irregular are subdivided further. This is analogous to the
idea in importance sampling of generating points in regions where the function
is more complex or takes higher values (Pinheiro and Chao,, 2006).

##  3 Illustration

We illustrate some of the key concepts and methods using simulations and a
real data analysis. The general approach to fitting SSFs presented in Section
2 requires implementing the (approximate) likelihood function, rather than
relying on existing CLR software. Writing custom code greatly increases the
flexibility of the model formulation and inference methods. To help readers
implement their own step selection analyses, we provide R code that can be
used as a starting point. We aimed for a trade-off between simplicity and
flexibility, and we provide basic functions that can be tailored to fit a wide
range of model formulations. The documented code and detailed examples are
provided in Appendix C.

###  3.1 Comparing sampling designs in simulations

Different methods of fitting SSFs can be viewed as different numerical
integration approaches for the same underlying model. When using Monte Carlo
and quadrature methods, the placement of integration points is known to affect
the precision of the results (Rizzo,, 2019), and this provides a rigorous
grounding for the intuition that some methods of generating locations perform
better than others. The closer the distribution of integration points is to
the true SSF, the more precise the approximation will be. In this section, we
use simulations to assess the effect of several design choices on our ability
to recover the parameters of an SSF. We do not intend these simulations to be
exhaustive, or to provide general guidelines to select the best sampling
design in every application, as the choice of method and model is study-
specific. Rather, our aim is to showcase the qualitative effect of different
design choices, encouraging researchers to explore the various options for
themselves when they perform step selection analysis .

We considered three designs for the integration points: (1) points sampled
uniformly at random over a disc of radius equal to the maximum observed step
length, (2) importance sampling where $h$ was based on a gamma distribution of
distances and a von Mises distribution of turning angles, and (3) uniform
quadrature where integration points were defined on a regular spatial grid
over a disc of radius equal to the maximum observed step length. In all cases,
the parameters of the importance function were chosen based on the empirical
distributions of step lengths and/or turning angles.

We first simulated a movement track of length 1000 from an SSF with known
parameters. For the habitat selection component, we used
$w(\bm{x}_{t},\bm{x}_{t+1})=\exp(\beta\times c_{1}(\bm{x}_{t+1}))$, where
$c_{1}$ was a simulated spatial covariate shown in Appendix B and $\beta=5$
was the corresponding selection parameter. The movement kernel was chosen as
$\phi(\bm{x}_{t+1}\mid\bm{x}_{1:t})=\exp(-4L_{t}+3\cos(\alpha_{t}))$, where
$L_{t}$ and $\alpha_{t}$ are the step length and turning angle at time $t$,
respectively. This formulation implies that step lengths followed a gamma
distribution with mean 0.5 and standard deviation 0.35, and turning angles
followed a von Mises distribution with mean zero and concentration 3. For each
scenario, we ran the following steps for increasing numbers of random points,
$K\in\\{5,10,\dots,195,200\\}$:

  1. 1.

Generate $K$ integration points for each observed step, based on a given
numerical integration method;

  2. 2.

Estimate the selection parameter $\beta$ using maximum likelihood estimation,
based on the integration points from step 1.

For Monte Carlo sampling methods, we repeated these steps 50 times to capture
estimator variability, yielding 50 estimated parameters $\widehat{\beta}$ for
each approach and each number of random points. For deterministic quadrature,
there is no sampling variability and so the procedure only needed to be run
once for any given $K$. Note that, for quadrature, the number of integration
points was constrained by the design of the grid, and so it does not always
exactly match the values of $K$ listed above. For each simulation, we
evaluated estimator bias due to the numerical integration, as the difference
between the estimated parameter and the asymptotic estimate obtained when $K$
is very large (here, $K=5000$). We did not compare the estimates to the true
parameter values used in simulation, because this would not make it possible
to separate bias due to integration error (which we are interested in) and
finite-sample bias inherent to maximum likelihood estimation (which does not
depend on integration method).

Figure 4 compares the bias in the selection parameter $\beta$ for the three
numerical integration approaches, over a range of numbers of random points.
The bias and variance decreased as the number of random points increased for
all tested methods. Overall, uniform Monte Carlo sampling had lower precision
(i.e., higher variation) than importance sampling; this is to be expected, as
importance sampling uses observed movement patterns to generate integration
points efficiently. Uniform quadrature was generally positively biased (i.e.,
the selection parameter was overestimated) and, although the bias decreased as
$K$ increased, both random sampling methods seemed to perform better in this
analysis. As described in Section 2.2, the bias is due to the log-
transformation of the likelihood for optimisation: the parameter value that
maximises the approximate log-likelihood can be biased even if the estimator
of the likelihood is unbiased (Gourieroux and Monfort,, 1996).

Figure 4: Results of simulation study. Bias (estimate $-$ truth) in estimated
selection parameter for one spatial covariate, as function of the number of
integration points. Results are compared for three methods for performing
numerical integration: uniform Monte Carlo (blue boxes), importance sampling
with gamma-distributed distances and von Mises-distributed angles (red boxes),
and uniform quadrature (black crosses). Each box represents 50 estimated
values. For quadrature, the number of integration points was constrained by
the design of the grid, so values of $K$ do not exactly match those used for
other methods.

These simulations show that the performance of any method of numerical
integration is dependent on the number of integration points. In practice,
there is no consensus about how high this number should be, and Thurfjell et
al., (2014) reported that different studies have used a wide range of numbers
of random points, between $K=2$ and $K=200$. The minimum number needed for a
given analysis depends on many factors, such as the length of the observed
times series, the sampling scheme used, and the complexity of the SSF model
formulation. For this reason, we recommend that practitioners try several
numbers of random points, until the parameter estimates stabilise, to ensure
that the approximation error in the results is small. This is consistent with
the results presented in Figure 2 of Fieberg et al., (2021), and similar to
the advice of Warton and Shepherd, (2010) and Northrup et al., (2013) in the
context of resource selection functions. Small samples can also lead to
numerical instability and failure to converge during model fitting.

###  3.2 Comparing models using AIC and the same set of integration points

To demonstrate how an understanding of numerical integration techniques and
their use in step selection analyses can facilitate model comparisons, we
considered location data from a red deer (Cervus elaphus) in Northern Germany,
automatically loaded with the amt R package as the data object deer (Signer et
al.,, 2019). The locations are on a regular 6-hour time grid, and the package
also provides a binary raster layer for forest cover (through the
get_sh_forest() function). We will use importance sampling with a single set
of integration points to fit multiple models and compare them using AIC.

We generated integration points using gamma-distributed step lengths and von
Mises-distributed turning angles, with parameters estimated from their
empirical distributions. We compared four SSF formulations, all with the same
habitat selection component $w$ (with forest cover as covariate), but with
different families of distributions of step lengths $L_{t}$ and turning angles
$\alpha_{t}$ in the movement kernel: (i) $L_{t}\sim\text{Exp}(\theta_{1})$ and
$\alpha_{t}\sim\text{uniform}(-\pi,\pi)$, (ii)
$L_{t}\sim\text{gamma}(\theta_{1},\theta_{2})$ and
$\alpha_{t}\sim\text{uniform}(-\pi,\pi)$, (iii)
$L_{t}\sim\text{gamma}(\theta_{1},\theta_{2})$ and $\alpha_{t}\sim\text{von
Mises}(0,\kappa)$, and (iv) $L_{t}\sim\text{Weibull}(\theta_{1},\theta_{2})$
and $\alpha_{t}\sim\text{wrapped Cauchy}(0,\kappa)$. Note that we make a
distinction between the distributions used to generate random locations (which
determine the importance function $h$) and the distributions used to specify
the movement kernel $\phi$. This allows us to fit all four SSFs using the same
integration points, i.e., the exact same data set, such that the models can be
compared using AIC.

The results are shown in Table 2. For this data set, AIC favoured the SSF
formulation (iv), where step lengths were modelled with a Weibull distribution
and turning angles with a wrapped Cauchy distribution. This example
illustrates two advantages of implementing maximum likelihood estimation using
numerical integration: the flexibility of choosing non-exponential models for
the movement kernel (the Weibull and wrapped Cauchy distributions are not in
the exponential family), and the ability to compare models with different
movement formulations.

| Step length | Turning angle | AIC |  $\Delta$AIC
---|---|---|---|---
(iv) | Weibull | wrapped Cauchy | 2500 | 0
(ii) | gamma | uniform | 2507 | 7
(iii) | gamma | von Mises | 2509 | 9
(i) | exponential | uniform | 2521 | 21
Table 2: Model comparison for deer analysis. The four SSFs are specified using
different distributions of step lengths and turning angles, and they are shown
in order of increasing AIC (i.e., the better model is at the top). $\Delta$AIC
is the difference in AIC between each model and the better model.

##  4 Rethinking step selection analysis

###  4.1 Beyond conditional logistic regression

It is common for step selection analysis to be presented as conditional
logistic regression (CLR), where each observed step is a “case” associated
with a set of “control” (random) steps. Indeed, in many important special
cases, the SSF likelihood is approximately equivalent to that of a CLR model
(Fortin et al.,, 2005; Forester et al.,, 2009), which allows parameters to be
estimated with CLR software (Signer et al.,, 2019; Therneau,, 2023). To
justify that approach, Forester et al., (2009) described an SSF as a discrete-
choice model, which is a popular way to present resource selection functions
(Cooper and Millspaugh,, 1999). In that framework, the animal is assumed to
have access to a finite number of discrete and mutually exclusive resource
units, and the model describes how it chooses one unit over the others. In a
step selection model, however, the animal has infinitely many movement options
(because it moves over a continuous space). To reduce the problem to a
discrete choice, Forester et al., (2009) conditioned each movement decision on
a set of random (“control”) points, and assumed that those represented the
animal’s options for that time step. This can be viewed as an approximation of
the target model, akin to the numerical integration approaches that we
presented in Section 2, and the two approaches lead to equivalent formulas.

The equivalence between CLR and the SSF likelihood has allowed ecologists to
leverage standard statistical software for quantifying drivers of movement and
habitat selection for many years. Recently, Muff et al., (2020) and Chatterjee
et al., (2023) demonstrated how a similar equivalence between CLR and Poisson
regression with fixed stratum-specific intercepts could be exploited to model
individual variability in habitat selection and movement parameters using
random effects. Despite these equivalencies, we emphasise that CLR is not the
model of interest, but rather a tool to fit the target SSF model, shown in
Equation 1. Notably, the equivalence between CLR and the SSF likelihood only
holds when both the habitat selection function $w$ and the movement kernel
$\phi$ have an exponential form. CLR is therefore limited to distributions
from the exponential family when modelling $\phi$, and most step selection
modelling research has therefore defaulted to using an exponential or gamma
distribution for step lengths and the von Mises distribution for turn angles
(Duchesne et al.,, 2015; Avgar et al.,, 2016). There is no such restriction in
the general model, however, which allows for the use of a much wider range of
distributions, such as the Weibull distribution for step length and the
wrapped Cauchy distribution for turning angle. As our applied example
illustrates (Section 3.2), these alternatives can potentially lead to improved
model fit. Likewise, it is possible to model habitat selection using functions
other than the exponential (e.g., Potts et al.,, 2014; Schlägel et al.,,
2017), although it is a natural choice for continuous variables since
selection functions should be positive and unbounded (McDonald,, 2013).

###  4.2 Separating model and inference

As we have already mentioned, it is important to distinguish between the
choice of model formulation and the method of parameter estimation. Here, the
model is of the form shown in Equation 1, and the main modelling decision is
defining the functional forms of $w$ and $\phi$. For a given model
formulation, many possible numerical integration methods can be used to
approximate the likelihood, and therefore to estimate the model parameters, as
described in Section 2. These methods, together with any implementation scheme
such as CLR or other likelihood maximisation, constitute the inference
procedure, and are not the model itself. Choices of model and inference
procedure are both important for the analysis: the model formulation should
capture important features of the data-generating process (i.e., animal
movement and habitat selection), whilst the inference procedure should be
chosen carefully to reduce the approximation error. We posit that much
confusion has arisen in the context of step selection analysis due to the
conflation of model and inference, especially when comparing techniques.

One particular area of confusion has been the interpretation of the
integration points $\\{\bm{z}_{1},\dots,\bm{z}_{K}\\}$ generated as part of
model fitting. Due to the historical influence of Fortin et al., (2005), who
assumed that the movement model $\phi$ was known and used it to generate
random locations (Section 2.2.1), the $\bm{z}_{k}$ are commonly assumed to be
a sample of “available” points (or, equivalently, steps connecting the
previous location and $\bm{z}_{k}$ are seen as possible movements). However,
this is crucially not the correct interpretation in other approaches, such as
uniform sampling (Section 2.2.2) and importance sampling (Section 2.2.3).
Indeed, with those methods, the sample of random locations does not have any
particular biological interpretation, and it merely constitutes a numerical
tool to approximate an integral over space. This is the case in most modern
step selection analyses, in which availability is estimated jointly with
habitat selection through the parameters of the function $\phi$, rather than
assumed known a priori (Rhodes et al.,, 2005; Forester et al.,, 2009; Avgar et
al.,, 2016).

Another important distinction is between the distributions used to model
movement, and the distributions used to generate random locations. The
modelled distributions are specified through the choice of $\phi$,
independently of the choice of sampling function (e.g., $h$ in importance
sampling). Forester et al., (2009) and Avgar et al., (2016) showed that it is
mathematically convenient to generate random locations from the same family of
distributions that is used in $\phi$, but this is not a necessity. As we
demonstrate in Section 3.2, it is for example possible to use a gamma
distribution of distances to generate random locations (i.e., in the
importance function $h$), but specify $\phi$ so as to model step lengths with
a Weibull distribution in the analysis. Furthermore, it is straightforward to
fit SSFs with different movement kernel formulations on the exact same data
set (including identical random locations). This makes it possible to select
the formulation for the movement kernel based on standard model selection
criteria such as AIC. This is not possible within the workflow outlined by
Avgar et al., (2016) and implemented by Signer et al., (2019), which requires
$h$ and $\phi$ to be from the same family. When the model and method of
inference are separate, it is also easier to determine how different
functional forms for the movement kernel $\phi$ can be implemented in
practice, and we discuss this in Appendix D.

In fact, the distribution from which the random locations are generated only
matters insofar as different distributions might require different numbers of
integration points to achieve low error (Section 3.1). For a large enough
number of integration points, the choice of distribution is inconsequential.
In practice, when the computational cost is high, it might not be an option to
increase the size $K$ of the random sample arbitrarily; in this case, it is
preferable to choose an importance function $h$ that reduces the estimation
variance for a given $K$. Importantly, the choice of $h$ does not reflect a
modelling assumption, and it is used merely as an inferential tool to reduce
the approximation error.

##  5 Conclusion

The general problem of evaluating complex integrals has been studied
extensively, and many different approaches could be adapted in the context of
step selection analysis. For example, the Laplace approximation is a versatile
method where it is not necessary to evaluate the function at integration
points. Instead, it is assumed that the integrand is well approximated by a
normal distribution, for which the integral is known. In step selection
analyses, the integrand typically combines the movement kernel and the habitat
selection function, resulting in a complex (possibly multimodal) function, and
it is not clear whether the Gaussian assumption would be reasonable. As
suggested in Section 2.3, another promising direction is to combine ideas from
quadrature and importance sampling, and choose the nodes of a deterministic
grid to improve the approximation in areas where the SSF is irregular. This is
similar to adaptive quadrature, which is already widely used in ecology for
non-Gaussian generalised linear mixed models (Bolker et al.,, 2009).

All model fitting approaches are approximately equivalent in the limit where
the number of points in the Monte Carlo sample or the quadrature grid is large
(i.e., as $K\to\infty$). It might therefore seem unnecessary to concern
ourselves with the design choices described in this review, because we can
always increase $K$ until the bias and/or variance in the estimation are
negligible. Because SSFs are relatively simple models, computation time is
often small for modest-sized data sets, and the additional cost of increasing
$K$ might be moderate. However, we have shown that using simplistic numerical
integration techniques can cause bias to persist even for quite large $K$
(e.g., uniform Monte Carlo and quadrature; Figure 4), so thinking carefully
about the numerical integration technique employed may often be preferable to
simply increasing $K$. Furthermore, step selection models are becoming more
sophisticated and complex, so computational effort might become the bottleneck
of many analyses (e.g., multistate models; Nicosia et al.,, 2017). We
anticipate that the sampling design will become increasingly critical in those
cases, as well as in studies with large data sets.

We are certainly not recommending that all biologists stop using CLR software
for step selection analyses, as it is a convenient, fast, and stable
implementation for many purposes. In cases where the habitat selection and
movement are modelled in exponential form, the CLR approach outlined by
Forester et al., (2009), Duchesne et al., (2015), and Avgar et al., (2016) can
safely be applied, e.g., with the amt R package (Signer et al.,, 2019).
However, even in that context, understanding the role numerical integration
plays in parameter estimation can shed light on several important technical
details. In particular, this approach clarifies the role of “control” points
in step selection analysis: although those locations are generally a numerical
tool rather than a biologically-relevant set of spatial locations, placing
them in areas where the animal is most likely to move decreases the required
number of points. In addition, importance sampling explains the post-hoc
parameter adjustment proposed by Avgar et al., (2016), which is needed to
correct for the choice of importance function (see also Appendix S3 of that
paper). Overall, we have shown that the interpretation of step selection
models ought to be separated from any specific model fitting approach, and
that understanding this dichotomy between model construction and model fitting
leads to much broader application of step selection functions than those
confined to a CLR approach.

### Acknowledgements

JF was supported by National Aeronautics and Space Administration award
80NSSC21K1182 and received partial salary support from the Minnesota
Agricultural Experimental Station.

## References

  * Arce Guillen et al., (2023) Arce Guillen, R., Lindgren, F., Muff, S., Glass, T. W., Breed, G. A., and Schlägel, U. E. (2023).  Accounting for unobserved spatial variation in step selection analyses of animal movement via spatial random effects.  bioRxiv, pages 2023–01.
  * Avgar et al., (2016) Avgar, T., Potts, J. R., Lewis, M. A., and Boyce, M. S. (2016).  Integrated step selection analysis: bridging the gap between resource selection and animal movement.  Methods in Ecology and Evolution, 7(5):619–630.
  * Bachl et al., (2019) Bachl, F. E., Lindgren, F., Borchers, D. L., and Illian, J. B. (2019).  inlabru: an R package for bayesian spatial modelling from ecological survey data.  Methods in Ecology and Evolution, 10:760–766.
  * Bolker, (2008) Bolker, B. M. (2008).  Ecological models and data in R.  Princeton University Press, Princeton.
  * Bolker et al., (2009) Bolker, B. M., Brooks, M. E., Clark, C. J., Geange, S. W., Poulsen, J. R., Stevens, M. H. H., and White, J.-S. S. (2009).  Generalized linear mixed models: a practical guide for ecology and evolution.  Trends in Ecology & Evolution, 24(3):127–135.
  * Boyce et al., (2003) Boyce, M. S., Mao, J. S., Merrill, E. H., Fortin, D., Turner, M. G., Fryxell, J., and Turchin, P. (2003).  Scale and heterogeneity in habitat selection by elk in Yellowstone National Park.  Ecoscience, 10(4):421–431.
  * Chatterjee et al., (2023) Chatterjee, N., Wolfson, D., Kim, D., Gomez, J. V., Freeman, S., Bacheler, N. M., Shertzer, K., Taylor, C., and Fieberg, J. (2023).  Modeling individual variability in habitat selection and movement using integrated step-selection analyses.  bioRxiv, pages 2023–07.
  * Cooper and Millspaugh, (1999) Cooper, A. B. and Millspaugh, J. J. (1999).  The application of discrete choice models to wildlife resource selection studies.  Ecology, 80(2):566–575.
  * Craiu et al., (2008) Craiu, R. V., Duchesne, T., and Fortin, D. (2008).  Inference methods for the conditional logistic regression model with longitudinal data.  Biometrical Journal, 50(1):97–109.
  * Duchesne et al., (2015) Duchesne, T., Fortin, D., and Rivest, L.-P. (2015).  Equivalence between step selection functions and biased correlated random walks for statistical inference on animal movement.  PloS One, 10(4):e0122947.
  * Fieberg, (2023) Fieberg, J. (2023).  Statistics for Ecologists: A Frequentist and Bayesian Treatment of Modern Regression Models.  An open-source online textbook.  https://statistics4ecologists-v1.netlify.app/.
  * Fieberg et al., (2021) Fieberg, J., Signer, J., Smith, B., and Avgar, T. (2021).  A ‘How to’ guide for interpreting parameters in habitat-selection analyses.  Journal of Animal Ecology, 90(5):1027–1043.
  * Forester et al., (2009) Forester, J. D., Im, H. K., and Rathouz, P. J. (2009).  Accounting for animal movement in estimation of resource selection functions: sampling and data analysis.  Ecology, 90(12):3554–3565.
  * Fortin et al., (2005) Fortin, D., Beyer, H. L., Boyce, M. S., Smith, D. W., Duchesne, T., and Mao, J. S. (2005).  Wolves influence elk movements: behavior shapes a trophic cascade in Yellowstone National Park.  Ecology, 86(5):1320–1330.
  * Gourieroux and Monfort, (1996) Gourieroux, C. and Monfort, A. (1996).  Simulation-based econometric methods.  Oxford University Press.
  * Johnson et al., (2008) Johnson, D. S., Thomas, D. L., Ver Hoef, J. M., and Christ, A. (2008).  A general framework for the analysis of animal resource selection from telemetry data.  Biometrics, 64(3):968–976.
  * Klappstein et al., (2022) Klappstein, N. J., Potts, J. R., Michelot, T., Börger, L., Pilfold, N. W., Lewis, M. A., and Derocher, A. E. (2022).  Energy-based step selection analysis: modelling the energetic drivers of animal movement and habitat use.  Journal of Animal Ecology, 91(5):946–957.
  * Klappstein et al., (2023) Klappstein, N. J., Thomas, L., and Michelot, T. (2023).  Flexible hidden Markov models for behaviour-dependent habitat selection.  Movement Ecology, 11(1):30.
  * Matthiopolous et al., (2023) Matthiopolous, J., Fieberg, J., and Aarts, G. (2023).  Species-Habitat Associations: Spatial data, predictive models, and ecological insights, 2nd edition.  University of Minnesota Libraries Publishing, University of Minnesota.
  * McDonald, (2013) McDonald, T. L. (2013).  The point process use-availability or presence-only likelihood and comments on analysis.  Journal of Animal Ecology, 82(6):1174–1182.
  * Muff et al., (2020) Muff, S., Signer, J., and Fieberg, J. (2020).  Accounting for individual-specific variation in habitat-selection studies: Efficient estimation of mixed-effects models using bayesian or frequentist computation.  Journal of Animal Ecology, 89(1):80–92.
  * Nicosia et al., (2017) Nicosia, A., Duchesne, T., Rivest, L.-P., and Fortin, D. (2017).  A multi-state conditional logistic regression model for the analysis of animal movement.  The Annals of Applied Statistics, 11(3):1537 – 1560.
  * Northrup et al., (2013) Northrup, J. M., Hooten, M. B., Anderson, C. R., and Wittemyer, G. (2013).  Practical guidance on characterizing availability in resource selection functions under a use–availability design.  Ecology, 94(7):1456–1463.
  * Pinheiro and Chao, (2006) Pinheiro, J. C. and Chao, E. C. (2006).  Efficient Laplacian and adaptive Gaussian quadrature algorithms for multilevel generalized linear mixed models.  Journal of Computational and Graphical Statistics, 15(1):58–81.
  * Pohle et al., (2023) Pohle, J., Signer, J., Eccard, J. A., Dammhahn, M., and Schlägel, U. E. (2023).  How to account for behavioural states in step-selection analysis: a model comparison.  arXiv preprint arXiv:2304.12964.
  * Potts et al., (2014) Potts, J. R., Bastille-Rousseau, G., Murray, D. L., Schaefer, J. A., and Lewis, M. A. (2014).  Predicting local and non-local effects of resources on animal space use using a mechanistic step selection model.  Methods in Ecology and Evolution, 5(3):253–262.
  * Rhodes et al., (2005) Rhodes, J. R., McAlpine, C. A., Lunney, D., and Possingham, H. P. (2005).  A spatially explicit habitat selection model incorporating home range behavior.  Ecology, 86(5):1199–1205.
  * Rizzo, (2019) Rizzo, M. L. (2019).  Statistical computing with R.  CRC Press.
  * Robert and Casella, (2010) Robert, C. P. and Casella, G. (2010).  Introducing Monte Carlo methods with R, volume 18.  Springer.
  * Schlägel and Lewis, (2014) Schlägel, U. E. and Lewis, M. A. (2014).  Detecting effects of spatial memory and dynamic information on animal movement decisions.  Methods in Ecology and Evolution, 5(11):1236–1246.
  * Schlägel et al., (2017) Schlägel, U. E., Merrill, E. H., and Lewis, M. A. (2017).  Territory surveillance and prey management: Wolves keep track of space and time.  Ecology and Evolution, 7(20):8388–8405.
  * Signer et al., (2019) Signer, J., Fieberg, J., and Avgar, T. (2019).  Animal movement tools (amt): R package for managing tracking data and conducting habitat selection analyses.  Ecology and evolution, 9(2):880–890.
  * Simpson et al., (2016) Simpson, D., Illian, J. B., Lindgren, F., Sørbye, S. H., and Rue, H. (2016).  Going off grid: Computationally efficient inference for log-Gaussian Cox processes.  Biometrika, 103(1):49–70.
  * Therneau, (2023) Therneau, T. M. (2023).  A Package for Survival Analysis in R.  R package version 3.5-5.
  * Thurfjell et al., (2014) Thurfjell, H., Ciuti, S., and Boyce, M. S. (2014).  Applications of step-selection functions in ecology and conservation.  Movement Ecology, 2:1–12.
  * Warton and Shepherd, (2010) Warton, D. I. and Shepherd, L. C. (2010).  Poisson point process models solve the “pseudo-absence problem” for presence-only data in ecology.  The Annals of Applied Statistics, pages 1383–1402.