amtThis 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.
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 rowsIn 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 hourThe 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)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_)) 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-09summary(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-13summary(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)
#AIC(m1$model)
#AIC(m2$model)To be done.
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-14devtools::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-04-21
##
## package * version date source
## amt * 0.0.4.0 2018-04-21 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.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.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.5 2018-02-18 CRAN (R 3.4.3)
## digest 0.6.15 2018-01-28 CRAN (R 3.4.3)
## dplyr 0.7.5 2018-04-20 local
## evaluate 0.10.1 2017-06-24 CRAN (R 3.4.3)
## fitdistrplus 1.0-9 2017-03-24 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
## 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)
## 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-49 2018-02-23 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
## mvtnorm 1.0-7 2018-01-26 CRAN (R 3.4.3)
## pillar 1.2.1 2018-02-27 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.1.2.9002 2018-02-07 Github (r-lib/progress@97f2c4e)
## 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.16 2018-03-13 CRAN (R 3.4.3)
## rematch2 2.0.1 2017-06-20 CRAN (R 3.4.3)
## rgdal 1.2-18 2018-03-17 cran (@1.2-18)
## rlang 0.2.0.9001 2018-04-13 Github (tidyverse/rlang@82b2727)
## rmarkdown 1.9 2018-03-01 CRAN (R 3.4.3)
## rprojroot 1.3-2 2018-01-03 CRAN (R 3.4.3)
## selectr 0.3-2 2018-03-05 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.7 2018-03-12 CRAN (R 3.4.3)
## stringr 1.3.0 2018-02-19 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)
## tidyselect 0.2.4 2018-02-26 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.2 2018-03-15 CRAN (R 3.4.3)
## xml2 1.2.0 2018-01-24 CRAN (R 3.4.3)
## yaml 2.1.18 2018-03-08 CRAN (R 3.4.3)