# Fitting Step-Selection Functions with `amt`

#### Johannes Signer

#### 2018-09-12

## 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.5350    0.5857   0.1087 -4.921 8.61e-07 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ##                  exp(coef) exp(-coef) lower .95 upper .95
    ## forestnon-forest    0.5857      1.708    0.4733    0.7247
    ##
    ## Rsquare= 0.002   (max possible= 0.293 )
    ## Likelihood ratio test= 23.42  on 1 df,   p=1e-06
    ## Wald test            = 24.22  on 1 df,   p=9e-07
    ## Score (logrank) test = 24.33  on 1 df,   p=8e-07

    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.86021   2.36365  0.32917  2.613 0.008969 **
    ## log_sl                   0.20068   1.22223  0.04968  4.039 5.36e-05 ***
    ## cos_ta                   0.14774   1.15922  0.20689  0.714 0.475152
    ## forestnon-forest:cos_ta -0.40580   0.66644  0.11674 -3.476 0.000509 ***
    ## forestnon-forest:log_sl -0.26805   0.76487  0.05830 -4.598 4.27e-06 ***
    ## log_sl:cos_ta            0.01538   1.01550  0.03490  0.441 0.659497
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ##                         exp(coef) exp(-coef) lower .95 upper .95
    ## forestnon-forest           2.3637     0.4231    1.2399    4.5058
    ## log_sl                     1.2222     0.8182    1.1088    1.3472
    ## cos_ta                     1.1592     0.8627    0.7728    1.7389
    ## forestnon-forest:cos_ta    0.6664     1.5005    0.5301    0.8378
    ## forestnon-forest:log_sl    0.7649     1.3074    0.6823    0.8575
    ## log_sl:cos_ta              1.0155     0.9847    0.9484    1.0874
    ##
    ## Rsquare= 0.005   (max possible= 0.293 )
    ## Likelihood ratio test= 59.18  on 6 df,   p=7e-11
    ## Wald test            = 62.01  on 6 df,   p=2e-11
    ## Score (logrank) test = 63.07  on 6 df,   p=1e-11

    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.86971   2.38622  0.32836  2.649 0.008082 **
    ## log_sl                   0.20166   1.22343  0.04965  4.061 4.88e-05 ***
    ## cos_ta                   0.22843   1.25663  0.09640  2.370 0.017805 *
    ## forestnon-forest:cos_ta -0.40863   0.66456  0.11658 -3.505 0.000456 ***
    ## forestnon-forest:log_sl -0.26977   0.76355  0.05814 -4.640 3.49e-06 ***
    ## ---
    ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    ##
    ##                         exp(coef) exp(-coef) lower .95 upper .95
    ## forestnon-forest           2.3862     0.4191    1.2537    4.5417
    ## log_sl                     1.2234     0.8174    1.1100    1.3485
    ## cos_ta                     1.2566     0.7958    1.0403    1.5180
    ## forestnon-forest:cos_ta    0.6646     1.5048    0.5288    0.8351
    ## forestnon-forest:log_sl    0.7636     1.3097    0.6813    0.8557
    ##
    ## Rsquare= 0.005   (max possible= 0.293 )
    ## Likelihood ratio test= 58.99  on 5 df,   p=2e-11
    ## Wald test            = 61.41  on 5 df,   p=6e-12
    ## Score (logrank) test = 62.09  on 5 df,   p=4e-12

    #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        -1.329e-01  8.756e-01  1.409e-01 -0.943  0.34545
    ## sl_                      7.177e-04  1.001e+00  1.549e-04  4.633 3.61e-06
    ## cos_ta                   2.169e-01  1.242e+00  1.102e-01  1.968  0.04910
    ## forestnon-forest:cos_ta -3.813e-01  6.830e-01  1.172e-01 -3.254  0.00114
    ## forestnon-forest:sl_    -9.820e-04  9.990e-01  1.998e-04 -4.915 8.89e-07
    ## sl_:cos_ta               6.048e-06  1.000e+00  1.324e-04  0.046  0.96357
    ##
    ## 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.8756     1.1421    0.6643    1.1539
    ## sl_                        1.0007     0.9993    1.0004    1.0010
    ## cos_ta                     1.2422     0.8050    1.0009    1.5418
    ## forestnon-forest:cos_ta    0.6830     1.4642    0.5428    0.8593
    ## forestnon-forest:sl_       0.9990     1.0010    0.9986    0.9994
    ## sl_:cos_ta                 1.0000     1.0000    0.9997    1.0003
    ##
    ## Rsquare= 0.005   (max possible= 0.293 )
    ## Likelihood ratio test= 58.79  on 6 df,   p=8e-11
    ## Wald test            = 65.01  on 6 df,   p=4e-12
    ## Score (logrank) test = 68.94  on 6 df,   p=7e-13

## Session

    devtools::session_info()

    ##  setting  value
    ##  version  R version 3.4.4 (2018-03-15)
    ##  system   x86_64, linux-gnu
    ##  ui       X11
    ##  language en_US
    ##  collate  C
    ##  tz       Europe/Berlin
    ##  date     2018-09-12
    ##
    ##  package      * version    date       source
    ##  amt          * 0.0.5.0    2018-09-12 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.4      2018-03-16 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.4      2018-03-16 local
    ##  crayon         1.3.4      2017-09-16 CRAN (R 3.4.3)
    ##  datasets     * 3.4.4      2018-03-16 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.6      2018-06-29 cran (@0.7.6)
    ##  evaluate       0.10.1     2017-06-24 CRAN (R 3.4.3)
    ##  fitdistrplus   1.0-11     2018-09-10 cran (@1.0-11)
    ##  glue           1.2.0      2017-10-29 CRAN (R 3.4.3)
    ##  graphics     * 3.4.4      2018-03-16 local
    ##  grDevices    * 3.4.4      2018-03-16 local
    ##  grid           3.4.4      2018-03-16 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)
    ##  lsei           1.2-0      2017-10-23 cran (@1.2-0)
    ##  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-50     2018-04-30 CRAN (R 3.4.3)
    ##  Matrix         1.2-14     2018-04-09 CRAN (R 3.4.3)
    ##  memoise        1.1.0      2017-04-21 CRAN (R 3.4.3)
    ##  methods      * 3.4.4      2018-03-16 local
    ##  mvtnorm        1.0-8      2018-05-31 CRAN (R 3.4.3)
    ##  npsurv         0.4-0      2017-10-14 cran (@0.4-0)
    ##  pillar         1.2.3      2018-05-25 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.2.0      2018-06-14 CRAN (R 3.4.3)
    ##  purrr          0.2.5      2018-05-29 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.18    2018-07-23 cran (@0.12.18)
    ##  rematch2       2.0.1      2017-06-20 CRAN (R 3.4.3)
    ##  rgdal          1.3-3      2018-06-22 CRAN (R 3.4.3)
    ##  rlang          0.2.2      2018-08-16 cran (@0.2.2)
    ##  rmarkdown      1.10       2018-06-11 CRAN (R 3.4.3)
    ##  rprojroot      1.3-2      2018-01-03 CRAN (R 3.4.3)
    ##  selectr        0.4-1      2018-04-06 CRAN (R 3.4.3)
    ##  sp           * 1.3-1      2018-06-05 cran (@1.3-1)
    ##  splines        3.4.4      2018-03-16 local
    ##  stats        * 3.4.4      2018-03-16 local
    ##  stringi        1.2.3      2018-06-12 CRAN (R 3.4.3)
    ##  stringr        1.3.1      2018-05-10 CRAN (R 3.4.3)
    ##  survival       2.42-3     2018-04-16 CRAN (R 3.4.3)
    ##  tibble         1.4.2      2018-01-22 CRAN (R 3.4.3)
    ##  tidyr          0.8.1      2018-05-18 cran (@0.8.1)
    ##  tidyselect     0.2.4      2018-02-26 CRAN (R 3.4.3)
    ##  tools          3.4.4      2018-03-16 local
    ##  utf8           1.1.4      2018-05-24 CRAN (R 3.4.3)
    ##  utils        * 3.4.4      2018-03-16 local
    ##  withr          2.1.2      2018-04-27 Github (jimhester/withr@79d7b0d)
    ##  xml2           1.2.0      2018-01-24 CRAN (R 3.4.3)
    ##  yaml           2.1.19     2018-05-01 CRAN (R 3.4.3)