In this appendix, we will walk the user through the process of fitting and interpreting parameters and output when conducting an *integrated step-selection analysis*. Step-selection functions model transitions or “steps” connecting sequential locations (\(\Delta t\) units apart) in geographical space:

\[\begin{equation} u(s,t+\Delta t) | u(s', t) = \frac{w(X(s); \beta(\Delta t))\phi(s, s'; \gamma(\Delta t))}{\int_{\tilde{s}\in G}w(X(\tilde{s}, s'); \beta(\Delta t))\phi(\tilde{s},s'; \gamma(\Delta t))d\tilde{s}} \end{equation}\]

where \(u(s, t+\Delta t) | u(s', t)\) gives the conditional probability of finding the individual at location \(s\) at time \(t+\Delta t\) given it was at location \(s'\) at time \(t\), \(w(X(s); \beta(\Delta t))\) is referred to as a step-selection function, and \(\phi(s,s'; \gamma(\Delta t))\) is a selection-independent movement kernel that describes how the animal would move in homogeneous habitat or in the absence of habitat selection (i.e., when \(w(X(s); \beta(\Delta t))=\) a constant for all \(s\)).

Integrated step-selection analyses are implemented using the following steps:

- Estimate a tentative selection-free movement kernel, \(\phi(s,s'; \gamma(\Delta t))\), using observed step-lengths and turn angles, giving \(\hat{\phi}(s,s'; \hat{\gamma}(\Delta t))\).
- Generate time-dependent available locations by simulating potential movements from the previously observed location: \(a(s, t+\Delta t) | u(s', t) = \phi(s,s'; \gamma(\Delta t))\).
- Estimate \(\beta\) using conditional logistic regression, with strata formed by combining time-dependent used and available locations.
- Re-estimate the movement parameters in \(\phi(s,s'; \gamma(\Delta t))\) using regression coefficients associated with movement characteristics (e.g., step-length, log step-length, and cosine of the turn angle). This step adjusts the parameters in \(\phi(s,s'; \gamma(\Delta t))\) to account for the effect of habitat selection when estimating the movement kernel (Avgar, Potts, Lewis, & Boyce, 2016).

The first two steps are accomplished in the **Format and Generate Random Steps** and **Tentative Movement Parameters** sections. The third step is accomplished in the **Fit iSSF** section. The last step is accomplished in the **Interpreting Movement Parameters** section.

We will begin by preparing our data to allow us to estimate parameters in the movement kernel, \(\phi(s,s';\gamma(\Delta t))\) and the habitat-selection function, \(w(X(s); \beta(\Delta t))\). This step will involve loading our location data and covariates, translating the locations into a series of regular steps (transitions between sequential locations spaced \(\Delta t\) units apart), sampling random available steps from a tentative statistical distribution formulated in terms of step-length distribution and turn-angle distribution, and associating the habitat covariates with the steps.

First, let’s read the fisher location and spatial covariate data. These data are available as objects in the `amt`

package (`amt_fisher`

and `env_covar`

), and we access them easily if the package is attached. The object `amt_fisher`

is already of class `track_xyt`

from the package `amt`

, which we have already seen in the RSF examples (see Supplementary Appendix A). The other three objects are all of class `RasterLayer`

, which can be created by reading most standard raster files with the function `raster::raster()`

. For ease, we have already loaded them, but your own data will likely need to be read in using functions from other packages.

```
fisher <- amt_fisher
landuse <- amt_fisher_covar$landuse
elevation <- amt_fisher_covar$elevation
popden <- amt_fisher_covar$popden
```

In this example, we will work with data from a single individual, “Lupe”, so let’s subset the `fisher`

data and sort the locations using the timestamp of the observations, `t_`

.

```
dat <- fisher %>%
filter(name == "Lupe") %>%
arrange(t_)
```

The raw land-use data we loaded here has numerical codes, which we will group into different landcover types using a custom-written function, `reclass_landuse()`

, below. We will apply this function when preparing the data.

```
reclass_landuse <- function(x) {
fct_collapse(factor(x),
forest = c("30","40","50","60", "70","80", "90","100"),
grass = c("120", "140"),
wet = c("160"))
}
```

Most step-selection analyses assume data are collected at regular time steps (for an exception, see Munden et al., 2020). Although your device may be programmed to collect data at regular time steps, there will typically be some small amount of variation in the time between locations. If we consider that variability to be negligible, we can round it off by resampling the track.

Check the actual sampling rate:

`summarize_sampling_rate(dat)`

```
## # A tibble: 1 x 9
## min q1 median mean q3 max sd n unit
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <chr>
## 1 0.400 1.97 2.03 8.93 3.98 1080. 48.1 3003 min
```

We can see summary statistics for the difference between locations in minutes (see the final column for the unit). These tags were programmed to take a location every two minutes, but there is some variation in the time between locations, and missed locations lead to some gaps in the data. We will therefore need to resample the track when we format the data for this analysis.

We will take advantage of `amt`

’s piped (`%>%`

) workflow here to prepare data for model fitting. This workflow will help make our code both readable and compact. Most of the functions from `amt`

have sensible defaults for each argument to get you started. You can easily see these defaults by checking the help file for that function (e.g., `?random_steps`

) or change them as you see fit.

In this workflow, we begin with Lupe’s location data (`dat`

), we resample the track to get regular time steps (`track_resample()`

), convert from points to steps (`steps_by_burst()`

), generate random available steps (`random_steps()`

), and attach the environmental covariates to both the observed and available steps (`extract_covariates()`

). Finally, we’ll use `dplyr::mutate()`

to scale and center our covariates and add the additional movement covariates that we will need to conduct an integrated step-selection analysis.

```
ssf_dat <- dat %>%
track_resample(rate = minutes(2), tolerance = seconds(20)) %>%
steps_by_burst() %>%
random_steps() %>%
extract_covariates(landuse, where = "both") %>%
extract_covariates(elevation, where = "both") %>%
extract_covariates(popden, where = "both") %>%
mutate(elevation_start = scale(elevation_start),
elevation_end = scale(elevation_end),
popden_start = scale(popden_start),
popden_end = scale(popden_end),
landuseC_start = reclass_landuse(landuse_start),
landuseC_end = reclass_landuse(landuse_end),
forest_start = landuseC_start == "forest",
forest_end = landuseC_end == "forest",
cos_ta_ = cos(ta_),
log_sl_ = log(sl_)
) %>%
filter(!is.na(ta_))
```

Before we proceed, we make a few notes about the code above:

- by default the
`random_steps()`

function fits a tentative gamma distribution to the observed step lengths and a tentative von Mises distribution to the observed turn angles. It then generates random available points by sampling step-lengths and turn angles from these fitted distributions and combining these random steps with the starting locations associated with each observed movement step. - by including the argument
`where = "both"`

to the`extract_covariates()`

function, we are able to extract covariates associated with both the starting and ending locations of each movement step. It is often sensible to use covariates associated with the starting location to model the effect of habitat on the movement kernel, \(\phi(s,s';\gamma(\Delta t)\), and covariates associated with the ending location to model the habitat-selection function, \(w(X(s); \beta(\Delta t))\) (Signer, Fieberg, & Avgar, 2019).

Here’s what the resulting data set looks like now. Each row is a single step with all of its associated attributes:

`burst_`

is a unique identifier associated with a set of consecutive steps that are spaced \(\Delta t\) units apart.`x1_`

and`y1_`

are the coordinates associated with the starting location of the step.`x2_`

and`y2_`

are the coordinates associated with the ending location of the step.`sl_`

and`ta_`

are the step-length and turn angle associated with the step.`t1_`

and`t2_`

are the timestamps associated with the start and end of the step.`dt_`

is the time duration associated with the step.`step_id_`

is a unique identifier associated with each step’s observed and random locations.`case_`

is an indicator variable equal to`TRUE`

for observed steps and`FALSE`

for random steps.- (
`landuse_start`

,`elevation_start`

,`popden_start`

) and (`landuse_end`

,`elevation_end`

,`popden_end`

) are covariate values for (landuse categories, elevation, and population density) associated with the start and end of the movement step, respectively. `cos_ta_`

and`log_sl`

capture the cosine of the turn angle and the log of the step-length. These covariates will allow us to adjust/refine the parameters of our tentative step-length and turn-angle distributions after fitting our integrated step-selection model.

`print(ssf_dat, n = 3, width = Inf)`

```
## # A tibble: 16,175 x 25
## burst_ x1_ x2_ y1_ y2_ sl_ direction_p ta_
## * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 6 1783660. 1783673. 2402636. 2402664. 31.2 2.90 -1.75
## 2 6 1783660. 1783659. 2402636. 2402633. 3.36 2.90 1.48
## 3 6 1783660. 1783654. 2402636. 2402631. 6.92 2.90 0.907
## t1_ t2_ dt_ step_id_ case_ landuse_start
## * <dttm> <dttm> <drtn> <int> <lgl> <dbl>
## 1 2010-12-17 09:05:14 2010-12-17 09:07:04 1.83 mins 1 FALSE 50
## 2 2010-12-17 09:05:14 2010-12-17 09:07:04 1.83 mins 1 FALSE 50
## 3 2010-12-17 09:05:14 2010-12-17 09:07:04 1.83 mins 1 FALSE 50
## landuse_end elevation_start[,1] elevation_end[,1] popden_start[,1]
## * <dbl> <dbl> <dbl> <dbl>
## 1 50 -0.275 -0.810 -0.754
## 2 50 -0.275 -0.230 -0.754
## 3 50 -0.275 -0.230 -0.754
## popden_end[,1] landuseC_start landuseC_end forest_start forest_end cos_ta_
## * <dbl> <fct> <fct> <lgl> <lgl> <dbl>
## 1 -0.763 forest forest TRUE TRUE -0.181
## 2 -0.763 forest forest TRUE TRUE 0.0862
## 3 -0.763 forest forest TRUE TRUE 0.616
## log_sl_
## * <dbl>
## 1 3.44
## 2 1.21
## 3 1.93
## # ... with 16,172 more rows
```

Let’s have a closer look at what is happening under the hood when applying the `random_steps()`

function. This function conveniently fits tentative step-length and turn-angle distributions and then samples from these distributions to generate random steps. The default arguments lead to `random_steps()`

fitting a gamma distribution for the step lengths and a von Mises distribution for the turn angles.

The tentative parameters in these statistical distributions (gamma, von Mises) are stored as attributes of the resulting object. We can view the parameters of the tentative step-length distribution using:

`sl_distr(ssf_dat)`

```
## $name
## [1] "gamma"
##
## $params
## $params$shape
## [1] 1.42
##
## $params$scale
## [1] 35.9
##
##
## attr(,"class")
## [1] "gamma_distr" "sl_distr" "amt_distr" "list"
```

Similarly, we can access the parameters of the tentative turn-angle distribution using:

`ta_distr(ssf_dat)`

```
## $name
## [1] "vonmises"
##
## $params
## $params$kappa
## [1] 0.484
##
## $params$mu
## [1] 0
##
##
## attr(,"class")
## [1] "vonmises_distr" "ta_distr" "amt_distr" "list"
```

After we fit the iSSF, we will see how we can update these distributions to adjust for the influence of habitat selection when estimating parameters of the movement kernel.

In this section, we will fit a basic iSSF, demonstrating how to interpret habitat effects and movement parameters. In the next section, we will fit more complex models and show how those models differ from the basic approach presented here.

We will fit a model with the same habitat covariates introduced in the main text. Suppose we hypothesize that habitat selection at the step scale depends on human population density (`popden`

), elevation (`elevation`

), and land-use category (`landuseC`

). Furthermore, suppose we hypothesize that, in the absence of habitat selection, step lengths and turn angles follow a constant distribution, i.e., the selection-free movement kernel does not depend on covariates. We can fit that model as follows:

```
m1 <- ssf_dat %>%
fit_issf(case_ ~ popden_end + elevation_end + landuseC_end +
sl_ + log_sl_ + cos_ta_ +
strata(step_id_), model = TRUE)
```

Notice that our model formula includes habitat covariates (`popden`

, `elevation`

, and `landuseC`

) and movement covariates (`sl_`

, `log_sl_`

, and `cos_ta_`

), as well as an indicator for the column that indicates the stratum, i.e., the grouping that matches each observed step to its corresponding available steps (`step_id_`

). The extra argument, `model = TRUE`

, is passed along to `survival::clogit()`

to tell it to store the original model data in the resulting object (`m1`

). We will need this extra information later when we use our fitted model to generate predictions.

We begin by exploring the coefficients in the fitted model using the `summary()`

function:

`summary(m1)`

```
## Call:
## coxph(formula = Surv(rep(1, 16175L), case_) ~ popden_end + elevation_end +
## landuseC_end + sl_ + log_sl_ + cos_ta_ + strata(step_id_),
## data = data, model = TRUE, method = "exact")
##
## n= 16175, number of events= 1155
##
## coef exp(coef) se(coef) z Pr(>|z|)
## popden_end -0.107806 0.897802 0.106930 -1.01 0.3134
## elevation_end 0.169540 1.184760 0.059918 2.83 0.0047 **
## landuseC_endgrass -1.195003 0.302703 0.510012 -2.34 0.0191 *
## landuseC_endwet -0.047474 0.953635 0.369706 -0.13 0.8978
## sl_ 0.000784 1.000784 0.001306 0.60 0.5483
## log_sl_ 0.023037 1.023305 0.057365 0.40 0.6880
## cos_ta_ 0.008637 1.008675 0.045452 0.19 0.8493
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## popden_end 0.898 1.114 0.728 1.107
## elevation_end 1.185 0.844 1.053 1.332
## landuseC_endgrass 0.303 3.304 0.111 0.823
## landuseC_endwet 0.954 1.049 0.462 1.968
## sl_ 1.001 0.999 0.998 1.003
## log_sl_ 1.023 0.977 0.914 1.145
## cos_ta_ 1.009 0.991 0.923 1.103
##
## Concordance= 0.522 (se = 0.01 )
## Likelihood ratio test= 18.3 on 7 df, p=0.01
## Wald test = 17.8 on 7 df, p=0.01
## Score (logrank) test = 18 on 7 df, p=0.01
```

For now, we will just focus on the habitat-selection coefficients, i.e., `popden`

, `elevation`

, and the two `landuseC`

coefficients. The interpretation of the parameters is very similar to the RSF interpretation (Supplementary Appendix A) but at a local scale with habitat availability determined by the selection-free movement kernel. We can exponentiate habitat-selection parameters in a fitted step-selection model to compare the relative rates of use of two locations that differ by 1 unit of the explanatory variable but are otherwise equivalent - i.e., they should be equally accessible and have identical values for all other explanatory variables.

Calculating the relative use of location \(s_1\) versus location \(s_2\) is fairly straightforward when \(s_1\) and \(s_2\) share the same values for all but one covariate (see e.g., Supplementary Appendix A). For more complex scenarios, we have implemented a function in `amt`

that will calculate the log-relative intensity (referred to as the log-Relative Selection Strength, or log-RSS, by Avgar, Lele, Keim, & Boyce, 2017) by leveraging the generic `predict()`

function available for many ordinary model classes (e.g., `lm()`

or `glm()`

) models in `R`

. Regardless of whether there are interactions or polynomial terms, the function will correctly calculate log-RSS for you. If you prefer to quantify relative-use (i.e., RSS), you can simply exponentiate the results.

Here we demonstrate the use of`log_rss()`

with the simple example presented in equation 8 of the main text. Suppose we have two locations with the following covariates:

- \(s_1\):
`elevation_end = 3`

,`popden_end = 1.5`

,`landuseC_end = "wet"`

- \(s_2\):
`elevation_end = 2`

,`popden_end = 1.5`

,`landuseC_end = "wet"`

We want to calculate the log-RSS for a step ending in \(s_1\) versus a step ending in \(s_2\), assuming \(s_1\) and \(s_2\) are equally accessible to the animal.

The function `log_rss()`

expects the two locations to be formatted as separate `data.frame`

objects. Each `data.frame`

must include all covariates used to fit the model. Furthermore, `factor`

variables should have the same `levels`

as the original data.

```
# data.frame for s1
s1 <- data.frame(
elevation_end = 3,
popden_end = 1.5,
landuseC_end = factor("wet", levels = levels(ssf_dat$landuseC_end)),
sl_ = 100,
log_sl_ = log(100),
cos_ta_ = 1)
# data.frame for s2
s2 <- data.frame(
elevation_end = 2,
popden_end = 1.5,
landuseC_end = factor("wet", levels = levels(ssf_dat$landuseC_end)),
sl_ = 100,
log_sl_ = log(100),
cos_ta_ = 1)
```

Now that we have specified each location as a `data.frame`

, we can pass them along to `log_rss()`

for the calculation. The function will return an object of class `log_rss`

, which is also more generally a `list`

. The `list`

element `"df"`

contains a `data.frame`

which contains the log-RSS calculation and could easily be used to make a plot when considering relative selection strength across a range of environmental characteristics.

```
lr1 <- log_rss(m1, x1 = s1, x2 = s2)
lr1$df
```

```
## elevation_end_x1 popden_end_x1 landuseC_end_x1 sl_x1 log_sl_x1 cos_ta_x1
## 1 3 1.5 wet 100 4.61 1
## log_rss
## 1 0.17
```

As we will see below, this function is designed to be able to consider several locations as `x1`

, relative to a single location in `x2`

. Because the `[["df"]]`

element of the `log_rss`

object is designed for plotting, it includes the value of all the covariates at location `x1`

by appending `_x1`

to the end of the predictor names. The `log_rss`

object also contains the original `x1`

and `x2`

as elements, along with the model formula.

We can see that the log-RSS for this example is 0.17, which is equivalent to \(\beta_{elevation}\) since we have considered the simple case in which \(s_1\) and \(s_2\) differ by only a single unit of elevation. Now that we know how the function works, we can consider more complex cases.

Suppose we are interested in making a figure showing the relative selection strength our fisher using a range of elevations. The function `log_rss()`

can accept a `data.frame`

with many rows as as its first argument, `x1`

, to go with a single comparison location, `x2`

(i.e., `x2`

must only have 1 row).

```
# Make a new data.frame for s1
s1 <- data.frame(
elevation_end = seq(from = -2, to = 2, length.out = 200),
popden_end = 1.5,
landuseC_end = factor("wet", levels = levels(ssf_dat$landuseC_end)),
sl_ = 100,
log_sl_ = log(100),
cos_ta_ = 1)
# data.frame for s2
s2 <- data.frame(
elevation_end = 0, # mean of elev, since we scaled and centered
popden_end = 1.5,
landuseC_end = factor("wet", levels = levels(ssf_dat$landuseC_end)),
sl_ = 100,
log_sl_ = log(100),
cos_ta_ = 1)
# Calculate log-RSS
lr2 <- log_rss(m1, s1, s2)
# Plot using ggplot2
ggplot(lr2$df, aes(x = elevation_end_x1, y = log_rss)) +
geom_line(size = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
xlab("Elevation (SD)") +
ylab("log-RSS vs Mean Elevation") +
theme_bw()
```

The line depicts the log-RSS between each value of \(s_1\) relative to \(s_2\). Note that these locations differ only in their values of `elevation`

, and the log-RSS is equal to 0 when the two locations are identical (i.e., when `elevation`

= 0). If you prefer to plot RSS to log-RSS, simply exponentiate the y-axis (and move your horizontal line for no selection to \(exp(0) = 1\)).

```
# Plot using ggplot2
ggplot(lr2$df, aes(x = elevation_end_x1, y = exp(log_rss))) +
geom_line(size = 1) +
geom_hline(yintercept = exp(0), linetype = "dashed", color = "gray30") +
xlab("Elevation (SD)") +
ylab("RSS vs Mean Elevation") +
theme_bw()
```

The RSS plots above do not give any indication of model uncertainty. One method for constructing confidence intervals associated with the RSS predictions is via non-parametric bootstrapping (Fieberg, Vitense, & Johnson, 2020). When using the argument `ci="boot"`

to the `log_rss()`

function, the function repeatedly resamples strata with replacement, refits the model, and calculates log-RSS. Quantiles of the bootstrap log-RSS statistics are used to construct, by default, a 95% confidence interval around the original log-RSS calculation. The size of the confidence interval and number of bootstrap iterations can be controlled by the user by changing the `ci_level`

and `n_boot`

arguments. Warning: this code can take some time to execute.

```
# Calculate log-RSS
lr2_ci_boot <- log_rss(m1, s1, s2, ci = "boot", ci_level = 0.95, n_boot = 1000)
```

```
## Generating bootstrapped confidence intervals...
## Expect this to take some time...
## ... finished bootstrapping.
```

```
# Check the header of the data.frame
head(lr2_ci_boot$df)
```

```
## elevation_end_x1 popden_end_x1 landuseC_end_x1 sl_x1 log_sl_x1 cos_ta_x1
## 1 -2.00 1.5 wet 100 4.61 1
## 2 -1.98 1.5 wet 100 4.61 1
## 3 -1.96 1.5 wet 100 4.61 1
## 4 -1.94 1.5 wet 100 4.61 1
## 5 -1.92 1.5 wet 100 4.61 1
## 6 -1.90 1.5 wet 100 4.61 1
## log_rss lwr upr
## 1 -0.339 -0.578 -0.0755
## 2 -0.336 -0.572 -0.0747
## 3 -0.332 -0.566 -0.0739
## 4 -0.329 -0.560 -0.0732
## 5 -0.325 -0.555 -0.0724
## 6 -0.322 -0.549 -0.0717
```

You can see that the resulting `lr2_ci_boot`

data.frame now has columns `lwr`

and `upr`

, indicating the lower and upper bounds of the confidence interval. We can now add this confidence envelope to our plot.

```
ggplot(lr2_ci_boot$df, aes(x = elevation_end_x1, y = log_rss)) +
geom_ribbon(aes(ymin = lwr, ymax = upr),
linetype = "dashed",
color = "black", fill = "gray80", alpha = 0.5) +
geom_line(size = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
xlab("Elevation (SD)") +
ylab("log-RSS vs Mean Elevation") +
theme_bw()
```

To avoid the computational time required to calculate bootstrap confidence intervals, you may decide that using a large-sample-based confidence interval is appropriate. We can use the standard errors from our fitted model to estimate these confidence intervals based on a normal approximation to the sampling distribution of \(\hat{\beta}\).

```
lr2_ci_se <- log_rss(m1, s1, s2, ci = "se", ci_level = 0.95)
# Check the header of the data.frame
head(lr2_ci_se$df)
```

```
## elevation_end_x1 popden_end_x1 landuseC_end_x1 sl_x1 log_sl_x1 cos_ta_x1
## 1 -2.00 1.5 wet 100 4.61 1
## 2 -1.98 1.5 wet 100 4.61 1
## 3 -1.96 1.5 wet 100 4.61 1
## 4 -1.94 1.5 wet 100 4.61 1
## 5 -1.92 1.5 wet 100 4.61 1
## 6 -1.90 1.5 wet 100 4.61 1
## log_rss lwr upr
## 1 -0.339 -0.574 -0.104
## 2 -0.336 -0.568 -0.103
## 3 -0.332 -0.562 -0.102
## 4 -0.329 -0.557 -0.101
## 5 -0.325 -0.551 -0.100
## 6 -0.322 -0.545 -0.099
```

We see, in this case, that the two confidence intervals (bootstrap, normal-based) lead to similar conclusions.

```
# Add "type" label to both predictions data.frames
lr2_ci_boot$df$type <- "Bootstrap"
lr2_ci_se$df$type <- "SE"
# Combine two prediction data.frames
lr2_ci_both <- rbind(lr2_ci_boot$df, lr2_ci_se$df)
# Plot
ggplot(lr2_ci_both, aes(x = elevation_end_x1, y = log_rss,
group = factor(type),
fill = factor(type))) +
geom_ribbon(aes(ymin = lwr, ymax = upr),
linetype = "dashed", color = "black", alpha = 0.25) +
geom_line(size = 1, color = "black") +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray30") +
scale_fill_manual(name = "CI Type", breaks = c("Bootstrap", "SE"),
values = c("blue", "orange")) +
xlab("Elevation (SD)") +
ylab("log-RSS vs Mean Elevation") +
theme_bw()
```

Previously, we viewed tentative estimates of parameters describing the distribution of step-lengths and turn-angles. When estimating these parameters, however, we did not account or adjust for the influence of habitat selection. Including movement characteristics (`sl_`

, `log_sl_`

, and `cos_ta_`

) in our iSSF allows us to update our estimates of these parameters, which describe selection-free movement distributions. Mathematical formulas required to perform these updates are given in Supplementary Appendix C. Here, we will see how to update the distributions using the functions from `amt`

.

In our model, we included parameters for the step length (`sl_`

) and its natural logarithm (`log_sl_`

). We included these parameters to update the scale and shape parameters of our tentative gamma distribution, respectively. Different terms will need to be included in the iSSF, depending on the assumed step-length distribution (for details, see Supplementary Appendix C).

For our basic model – where the movement parameters are *not* interacting with any covariates – `amt`

has a simple function to update the step-length distribution.

```
# Update the distribution
updated_sl <- update_sl_distr(m1)
# Check the class
class(updated_sl)
```

`## [1] "gamma_distr" "sl_distr" "amt_distr" "list"`

```
# Print
print(updated_sl)
```

```
## $name
## [1] "gamma"
##
## $params
## $params$shape
## [1] 1.44
##
## $params$scale
## [1] 37
##
##
## attr(,"class")
## [1] "gamma_distr" "sl_distr" "amt_distr" "list"
```

We can see that the updated distribution object has classes that are designed for working with step-length/turn-angle distributions in `amt`

. In the most general sense, it is a `list`

, so you can access its named components either with `$name`

or `[["name"]]`

, as you would with any other list. But it is also of class `amt_distr`

, a general class for working with distributions in `amt`

; class `sl_distr`

, a more specific class for working with step-length distributions; and class `gamma_distr`

, an even more specific class for working with gamma distributions. We can see that the object has a name (`"gamma"`

) and the appropriate parameters for the gamma distribution – shape (`$params$shape`

) and scale (`$params$scale`

).

We can compare the tentative and updated distributions to see how accounting for habitat-selection changes our estimates of the selection-free step-length distribution.

```
# data.frame for plotting
plot_sl <- data.frame(x = rep(NA, 100))
# x-axis is sequence of possible step lengths
plot_sl$x <- seq(from = 0, to = 400, length.out = 100)
# y-axis is the probability density under the given gamma distribution
# For the tentative distribution
plot_sl$tentative <- dgamma(
x = plot_sl$x,
shape = m1$sl_$params$shape,
scale = m1$sl_$params$scale)
# For the updated distribution
plot_sl$updated <- dgamma(
x = plot_sl$x,
shape = updated_sl$params$shape,
scale = updated_sl$params$scale)
# Pivot from wide data to long data
plot_sl <- plot_sl %>%
pivot_longer(cols = -x)
# Plot
ggplot(plot_sl, aes(x = x, y = value, color = factor(name))) +
geom_line(size = 1) +
xlab("Step Length (m)") +
ylab("Probability Density") +
scale_color_manual(name = "Distribution",
breaks = c("tentative", "updated"),
values = c("blue", "orange")) +
theme_bw()
```