#'---
#'title: SSF analysis of fisher using Stan 
#'author: "S. Muff, J. Signer, J. Fieberg"
#'date: "`r format(Sys.time(), '%d %B, %Y')`"
#'output:
#'  html_document:
#'    toc: yes
#'---

#+ include = FALSE
knitr::opts_chunk$set(warning = FALSE, message = FALSE, cache = TRUE)

#' ## 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
#+ echo=TRUE, message=FALSE, warning=FALSE
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()


stan_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()

summary(res_stan)$summary[1:2, ]
tail(summary(res_stan)$summary)

traceplot(res_stan,pars=c("beta[1]","a_re[1]","a_re[2]"))

#' ## Session Info
#'
devtools::session_info()