---
title: 'Lab 05: Analysis in RMark'
author: "WILD 502 - Jay Rotella"
output:
  html_document: default
---

## Bring in the Data and Process it

Here, the data are for imaginary swamp squirrels. There are data for 6 sampling occasions and for 2 groups: squirrels first captured when they were juveniles or when they were older. There are also 2 continuous covariates: birth (birth date) and tail (tail length). The code creates tail length squared and provides information on precipitation levels during capture work each year, which might have affected capture probability (*p*).

```{r}
library(knitr)
library(ggplot2)
library(RMark)
sq <- convert.inp("http://www.montana.edu/rotella/documents/502/SwampSquirrels.inp", 
                  group.df = data.frame(age = c("juv","older")),
                  covariates = c("birth", "tail"))
sq$tail.sq <- sq$tail^2

sq.pr <- process.data(sq, model = "CJS",
                      age.var = 1,
                      initial.ages = c(0, 1),
                      groups = ("age"))

# Create default design data  
# For Phi, values are 0 or older; for p the values are 1 or older
# This provides an individual covariate for each parameter type
# that changes as the animal ages
sq.ddl = make.design.data(sq.pr,
                          parameters = list(Phi = list(age.bins = c(0, 1, 6)),
                                            p = list(age.bins = c(1, 2, 6))),
                          right = FALSE)

# Create precipitation variable
sq.ddl$p$precip = 0 # Start with 0 for all and then over-write
sq.ddl$p$precip[sq.ddl$p$time == 2] = 0.11
sq.ddl$p$precip[sq.ddl$p$time == 3] = 1.38
sq.ddl$p$precip[sq.ddl$p$time == 4] = 1.12
sq.ddl$p$precip[sq.ddl$p$time == 5] = 1.58
sq.ddl$p$precip[sq.ddl$p$time == 6] = 1.87
```

## Set up Function for Competing Models

```{r}
run.sq = function() {
  #  Define range of models for Phi
  #
  Phi.dot = list(formula =  ~ 1)
  Phi.age = list(formula =  ~ age)
  Phi.age.birth = list(formula =  ~ age + birth)
  Phi.age.tail = list(formula =  ~ age + tail)
  Phi.age.tail.sq = list(formula =  ~ age + tail + tail.sq)
  Phi.age.birth.tail = list(formula =  ~ age + birth + tail)
  Phi.age.birth.tail.sq = list(formula =  ~ age + birth + tail + tail.sq) 
  
  #  Define range of models for p
  p.dot = list(formula =  ~ 1)
  p.time = list(formula =  ~ -1 + time)
  p.precip = list(formula = ~ precip)
  
  # Create models for all combinations of phi & p
  sq.model.list = create.model.list("CJS")
  
  # NOTE: to avoid having all the output for each model appear when you
  # call the function, add ', output=FALSE' after 'ddl=sq.ddl' below.
  # Here, I don't do that so you can see the output for each model,
  # but this might not be desired if you have many models.
  sq.results = mark.wrapper(sq.model.list,
                            data = sq.pr, ddl = sq.ddl)
  #
  # Return model table and list of models
  #
  return(sq.results)
}
```

## Run the Models 

The code below runs all combinations of structures in the function for $\phi$ and $p$. The code also produces a model-selection table.

```{r}
sq.res <- run.sq()
sq.res
names(sq.res)
```

## Examine Output for Models of Interest

```{r}
top <- sq.res$Phi.age.birth.tail.sq.p.time
top$results$beta
```

## Make Predictions for $\phi$ from Top Model

```{r}
min.tail <- min(sq$tail)
max.tail <- max(sq$tail)
tail.values <- seq(from = min.tail, to = max.tail, by = 0.25)
birth.values <- c(quantile(sq$birth, 0.05), 
                  mean(sq$birth),
                  quantile(sq$birth, 0.95))
pred.dat <- expand.grid(birth = birth.values,
                       tail = tail.values)
pred.dat$tail.sq <- pred.dat$tail^2

# make predictions for rows of 'sq.ddl' associated with juveniles, 
# i.e., (par.index=1, 2) = a juvenile and an adult 
pred.top <- covariate.predictions(top, data = pred.dat, indices = c(1, 2))

# view head of the prediction data.frame
head(pred.top$estimates)
```

### Plot phi versus tail length: make the plot

Here, you can see how to make the plot using the 'ggplot2' package

```{r}
pred = pred.top$estimates
# use parameter index to create age_class variable
pred$age_class = ifelse(pred$par.index == 1, "juv", "older")
# Store information on birth date values
pred$birthdate = ifelse(pred$birth == mean(sq$birth), "b - bd_avg",
    ifelse(pred$birth == quantile(sq$birth, 0.05), "a - bd_.05",
           "c - bd_.95"))
# build and store the plot in object 'p'
ggplot(pred, aes(x = tail, y = estimate, group = age_class)) + 
  geom_line(aes(color = age_class), size = 1.5) + 
  geom_ribbon(aes(ymin = lcl, ymax = ucl), alpha = 0.2) + 
    xlab("Tail Length") + ylab("Apparent Survival Rate") +
    ylim(0, 1) + facet_wrap( ~ birthdate) +
  theme(legend.position="top")
```

### Clean up files created

```{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)
```