SSF analysis of fisher using Stan
S. Muff, J. Signer, J. Fieberg
06 Juli, 2019
Purpose
The purpose of this document is to illustrate how a simple SSF with random effects can be fitted to tracking data using Stan. We use a data set containing fisher locations from:
- LaPoint, S., Gallery, P., Wikelski, M. and Kays, R. (2013). Animal behavior, cost-based corridor models, and real corridors. Landscape Ecology, 28, 1615-1630.
Load libraries and prepare data
library(rstan)
library(tidyverse)
library(amt)
library(here)
library(tictoc)
dat_ssf <- read_rds("../Fisher_analysis/ssf_dat.rds")
model <- "
data {
int<lower=1> N; // no data points
int<lower=1> I; // no steps (over all individuals)
int<lower=1> J; // no individuals
int<lower=0> y[N]; // response
real forest[N]; // covariate
int<lower=1, upper=I> stepid[N]; // step id
int<lower=1, upper=J> indid[N]; // individual id
}
parameters {
vector[1] beta; // fixed effects
vector[I] a_re; // RE for steps
vector[J] i_re; // RE effects for each individual
real<lower = 0> sigmaind;
}
model {
real mu;
// priors
// maybe add a prior for sd here
sigmaind ~ inv_gamma(0.01, 0.01);
a_re ~ normal(0, 1000000);
i_re ~ normal(0, sqrt(sigmaind));
//likelihood
for (i in 1:N) {
mu = a_re[stepid[i]] + (beta[1] + i_re[indid[i]]) * forest[i];
y[i] ~ poisson_log(mu);
}
}
"
# compie the model
tic()
stan_mod <- stan_model(model_code = model)
toc()## 67.174 sec elapsedstan_dat <- list(N = nrow(dat_ssf), I = length(unique(dat_ssf$step_id)),
J = length(unique(dat_ssf$id)),
y = dat_ssf$y,
forest = dat_ssf$forest,
stepid = as.numeric(factor(dat_ssf$step_id)),
indid = dat_ssf$id)
tic()
res_stan <- sampling(stan_mod, stan_dat, cores = 2, chains = 2, iter = 2000)
toc()## 885.403 sec elapsedsummary(res_stan)$summary[1:2, ]## mean se_mean sd 2.5% 25% 50%
## beta[1] 0.4861758 0.01034900 0.1691656 0.1399171 0.3848552 0.4847778
## a_re[1] -2.9380057 0.04253584 1.2452178 -5.8908807 -3.5634641 -2.7429864
## 75% 97.5% n_eff Rhat
## beta[1] 0.5865577 0.842705 267.1944 1.001781
## a_re[1] -2.0478205 -1.137857 856.9996 1.001217tail(summary(res_stan)$summary)## mean se_mean sd 2.5% 25%
## i_re[3] -2.734904e-01 0.010552518 0.1844857 -6.605803e-01 -3.812069e-01
## i_re[4] 1.651602e-01 0.011451883 0.1991103 -2.565372e-01 4.492330e-02
## i_re[5] -5.346477e-02 0.011362253 0.2061377 -5.003915e-01 -1.744819e-01
## i_re[6] 4.960285e-01 0.010343257 0.1968469 1.346683e-01 3.682886e-01
## sigmaind 1.816949e-01 0.009768199 0.1976612 3.566320e-02 7.787830e-02
## lp__ -2.110145e+04 2.350206492 57.3884365 -2.121816e+04 -2.113894e+04
## 50% 75% 97.5% n_eff Rhat
## i_re[3] -2.662401e-01 -1.603033e-01 7.192235e-02 305.6421 1.0006951
## i_re[4] 1.672997e-01 2.802348e-01 5.841791e-01 302.2968 0.9990262
## i_re[5] -4.978969e-02 7.792184e-02 3.325547e-01 329.1441 1.0056600
## i_re[6] 4.862030e-01 6.159896e-01 9.012251e-01 362.1950 1.0018425
## sigmaind 1.252237e-01 2.140179e-01 6.888170e-01 409.4620 1.0024960
## lp__ -2.110166e+04 -2.106237e+04 -2.099002e+04 596.2615 0.9991913traceplot(res_stan,pars=c("beta[1]","a_re[1]","a_re[2]"))
Session Info
devtools::session_info()## ─ Session info ──────────────────────────────────────────────────────────
## setting value
## version R version 3.6.0 (2019-04-26)
## os Ubuntu 18.04.2 LTS
## system x86_64, linux-gnu
## ui X11
## language de_CH:de
## collate de_CH.UTF-8
## ctype de_CH.UTF-8
## tz Europe/Zurich
## date 2019-07-06
##
## ─ Packages ──────────────────────────────────────────────────────────────
## package * version date lib
## amt * 0.0.6 2019-03-19 [1]
## assertthat 0.2.1 2019-03-21 [2]
## backports 1.1.4 2019-04-10 [1]
## broom 0.4.2 2017-02-13 [2]
## callr 3.2.0 2019-03-15 [2]
## cellranger 1.1.0 2016-07-27 [2]
## cli 1.1.0 2019-03-19 [2]
## codetools 0.2-16 2018-12-24 [4]
## colorspace 1.4-1 2019-03-18 [1]
## crayon 1.3.4 2017-09-16 [2]
## desc 1.2.0 2018-05-01 [2]
## devtools 2.1.0 2019-07-03 [2]
## digest 0.6.19 2019-05-20 [1]
## dplyr * 0.8.3 2019-07-04 [1]
## evaluate 0.10.1 2017-06-24 [2]
## forcats * 0.2.0 2017-01-23 [2]
## foreign 0.8-71 2018-07-20 [1]
## fs 1.3.1 2019-05-06 [1]
## ggplot2 * 3.2.0 2019-06-16 [1]
## glue 1.3.1 2019-03-12 [1]
## gridExtra 2.3 2017-09-09 [2]
## gtable 0.3.0 2019-03-25 [2]
## haven 2.1.1 2019-07-04 [1]
## here * 0.1 2017-05-28 [1]
## highr 0.6 2016-05-09 [2]
## hms 0.4.2 2018-03-10 [2]
## htmltools 0.3.6 2017-04-28 [1]
## httr 1.4.0 2018-12-11 [2]
## inline 0.3.15 2018-05-18 [2]
## jsonlite 1.6 2018-12-07 [1]
## knitr 1.23 2019-05-18 [1]
## labeling 0.3 2014-08-23 [2]
## lattice 0.20-38 2018-11-04 [2]
## lazyeval 0.2.2 2019-03-15 [1]
## loo 2.1.0 2019-03-13 [2]
## lubridate 1.7.4 2018-04-11 [1]
## magrittr 1.5 2014-11-22 [2]
## Matrix 1.2-17 2019-03-22 [2]
## matrixStats 0.54.0 2018-07-23 [1]
## memoise 1.1.0 2017-04-21 [2]
## mnormt 1.5-5 2016-10-15 [1]
## modelr 0.1.1 2017-07-24 [2]
## munsell 0.5.0 2018-06-12 [2]
## nlme 3.1-140 2019-05-12 [4]
## pillar 1.4.0 2019-05-11 [2]
## pkgbuild 1.0.3 2019-03-20 [2]
## pkgconfig 2.0.2 2018-08-16 [2]
## pkgload 1.0.2 2018-10-29 [1]
## plyr 1.8.4 2016-06-08 [1]
## prettyunits 1.0.2 2015-07-13 [2]
## processx 3.4.0 2019-07-03 [1]
## ps 1.3.0 2018-12-21 [1]
## psych 1.8.4 2018-05-06 [2]
## purrr * 0.3.2 2019-03-15 [1]
## R6 2.4.0 2019-02-14 [2]
## raster 2.9-5 2019-05-14 [1]
## Rcpp 1.0.1 2019-03-17 [2]
## readr * 1.3.1 2018-12-21 [1]
## readxl 1.3.1 2019-03-13 [1]
## remotes 2.1.0 2019-06-24 [2]
## reshape2 1.4.3 2017-12-11 [1]
## rlang 0.4.0 2019-06-25 [1]
## rmarkdown 1.7 2017-11-10 [2]
## rprojroot 1.3-2 2018-01-03 [2]
## rstan * 2.18.2 2018-11-07 [1]
## rstudioapi 0.10 2019-03-19 [2]
## rvest 0.3.2 2016-06-17 [2]
## scales 1.0.0 2018-08-09 [1]
## sessioninfo 1.1.1 2018-11-05 [2]
## sp 1.3-1 2018-06-05 [1]
## StanHeaders * 2.18.1-10 2019-06-14 [1]
## stringi 1.4.3 2019-03-12 [1]
## stringr * 1.4.0 2019-02-10 [2]
## survival 2.44-1.1 2019-04-01 [2]
## testthat 2.1.1 2019-04-23 [1]
## tibble * 2.1.3 2019-06-06 [1]
## tictoc * 1.0 2019-05-20 [2]
## tidyr * 0.8.3 2019-03-01 [1]
## tidyselect 0.2.5 2018-10-11 [1]
## tidyverse * 1.2.1 2017-11-14 [2]
## usethis 1.5.0 2019-04-07 [2]
## withr 2.1.2 2018-03-15 [2]
## xfun 0.4 2018-10-23 [2]
## xml2 1.2.0 2018-01-24 [1]
## yaml 2.2.0 2018-07-25 [2]
## source
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.4.0)
## CRAN (R 3.4.4)
## CRAN (R 3.4.0)
## CRAN (R 3.4.4)
## CRAN (R 3.5.2)
## CRAN (R 3.6.0)
## CRAN (R 3.4.2)
## CRAN (R 3.4.4)
## Github (hadley/devtools@4e08802)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.1)
## CRAN (R 3.4.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.2)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.1)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.4.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.2)
## CRAN (R 3.6.0)
## CRAN (R 3.4.2)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.4.4)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.2)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.4.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## Github (collectivemedia/tictoc@211fe59)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
## CRAN (R 3.4.3)
## CRAN (R 3.4.4)
## CRAN (R 3.4.4)
## CRAN (R 3.4.4)
## CRAN (R 3.6.0)
## CRAN (R 3.6.0)
##
## [1] /home/steffi/R/x86_64-pc-linux-gnu-library/3.6
## [2] /usr/local/lib/R/site-library
## [3] /usr/lib/R/site-library
## [4] /usr/lib/R/library