--- title: "Lab 02 - Analysis in RMark" author: "WILD 502 - Jay Rotella" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk\$set(echo = TRUE) ``` ### Introduction In lab this week, we have known fate data for 1 interval of time and the input file also contains 4 individual covariates (Area, Sex, Mass, & Length). Here, we'll use RMark and work with the same models that you evaluated directly in MARK for lab 2. ### Import and process the data ```{r} library(RMark) library(ggplot2) # bring in the data file fawns <- convert.inp("http://www.montana.edu/rotella/documents/502/lab02-fawns.inp", covariates = c("area", "sex", "mass", "length")) # add "length.sq" variable for use in quadratic models fawns\$length.sq <- (fawns\$length)^2 head(fawns) # convert area & sex variables to factors so that we can treat them as # grouping variables in the analysis to ease working with output later fawns\$area <- factor(fawns\$area, levels = c(0, 1), labels = c("control", "treatment")) fawns\$sex <- factor(fawns\$sex, levels = c(0, 1), labels = c("female", "male")) # Store versions of fate for use in final plotting # first a factor version that has nice labels fawns\$fate <- factor(as.numeric(fawns\$ch), levels = c(11, 10), labels = c("died", "lived")) # second, a numeric version that's 0 = dead, 1 = lived fawns\$Fate <- 2 - as.numeric(fawns\$fate) summary(fawns) # we know from other work on this dataset that there is 1 record with # a very high value for mass, which we want to drop fawns2 = fawns[fawns\$mass < 70, ] ``` ### Process the data for use in RMark ```{r} # Process data fawns.processed = process.data(fawns2, model = "Known", groups = c("area", "sex")) # # Create default design data fawns.ddl = make.design.data(fawns.processed) ``` ### Build function for creating models Here, we write a function that establishes which models are to be run in MARK. ```{r} # setup a function run.fawns <- function() { # Define range of models for S S.dot = list(formula = ~ 1) S.area = list(formula = ~ area) S.mass = list(formula = ~ mass) S.length = list(formula = ~ length) S.length.sq = list(formula = ~ length + length.sq) S.sex = list(formula = ~ sex) S.area.length = list(formula = ~ area + length) S.length.mass = list(formula = ~ length + mass) # sex-specific intercepts & common slope S.sex.length = list(formula = ~ sex + length) # 1 intercept & sex-specific slopes S.length.by.sex = list(formula = ~ sex:length) # sex-specific intercepts & slopes S.sex.x.length = list(formula = ~ sex * length) # Create model list model.list = create.model.list("Known") # NOTE: to avoid having all the output for each model appear when you # call the function, add ', output=FALSE' after 'ddl=fawns.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. fawns.results = mark.wrapper(model.list, data = fawns.processed, ddl = fawns.ddl, invisible = TRUE, threads = 2) # Return model table and list of models return(fawns.results) } ``` ### 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. 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} fawns.results = run.fawns() ``` ### Look at the model-selection table ```{r} fawns.results ``` ### Work with the top model You would not need to re-run the top model since we've already run it. But, by doing it this way, you can see how you can to run a model of interest without using a function as was done earlier. ```{r} # here's how to run a model outside the function top <- mark(fawns.processed,ddl = fawns.ddl, model = "Known", model.parameters = list("S" = list(formula = ~ sex * length))) # look at beta-hats top\$results\$beta # create values of length to use for predictions min.length = min(fawns2\$length) max.length = max(fawns2\$length) length.values = seq(from = min.length, to = max.length, length = 100) # determine which parameter indices go with males and females fawns.ddl # make predictions for rows of 'fawns.ddl' associated with # females (par.index=1) and males (par.index=3) pred.top <- covariate.predictions(top, data = data.frame(length = length.values, length.sq = length.values ^ 2), indices = c(1, 3)) # store values of sex in pred.top female.rows <- which(pred.top\$estimates\$par.index == 1) male.rows <- which(pred.top\$estimates\$par.index == 3) pred.top\$estimates\$sex <- NA pred.top\$estimates\$sex[female.rows] <- "female" pred.top\$estimates\$sex[male.rows] <- "male" head(pred.top\$estimates) ``` ### Working with output: Plotting Here, we plot estimates of "S" for male and female fawns of different lengths. We also want to plot CI's for the values. To do so, we use the predicted values we created in the last step. Here, I use the ggplot2 package, which makes nice graphs. ```{r fig.width=7, fig.height=6, warning=FALSE} pred <- pred.top\$estimates # build and store the plot in object 'p' p <- ggplot(pred, aes(x = length, y = estimate, group = sex)) + geom_line(size = 1.5, aes(colour = sex)) + geom_ribbon(aes(ymin = lcl, ymax = ucl), alpha = 0.2) + scale_colour_brewer(palette = "Set1") + theme(legend.position = c(0.15, 0.8), legend.justification = c(1, 0)) + xlab("Body Length (cm)") + ylab("Estimated Survival Rate") + geom_jitter(data = fawns2, aes(y = Fate, color = sex), height = 0.025) # print the plot p # 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) ```