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.
library(knitr)
library(RMark)## This is RMark 2.2.4
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zipedm <- 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.9231757Phi.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.9194581Phi.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)
##
## Npar : 7 (unadjusted=5)
## -2lnL: 316.1166
## AICc : 330.6794 (unadjusted=326.41507)
##
## 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 1Phi.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)
##
## Npar : 12 (unadjusted=10)
## -2lnL: 313.0805
## AICc : 338.6887 (unadjusted=334.20293)
##
## 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.763752Phi.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.20293It’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)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 |
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.4899840Phi.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.9704014if 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)