---
title: "Lab03 - Analysis in RMark"
author: "WILD 502 - Jay Rotella"
output: html_document
---

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

## CJS model for data on male European dippers

Here, we'll run the 4 models for male European dippers where apparent survival rate (*phi*) and detection rate (*p*) are each either allowed to vary by time (t) or constrained to be constant across years (. or dot). Here, there are only 2 structures for each model, so I don't write a function to create all possible combinations. Instead, the models are coded individually, which also allows me to easily adjust the parameter count to the correct value of 11 for the "Phi(time), p(time)" model. If I don't do that, RMark will adjust it to 12. Notice that all CJS models here are run using the logit link. With the logit link, the parameter count will be too low for the "Phi(dot), p(time)" model because some of the estimates of *p* are 1.0 for the model, which causes difficulty with the routine MARK uses for counting parameters. RMark's default behavior is to adjust the number of parameters to the number of columns in the design matrix, which is what we want for the "Phi(dot), p(time)" model.

### Set up the data & run some models

```{r}
library(knitr)
library(RMark)

edm <- convert.inp("http://www.montana.edu/rotella/documents/502/ed-males.inp",
      group.df = NULL)
head(edm)

# process the data
edm.pr <- process.data(edm, model = "CJS")

# Setup model structures for each parameter
Phi.dot = list(formula = ~ 1)
# force use of an identity matrix by putting '-1' in formula
Phi.time = list(formula = ~ -1 + time) 

p.dot = list(formula = ~ 1)
# force use of an identity matrix by putting '-1' in formula
p.time = list(formula = ~ -1 + time) 


Phi.dot.p.dot = mark(edm.pr, 
                     model.parameters = list(Phi = Phi.dot, p = p.dot))
Phi.time.p.dot = mark(edm.pr, 
                      model.parameters = list(Phi = Phi.time, p = p.dot))
Phi.dot.p.time = mark(edm.pr, 
                      model.parameters = list(Phi = Phi.dot, p = p.time))
Phi.time.p.time = mark(edm.pr, 
                      model.parameters = list(Phi = Phi.time, p = p.time))
Phi.time.p.time = adjust.parameter.count(Phi.time.p.time, 11)
```

### Run the models

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. 

```{r}
edm.results <- model.table(model.list = c("Phi.dot.p.dot", 
                                          "Phi.dot.p.time",
                                          "Phi.time.p.dot",
                                          "Phi.time.p.time"),
                           model.name = FALSE)
```

### Look at the model-selection table
```{r}
kable(edm.results, digits = 3)
```

### Look at output from top model

Although you've already seen the output above, you might like to see how to get the real parameter estimates from the top model. Here's a way that works.

```{r}
Phi.dot.p.dot$results$beta

Phi.dot.p.dot$results$real
```

### Clean up the files

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

```{r}
#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)
```