Fitting Step-Selection Functions with amt

Johannes Signer

2019-03-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 
##   <S3: tab> <S3: tab> <S3: ta> <S3: ta> <S3: ta> <S3: ta> <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)
## 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.5318    0.5875   0.1100 -4.834 1.34e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## forestnon-forest    0.5875      1.702    0.4736    0.7289
## 
## Concordance= 0.53  (se = 0.007 )
## Rsquare= 0.002   (max possible= 0.293 )
## Likelihood ratio test= 22.71  on 1 df,   p=2e-06
## Wald test            = 23.37  on 1 df,   p=1e-06
## Score (logrank) test = 23.48  on 1 df,   p=1e-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.79024   2.20392  0.32997  2.395 0.016626 *  
## log_sl                   0.19471   1.21496  0.05015  3.883 0.000103 ***
## cos_ta                  -0.37366   0.68821  0.20722 -1.803 0.071356 .  
## forestnon-forest:cos_ta -0.24225   0.78486  0.11803 -2.053 0.040119 *  
## forestnon-forest:log_sl -0.25370   0.77592  0.05833 -4.350 1.36e-05 ***
## log_sl:cos_ta            0.04330   1.04425  0.03495  1.239 0.215459    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## forestnon-forest           2.2039     0.4537    1.1543    4.2079
## log_sl                     1.2150     0.8231    1.1012    1.3404
## cos_ta                     0.6882     1.4530    0.4585    1.0330
## forestnon-forest:cos_ta    0.7849     1.2741    0.6228    0.9891
## forestnon-forest:log_sl    0.7759     1.2888    0.6921    0.8699
## log_sl:cos_ta              1.0442     0.9576    0.9751    1.1183
## 
## Concordance= 0.598  (se = 0.012 )
## Rsquare= 0.007   (max possible= 0.293 )
## Likelihood ratio test= 86.49  on 6 df,   p=<2e-16
## Wald test            = 85.65  on 6 df,   p=2e-16
## Score (logrank) test = 87.28  on 6 df,   p=<2e-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.80337   2.23305  0.32987  2.435 0.014874 *  
## log_sl                   0.19273   1.21255  0.05022  3.838 0.000124 ***
## cos_ta                  -0.14731   0.86303  0.09768 -1.508 0.131512    
## forestnon-forest:cos_ta -0.24980   0.77896  0.11792 -2.118 0.034139 *  
## forestnon-forest:log_sl -0.25583   0.77428  0.05837 -4.383 1.17e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                         exp(coef) exp(-coef) lower .95 upper .95
## forestnon-forest           2.2330     0.4478    1.1698    4.2627
## log_sl                     1.2126     0.8247    1.0989    1.3380
## cos_ta                     0.8630     1.1587    0.7127    1.0451
## forestnon-forest:cos_ta    0.7790     1.2838    0.6182    0.9815
## forestnon-forest:log_sl    0.7743     1.2915    0.6906    0.8681
## 
## Concordance= 0.597  (se = 0.012 )
## Rsquare= 0.007   (max possible= 0.293 )
## Likelihood ratio test= 84.95  on 5 df,   p=<2e-16
## Wald test            = 83.34  on 5 df,   p=<2e-16
## Score (logrank) test = 84.94  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.0891655  0.9146942  0.1419172 -0.628 0.529812
## sl_                      0.0007624  1.0007627  0.0001543  4.941 7.77e-07
## cos_ta                  -0.3799775  0.6838768  0.1102405 -3.447 0.000567
## forestnon-forest:cos_ta -0.1750264  0.8394349  0.1181076 -1.482 0.138361
## forestnon-forest:sl_    -0.0010231  0.9989774  0.0001978 -5.173 2.30e-07
## sl_:cos_ta               0.0004084  1.0004085  0.0001313  3.111 0.001865
##                            
## 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.9147     1.0933    0.6926    1.2080
## sl_                        1.0008     0.9992    1.0005    1.0011
## cos_ta                     0.6839     1.4623    0.5510    0.8488
## forestnon-forest:cos_ta    0.8394     1.1913    0.6660    1.0581
## forestnon-forest:sl_       0.9990     1.0010    0.9986    0.9994
## sl_:cos_ta                 1.0004     0.9996    1.0002    1.0007
## 
## Concordance= 0.6  (se = 0.013 )
## Rsquare= 0.009   (max possible= 0.293 )
## Likelihood ratio test= 106.7  on 6 df,   p=<2e-16
## Wald test            = 111.2  on 6 df,   p=<2e-16
## Score (logrank) test = 117.8  on 6 df,   p=<2e-16

Session

devtools::session_info()
## ─ Session info ──────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.5.3 (2019-03-11)
##  os       Ubuntu 18.04.2 LTS          
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language en_US                       
##  collate  C                           
##  ctype    en_US.UTF-8                 
##  tz       Europe/Berlin               
##  date     2019-03-19                  
## 
## ─ Packages ──────────────────────────────────────────────────────────────
##  package      * version  date       lib source        
##  amt          * 0.0.6    2019-03-19 [1] local         
##  assertthat     0.2.0    2017-04-11 [3] CRAN (R 3.5.3)
##  backports      1.1.3    2018-12-14 [3] CRAN (R 3.5.3)
##  boot           1.3-20   2017-07-30 [6] CRAN (R 3.5.0)
##  callr          3.1.1    2018-12-21 [3] CRAN (R 3.5.3)
##  circular       0.4-93   2017-06-29 [3] CRAN (R 3.5.3)
##  cli            1.0.1    2018-09-25 [3] CRAN (R 3.5.3)
##  codetools      0.2-16   2018-12-24 [6] CRAN (R 3.5.2)
##  colorspace     1.4-0    2019-01-13 [3] CRAN (R 3.5.3)
##  crayon         1.3.4    2017-09-16 [3] CRAN (R 3.5.3)
##  desc           1.2.0    2018-05-01 [3] CRAN (R 3.5.3)
##  devtools       2.0.1    2018-10-26 [3] CRAN (R 3.5.3)
##  digest         0.6.18   2018-10-10 [3] CRAN (R 3.5.3)
##  dplyr        * 0.8.0.1  2019-02-15 [3] CRAN (R 3.5.3)
##  evaluate       0.13     2019-02-12 [3] CRAN (R 3.5.3)
##  fansi          0.4.0    2018-10-05 [3] CRAN (R 3.5.3)
##  fitdistrplus   1.0-14   2019-01-23 [3] CRAN (R 3.5.3)
##  fs             1.2.6    2018-08-23 [3] CRAN (R 3.5.3)
##  ggplot2      * 3.1.0    2018-10-25 [3] CRAN (R 3.5.3)
##  glue           1.3.1    2019-03-12 [3] CRAN (R 3.5.3)
##  gtable         0.2.0    2016-02-26 [3] CRAN (R 3.5.3)
##  htmltools      0.3.6    2017-04-28 [3] CRAN (R 3.5.3)
##  knitr        * 1.22     2019-03-08 [3] CRAN (R 3.5.3)
##  labeling       0.3      2014-08-23 [3] CRAN (R 3.5.3)
##  lattice        0.20-38  2018-11-04 [6] CRAN (R 3.5.1)
##  lazyeval       0.2.1    2017-10-29 [3] CRAN (R 3.5.3)
##  lsei           1.2-0    2017-10-23 [3] CRAN (R 3.5.3)
##  lubridate    * 1.7.4    2018-04-11 [3] CRAN (R 3.5.3)
##  magrittr       1.5      2014-11-22 [3] CRAN (R 3.5.3)
##  MASS           7.3-51.1 2018-11-01 [6] CRAN (R 3.5.1)
##  Matrix         1.2-16   2019-03-08 [6] CRAN (R 3.5.2)
##  memoise        1.1.0    2017-04-21 [3] CRAN (R 3.5.3)
##  munsell        0.5.0    2018-06-12 [3] CRAN (R 3.5.3)
##  mvtnorm        1.0-10   2019-03-05 [3] CRAN (R 3.5.3)
##  npsurv         0.4-0    2017-10-14 [3] CRAN (R 3.5.3)
##  pillar         1.3.1    2018-12-15 [3] CRAN (R 3.5.3)
##  pkgbuild       1.0.2    2018-10-16 [3] CRAN (R 3.5.3)
##  pkgconfig      2.0.2    2018-08-16 [3] CRAN (R 3.5.3)
##  pkgload        1.0.2    2018-10-29 [3] CRAN (R 3.5.3)
##  plyr           1.8.4    2016-06-08 [3] CRAN (R 3.5.3)
##  prettyunits    1.0.2    2015-07-13 [3] CRAN (R 3.5.3)
##  processx       3.3.0    2019-03-10 [3] CRAN (R 3.5.3)
##  ps             1.3.0    2018-12-21 [3] CRAN (R 3.5.3)
##  purrr          0.3.1    2019-03-03 [3] CRAN (R 3.5.3)
##  R6             2.4.0    2019-02-14 [3] CRAN (R 3.5.3)
##  raster       * 2.8-19   2019-01-30 [3] CRAN (R 3.5.3)
##  Rcpp           1.0.0    2018-11-07 [3] CRAN (R 3.5.3)
##  remotes        2.0.2    2018-10-30 [3] CRAN (R 3.5.3)
##  rgdal          1.4-3    2019-03-14 [3] CRAN (R 3.5.3)
##  rlang          0.3.1    2019-01-08 [3] CRAN (R 3.5.3)
##  rmarkdown      1.11     2018-12-08 [3] CRAN (R 3.5.3)
##  rprojroot      1.3-2    2018-01-03 [3] CRAN (R 3.5.3)
##  scales         1.0.0    2018-08-09 [3] CRAN (R 3.5.3)
##  sessioninfo    1.1.1    2018-11-05 [3] CRAN (R 3.5.3)
##  sp           * 1.3-1    2018-06-05 [3] CRAN (R 3.5.3)
##  stringi        1.4.3    2019-03-12 [3] CRAN (R 3.5.3)
##  stringr        1.4.0    2019-02-10 [3] CRAN (R 3.5.3)
##  survival       2.43-3   2018-11-26 [6] CRAN (R 3.5.1)
##  testthat       2.0.1    2018-10-13 [3] CRAN (R 3.5.3)
##  tibble         2.0.1    2019-01-12 [3] CRAN (R 3.5.3)
##  tidyr          0.8.3    2019-03-01 [3] CRAN (R 3.5.3)
##  tidyselect     0.2.5    2018-10-11 [3] CRAN (R 3.5.3)
##  usethis        1.4.0    2018-08-14 [3] CRAN (R 3.5.3)
##  utf8           1.1.4    2018-05-24 [3] CRAN (R 3.5.3)
##  withr          2.1.2    2018-03-15 [3] CRAN (R 3.5.3)
##  xfun           0.5      2019-02-20 [3] CRAN (R 3.5.3)
##  yaml           2.2.0    2018-07-25 [3] CRAN (R 3.5.3)
## 
## [1] /tmp/RtmpFyX5Ry/Rinst3ee77a7b10c5
## [2] /tmp/RtmppK4TXN/temp_libpath3e6ca82d4c
## [3] /home/jsigner/R/x86_64-pc-linux-gnu-library/3.5
## [4] /usr/local/lib/R/site-library
## [5] /usr/lib/R/site-library
## [6] /usr/lib/R/library