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