Fitting Step-Selection Functions with `amt`

Johannes Signer

2025年08月23日

About

This vignette briefly introduces how one can fit a Step-Selection Function (SSF) with the amt package. We will be using the example data of one red deer from northern Germany and one covariate: a forest cover map. For a more through discussion see also Fieberg et al. 2020¹ and supplement B.

Getting the data ready

First we load the required libraries and the relocation data (called deer)

 library(lubridate)
 library(amt)
 data("deer")
deer

## # A tibble: 826 ×ばつ 4
## x_ y_ t_ burst_
## * <dbl> <dbl> <dttm> <dbl>
## 1 4314068. 3445807. 2008年03月30日 00:01:47 1
## 2 4314053. 3445768. 2008年03月30日 06:00:54 1
## 3 4314105. 3445859. 2008年03月30日 12:01:47 1
## 4 4314044. 3445785. 2008年03月30日 18:01:24 1
## 5 4313015. 3445858. 2008年03月31日 00:01:23 1
## 6 4312860. 3445857. 2008年03月31日 06:01:45 1
## 7 4312854. 3445856. 2008年03月31日 12:01:11 1
## 8 4312858. 3445858. 2008年03月31日 18:01:55 1
## 9 4312745. 3445862. 2008年04月01日 00:01:24 1
## 10 4312651. 3446024. 2008年04月01日 06:00:54 1
## # i 816 more rows

In order to continue, we need a regular sampling rate. To check the current sampling rate, we use summarize_sampling_rate:

 summarize_sampling_rate(deer)

## # A tibble: 1 ×ばつ 9
## min q1 median mean q3 max sd n unit 
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 5.96 6.00 6.01 11.5 6.02 3924. 137. 825 hour

The median sampling rate is 6h, which is what we aimed for.

Next, we have to get the environmental covariates. A forest layer is included in the package. Note, that this a regular SpatRast.

sh_forest <- get_sh_forest()
sh_forest

## class : SpatRaster 
## size : 720, 751, 1 (nrow, ncol, nlyr)
## resolution : 25, 25 (x, y)
## extent : 4304725, 4323500, 3437725, 3455725 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## source(s) : memory
## name : forest 
## min value : 0 
## max value : 1

Prepare Data for SSF

Steps

Before fitting a SSF we have to do some data preparation. First, we change from a point representation to a step representation, using the function steps_by_burst, which in contrast to the steps function accounts for bursts.

ssf1 <- deer |> steps_by_burst()

Control/random steps

The generic function random_steps provides a methods for a track_xy*, where each observed step is paired with n_control control steps (i.e., steps that share the same starting location but have different turn angles and step lengths). The distributions for drawing step lengths and turning angles are usually obtained by fitting known parametric distribution to the observed step length and turn angles.

The function random_steps has seven arguments. For most use cases the defaults are just fine, but there might situation where the user wants to adjust some of the arguments. The arguments are:

x: This is the track_xy* for which the random steps are created. That is, for each step in x n_control random steps are created.
n_control: The number of random steps that should be created for each observed step.
sl_distr: This is the distribution of the step lengths. By default a gamma distribution is fit to the observed step lengths of the x. But any amt_distr is suitable here. ²
ta_distr: This is the turn angle distribution, with the default being a von Mises distribution.
rand_sl: These are the random step lengths, by default 1e5 random numbers from the distribution fitted in 3³.
rand_ta: These are the random turn angles, by default 1e4 random numbers from the distribution fitted in 4.
include_observed: This argument is by default TRUE and indicates if the observed steps should be included or not.

The default situation

In most situations the following code snippet should work⁴.

ssf1 <- ssf1 |> random_steps(n_control = 15)

A exponential distribution for step lengths

todo

Extract covariates

As a last step, we have to extract the covariates at the end point of each step. We can do this with extract_covariates.

ssf1 <- ssf1 |> extract_covariates(sh_forest)

Since the forest layers is coded as 1 = forest and 2 != forest, we create a factor with appropriate levels. We also calculate the log of the step length and the cosine of the turn angle, which we may use later for a integrated step selection function.

ssf1 <- ssf1 |> 
 mutate(forest = factor(forest, levels = 1:0, labels = c("forest", "non-forest")), 
 cos_ta = cos(ta_), 
 log_sl = log(sl_))

Fitting SSF

Now all pieces are there to fit a SSF. We will use fit_clogit, which is a wrapper around survival::clogit.

m0 <- ssf1 |> fit_clogit(case_ ~ forest + strata(step_id_))
m1 <- ssf1 |> fit_clogit(case_ ~ forest + forest:cos_ta + forest:log_sl + log_sl * cos_ta + strata(step_id_))
m2 <- ssf1 |> fit_clogit(case_ ~ forest + forest:cos_ta + forest:log_sl + log_sl + cos_ta + strata(step_id_))
 summary(m0)

## Call:
## coxph(formula = Surv(rep(1, 12096L), case_) ~ forest + strata(step_id_), 
## data = data, method = "exact")
## 
## n= 12096, number of events= 756 
## 
## coef exp(coef) se(coef) z Pr(>|z|) 
## forestnon-forest -0.5145 0.5978 0.1088 -4.727 2.28e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## exp(coef) exp(-coef) lower .95 upper .95
## forestnon-forest 0.5978 1.673 0.483 0.74
## 
## Concordance= 0.529 (se = 0.007 )
## Likelihood ratio test= 21.65 on 1 df, p=3e-06
## Wald test = 22.34 on 1 df, p=2e-06
## Score (logrank) test = 22.44 on 1 df, p=2e-06

 summary(m1)

## Call:
## coxph(formula = Surv(rep(1, 12096L), case_) ~ forest + forest:cos_ta + 
## forest:log_sl + log_sl * cos_ta + strata(step_id_), data = data, 
## method = "exact")
## 
## n= 12096, number of events= 756 
## 
## coef exp(coef) se(coef) z Pr(>|z|) 
## forestnon-forest 0.83934 2.31483 0.32799 2.559 0.010497 * 
## log_sl 0.17416 1.19025 0.04958 3.513 0.000443 ***
## cos_ta -0.20329 0.81604 0.20531 -0.990 0.322092 
## forestnon-forest:cos_ta -0.31159 0.73228 0.11769 -2.648 0.008106 ** 
## forestnon-forest:log_sl -0.25554 0.77450 0.05786 -4.416 1e-05 ***
## log_sl:cos_ta 0.01656 1.01669 0.03450 0.480 0.631313 
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## exp(coef) exp(-coef) lower .95 upper .95
## forestnon-forest 2.3148 0.4320 1.2171 4.4026
## log_sl 1.1902 0.8402 1.0800 1.3117
## cos_ta 0.8160 1.2254 0.5457 1.2203
## forestnon-forest:cos_ta 0.7323 1.3656 0.5814 0.9223
## forestnon-forest:log_sl 0.7745 1.2912 0.6915 0.8675
## log_sl:cos_ta 1.0167 0.9836 0.9502 1.0878
## 
## Concordance= 0.609 (se = 0.013 )
## Likelihood ratio test= 90.25 on 6 df, p=<2e-16
## Wald test = 88.38 on 6 df, p=<2e-16
## Score (logrank) test = 90.65 on 6 df, p=<2e-16

 summary(m2)

## Call:
## coxph(formula = Surv(rep(1, 12096L), case_) ~ forest + forest:cos_ta + 
## forest:log_sl + log_sl + cos_ta + strata(step_id_), data = data, 
## method = "exact")
## 
## n= 12096, number of events= 756 
## 
## coef exp(coef) se(coef) z Pr(>|z|) 
## forestnon-forest 0.84859 2.33635 0.32724 2.593 0.009510 ** 
## log_sl 0.17369 1.18969 0.04956 3.505 0.000457 ***
## cos_ta -0.11644 0.89008 0.09693 -1.201 0.229614 
## forestnon-forest:cos_ta -0.31497 0.72981 0.11750 -2.681 0.007348 ** 
## forestnon-forest:log_sl -0.25707 0.77331 0.05776 -4.451 8.56e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## exp(coef) exp(-coef) lower .95 upper .95
## forestnon-forest 2.3363 0.4280 1.2302 4.4370
## log_sl 1.1897 0.8406 1.0796 1.3111
## cos_ta 0.8901 1.1235 0.7361 1.0763
## forestnon-forest:cos_ta 0.7298 1.3702 0.5797 0.9188
## forestnon-forest:log_sl 0.7733 1.2931 0.6905 0.8660
## 
## Concordance= 0.609 (se = 0.013 )
## Likelihood ratio test= 90.02 on 5 df, p=<2e-16
## Wald test = 87.6 on 5 df, p=<2e-16
## Score (logrank) test = 89.63 on 5 df, p=<2e-16

Interpretation of coefficients

See the discussion in Fieberg et al 2021.

A note on piping

All steps described above, could easily be wrapped into one piped workflow:

m1 <- deer |> 
 steps_by_burst() |> random_steps(n = 15) |> 
 extract_covariates(sh_forest) |> 
 mutate(forest = factor(forest, levels = 1:0, labels = c("forest", "non-forest")), 
 cos_ta = cos(ta_), 
 log_sl = log(sl_)) |> 
 fit_clogit(case_ ~ forest + forest:cos_ta + forest:sl_ + sl_ * cos_ta + strata(step_id_))

 summary(m1)

## Call:
## coxph(formula = Surv(rep(1, 12096L), case_) ~ forest + forest:cos_ta + 
## forest:sl_ + sl_ * cos_ta + strata(step_id_), data = data, 
## method = "exact")
## 
## n= 12096, number of events= 756 
## 
## coef exp(coef) se(coef) z Pr(>|z|) 
## forestnon-forest -0.1036793 0.9015144 0.1407961 -0.736 0.46150 
## sl_ 0.0006758 1.0006760 0.0001538 4.394 1.11e-05 ***
## cos_ta -0.3355283 0.7149603 0.1098064 -3.056 0.00225 ** 
## forestnon-forest:cos_ta -0.2134928 0.8077580 0.1183237 -1.804 0.07118 . 
## forestnon-forest:sl_ -0.0009751 0.9990253 0.0001986 -4.909 9.14e-07 ***
## sl_:cos_ta 0.0003849 1.0003850 0.0001271 3.027 0.00247 ** 
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## exp(coef) exp(-coef) lower .95 upper .95
## forestnon-forest 0.9015 1.1092 0.6841 1.1880
## sl_ 1.0007 0.9993 1.0004 1.0010
## cos_ta 0.7150 1.3987 0.5765 0.8866
## forestnon-forest:cos_ta 0.8078 1.2380 0.6406 1.0186
## forestnon-forest:sl_ 0.9990 1.0010 0.9986 0.9994
## sl_:cos_ta 1.0004 0.9996 1.0001 1.0006
## 
## Concordance= 0.603 (se = 0.013 )
## Likelihood ratio test= 105.2 on 6 df, p=<2e-16
## Wald test = 107.8 on 6 df, p=<2e-16
## Score (logrank) test = 116.5 on 6 df, p=<2e-16

Session

sessioninfo::session_info()

## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.5.1 (2025年06月13日)
## os macOS Sequoia 15.6
## system aarch64, darwin20
## ui X11
## language (EN)
## collate C
## ctype en_US.UTF-8
## tz Europe/Berlin
## date 2025年08月23日
## pandoc 3.4 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
## quarto 1.6.40 @ /usr/local/bin/quarto
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## amt * 0.3.0.0 2025年08月23日 [1] local
## backports 1.5.0 2024年05月23日 [3] CRAN (R 4.5.0)
## boot 1.3-31 2024年08月28日 [3] CRAN (R 4.5.1)
## bslib 0.9.0 2025年01月30日 [3] CRAN (R 4.5.0)
## cachem 1.1.0 2024年05月16日 [3] CRAN (R 4.5.0)
## checkmate 2.3.3 2025年08月18日 [3] CRAN (R 4.5.0)
## circular 0.5-1 2024年08月29日 [3] CRAN (R 4.5.0)
## class 7.3-23 2025年01月01日 [3] CRAN (R 4.5.1)
## classInt 0.4-11 2025年01月08日 [3] CRAN (R 4.5.0)
## cli 3.6.5 2025年04月23日 [3] CRAN (R 4.5.0)
## codetools 0.2-20 2024年03月31日 [3] CRAN (R 4.5.1)
## data.table 1.17.8 2025年07月10日 [3] CRAN (R 4.5.0)
## DBI 1.2.3 2024年06月02日 [3] CRAN (R 4.5.0)
## digest 0.6.37 2024年08月19日 [3] CRAN (R 4.5.0)
## dplyr * 1.1.4 2023年11月17日 [3] CRAN (R 4.5.0)
## e1071 1.7-16 2024年09月16日 [3] CRAN (R 4.5.0)
## evaluate 1.0.4 2025年06月18日 [3] CRAN (R 4.5.0)
## farver 2.1.2 2024年05月13日 [3] CRAN (R 4.5.0)
## fastmap 1.2.0 2024年05月15日 [3] CRAN (R 4.5.0)
## fitdistrplus 1.2-4 2025年07月03日 [3] CRAN (R 4.5.0)
## generics 0.1.4 2025年05月09日 [3] CRAN (R 4.5.0)
## ggforce 0.5.0 2025年06月18日 [3] CRAN (R 4.5.0)
## ggplot2 * 3.5.2 2025年04月09日 [3] CRAN (R 4.5.0)
## ggraph * 2.2.1 2024年03月07日 [3] CRAN (R 4.5.0)
## ggrepel 0.9.6 2024年09月07日 [3] CRAN (R 4.5.0)
## glue 1.8.0 2024年09月30日 [3] CRAN (R 4.5.0)
## graphlayouts 1.2.2 2025年01月23日 [3] CRAN (R 4.5.0)
## gridExtra 2.3 2017年09月09日 [3] CRAN (R 4.5.0)
## gtable 0.3.6 2024年10月25日 [3] CRAN (R 4.5.0)
## htmltools 0.5.8.1 2024年04月04日 [3] CRAN (R 4.5.0)
## igraph 2.1.4 2025年01月23日 [3] CRAN (R 4.5.0)
## jquerylib 0.1.4 2021年04月26日 [3] CRAN (R 4.5.0)
## jsonlite 2.0.0 2025年03月27日 [3] CRAN (R 4.5.0)
## KernSmooth 2.23-26 2025年01月01日 [3] CRAN (R 4.5.1)
## knitr * 1.50 2025年03月16日 [3] CRAN (R 4.5.0)
## labeling 0.4.3 2023年08月29日 [3] CRAN (R 4.5.0)
## lattice 0.22-7 2025年04月02日 [3] CRAN (R 4.5.1)
## lifecycle 1.0.4 2023年11月07日 [3] CRAN (R 4.5.0)
## lubridate * 1.9.4 2024年12月08日 [3] CRAN (R 4.5.0)
## magrittr 2.0.3 2022年03月30日 [3] CRAN (R 4.5.0)
## MASS 7.3-65 2025年02月28日 [3] CRAN (R 4.5.1)
## Matrix 1.7-3 2025年03月11日 [3] CRAN (R 4.5.1)
## memoise 2.0.1 2021年11月26日 [3] CRAN (R 4.5.0)
## mvtnorm 1.3-3 2025年01月10日 [3] CRAN (R 4.5.0)
## pillar 1.11.0 2025年07月04日 [3] CRAN (R 4.5.0)
## pkgconfig 2.0.3 2019年09月22日 [3] CRAN (R 4.5.0)
## polyclip 1.10-7 2024年07月23日 [3] CRAN (R 4.5.0)
## proxy 0.4-27 2022年06月09日 [3] CRAN (R 4.5.0)
## purrr 1.1.0 2025年07月10日 [3] CRAN (R 4.5.0)
## R6 2.6.1 2025年02月15日 [3] CRAN (R 4.5.0)
## rbibutils 2.3 2024年10月04日 [3] CRAN (R 4.5.0)
## RColorBrewer 1.1-3 2022年04月03日 [3] CRAN (R 4.5.0)
## Rcpp 1.1.0 2025年07月02日 [3] CRAN (R 4.5.0)
## Rdpack 2.6.4 2025年04月09日 [3] CRAN (R 4.5.0)
## rlang 1.1.6 2025年04月11日 [3] CRAN (R 4.5.0)
## rmarkdown 2.29 2024年11月04日 [3] CRAN (R 4.5.0)
## rstudioapi 0.17.1 2024年10月22日 [3] CRAN (R 4.5.0)
## sass 0.4.10 2025年04月11日 [3] CRAN (R 4.5.0)
## scales 1.4.0 2025年04月24日 [3] CRAN (R 4.5.0)
## sessioninfo 1.2.3 2025年02月05日 [3] CRAN (R 4.5.0)
## sf 1.0-21 2025年05月15日 [3] CRAN (R 4.5.0)
## survival 3.8-3 2024年12月17日 [3] CRAN (R 4.5.1)
## terra 1.8-60 2025年07月21日 [3] CRAN (R 4.5.0)
## tibble 3.3.0 2025年06月08日 [3] CRAN (R 4.5.0)
## tidygraph * 1.3.1 2024年01月30日 [3] CRAN (R 4.5.0)
## tidyr 1.3.1 2024年01月24日 [3] CRAN (R 4.5.0)
## tidyselect 1.2.1 2024年03月11日 [3] CRAN (R 4.5.0)
## timechange 0.3.0 2024年01月18日 [3] CRAN (R 4.5.0)
## tweenr 2.0.3 2024年02月26日 [3] CRAN (R 4.5.0)
## units 0.8-7 2025年03月11日 [3] CRAN (R 4.5.0)
## utf8 1.2.6 2025年06月08日 [3] CRAN (R 4.5.0)
## vctrs 0.6.5 2023年12月01日 [3] CRAN (R 4.5.0)
## viridis 0.6.5 2024年01月29日 [3] CRAN (R 4.5.0)
## viridisLite 0.4.2 2023年05月02日 [3] CRAN (R 4.5.0)
## withr 3.0.2 2024年10月28日 [3] CRAN (R 4.5.0)
## xfun 0.53 2025年08月19日 [3] CRAN (R 4.5.0)
## yaml 2.3.10 2024年07月26日 [3] CRAN (R 4.5.0)
## 
## [1] /private/var/folders/ln/h3zng0fs2pq7mhn_hzn0d8x00000gn/T/RtmpvmiLWg/Rinst1761f5e417fe9
## [2] /private/var/folders/ln/h3zng0fs2pq7mhn_hzn0d8x00000gn/T/RtmpfKQTmf/temp_libpath16445556b59c1
## [3] /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library
## * ── Packages attached to the search path.
## 
## ──────────────────────────────────────────────────────────────────────────────

https://www.biorxiv.org/content/10.1101/2020.11.12.379834v4 ↩︎
See also ?fit_distr.↩︎
Note, this possible because of the Glivenko-Cantelli theorem and works as long as the sample from the original distribution (the sample you provide here) is large enough.↩︎
And how it was implemented in amt up to version 0.0.6. This should be backward compatible and not break existing code.↩︎

Fitting Step-Selection Functions with amt