--- title: "Lab 10 - Analyses in RMark" author: "WILD 502 - Jay Rotella" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk\$set(echo = TRUE) ``` ## Single-Species Single-Season Occupancy Models Here, we'll run the single-season occupancy models for data on Weta that you analyzed directly in MARK in lab 10. Here, we're dealing with 2 types of parameters: 'Psi' & 'p'. ### Bring in the data ```{r} library(RMark) # you can ignore warning about final line issue that might appear weta <- convert.inp("http://www.montana.edu/rotella/documents/502/wetaTwoGroups.inp", group.df = data.frame(Browse = c("Browsed", "Unbrowsed")), covariates = c( "obs11", "obs12", "obs13", "obs14", "obs15", "obs21", "obs22", "obs23", "obs24", "obs25")) # store the number of sites n.sites <- length(weta\$ch) # store the maximum number of visits n.visits <- max(nchar(weta\$ch)) # make a factor that records whether observer 1, 2, or 3 visited a site # on each visit. A '9' was used to indicate that no observer visited the site. # So, if there's a 9 for either observer 1 or 2 on a given date, we # assign a '.' for that day obs <- matrix(NA, n.sites, n.visits) for (i in 1:n.sites) { obs[i, 1] = ifelse(weta\$obs11[i] == 1, 1, ifelse(weta\$obs21[i] == 1, 2, ifelse(weta\$obs11[i] == 9, ".", 3))) obs[i, 2] = ifelse(weta\$obs12[i] == 1, 1, ifelse(weta\$obs22[i] == 1, 2, ifelse(weta\$obs12[i] == 9, ".", 3))) obs[i, 3] = ifelse(weta\$obs13[i] == 1, 1, ifelse(weta\$obs23[i] == 1, 2, ifelse(weta\$obs13[i] == 9, ".", 3))) obs[i, 4] = ifelse(weta\$obs14[i] == 1, 1, ifelse(weta\$obs24[i] == 1, 2, ifelse(weta\$obs14[i] == 9, ".", 3))) obs[i, 5] = ifelse(weta\$obs15[i] == 1, 1, ifelse(weta\$obs25[i] == 1, 2, ifelse(weta\$obs15[i] == 9, ".", 3))) } head(weta) head(obs) # can drop obs11, obs12, ..., obs25 from weta since not used weta <- weta[,-c(4:13)] # add 'Obs' values for each day to 'weta' weta\$Obs1 = as.factor(obs[, 1]) weta\$Obs2 = as.factor(obs[, 2]) weta\$Obs3 = as.factor(obs[, 3]) weta\$Obs4 = as.factor(obs[, 4]) weta\$Obs5 = as.factor(obs[, 5]) head(weta) ``` ### Build function to process data & create models First, the function puts the observer and time covariates into proper form for use in modeling and processes the data. Second, thefunction decleares the structures for 'p' and 'Psi' that will be used in all-possible combinations of the structures for each parameter type when the function is executed. ```{r} # Create function to fit the 18 models in the book fit.weta.models <- function() { # use 'make.time.factor' to create time-varying dummy variables Obs1 and Obs2 # observer 3 is used as the intercept to match what we did in MARK weta = make.time.factor(weta, "Obs", 1:5, intercept = 3) # Process data weta.process = process.data(weta, model = "Occupancy", groups = "Browse") weta.ddl = make.design.data(weta.process) # time factor variable copied to 'Day' to match name used in the lab weta.ddl\$p\$Day = weta.ddl\$p\$time # Define p models p.dot = list(formula = ~ 1) p.browse = list(formula = ~ Browse) p.day = list(formula = ~ Day) p.obs = list(formula = ~ Obs1 + Obs2) p.day.obs = list(formula = ~ Day + Obs1 + Obs2) p.day.browse = list(formula = ~ Day + Browse) p.obs.browse = list(formula = ~ Obs1 + Obs2 + Browse) p.day.obs.browse = list(formula = ~ Day + Obs1 + Obs2 + Browse) # Define Psi models Psi.dot = list(formula = ~ 1) Psi.browse = list(formula = ~ Browse) # Create model list cml = create.model.list("Occupancy") # Run and return marklist of models return(mark.wrapper(cml, data = weta.process, ddl = weta.ddl)) } ``` ### Run the models and examine the output It's very simple to then run the function and each of the models. As implemented here, the output from each of the models appears in the console where it can be reviewed when the function written above is called. If you don't want the output from every model to appear automatically, see the note above the "mark.wrapper" command for how to set the option of suppressing the output. One can also examine model-specific output in other ways as shown below. ```{r} weta.models <- fit.weta.models() ``` ### Examine model-selection table Once the models are run, we can examine the model-selection table. We can also examine model-specific output. ```{r} options(width = 160) weta.models ``` ### Use AIC instead of AICc for model selection Here, we use AIC so that we can compare direclty with results presented in the Occupancy book. You can see that things match (as they should) for most models but that something appears amiss in the book with the Psi(.)p(Obs+Browse) model as it should have 5 parameters but is only listed as having 4 in the book (it is assigned k=5 here). ```{r} options(width = 160) # Store and view a model table that uses AIC rather than AICc AIC.table = weta.models AIC.table\$model.table = model.table(weta.models, use.AIC = TRUE) AIC.table ``` ### Look at output from top model ```{r} # View the model names names(weta.models) # examine the output from top-ranked model (#7) and # store top model in object with short name top <- weta.models\$p.day.obs.Psi.browse # look at summary of top model's output summary(top, showall = FALSE) ``` ### Model Average Psi Here, model-averaging is done for Psi in each habitat setting. ```{r} Mod.Avg.Psi <- model.average(weta.models, "Psi") Mod.Avg.Psi ``` ### Additional Work There is much more that can be done in RMark. And, in the help for the RMark package, the weta example has been worked out in detail. You can see the details by typing "?weta" after loading the RMark package. If you read through the example, you'll see that much of what's provided on pages 116-122 of the Occupancy book is worked out in the RMark 'weta' example.