WILD 502: Lab 04 - CJS models part II
========================================================
#### Jay Rotella, Ecology Dept., Montana State University

Here, you'll find code for running each of the models used in lab this week with the data for swifts found in the "swift.inp" data set. You'll also find the model-selection table, model-specific results, and results based on some follow-up work with the output from the top model.

### Bring in the data
```{r, warning=FALSE}
library(RMark)

swift=convert.inp("http://www.montana.edu/rotella/documents/502/AA.inp",
      group.df=data.frame(colony=c("Poor","Good")), covariates=NULL)
head(swift)
```

### Process the data

```{r}
options(width = 150)
# Process data
swift.pr=process.data(swift, begin.time = 1, model = "CJS", groups=("colony"))
#  
# Create default design data
swift.ddl=make.design.data(swift.pr)

# Examine the design data. Because the output is so long, 
# only the head of each piece of design data is shown
head(swift.ddl$Phi)
head(swift.ddl$p)
```

### Build function for creating models

Here, we set up a function that contains the structures for 'Phi' and 'p'. It will run all combinations of the structures for each parameter type when executed.


```{r}
run.swift=function()
{
#
# Process data
swift.pr=process.data(swift, begin.time = 1, model = "CJS", groups=("colony"))
#  
# Create default design data
#
swift.ddl=make.design.data(swift.pr)
#
#  Define range of models for Phi
#
Phi.dot=list(formula=~1)
Phi.time=list(formula=~time)
Phi.colony=list(formula=~colony)
Phi.colonyXtime=list(formula=~colony*time)
#
#  Define range of models for p
#
p.dot=list(formula=~1)
p.time=list(formula=~time)
#
# Run all pairings of models for phi & p.

# NOTE: if you do not want to see the output for each model, you can add the text  
# ', output=FALSE' after 'ddl=swift.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 a
# lot of models! 
# I added 'adjust=FALSE' below to prevent RMark from adjusting the number of 
# parameters, which would be problematic for models like phi(t),p(t) or
# phi(c*t),p(t) in which not all parameters can be estimated.
swift.model.list=create.model.list("CJS")
swift.results=mark.wrapper(swift.model.list,
      data=swift.pr,ddl=swift.ddl, adjust=FALSE)
#
# Return model table and list of models
#
return(swift.results)
}
```

### 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. 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}
swift.results=run.swift()
```


### 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 = 150)
swift.results
names(swift.results)

# examine the output from top-ranked model (model #2) and
# store top model in object with short name 
top=swift.results$Phi.colony.p.time

# look at output
top$results$beta
top$results$real[c(1,2),]
```

### Adjust model-selection results for estimated overdispersion

Here, we adjust using an estimate of 1.035414. To estimate this quantity, you'd need to work in MARK directly. Notice that these SEs have been inflated by the square root of 1.035414.
```{r}
options(width = 150)
QAICc.table=adjust.chat(1.035414, swift.results)
QAICc.table

# store top model in object with short name 
top.qaic=QAICc.table$Phi.colony.p.time

# look at estimates of real parameters and notice that the SEs have been adjusted
phi.hat.qaic=get.real(top.qaic, "Phi", se=TRUE)
# only need to print out rows 1 and 29 to get estimate for each colony
phi.hat.qaic[c(1,29),c(1:6)]

# compare with results obtained from AICc
top$results$real[c(1,2),1:6]
```

### Work with the results: estimate difference in survival

In lab, we were interested in estimating the difference in survival rates for the 2 colonies and putting a confidence interval on that estimated difference. Here, we use the output from the model and the delta method to develop the SE and 95% CI for expressing the uncertainty for the estimated difference. First, the code uses the 'deltamethod' function in the 'msm' package to calculate the SE for the estimated difference. If you don't have the package 'msm' and you're trying to use the code below, you'll need to install the package and then load it before executing the code below. The code also shows the actual (& algebraically simple) calculation for estimating that SE  based on the variance and covariance of the estimated survival rates for the 2 colonies. For the variance of a difference, we sum the variances and subtract 2 times the covariance.

```{r}
# store var-cov for real estimates
rr=get.real(top.qaic, "Phi", se=TRUE,vcv=TRUE)
sigma=rr$vcv.real # store var-cov in "sigma"

# load package 'msm' that implements the delta method
# then run delta method on difference in survival
library(msm)
s.good = rr$estimates$estimate[1]  # first row w/ output for good colony
s.poor = rr$estimates$estimate[29] # first row w/ output for poor colony
Diff= s.good - s.poor 
seDiff=deltamethod(~x1-x2,c(s.good,s.poor),sigma)
lclDiff=Diff-1.96*seDiff
uclDiff=Diff+1.96*seDiff
round(c(Diff,seDiff,lclDiff,uclDiff),3)
# compare SE to that obtained by direct calculation
# where: variance of difference = var(A)+var(B)-2*cov(A,B)
(se.diff.direct=sqrt(sum(diag(sigma))-2*sigma[1,2]))

# cleanup files created by MARK with the 2 lines below
#rm(list=ls(all=TRUE))
#cleanup(ask = FALSE)
```