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: 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 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.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-06summary(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-16summary(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)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 -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-16devtools::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