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