Fitting Step-Selection Functions with amt

Johannes Signer

2018-02-08

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.00
##  2 4314053 3445768 2008-03-30 06:00:54   1.00
##  3 4314105 3445859 2008-03-30 12:01:47   1.00
##  4 4314044 3445785 2008-03-30 18:01:24   1.00
##  5 4313015 3445858 2008-03-31 00:01:23   1.00
##  6 4312860 3445857 2008-03-31 06:01:45   1.00
##  7 4312854 3445856 2008-03-31 12:01:11   1.00
##  8 4312858 3445858 2008-03-31 18:01:55   1.00
##  9 4312745 3445862 2008-04-01 00:01:24   1.00
## 10 4312651 3446024 2008-04-01 06:00:54   1.00
## # ... 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)
## [1] 4177.111
AIC(m1$model)
## [1] 4157.206
AIC(m2$model)
## [1] 4155.369

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-02-08                  
## 
##  package      * version    date       source                            
##  amt          * 0.0.2.0    2018-02-08 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        2016-11-13 CRAN (R 3.4.3)                    
##  bindrcpp     * 0.2        2017-06-17 CRAN (R 3.4.3)                    
##  boot           1.3-20     2017-07-30 CRAN (R 3.4.1)                    
##  broom          0.4.3      2017-11-20 CRAN (R 3.4.3)                    
##  cellranger     1.1.0      2016-07-27 CRAN (R 3.4.3)                    
##  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)        
##  colorspace     1.3-2      2016-12-14 CRAN (R 3.4.3)                    
##  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.4     2017-11-09 CRAN (R 3.4.3)                    
##  digest         0.6.15     2018-01-28 CRAN (R 3.4.3)                    
##  dplyr        * 0.7.4      2017-09-28 CRAN (R 3.4.3)                    
##  evaluate       0.10.1     2017-06-24 CRAN (R 3.4.3)                    
##  fitdistrplus   1.0-9      2017-03-24 CRAN (R 3.4.3)                    
##  forcats      * 0.2.0      2017-01-23 CRAN (R 3.4.3)                    
##  foreign        0.8-69     2017-06-21 CRAN (R 3.4.0)                    
##  ggplot2      * 2.2.1      2016-12-30 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                             
##  gtable         0.2.0      2016-02-26 CRAN (R 3.4.3)                    
##  haven          1.1.1      2018-01-18 CRAN (R 3.4.3)                    
##  hms            0.4.1      2018-01-24 CRAN (R 3.4.3)                    
##  htmltools      0.3.6      2017-04-28 CRAN (R 3.4.3)                    
##  httr           1.3.1      2017-08-20 CRAN (R 3.4.3)                    
##  jsonlite       1.5        2017-06-01 CRAN (R 3.4.3)                    
##  knitr        * 1.19       2018-01-29 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.2      2018-02-06 CRAN (R 3.4.3)                    
##  magrittr       1.5        2014-11-22 CRAN (R 3.4.3)                    
##  MASS           7.3-48     2017-12-24 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                             
##  mnormt         1.5-5      2016-10-15 CRAN (R 3.4.3)                    
##  modelr         0.1.1      2017-07-24 CRAN (R 3.4.3)                    
##  munsell        0.4.3      2016-02-13 CRAN (R 3.4.3)                    
##  mvtnorm        1.0-7      2018-01-26 CRAN (R 3.4.3)                    
##  nlme           3.1-131    2017-02-06 CRAN (R 3.4.0)                    
##  parallel       3.4.3      2017-12-01 local                             
##  pillar         1.1.0      2018-01-14 CRAN (R 3.4.3)                    
##  pkgconfig      2.0.1      2017-03-21 CRAN (R 3.4.3)                    
##  plyr           1.8.4      2016-06-08 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)   
##  psych          1.7.8      2017-09-09 CRAN (R 3.4.3)                    
##  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.15    2018-01-20 CRAN (R 3.4.3)                    
##  readr        * 1.1.1      2017-05-16 CRAN (R 3.4.3)                    
##  readxl         1.0.0      2017-04-18 CRAN (R 3.4.3)                    
##  rematch2       2.0.1      2017-06-20 CRAN (R 3.4.3)                    
##  reshape2       1.4.3      2017-12-11 CRAN (R 3.4.3)                    
##  rgdal          1.2-16     2017-11-21 CRAN (R 3.4.3)                    
##  rlang          0.1.6.9003 2018-02-07 Github (tidyverse/rlang@4b487f0)  
##  rmarkdown      1.8        2017-11-17 CRAN (R 3.4.3)                    
##  rprojroot      1.3-2      2018-01-03 CRAN (R 3.4.3)                    
##  rstudioapi     0.7        2017-09-07 CRAN (R 3.4.3)                    
##  rvest          0.3.2      2016-06-17 CRAN (R 3.4.3)                    
##  scales         0.5.0      2017-08-24 CRAN (R 3.4.3)                    
##  selectr        0.3-1      2016-12-19 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.6      2017-11-17 CRAN (R 3.4.3)                    
##  stringr      * 1.2.0      2017-02-18 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)                    
##  tidyverse    * 1.2.1      2017-11-14 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.1      2017-12-19 CRAN (R 3.4.3)                    
##  XML            3.98-1.9   2017-06-19 CRAN (R 3.4.3)                    
##  xml2           1.2.0      2018-01-24 CRAN (R 3.4.3)                    
##  yaml           2.1.16     2017-12-12 CRAN (R 3.4.3)