## CJS model for data on male European dippers

Here, we’ll run the 4 models for male European dippers where apparent survival rate (phi) and detection rate (p) are each either allowed to vary by time (t) or constrained to be constant across years (. or dot). Here, there are only 2 structures for each model, so I don’t write a function to create all possible combinations. Instead, the models are coded individually, which also allows me to easily adjust the parameter count to the correct value of 11 for the “Phi(time), p(time)” model. If I don’t do that, RMark will adjust it to 12. Notice that all CJS models here are run using the logit link. With the logit link, the parameter count will be too low for the “Phi(dot), p(time)” model because some of the estimates of p are 1.0 for the model, which causes difficulty with the routine MARK uses for counting parameters. RMark’s default behavior is to adjust the number of parameters to the number of columns in the design matrix, which is what we want for the “Phi(dot), p(time)” model.

### Set up the data & run some models

library(knitr)
library(RMark)
## This is RMark 2.2.4
##  Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
edm <- convert.inp("http://www.montana.edu/rotella/documents/502/ed-males.inp",
group.df = NULL)
head(edm)
##        ch freq
## 1 1111110    1
## 2 1111000    1
## 3 1100000    1
## 4 1100000    1
## 5 1100000    1
## 6 1100000    1
# process the data
edm.pr <- process.data(edm, model = "CJS")

# Setup model structures for each parameter
Phi.dot = list(formula = ~ 1)
# force use of an identity matrix by putting '-1' in formula
Phi.time = list(formula = ~ -1 + time)

p.dot = list(formula = ~ 1)
# force use of an identity matrix by putting '-1' in formula
p.time = list(formula = ~ -1 + time)

Phi.dot.p.dot = mark(edm.pr,
model.parameters = list(Phi = Phi.dot, p = p.dot))
##
## Output summary for CJS model
## Name : Phi(~1)p(~1)
##
## Npar :  2
## -2lnL:  318.4938
## AICc :  322.5527
##
## Beta
##                  estimate        se        lcl       ucl
## Phi:(Intercept) 0.2648275 0.1446688 -0.0187234 0.5483784
## p:(Intercept)   2.4862984 0.5120845  1.4826127 3.4899840
##
##
## Real Parameter Phi
##
##           1         2         3         4         5         6
## 1 0.5658226 0.5658226 0.5658226 0.5658226 0.5658226 0.5658226
## 2           0.5658226 0.5658226 0.5658226 0.5658226 0.5658226
## 3                     0.5658226 0.5658226 0.5658226 0.5658226
## 4                               0.5658226 0.5658226 0.5658226
## 5                                         0.5658226 0.5658226
## 6                                                   0.5658226
##
##
## Real Parameter p
##
##           2         3         4         5         6         7
## 1 0.9231757 0.9231757 0.9231757 0.9231757 0.9231757 0.9231757
## 2           0.9231757 0.9231757 0.9231757 0.9231757 0.9231757
## 3                     0.9231757 0.9231757 0.9231757 0.9231757
## 4                               0.9231757 0.9231757 0.9231757
## 5                                         0.9231757 0.9231757
## 6                                                   0.9231757
Phi.time.p.dot = mark(edm.pr,
model.parameters = list(Phi = Phi.time, p = p.dot))
##
## Output summary for CJS model
## Name : Phi(~-1 + time)p(~1)
##
## Npar :  7
## -2lnL:  315.4939
## AICc :  330.0567
##
## Beta
##                 estimate        se        lcl       ucl
## Phi:time1      0.4512439 0.6301250 -0.7838010 1.6862889
## Phi:time2     -0.1673376 0.4022469 -0.9557415 0.6210663
## Phi:time3     -0.0159046 0.3404187 -0.6831252 0.6513160
## Phi:time4      0.4596695 0.3426110 -0.2118482 1.1311871
## Phi:time5      0.2974186 0.3117960 -0.3137017 0.9085388
## Phi:time6      0.5406789 0.3424413 -0.1305061 1.2118640
## p:(Intercept)  2.4350067 0.5162517  1.4231532 3.4468601
##
##
## Real Parameter Phi
##
##           1         2         3         4         5         6
## 1 0.6109349 0.4582629 0.4960239 0.6129358 0.5738113 0.6319703
## 2           0.4582629 0.4960239 0.6129358 0.5738113 0.6319703
## 3                     0.4960239 0.6129358 0.5738113 0.6319703
## 4                               0.6129358 0.5738113 0.6319703
## 5                                         0.5738113 0.6319703
## 6                                                   0.6319703
##
##
## Real Parameter p
##
##           2         3         4         5         6         7
## 1 0.9194581 0.9194581 0.9194581 0.9194581 0.9194581 0.9194581
## 2           0.9194581 0.9194581 0.9194581 0.9194581 0.9194581
## 3                     0.9194581 0.9194581 0.9194581 0.9194581
## 4                               0.9194581 0.9194581 0.9194581
## 5                                         0.9194581 0.9194581
## 6                                                   0.9194581
Phi.dot.p.time = mark(edm.pr,
model.parameters = list(Phi = Phi.dot, p = p.time))
##
## Note: only 5 parameters counted of 7 specified parameters
## AICc and parameter count have been adjusted upward
##
## Output summary for CJS model
## Name : Phi(~1)p(~-1 + time)
##
## -2lnL:  316.1166
##
## Beta
##                  estimate           se           lcl          ucl
## Phi:(Intercept)  0.225255    0.1388719 -4.693380e-02 4.974439e-01
## p:time2          1.385924    1.0473518 -6.668858e-01 3.438734e+00
## p:time3         19.815873 9119.0338000 -1.785349e+04 1.789312e+04
## p:time4          2.062320    0.9739255  1.534258e-01 3.971214e+00
## p:time5          2.646892    1.0057624  6.755977e-01 4.618187e+00
## p:time6          2.683137    1.0014551  7.202854e-01 4.645989e+00
## p:time7         19.574117 6166.4869000 -1.206674e+04 1.210589e+04
##
##
## Real Parameter Phi
##
##           1         2         3         4         5         6
## 1 0.5560768 0.5560768 0.5560768 0.5560768 0.5560768 0.5560768
## 2           0.5560768 0.5560768 0.5560768 0.5560768 0.5560768
## 3                     0.5560768 0.5560768 0.5560768 0.5560768
## 4                               0.5560768 0.5560768 0.5560768
## 5                                         0.5560768 0.5560768
## 6                                                   0.5560768
##
##
## Real Parameter p
##
##           2 3         4         5         6 7
## 1 0.7999407 1 0.8871866 0.9338192 0.9360243 1
## 2           1 0.8871866 0.9338192 0.9360243 1
## 3             0.8871866 0.9338192 0.9360243 1
## 4                       0.9338192 0.9360243 1
## 5                                 0.9360243 1
## 6                                           1
Phi.time.p.time = mark(edm.pr,
model.parameters = list(Phi = Phi.time, p = p.time))
##
## Note: only 10 parameters counted of 12 specified parameters
##
## AICc and parameter count have been adjusted upward
##
## Output summary for CJS model
## Name : Phi(~-1 + time)p(~-1 + time)
##
## -2lnL:  313.0805
##
## Beta
##             estimate        se        lcl        ucl
## Phi:time1  0.8329093 0.9705503 -1.0693693  2.7351879
## Phi:time2 -0.3101550 0.3969582 -1.0881930  0.4678830
## Phi:time3  0.0211523 0.3500130 -0.6648731  0.7071777
## Phi:time4  0.4447981 0.3522276 -0.2455680  1.1351641
## Phi:time5  0.2851896 0.3171429 -0.3364105  0.9067898
## Phi:time6  1.1734788 0.0000000  1.1734788  1.1734788
## p:time2    0.9315578 1.1041357 -1.2325481  3.0956638
## p:time3   23.2177800 0.0000000 23.2177800 23.2177800
## p:time4    2.3051463 1.0376230  0.2714053  4.3388873
## p:time5    2.5477074 1.0301727  0.5285689  4.5668459
## p:time6    2.6798792 1.0270738  0.6668145  4.6929439
## p:time7    1.1733609 0.0000000  1.1733609  1.1733609
##
##
## Real Parameter Phi
##
##           1         2         3         4         5         6
## 1 0.6969697 0.4230769 0.5052879 0.6094017 0.5708181 0.7637733
## 2           0.4230769 0.5052879 0.6094017 0.5708181 0.7637733
## 3                     0.5052879 0.6094017 0.5708181 0.7637733
## 4                               0.6094017 0.5708181 0.7637733
## 5                                         0.5708181 0.7637733
## 6                                                   0.7637733
##
##
## Real Parameter p
##
##           2 3         4         5         6        7
## 1 0.7173912 1 0.9093024 0.9274193 0.9358289 0.763752
## 2           1 0.9093024 0.9274193 0.9358289 0.763752
## 3             0.9093024 0.9274193 0.9358289 0.763752
## 4                       0.9274193 0.9358289 0.763752
## 5                                 0.9358289 0.763752
## 6                                           0.763752
Phi.time.p.time = adjust.parameter.count(Phi.time.p.time, 11)
##
## Number of parameters adjusted from 12 to 11
## Adjusted AICc = 336.434326153846
## Unadjusted AICc = 334.20293

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

edm.results <- model.table(model.list = c("Phi.dot.p.dot",
"Phi.dot.p.time",
"Phi.time.p.dot",
"Phi.time.p.time"),
model.name = FALSE)

### Look at the model-selection table

kable(edm.results, digits = 3)
Phi p model npar AICc DeltaAICc weight Deviance
1 ~1 ~1 Phi.dot.p.dot 2 322.553 0.000 0.960 41.815
3 ~-1 + time ~1 Phi.time.p.dot 7 330.057 7.504 0.023 38.815
2 ~1 ~-1 + time Phi.dot.p.time 7 330.679 8.127 0.017 39.437
4 ~-1 + time ~-1 + time Phi.time.p.time 11 336.434 13.882 0.001 36.401

### Look at output from top model

Although you’ve already seen the output above, you might like to see how to get the real parameter estimates from the top model. Here’s a way that works.

Phi.dot.p.dot$results$beta
##                  estimate        se        lcl       ucl
## Phi:(Intercept) 0.2648275 0.1446688 -0.0187234 0.5483784
## p:(Intercept)   2.4862984 0.5120845  1.4826127 3.4899840
Phi.dot.p.dot$results$real
##                  estimate        se       lcl       ucl fixed    note
## Phi g1 c1 a0 t1 0.5658226 0.0355404 0.4953193 0.6337593
## p g1 c1 a1 t2   0.9231757 0.0363182 0.8149669 0.9704014

### Clean up the files

if you’re running this code, it’s wise to cleanup files you can use the next 2 lines without comments to do that

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