---
title: "Lab 09 - Analyses in RMark"
author: "WILD 502 - Jay Rotella"
output: html_document
---

This markdown file runs 3 robust-design models of interest in RMark. The models are described in section 15.6.1 of Cooch & White' "Gentle Introduction to MARK". They represent models with either (i) no temporary emigration out of the study population, (ii) random temporary emigration, or (iii) Markovian temporary
emigration. 

### Bring in the data

```{r, warning=FALSE}
library(RMark)
rd.inp <- convert.inp("http://www.montana.edu/rotella/documents/502/rd_simple1.inp")

# set up time intervals for 5 primary occasions each with 3 secondary occasions
# these values are the interval lengths between occasions
# e.g, 0 time between occ 1 & 2, 0 between 2 & 3, 1 year between 3 & 4
time.intervals <- c(0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0)

rd <- process.data(data = rd.inp,
                  model = "Robust",
                  time.intervals = time.intervals)
                  
names(rd)

# take a look at the data
head(rd$data)

# check the number of primary occasions 
rd$nocc

# check the number of secondary occasions 
rd$nocc.secondary
```

### Set up models to run

For this week's lab, we'll set up specific models in the RMark function. This is done because here we don't just want all possible combinations of the parameterizations for S, p, and gamma. Rather, we want only specific combinations and so need to dictate which ones we want one model at a time.

```{r}
run.robust <- function() {

  # In all 3 models, S varies by year, p varies by primary occasion but not secondary,
  # So we set up structures for S and p before building competing models of TE
  
  # Apparent survival varies by year
  S.time = list(formula =  ~ time)
  
  # p varies by primary session but not among secondaries within primary
  # session = varies by primary, time=varies by secondary
  # session:time = varies by both
  # p=c due to use of "share=TRUE"
  p.session = list(formula =  ~ session, share = TRUE)
  
  # set up Gamma structure for Markovian  time-varying temporary emigration
  # With time-varying S,p, and TE, we need to constrain some gamma's
  # A conventional way to do it is to constrain
  #    gamma"_(k) = gamm"_(k-1) & gamma'_(k) = gamm'_(k-1)
  # so that all S and remaining gamma's are identifiable.
  GammaDoublePrime.markov = list(formula =  ~ time)
  GammaPrime.time = list(formula =  ~ time)
  
  # Model 1: TE is Markovian & time-varying, S(t), p(varies by primary)
  # use design parameters to constrain some of the gammas
  # specifically, this uses time-binning to constrain gamma" and gamma'
  # bins need to be same for each type of gamma to make this happen
  mod.markov = mark(data = rd, 
                    design.parameters = list(
                      GammaDoublePrime = list(time.bins = c(1, 2, 3, 5)),
                      GammaPrime = list(time.bins = c(2, 3, 5))),
                    right = FALSE,
                    model.parameters = list(
                      S = S.time,
                      p = p.session,
                      GammaDoublePrime = GammaDoublePrime.markov,
                      GammaPrime = GammaPrime.time))
  
  # Model 2: random, time-varying TE, S(t), p(varies by primary)
  # set up Gamma structure of random temporary emigration
  # by using "share=TRUE"
  # TE rate varies by year
  # Even though Gamma' is shared with Gamma" such that don't
  # provide model.parameters explicitly for GammaPrime,
  # still need to provide time bins for GammaPrime and need to
  # make them the same as those used for GammaDoublePrime
  GammaDoublePrime.random = list(formula =  ~ time, share = TRUE)
  
  mod.random = mark(data = rd,
                    design.parameters = list(
                      GammaDoublePrime = list(time.bins = c(1, 2, 3, 5)),
                      GammaPrime = list(time.bins = c(2, 3, 5))),
                    right = FALSE,
                    model.parameters = list(
                      S = S.time,
                      p = p.session,
                      GammaDoublePrime = GammaDoublePrime.random))
            
  # Model 3 = no TE
  # fix Gamma's = 0
  GammaDoublePrime.zeroTE = list(formula =  ~ 1,
                                 fixed = 0)
  
  GammaPrime.zero = list(formula = ~ 1,
                         fixed = 0)
  
  mod.ZeroTE = mark(data = rd, 
                    model.parameters = list(
                      S = S.time,
                      p = p.session,
                      GammaDoublePrime = GammaDoublePrime.zeroTE,
                      GammaPrime = GammaPrime.zero))


# Return model table and list of models
return(collect.models())
}
```

### Run the models and examine the output

It's very simple to then use the function to run each of the models. As implemented here, the output from each of the models appears in the console where it can be reviewed. One can also examine model-specific output in other ways as shown below. 

```{r}
robust.results=run.robust()
```

### Examine model-selection table

```{r}
options(width = 90)
robust.results
names(robust.results)
```

### Examine output from each model

```{r}
round(robust.results$mod.markov$results$real[, 1:4], 3)
round(robust.results$mod.markov$results$derived$`N Population Size`, 3)

round(robust.results$mod.random$results$real[, 1:4], 3)
round(robust.results$mod.random$results$derived$`N Population Size`, 3)

round(robust.results$mod.ZeroTE$results$real[, 1:4], 3)
round(robust.results$mod.ZeroTE$results$derived$`N Population Size`, 3)
```
