# Fitting Step-Selection Functions with `amt`

#### Johannes Signer

#### 2018-04-21

## 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.

## Getting the data ready

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

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

    ## # A tibble: 826 x 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.
    ## # ... with 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 x 9
    ##   min              q1      median   mean    q3     max      sd     n unit
    ##   <S3: table>      <S3: t> <S3: ta> <S3: t> <S3: > <S3:> <dbl> <int> <chr>
    ## 1 5.96444444444444 5.9966… 6.00583… 11.461… 6.016… 3923…  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 `RasterLayer`.

    data("sh_forest")
    sh_forest

    ## class       : RasterLayer
    ## dimensions  : 720, 751, 540720  (nrow, ncol, ncell)
    ## 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
    ## data source : in memory
    ## names       : sh.forest
    ## values      : 1, 2  (min, max)

## Prepare Data for SSF

Before fitting a step selection, the data well need to prepared. 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()

Next, we generate random steps with the function `random_steps`. This function fits by default a Gamma distribution to the step lengths and a von Mises distribution to the turn angles, and then pairs each observed step with `n` random steps.

    ssf1 <- ssf1 %>% random_steps(n = 15)

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(sh.forest, levels = 1:2, 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, 12144L), case_) ~ forest + strata(step_id_),
    ##     data = data, method = "exact")
    ##
    ##   n= 12144, number of events= 759
    ##
    ##                     coef exp(coef) se(coef)      z Pr(>|z|)
    ## forestnon-forest -0.6473    0.5235   0.1095 -5.913 3.37e-09 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ##                  exp(coef) exp(-coef) lower .95 upper .95
    ## forestnon-forest    0.5235       1.91    0.4224    0.6487
    ##
    ## Rsquare= 0.003   (max possible= 0.293 )
    ## Likelihood ratio test= 33.68  on 1 df,   p=6.502e-09
    ## Wald test            = 34.96  on 1 df,   p=3.369e-09
    ## Score (logrank) test = 35.19  on 1 df,   p=2.997e-09

    summary(m1)

    ## Call:
    ## coxph(formula = Surv(rep(1, 12144L), case_) ~ forest + forest:cos_ta +
    ##     forest:log_sl + log_sl * cos_ta + strata(step_id_), data = data,
    ##     method = "exact")
    ##
    ##   n= 12144, number of events= 759
    ##
    ##                             coef exp(coef) se(coef)      z Pr(>|z|)
    ## forestnon-forest         0.60922   1.83900  0.33260  1.832 0.067000 .
    ## log_sl                   0.17732   1.19402  0.04949  3.583 0.000340 ***
    ## cos_ta                   0.15503   1.16770  0.20813  0.745 0.456354
    ## forestnon-forest:cos_ta -0.40263   0.66856  0.11816 -3.407 0.000656 ***
    ## forestnon-forest:log_sl -0.23723   0.78881  0.05830 -4.069 4.71e-05 ***
    ## log_sl:cos_ta            0.01408   1.01418  0.03491  0.403 0.686729
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ##                         exp(coef) exp(-coef) lower .95 upper .95
    ## forestnon-forest           1.8390     0.5438    0.9582    3.5294
    ## log_sl                     1.1940     0.8375    1.0836    1.3156
    ## cos_ta                     1.1677     0.8564    0.7765    1.7559
    ## forestnon-forest:cos_ta    0.6686     1.4958    0.5303    0.8428
    ## forestnon-forest:log_sl    0.7888     1.2677    0.7036    0.8843
    ## log_sl:cos_ta              1.0142     0.9860    0.9471    1.0860
    ##
    ## Rsquare= 0.005   (max possible= 0.293 )
    ## Likelihood ratio test= 63.58  on 6 df,   p=8.392e-12
    ## Wald test            = 67.15  on 6 df,   p=1.572e-12
    ## Score (logrank) test = 68.59  on 6 df,   p=7.973e-13

    summary(m2)

    ## Call:
    ## coxph(formula = Surv(rep(1, 12144L), case_) ~ forest + forest:cos_ta +
    ##     forest:log_sl + log_sl + cos_ta + strata(step_id_), data = data,
    ##     method = "exact")
    ##
    ##   n= 12144, number of events= 759
    ##
    ##                             coef exp(coef) se(coef)      z Pr(>|z|)
    ## forestnon-forest         0.62002   1.85896  0.33146  1.871 0.061406 .
    ## log_sl                   0.17834   1.19524  0.04944  3.607 0.000310 ***
    ## cos_ta                   0.22910   1.25747  0.09806  2.336 0.019478 *
    ## forestnon-forest:cos_ta -0.40553   0.66662  0.11798 -3.437 0.000587 ***
    ## forestnon-forest:log_sl -0.23911   0.78733  0.05809 -4.116 3.85e-05 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ##                         exp(coef) exp(-coef) lower .95 upper .95
    ## forestnon-forest           1.8590     0.5379    0.9708    3.5597
    ## log_sl                     1.1952     0.8367    1.0848    1.3169
    ## cos_ta                     1.2575     0.7952    1.0376    1.5239
    ## forestnon-forest:cos_ta    0.6666     1.5001    0.5290    0.8400
    ## forestnon-forest:log_sl    0.7873     1.2701    0.7026    0.8823
    ##
    ## Rsquare= 0.005   (max possible= 0.293 )
    ## Likelihood ratio test= 63.42  on 5 df,   p=2.382e-12
    ## Wald test            = 66.66  on 5 df,   p=5.081e-13
    ## Score (logrank) test = 67.81  on 5 df,   p=2.918e-13

    #AIC(m0$model)
    #AIC(m1$model)
    #AIC(m2$model)

### Interpretation of coefficients

To be done.

## 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(sh.forest, levels = 1:2, 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, 12144L), case_) ~ forest + forest:cos_ta +
    ##     forest:sl_ + sl_ * cos_ta + strata(step_id_), data = data,
    ##     method = "exact")
    ##
    ##   n= 12144, number of events= 759
    ##
    ##                               coef  exp(coef)   se(coef)      z Pr(>|z|)
    ## forestnon-forest        -2.617e-01  7.698e-01  1.411e-01 -1.855 0.063636
    ## sl_                      6.729e-04  1.001e+00  1.576e-04  4.270 1.96e-05
    ## cos_ta                   2.389e-01  1.270e+00  1.107e-01  2.158 0.030958
    ## forestnon-forest:cos_ta -4.130e-01  6.616e-01  1.174e-01 -3.519 0.000433
    ## forestnon-forest:sl_    -9.060e-04  9.991e-01  2.008e-04 -4.513 6.40e-06
    ## sl_:cos_ta               2.295e-05  1.000e+00  1.329e-04  0.173 0.862878
    ##
    ## forestnon-forest        .
    ## sl_                     ***
    ## cos_ta                  *
    ## forestnon-forest:cos_ta ***
    ## forestnon-forest:sl_    ***
    ## sl_:cos_ta
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ##                         exp(coef) exp(-coef) lower .95 upper .95
    ## forestnon-forest           0.7698     1.2991    0.5838    1.0150
    ## sl_                        1.0007     0.9993    1.0004    1.0010
    ## cos_ta                     1.2698     0.7875    1.0221    1.5775
    ## forestnon-forest:cos_ta    0.6616     1.5114    0.5257    0.8328
    ## forestnon-forest:sl_       0.9991     1.0009    0.9987    0.9995
    ## sl_:cos_ta                 1.0000     1.0000    0.9998    1.0003
    ##
    ## Rsquare= 0.005   (max possible= 0.293 )
    ## Likelihood ratio test= 65.38  on 6 df,   p=3.616e-12
    ## Wald test            = 71.61  on 6 df,   p=1.914e-13
    ## Score (logrank) test = 75.79  on 6 df,   p=2.642e-14

## Session

    devtools::session_info()

    ##  setting  value
    ##  version  R version 3.4.3 (2017-11-30)
    ##  system   x86_64, linux-gnu
    ##  ui       X11
    ##  language en_US
    ##  collate  C
    ##  tz       Europe/Berlin
    ##  date     2018-04-21
    ##
    ##  package      * version    date       source
    ##  amt          * 0.0.4.0    2018-04-21 local
    ##  ansistrings    1.0.0.9000 2018-02-07 Github (r-lib/ansistrings@f27619b)
    ##  assertthat     0.2.0      2017-04-11 CRAN (R 3.4.3)
    ##  backports      1.1.2      2017-12-13 CRAN (R 3.4.3)
    ##  base         * 3.4.3      2017-12-01 local
    ##  bindr          0.1.1      2018-03-13 CRAN (R 3.4.3)
    ##  bindrcpp     * 0.2.2      2018-03-29 cran (@0.2.2)
    ##  boot           1.3-20     2017-07-30 CRAN (R 3.4.1)
    ##  circular       0.4-93     2017-06-29 CRAN (R 3.4.3)
    ##  cli            1.0.0.9001 2018-02-07 Github (r-lib/cli@1b58269)
    ##  compiler       3.4.3      2017-12-01 local
    ##  crayon         1.3.4      2017-09-16 CRAN (R 3.4.3)
    ##  datasets     * 3.4.3      2017-12-01 local
    ##  devtools       1.13.5     2018-02-18 CRAN (R 3.4.3)
    ##  digest         0.6.15     2018-01-28 CRAN (R 3.4.3)
    ##  dplyr          0.7.5      2018-04-20 local
    ##  evaluate       0.10.1     2017-06-24 CRAN (R 3.4.3)
    ##  fitdistrplus   1.0-9      2017-03-24 CRAN (R 3.4.3)
    ##  glue           1.2.0      2017-10-29 CRAN (R 3.4.3)
    ##  graphics     * 3.4.3      2017-12-01 local
    ##  grDevices    * 3.4.3      2017-12-01 local
    ##  grid           3.4.3      2017-12-01 local
    ##  hms            0.4.2      2018-03-10 CRAN (R 3.4.3)
    ##  htmltools      0.3.6      2017-04-28 CRAN (R 3.4.3)
    ##  knitr        * 1.20       2018-02-20 CRAN (R 3.4.3)
    ##  lattice        0.20-35    2017-03-25 CRAN (R 3.4.3)
    ##  lazyeval       0.2.1      2017-10-29 CRAN (R 3.4.3)
    ##  lubridate    * 1.7.4      2018-04-11 cran (@1.7.4)
    ##  magrittr       1.5        2014-11-22 CRAN (R 3.4.3)
    ##  MASS           7.3-49     2018-02-23 CRAN (R 3.4.3)
    ##  Matrix         1.2-12     2017-11-16 CRAN (R 3.4.3)
    ##  memoise        1.1.0      2017-04-21 CRAN (R 3.4.3)
    ##  methods      * 3.4.3      2017-12-01 local
    ##  mvtnorm        1.0-7      2018-01-26 CRAN (R 3.4.3)
    ##  pillar         1.2.1      2018-02-27 CRAN (R 3.4.3)
    ##  pkgconfig      2.0.1      2017-03-21 CRAN (R 3.4.3)
    ##  prettyunits    1.0.2      2015-07-13 cran (@1.0.2)
    ##  progress       1.1.2.9002 2018-02-07 Github (r-lib/progress@97f2c4e)
    ##  purrr          0.2.4      2017-10-18 CRAN (R 3.4.3)
    ##  R6             2.2.2      2017-06-17 CRAN (R 3.4.3)
    ##  raster       * 2.6-7      2017-11-13 CRAN (R 3.4.3)
    ##  Rcpp           0.12.16    2018-03-13 CRAN (R 3.4.3)
    ##  rematch2       2.0.1      2017-06-20 CRAN (R 3.4.3)
    ##  rgdal          1.2-18     2018-03-17 cran (@1.2-18)
    ##  rlang          0.2.0.9001 2018-04-13 Github (tidyverse/rlang@82b2727)
    ##  rmarkdown      1.9        2018-03-01 CRAN (R 3.4.3)
    ##  rprojroot      1.3-2      2018-01-03 CRAN (R 3.4.3)
    ##  selectr        0.3-2      2018-03-05 CRAN (R 3.4.3)
    ##  sp           * 1.2-7      2018-01-19 CRAN (R 3.4.3)
    ##  splines        3.4.3      2017-12-01 local
    ##  stats        * 3.4.3      2017-12-01 local
    ##  stringi        1.1.7      2018-03-12 CRAN (R 3.4.3)
    ##  stringr        1.3.0      2018-02-19 CRAN (R 3.4.3)
    ##  survival       2.41-3     2017-04-04 CRAN (R 3.4.0)
    ##  tibble         1.4.2      2018-01-22 CRAN (R 3.4.3)
    ##  tidyr          0.8.0      2018-01-29 CRAN (R 3.4.3)
    ##  tidyselect     0.2.4      2018-02-26 CRAN (R 3.4.3)
    ##  tools          3.4.3      2017-12-01 local
    ##  utf8           1.1.3      2018-01-03 CRAN (R 3.4.3)
    ##  utils        * 3.4.3      2017-12-01 local
    ##  withr          2.1.2      2018-03-15 CRAN (R 3.4.3)
    ##  xml2           1.2.0      2018-01-24 CRAN (R 3.4.3)
    ##  yaml           2.1.18     2018-03-08 CRAN (R 3.4.3)