Fitting Step-Selection Functions with amt

Johannes Signer

2019-09-19

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 
##   <table>  <table>  <table>  <table>  <table>  <table>  <dbl> <int> <chr>
## 1 5.964444 5.996667 6.005833 11.46179 6.016111 3923.997  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)
## crs        : +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     : 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.5200    0.5945   0.1089 -4.776 1.79e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## forestnon-forest    0.5945      1.682    0.4803    0.7359
## 
## Concordance= 0.529  (se = 0.007 )
## Likelihood ratio test= 22.06  on 1 df,   p=3e-06
## Wald test            = 22.81  on 1 df,   p=2e-06
## Score (logrank) test = 22.9  on 1 df,   p=2e-06
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.77414   2.16873  0.33061  2.342 0.019202 *  
## log_sl                   0.18971   1.20890  0.04993  3.800 0.000145 ***
## cos_ta                  -0.26631   0.76620  0.20831 -1.278 0.201102    
## forestnon-forest:cos_ta -0.25404   0.77566  0.11761 -2.160 0.030767 *  
## forestnon-forest:log_sl -0.24688   0.78124  0.05839 -4.228 2.36e-05 ***
## log_sl:cos_ta            0.02454   1.02484  0.03500  0.701 0.483242    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## forestnon-forest           2.1687     0.4611    1.1345    4.1459
## log_sl                     1.2089     0.8272    1.0962    1.3332
## cos_ta                     0.7662     1.3051    0.5094    1.1525
## forestnon-forest:cos_ta    0.7757     1.2892    0.6160    0.9767
## forestnon-forest:log_sl    0.7812     1.2800    0.6968    0.8760
## log_sl:cos_ta              1.0248     0.9758    0.9569    1.0976
## 
## Concordance= 0.601  (se = 0.013 )
## Likelihood ratio test= 83.97  on 6 df,   p=5e-16
## Wald test            = 82.65  on 6 df,   p=1e-15
## Score (logrank) test = 83.99  on 6 df,   p=5e-16
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.78755   2.19801  0.32932  2.391 0.016783 *  
## log_sl                   0.18953   1.20868  0.04984  3.803 0.000143 ***
## cos_ta                  -0.13722   0.87178  0.09740 -1.409 0.158885    
## forestnon-forest:cos_ta -0.25885   0.77194  0.11740 -2.205 0.027471 *  
## forestnon-forest:log_sl -0.24921   0.77941  0.05819 -4.283 1.85e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## forestnon-forest           2.1980     0.4550    1.1527    4.1913
## log_sl                     1.2087     0.8273    1.0962    1.3327
## cos_ta                     0.8718     1.1471    0.7203    1.0551
## forestnon-forest:cos_ta    0.7719     1.2954    0.6133    0.9717
## forestnon-forest:log_sl    0.7794     1.2830    0.6954    0.8736
## 
## Concordance= 0.598  (se = 0.013 )
## Likelihood ratio test= 83.47  on 5 df,   p=<2e-16
## Wald test            = 81.76  on 5 df,   p=4e-16
## Score (logrank) test = 83.22  on 5 df,   p=<2e-16
#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        -0.1209311  0.8860950  0.1400205 -0.864  0.38777
## sl_                      0.0006510  1.0006512  0.0001568  4.151 3.31e-05
## cos_ta                  -0.2893298  0.7487652  0.1095607 -2.641  0.00827
## forestnon-forest:cos_ta -0.2178123  0.8042764  0.1181950 -1.843  0.06536
## forestnon-forest:sl_    -0.0009040  0.9990964  0.0002001 -4.517 6.27e-06
## sl_:cos_ta               0.0002409  1.0002410  0.0001317  1.829  0.06737
##                            
## 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.8861     1.1285    0.6734    1.1659
## sl_                        1.0007     0.9993    1.0003    1.0010
## cos_ta                     0.7488     1.3355    0.6041    0.9281
## forestnon-forest:cos_ta    0.8043     1.2434    0.6380    1.0139
## forestnon-forest:sl_       0.9991     1.0009    0.9987    0.9995
## sl_:cos_ta                 1.0002     0.9998    1.0000    1.0005
## 
## Concordance= 0.603  (se = 0.013 )
## Likelihood ratio test= 98.26  on 6 df,   p=<2e-16
## Wald test            = 97.25  on 6 df,   p=<2e-16
## Score (logrank) test = 111.3  on 6 df,   p=<2e-16

Session

devtools::session_info()
## ─ Session info ──────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.1 (2019-07-05)
##  os       Ubuntu 18.04.3 LTS          
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language en_US:en                    
##  collate  C                           
##  ctype    en_US.UTF-8                 
##  tz       Europe/Berlin               
##  date     2019-09-19                  
## 
## ─ Packages ──────────────────────────────────────────────────────────────
##  package      * version  date       lib source        
##  amt          * 0.0.7    2019-09-19 [1] local         
##  assertthat     0.2.1    2019-03-21 [3] CRAN (R 3.6.1)
##  backports      1.1.4    2019-04-10 [3] CRAN (R 3.6.1)
##  boot           1.3-23   2019-07-05 [6] CRAN (R 3.6.1)
##  callr          3.3.1    2019-07-18 [3] CRAN (R 3.6.1)
##  circular       0.4-93   2017-06-29 [3] CRAN (R 3.6.1)
##  cli            1.1.0    2019-03-19 [3] CRAN (R 3.6.1)
##  codetools      0.2-16   2018-12-24 [6] CRAN (R 3.5.2)
##  colorspace     1.4-1    2019-03-18 [3] CRAN (R 3.6.1)
##  crayon         1.3.4    2017-09-16 [3] CRAN (R 3.6.1)
##  desc           1.2.0    2018-05-01 [3] CRAN (R 3.6.1)
##  devtools       2.2.0    2019-09-07 [3] CRAN (R 3.6.1)
##  digest         0.6.20   2019-07-04 [3] CRAN (R 3.6.1)
##  dplyr        * 0.8.3    2019-07-04 [3] CRAN (R 3.6.1)
##  DT             0.9      2019-09-17 [3] CRAN (R 3.6.1)
##  ellipsis       0.2.0.1  2019-07-02 [3] CRAN (R 3.6.1)
##  evaluate       0.14     2019-05-28 [3] CRAN (R 3.6.1)
##  fansi          0.4.0    2018-10-05 [3] CRAN (R 3.6.1)
##  fitdistrplus   1.0-14   2019-01-23 [3] CRAN (R 3.6.1)
##  fs             1.3.1    2019-05-06 [3] CRAN (R 3.6.1)
##  ggplot2      * 3.2.1    2019-08-10 [3] CRAN (R 3.6.1)
##  glue           1.3.1    2019-03-12 [3] CRAN (R 3.6.1)
##  gtable         0.3.0    2019-03-25 [3] CRAN (R 3.6.1)
##  htmltools      0.3.6    2017-04-28 [3] CRAN (R 3.6.1)
##  htmlwidgets    1.3      2018-09-30 [3] CRAN (R 3.6.1)
##  knitr        * 1.25     2019-09-18 [3] CRAN (R 3.6.1)
##  labeling       0.3      2014-08-23 [3] CRAN (R 3.6.1)
##  lattice        0.20-38  2018-11-04 [6] CRAN (R 3.5.1)
##  lazyeval       0.2.2    2019-03-15 [3] CRAN (R 3.6.1)
##  lifecycle      0.1.0    2019-08-01 [3] CRAN (R 3.6.1)
##  lsei           1.2-0    2017-10-23 [3] CRAN (R 3.6.1)
##  lubridate    * 1.7.4    2018-04-11 [3] CRAN (R 3.6.1)
##  magrittr       1.5      2014-11-22 [3] CRAN (R 3.6.1)
##  MASS           7.3-51.4 2019-04-26 [6] CRAN (R 3.6.1)
##  Matrix         1.2-17   2019-03-22 [6] CRAN (R 3.6.1)
##  memoise        1.1.0    2017-04-21 [3] CRAN (R 3.6.1)
##  munsell        0.5.0    2018-06-12 [3] CRAN (R 3.6.1)
##  mvtnorm        1.0-11   2019-06-19 [3] CRAN (R 3.6.1)
##  npsurv         0.4-0    2017-10-14 [3] CRAN (R 3.6.1)
##  pillar         1.4.2    2019-06-29 [3] CRAN (R 3.6.1)
##  pkgbuild       1.0.5    2019-08-26 [3] CRAN (R 3.6.1)
##  pkgconfig      2.0.2    2018-08-16 [3] CRAN (R 3.6.1)
##  pkgload        1.0.2    2018-10-29 [3] CRAN (R 3.6.1)
##  prettyunits    1.0.2    2015-07-13 [3] CRAN (R 3.6.1)
##  processx       3.4.1    2019-07-18 [3] CRAN (R 3.6.1)
##  ps             1.3.0    2018-12-21 [3] CRAN (R 3.6.1)
##  purrr          0.3.2    2019-03-15 [3] CRAN (R 3.6.1)
##  R6             2.4.0    2019-02-14 [3] CRAN (R 3.6.1)
##  raster       * 3.0-2    2019-08-22 [3] CRAN (R 3.6.1)
##  Rcpp           1.0.2    2019-07-25 [3] CRAN (R 3.6.1)
##  remotes        2.1.0    2019-06-24 [3] CRAN (R 3.6.1)
##  rgdal          1.4-4    2019-05-29 [3] CRAN (R 3.6.1)
##  rlang          0.4.0    2019-06-25 [3] CRAN (R 3.6.1)
##  rmarkdown      1.15     2019-08-21 [3] CRAN (R 3.6.1)
##  rprojroot      1.3-2    2018-01-03 [3] CRAN (R 3.6.1)
##  scales         1.0.0    2018-08-09 [3] CRAN (R 3.6.1)
##  sessioninfo    1.1.1    2018-11-05 [3] CRAN (R 3.6.1)
##  sp           * 1.3-1    2018-06-05 [3] CRAN (R 3.6.1)
##  stringi        1.4.3    2019-03-12 [3] CRAN (R 3.6.1)
##  stringr        1.4.0    2019-02-10 [3] CRAN (R 3.6.1)
##  survival       2.44-1.1 2019-04-01 [6] CRAN (R 3.6.1)
##  testthat       2.2.1    2019-07-25 [3] CRAN (R 3.6.1)
##  tibble         2.1.3    2019-06-06 [3] CRAN (R 3.6.1)
##  tidyr          1.0.0    2019-09-11 [3] CRAN (R 3.6.1)
##  tidyselect     0.2.5    2018-10-11 [3] CRAN (R 3.6.1)
##  usethis        1.5.1    2019-07-04 [3] CRAN (R 3.6.1)
##  utf8           1.1.4    2018-05-24 [3] CRAN (R 3.6.1)
##  vctrs          0.2.0    2019-07-05 [3] CRAN (R 3.6.1)
##  withr          2.1.2    2018-03-15 [3] CRAN (R 3.6.1)
##  xfun           0.9      2019-08-21 [3] CRAN (R 3.6.1)
##  yaml           2.2.0    2018-07-25 [3] CRAN (R 3.6.1)
##  zeallot        0.1.0    2018-01-28 [3] CRAN (R 3.6.1)
## 
## [1] /tmp/Rtmp6e2sjN/Rinst545c195c972e
## [2] /tmp/Rtmp3D0rnU/temp_libpath95d4161a4fb
## [3] /home/jsigner/R/x86_64-pc-linux-gnu-library/3.6
## [4] /usr/local/lib/R/site-library
## [5] /usr/lib/R/site-library
## [6] /usr/lib/R/library