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