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

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

## Estimating Abundance for a Closed Population

Here, you'll find code for running the models used in lab this week with the "ClosedSim.inp" data set. These models are closed models, and the code below gives you a chance to see how to run closed models in RMark. It should not be too hard for you to modify the code below to run all the models for Lab 7 on "capture.inp". If you try that,  you'll see just how convenient it can be to apply code developed for 1 data set to another as long as the model sets for the 2 projects are similar enough.

### Bring in the data

```{r}
# Lab 07 for WILD 502 at Montana State University
library(RMark)
caps <- convert.inp(
  "http://www.montana.edu/rotella/documents/502/ClosedSim.inp",
  group.df = NULL,
  covariates = NULL)
head(caps)
tail(caps)
```

### Manipulate the data

Before we can run all of the models for this lab, we need to add a few covariates that will be used to constrain the estimated capture probabilities in various ways.

```{r}
# Process data
caps.pr <- process.data(caps, begin.time = 1, model = "Closed")

# Create default design data
caps.ddl <- make.design.data(caps.pr)

# add columns to ddl data for p for M(bh) models
# 1. add 2 categories of time (1st occ.[0] or later[1])
caps.ddl$p$t2 <- 0
caps.ddl$p$t2 = ifelse(caps.ddl$p$Time == 0, 0, 1)
# 2. add 2nd category of time (1st or 2nd occ.[0] or later[1])
caps.ddl$p$t3 <- 0
caps.ddl$p$t3 <- ifelse(caps.ddl$p$Time < 2, 0, 1)

# add columns to ddl data for c for M(bh) models
# Not used here but must be added for things to work
# 1. add 2 categories of time (1st occ.[0] or later[1])
caps.ddl$c$t2 <- 0
caps.ddl$c$t2 <- ifelse(caps.ddl$c$Time == 0, 0, 1)
# 2. add 2nd category of time (1st or 2nd occ.[0] or later[1])
caps.ddl$c$t3 <- 0
caps.ddl$c$t3 <- ifelse(caps.ddl$c$Time < 2, 0, 1)

# Define weekly temperature data
caps.ddl$p$temp <- c(-1.5, 0.2, 0.1, -0.4, 1.5, 1.1)
caps.ddl$c$temp <- c(0.2, 0.1, -0.4, 1.5, 1.1)
# examine the data
caps.ddl
```

### Build function for creating models

Here, we set up the structures for 'p' and 'c'. I use the "share=TRUE" or "share=FALSE" options in each of the structures to indicate whether 'p' & 'c' should share the same columns of the design matrix or not. Although this is not necessary for all of the structures below, it does add a covariate "c" to the design data with c=0 for rows pertaining to parameter 'p' and c=1 for rows pertaining to parameter 'c. This is nice as it gives us the opportunity to build some of the additive structures we're interested in. We can then use the covariate "c" in formula statements if we want to. But, we don't have to include that covariate if we don't want to (e.g., see the p.dot structure).

```{r}
# function for running set of models for phi and for p
run.caps <- function() {

  # Define parameters for p and c 
  # Note: "share=TRUE" indicates that 'p' & 'c'
  #  share the same columns in the design matrix
  p.dot = list(formula =  ~ 1, share = TRUE)
  p.time = list(formula =  ~ time, share = TRUE)
  p.time.behav.add = list(formula =  ~ time + c, share = TRUE)
  p.dot.behav = list(formula =  ~ 1, share = FALSE)
  p.bh.2p = list(formula =  ~ t2, share = FALSE)
  p.bh.3p = list(formula =  ~ t2 + t3, share = FALSE)
  p.temp = list(formula =  ~ temp, share = TRUE)
  
  # Create competing models based on strcutures for 'p' & 'c'
  caps.model.list = create.model.list("Closed")
  
  # NOTE: if you do not want to see the output for each model, add  
  # ', output=FALSE' after 'ddl=caps.ddl' below. 
  caps.results = mark.wrapper(caps.model.list,
                              data = caps.pr, 
                              ddl = caps.ddl)
                              
  # Return model table and list of models
  return(caps.results)
}
```

### Run the models and examine the output

It's worth noting that 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. And, if you don't want to see the output from every model appear in the output here, see the note above for the use of the "mark.wrapper" command where one has the option of suppressing the output.

```{r}
caps.results <- run.caps()
```

### Examine model-selection table

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

```{r}
caps.results
```

### Examine output from top model

```{r}
# examine model names and find the name of the top model
names(caps.results)

# examine the output from top-ranked model (#7)
caps.results$p.time.behav.add$results

# if you're running this code, it's wise to cleanup files
# you can use the next 2 lines without comments to do that

#rm(list=ls(all=TRUE))
#cleanup(ask = FALSE)
```
