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)