amt
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.
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)
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.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)
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 -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
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)