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

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Single-Species Multi-Season Occupancy Models

Here, we'll run the single-species, multi-season occupancy models for data on northern spotted owls that you analyzed directly in MARK in lab 11. Here, we're dealing with 5 types of parameters: 'Psi', 'Epsilon', 'Gamma' & 'p'. There are different parameterizations available. In 5 of the 6 models, we'll directly model Psi in year 1 as well as Epsilon and Gamma for each year, and p for each sampling occasion. For those models, Psi in years 2-5 will be derived from Psi in year 1 and the relevant Epsilon and Gamma estimates as reviewed in class. In the 6th model, we'll model Psi in each year directly along with Gamma and p and we'll derive Epsilon for each year from the estimates of annual estimates of Psi and Gamma.

### Bring in the data & process it
```{r}
library(RMark)

# you can ignore warning about final line issue that might appear
nso <- convert.inp("http://www.montana.edu/rotella/documents/502/nso.inp")

# Process data and set up the time intervals: 0 for intervals that
# occur between secondary occasions in the same year and 1 for
# intervals that occur between years (8 surveys )
# Here, we use the 'Psi(yr1), Epsilon(t), Gamma(t)' parameterization
nso.pr.EG <- process.data(nso,model = "RDOccupEG", 
                          time.intervals = c(0, 0, 0, 0, 0, 0, 0, 1,
                                             0, 0, 0, 0, 0, 0, 0, 1,
                                             0, 0, 0, 0, 0, 0, 0, 1,
                                             0, 0, 0, 0, 0, 0, 0, 1,
                                             0, 0, 0, 0, 0, 0, 0))

nso.ddl = make.design.data(nso.pr.EG)
# add design data for models that set Epsilon to -1 and Gamma to +1
# and that hold Epsilon and Gamma constant across years
nso.ddl$Epsilon$eps = -1
nso.ddl$Gamma$eps = 1

# Make a different version of processed data for the Psi & Gamma 
# parameterization that is used in model 6, i.e., 'Psi(t),Gamma(t)'
nso.pr.PG = process.data(nso, model = "RDOccupPG",
                         time.intervals = c(0, 0, 0, 0, 0, 0, 0, 1,
                                            0, 0, 0, 0, 0, 0, 0, 1, 
                                            0, 0, 0, 0, 0, 0, 0, 1,
                                            0, 0, 0, 0, 0, 0, 0, 1,
                                            0, 0, 0, 0, 0, 0, 0))

head(nso)
nso.ddl
```

### Build function to create models

First, the function puts the observer and time covariates into proper form for use in modeling and processes the data. Second, the function declares the structures for 'p' and 'Psi' that will be used in all-possible combinations of the structures for each parameter type when the function is executed.

```{r}
# Create function to fit the 18 models in the book
fit.nso.models <- function() {
  
  # Define formulas for constructing models
  # Note that in the design data for 'p', 'session' is a 
  # factor based on primary occasions
  p.year = list(formula =  ~ session)
  Epsilon.dot = list(formula =  ~ 1)
  Epsilon.zero = list(formula =  ~ 1, 
                      fixed = list(time = 1:8, value = 0))
  Epsilon.random.shared = list(formula =  ~ -1 + eps, 
                               share = TRUE)
  Epsilon.random.yr = list(formula =  ~ time:eps,
                           share = TRUE,
                           remove.intercept = TRUE)
  Epsilon.year = list(formula =  ~ time)
  Gamma.dot = list(formula =  ~ 1)
  Gamma.zero = list(formula =  ~ 1, 
                    fixed = list(time = 1:8, value = 0))
  Gamma.year = list(formula =  ~ time)
  
  Psi.dot = list(formula =  ~ 1)
  
  # Define the models
  
  # Model 1 - no change in occupancy (epsilon & gamma=0)
  mod.1 = mark(nso.pr.EG, nso.ddl,
               model.parameters = list(
                 p = p.year, 
                 Psi = Psi.dot,
                 Epsilon = Epsilon.zero,
                 Gamma = Gamma.zero))
  
  # Model 2 - random changes in occupancy (epsilon=1-gamma)
  mod.2 = mark(nso.pr.EG, nso.ddl,
               model.parameters = list(
                 p = p.year,
                 Psi = Psi.dot,
                 Epsilon = Epsilon.random.shared))
  
  # Model 3 - random changes in occupancy (epsilon=1-gamma)
  mod.3 = mark(nso.pr.EG, nso.ddl,
               model.parameters = list(
                 p = p.year,
                 Psi = Psi.dot,
                 Epsilon = Epsilon.random.yr))
  
  # Model 4 - annual variation in epsilon & gamma
  mod.4 = mark(nso.pr.EG, nso.ddl,
               model.parameters = list(
                 p = p.year, 
                 Psi = Psi.dot,
                 Epsilon = Epsilon.dot,
                 Gamma = Gamma.dot))
  
  # Model 5 - annual variation in epsilon & gamma
  mod.5=mark(nso.pr.EG, nso.ddl,
             model.parameters = list(
               p = p.year,
               Psi = Psi.dot,
               Epsilon = Epsilon.year,
               Gamma = Gamma.year))
  
  # Model 6 - Work with different parameterization for which model psi by year.
  #  For this model, psi & gamma are held constant across years
  mod.6 = mark(nso.pr.PG, 
               model.parameters = list(
                 p = p.year,
                 Psi = Psi.dot,
                 Gamma = Gamma.dot))
  
  # Collect the models to run when the function is called
  return(collect.models() )
}
```

### Run the models and examine the output

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

```{r}
nso.models <- fit.nso.models()
```

### Examine model-selection table

Once the models are run, we can examine the model-selection table. We can also examine model-specific output as shown below.

```{r}
options(width = 160)
nso.models
```

### Look at output from top model

```{r}
# View the model names; here they are just mod.1, mod.2, ..., mod.6,
# so we don't really need to rename the model like we did in some labs
names(nso.models)

# examine the output from top-ranked model (#6) and
nso.models$mod.6$results$real

# for mod.6, estimates of Epsilon are derived
nso.models$mod.6$results$derived
```

### Model-Specific Estimates of Psi for each year

Here, we look at model-specific estimates of psi in year 1 and for each of the subsequent years. Remember, the latter are derived parameters in models 1-5.


```{r}
# store & view model-specific estimates of Psi
(psi.mod.1 = nso.models$mod.1$results$derived$`psi Probability Occupied`)
(psi.mod.2 = nso.models$mod.2$results$derived$`psi Probability Occupied`)
(psi.mod.3 = nso.models$mod.3$results$derived$`psi Probability Occupied`)
(psi.mod.4 = nso.models$mod.4$results$derived$`psi Probability Occupied`)
(psi.mod.5 = nso.models$mod.5$results$derived$`psi Probability Occupied`)
# for mod.6, Psi-hat is held constant across all years, so
# we can just repeat that row 6 times (which takes a bit of work
# to get the results to look like the results for other models) 
psi.mod.6 <- do.call("rbind", 
                     replicate(5, 
                               nso.models$mod.6$results$real[1, -c(2, 5:6)],
                               simplify = FALSE))
row.names(psi.mod.6) <- 1:5
psi.mod.6
```

### Graph the estimates of psi for each year from the various models

First, we'll store the estimates in a data.frame that records which model they came from. Then, we'll use the 'ggplot2' package to plot the estimates.

```{r, fig.height=6, fig.width=8}
psi.hats <- rbind(psi.mod.1,
                 psi.mod.2,
                 psi.mod.3,
                 psi.mod.4,
                 psi.mod.5,
                 psi.mod.6)
# create a variable with the model number
psi.hats$model = rep(1:6, each = 5)
# create a variable for year
psi.hats$year = rep(1997:2001, 6)
head(psi.hats)

# make the graph
library(ggplot2)
ggplot(psi.hats, aes(x = year, y = estimate, group = factor(model))) + 
  geom_line(size = 1.5, aes(colour = factor(model))) + 
  geom_ribbon(aes(ymin = lcl, ymax = ucl), alpha = 0.2) + 
  scale_colour_brewer(palette = "Set1", name = "Model #") + 
  ylim(0, 1) +
  theme(legend.position = c(0.2, 0.02),  
        legend.justification = c(1, 0)) + 
  xlab("Year") + 
  ylab("Estimated Occupancy Rate") 
```
