#'---
#'title: RSF analysis of mountain goats (Section 4.1)
#'author: "S. Muff, J. Signer, J. Fieberg"
#'date: "`r format(Sys.time(), '%d %B, %Y')`"
#'output:
#'  html_document:
#'    toc: yes
#'---

#' ## Load libraries and data
#+ echo=TRUE, message=FALSE, warning=FALSE
library(ResourceSelection)
library(glmmTMB)
library(INLA)
data(goats)
str(goats)

#' Scale and center variables
#+ echo=TRUE, message=FALSE
goats$ELEVATION <- scale(goats$ELEVATION)
goats$TASP <- scale(goats$TASP)

#' Use and available data by animal
with(goats, prop.table(table(ID, STATUS), 1))


#' ## glmmTMB()
#'
#' ### Model M1:  
#' 
#' - Fixed Effects: elevation
#' - Random Effects: intercept only
#' 
#+ echo=TRUE, message=FALSE,cache=TRUE
goats.M1 <- glmmTMB(STATUS ~ ELEVATION  + (1|ID), family=binomial(), data = goats)
summary(goats.M1)

#' ### Model M2 
#' 
#' - Fixed Effects: elevation, aspect
#' - Random Effects: intercept only
#' 
#+ echo=TRUE, message=FALSE,cache=TRUE
goats.M2 <- glmmTMB(STATUS ~ ELEVATION + TASP + (1|ID), family=binomial(), data = goats)
summary(goats.M2)


#' ### Model M3 
#'
#' - Fixed Effects: elevation, aspect
#' - Random Effects: intercepts and slopes (elevation, aspect)
#' 
#+ echo=TRUE, message=FALSE,cache=TRUE
goats.M3 <- glmmTMB(STATUS ~ TASP +  ELEVATION + (1|ID) + (0+ELEVATION |ID) + (0+TASP|ID), family=binomial(), data = goats)
summary(goats.M3)



#' ### Model M4 (with fixed intercept variance)
#'
#' - Fixed Effects: elevation
#' - Random Effects: intercept (with large fixed variance), elevation, aspect
#' 
#' Here, we also use a weighted likelihood. To this end, we need to create a variable for the weights, where used points (`STATUS=1`) keep weight 1, and available points (`STATUS=0`) obtain a large weight $W$ (here $W=1000$):
#+ echo=TRUE, message=FALSE
goats$weight <- 1000^(1-goats$STATUS)


#' We fit the same model as under M3, again using `glmmTMB`. Note that we have to manually fix the variance of the intercept first. Start by setting up the model, but do not yet fit it:
#+ echo=TRUE, message=FALSE,cache=TRUE
goats.M4.tmp <- glmmTMB(STATUS ~   ELEVATION + TASP + (1|ID) + (0+ELEVATION |ID) + (0+TASP|ID), family=binomial(), data = goats,doFit=F, weights = weight)


#' Then fix the standard deviation of the first random term, which is the `(1|ID)` component  in the above model equation. We use $\sigma=10^3$, which corresponds to a variance of $10^6$:
#+ echo=TRUE, message=FALSE,cache=TRUE
goats.M4.tmp$parameters$theta[1] = log(1e3)


#' We need to tell `glmmTMB` not to change the first entry of the vector of variances, and give all other variances another indicator to make sure they can be freely estimated:
#+ echo=TRUE, message=FALSE,cache=TRUE
goats.M4.tmp$mapArg = list(theta=factor(c(NA,1:2)))

#' Then fit the model and look at the results:
#+ echo=TRUE, message=FALSE, cache=TRUE 
goats.M4 <- glmmTMB:::fitTMB(goats.M4.tmp)
summary(goats.M4)


#' ### Model M4 (with intercept variance estimated)
#'
#' For comparison, we again fit model M4, but without fixing the intercept variance, letting it be estimated instead. Importantly, estimating the intercept variance is the current standard procedure. For this particular RSF case, it does not lead to a real difference, as expected due to the many observations per individual. This confirms that the decision to fix or estimate the intercept variance is not critical for RSFs, in contrast to SSFs (see Discussion in the paper).
#+ echo=TRUE, message=FALSE, cache=TRUE
goats.M4.2 <- glmmTMB(STATUS ~   ELEVATION + TASP + (1|ID) + (0+ELEVATION |ID) + (0+TASP|ID), family=binomial(), data = goats, weights = weight)
summary(goats.M4.2)


#' ##  INLA (only model M4)
#'
#' Let us now carry the analysis of model M4 with random intercept $\mathsf{N}(0,\sigma_{ID}^2)$ and fixed variance $\sigma_{ID}^2=10^6$ using INLA. A peculiarity of INLA is that the same variable cannot be used more than once. So for ID we need to generate two new (but identical) variables:
#+ echo=TRUE, message=FALSE
goats$ID2 <- goats$ID3 <- goats$ID

#' For the fixed effects we use the INLA (default) priors $\beta \sim \mathsf{N}(0,\sigma_\beta^2)$ with $\sigma_\beta^2=10^4$. The precisions of the priors are thus set to:
#+ echo=TRUE, message=FALSE
prec.beta.TASP  <- 1e-4  
prec.beta.ELEVATION  <- 1e-4  

#' The INLA formula with the fixed effects `TASP` and `ELEVATION`, plus three random effects: one for the individual-specific intercept, and two random slopes for `TASP` and `ELEVATION`. Note that the precision (thus $1/\sigma^2$) for `ID` is fixed (`fixed=T`) at the value of $10^{-6}$ (thus the variance is fixed at $10^6$). The precisions for the random slopes for `TASP` and `ELEVATION` are given PC(1,0.05) priors:
#+ echo=TRUE, message=FALSE
formula.inla <-STATUS ~  TASP  + ELEVATION + 
  f(ID,model="iid",hyper=list(theta = list(initial=log(1e-6),fixed=T))) +
  f(ID2,TASP,values=1:10,model="iid",
    hyper=list(theta=list(initial=log(1),fixed=F,prior="pc.prec",param=c(1,0.05)))) + 
  f(ID3,ELEVATION,values=1:10,model="iid",
    hyper=list(theta=list(initial=log(1),fixed=F,prior="pc.prec",param=c(1,0.05)))) 


#' The actual INLA call is then given as follows:
#+  echo=TRUE, message=FALSE, cache=TRUE
inla.setOption(enable.inla.argument.weights=TRUE)
goats.M4.inla  <- inla(formula.inla, family ="binomial", data=goats, weights=goats$weight,
                control.fixed = list(
                  mean = 0,
                  prec = list(TASP = prec.beta.TASP,
                              ELEVATION = prec.beta.ELEVATION)
                )
)

 

#' The summary for the posterior distribution of the fixed effects is given as follows:
#+ echo=TRUE
goats.M4.inla$summary.fixed


#' Since variances are parameterized and treated as precisions, the summary of the respective posterior distributions is given for the precisions:
#+ echo=TRUE
goats.M4.inla$summary.hyperpar


#' Source R functions for calculating posterior means 
#' and medians of the precisions.
source("inla_emarginal.R")
source("inla_mmarginal.R")
inla_emarginal(goats.M4.inla)
inla_mmarginal(goats.M4.inla)

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