Here, you’ll find code for running the models used in lab this week with the “ClosedSim.inp” data set. These models are closed models, and the code below gives you a chance to see how to run closed models in RMark. It should not be too hard for you to modify the code below to run all the models for Lab 7 on “capture.inp”. If you try that, you’ll see just how convenient it can be to apply code developed for 1 data set to another as long as the model sets for the 2 projects are similar enough.
# Lab 07 for WILD 502 at Montana State University
library(RMark)## This is RMark 2.2.4
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zipcaps <- convert.inp(
"http://www.montana.edu/rotella/documents/502/ClosedSim.inp",
group.df = NULL,
covariates = NULL)
head(caps)## ch freq
## 1 000001 1
## 2 010001 1
## 3 000110 1
## 4 000100 1
## 5 110011 1
## 6 011111 1tail(caps)## ch freq
## 346 100000 1
## 347 000101 1
## 348 110000 1
## 349 000100 1
## 350 100011 1
## 351 001000 1Before we can run all of the models for this lab, we need to add a few covariates that will be used to constrain the estimated capture probabilities in various ways.
# Process data
caps.pr <- process.data(caps, begin.time = 1, model = "Closed")
# Create default design data
caps.ddl <- make.design.data(caps.pr)
# add columns to ddl data for p for M(bh) models
# 1. add 2 categories of time (1st occ.[0] or later[1])
caps.ddl$p$t2 <- 0
caps.ddl$p$t2 = ifelse(caps.ddl$p$Time == 0, 0, 1)
# 2. add 2nd category of time (1st or 2nd occ.[0] or later[1])
caps.ddl$p$t3 <- 0
caps.ddl$p$t3 <- ifelse(caps.ddl$p$Time < 2, 0, 1)
# add columns to ddl data for c for M(bh) models
# Not used here but must be added for things to work
# 1. add 2 categories of time (1st occ.[0] or later[1])
caps.ddl$c$t2 <- 0
caps.ddl$c$t2 <- ifelse(caps.ddl$c$Time == 0, 0, 1)
# 2. add 2nd category of time (1st or 2nd occ.[0] or later[1])
caps.ddl$c$t3 <- 0
caps.ddl$c$t3 <- ifelse(caps.ddl$c$Time < 2, 0, 1)
# Define weekly temperature data
caps.ddl$p$temp <- c(-1.5, 0.2, 0.1, -0.4, 1.5, 1.1)
caps.ddl$c$temp <- c(0.2, 0.1, -0.4, 1.5, 1.1)
# examine the data
caps.ddl## $p
## par.index model.index group time Time c t2 t3 temp
## 1 1 1 1 1 0 0 0 0 -1.5
## 2 2 2 1 2 1 0 1 0 0.2
## 3 3 3 1 3 2 0 1 1 0.1
## 4 4 4 1 4 3 0 1 1 -0.4
## 5 5 5 1 5 4 0 1 1 1.5
## 6 6 6 1 6 5 0 1 1 1.1
##
## $c
## par.index model.index group time Time c t2 t3 temp
## 1 1 7 1 2 0 1 0 0 0.2
## 2 2 8 1 3 1 1 1 0 0.1
## 3 3 9 1 4 2 1 1 1 -0.4
## 4 4 10 1 5 3 1 1 1 1.5
## 5 5 11 1 6 4 1 1 1 1.1
##
## $f0
## par.index model.index group age time Age Time
## 1 1 12 1 0 1 0 0
##
## $pimtypes
## $pimtypes$p
## $pimtypes$p$pim.type
## [1] "all"
##
##
## $pimtypes$c
## $pimtypes$c$pim.type
## [1] "all"
##
##
## $pimtypes$f0
## $pimtypes$f0$pim.type
## [1] "all"Here, we set up the structures for ‘p’ and ‘c’. I use the “share=TRUE” or “share=FALSE” options in each of the structures to indicate whether ‘p’ & ‘c’ should share the same columns of the design matrix or not. Although this is not necessary for all of the structures below, it does add a covariate “c” to the design data with c=0 for rows pertaining to parameter ‘p’ and c=1 for rows pertaining to parameter ’c. This is nice as it gives us the opportunity to build some of the additive structures we’re interested in. We can then use the covariate “c” in formula statements if we want to. But, we don’t have to include that covariate if we don’t want to (e.g., see the p.dot structure).
# function for running set of models for phi and for p
run.caps <- function() {
# Define parameters for p and c
# Note: "share=TRUE" indicates that 'p' & 'c'
# share the same columns in the design matrix
p.dot = list(formula = ~ 1, share = TRUE)
p.time = list(formula = ~ time, share = TRUE)
p.time.behav.add = list(formula = ~ time + c, share = TRUE)
p.dot.behav = list(formula = ~ 1, share = FALSE)
p.bh.2p = list(formula = ~ t2, share = FALSE)
p.bh.3p = list(formula = ~ t2 + t3, share = FALSE)
p.temp = list(formula = ~ temp, share = TRUE)
# Create competing models based on strcutures for 'p' & 'c'
caps.model.list = create.model.list("Closed")
# NOTE: if you do not want to see the output for each model, add
# ', output=FALSE' after 'ddl=caps.ddl' below.
caps.results = mark.wrapper(caps.model.list,
data = caps.pr,
ddl = caps.ddl)
# Return model table and list of models
return(caps.results)
}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. 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.
caps.results <- run.caps()## Using default formula for c## Using default formula for f0##
## p.bh.2p##
## Output summary for Closed model
## Name : p(~t2)c(~1)f0(~1)
##
## Npar : 4
## -2lnL: -912.5
## AICc : -904.4809
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.6393628 0.2839006 -2.1958080 -1.0829176
## p:t2 -0.3702683 0.1707744 -0.7049861 -0.0355505
## c:(Intercept) -0.9058468 0.0675210 -1.0381880 -0.7735056
## f0:(Intercept) 5.6441882 0.4931827 4.6775500 6.6108264
##
##
## Real Parameter p
## 1 2 3 4 5 6
## 0.1625518 0.1181954 0.1181954 0.1181954 0.1181954 0.1181954
##
##
## Real Parameter c
## 2 3 4 5 6
## 0.2878505 0.2878505 0.2878505 0.2878505 0.2878505
##
##
## Real Parameter f0
## 1
## 282.644##
## p.bh.3p##
## Output summary for Closed model
## Name : p(~t2 + t3)c(~1)f0(~1)
##
## Npar : 5
## -2lnL: -912.6323
## AICc : -902.6037
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.5591800 0.3201461 -2.1866664 -0.9316936
## p:t2 -0.3908367 0.1797584 -0.7431632 -0.0385103
## p:t3 0.0769582 0.2114710 -0.3375249 0.4914414
## c:(Intercept) -0.9058467 0.0675210 -1.0381880 -0.7735055
## f0:(Intercept) 5.4879331 0.6102503 4.2918424 6.6840237
##
##
## Real Parameter p
## 1 2 3 4 5 6
## 0.1737643 0.1245515 0.1331882 0.1331882 0.1331882 0.1331882
##
##
## Real Parameter c
## 2 3 4 5 6
## 0.2878505 0.2878505 0.2878505 0.2878505 0.2878505
##
##
## Real Parameter f0
## 1
## 241.757##
## p.dot##
## Output summary for Closed model
## Name : p(~1)c()f0(~1)
##
## Npar : 2
## -2lnL: -893.2344
## AICc : -889.2287
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.029691 0.0590212 -1.145372 -0.9140093
## f0:(Intercept) 4.195525 0.1748259 3.852867 4.5381842
##
##
## Real Parameter p
## 1 2 3 4 5 6
## 0.2631441 0.2631441 0.2631441 0.2631441 0.2631441 0.2631441
##
##
## Real Parameter c
## 2 3 4 5 6
## 0.2631441 0.2631441 0.2631441 0.2631441 0.2631441
##
##
## Real Parameter f0
## 1
## 66.3886##
## p.dot.behav##
## Output summary for Closed model
## Name : p(~1)c(~1)f0(~1)
##
## Npar : 3
## -2lnL: -907.7592
## AICc : -901.7477
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.5332857 0.1822754 -1.890545 -1.1760258
## c:(Intercept) -0.9058468 0.0675210 -1.038188 -0.7735056
## f0:(Intercept) 5.0555181 0.2980502 4.471340 5.6396965
##
##
## Real Parameter p
## 1 2 3 4 5 6
## 0.1775135 0.1775135 0.1775135 0.1775135 0.1775135 0.1775135
##
##
## Real Parameter c
## 2 3 4 5 6
## 0.2878505 0.2878505 0.2878505 0.2878505 0.2878505
##
##
## Real Parameter f0
## 1
## 156.8858##
## p.temp##
## Output summary for Closed model
## Name : p(~temp)c()f0(~1)
##
## Npar : 3
## -2lnL: -894.6979
## AICc : -888.6865
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.0207994 0.0594446 -1.1373109 -0.9042880
## p:temp -0.0559484 0.0462106 -0.1465211 0.0346243
## f0:(Intercept) 4.1942469 0.1749091 3.8514251 4.5370687
##
##
## Real Parameter p
## 1 2 3 4 5 6
## 0.2815317 0.2626987 0.2637837 0.2692522 0.2488562 0.2530629
##
##
## Real Parameter c
## 2 3 4 5 6
## 0.2626987 0.2637837 0.2692522 0.2488562 0.2530629
##
##
## Real Parameter f0
## 1
## 66.30378##
## p.time##
## Output summary for Closed model
## Name : p(~time)c()f0(~1)
##
## Npar : 7
## -2lnL: -965.7636
## AICc : -951.7102
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.1014370 0.1192541 -1.3351749 -0.8676990
## p:time2 -0.1200648 0.1634356 -0.4403985 0.2002690
## p:time3 -0.4874311 0.1736375 -0.8277606 -0.1471016
## p:time4 0.7646570 0.1516801 0.4673641 1.0619500
## p:time5 -0.0657684 0.1622223 -0.3837241 0.2521874
## p:time6 0.2325675 0.1567672 -0.0746962 0.5398311
## f0:(Intercept) 4.1251014 0.1795750 3.7731344 4.4770683
##
##
## Real Parameter p
## 1 2 3 4 5 6
## 0.2494707 0.2276723 0.1695432 0.4165919 0.2373605 0.2954896
##
##
## Real Parameter c
## 2 3 4 5 6
## 0.2276723 0.1695432 0.4165919 0.2373605 0.2954896
##
##
## Real Parameter f0
## 1
## 61.87408##
## p.time.behav.add##
## Output summary for Closed model
## Name : p(~time + c)c()f0(~1)
##
## Npar : 8
## -2lnL: -976.4347
## AICc : -960.366
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.7195022 0.4628964 -2.6267791 -0.8122254
## p:time2 -0.3843736 0.1860089 -0.7489510 -0.0197962
## p:time3 -0.8840206 0.2245540 -1.3241464 -0.4438947
## p:time4 0.1847690 0.2530719 -0.3112518 0.6807898
## p:time5 -0.7517652 0.2840413 -1.3084862 -0.1950442
## p:time6 -0.5226170 0.2975576 -1.1058299 0.0605960
## p:c 1.2839083 0.6630643 -0.0156978 2.5835145
## f0:(Intercept) 5.7897150 0.7920030 4.2373891 7.3420408
##
##
## Real Parameter p
## 1 2 3 4 5 6
## 0.1519353 0.1087207 0.068912 0.1773022 0.0778971 0.0960314
##
##
## Real Parameter c
## 2 3 4 5 6
## 0.3057706 0.2108824 0.4376205 0.2337316 0.2772366
##
##
## Real Parameter f0
## 1
## 326.9198Once the models are run, we can examine the model-selection table. We can also examine model-specific output.
caps.results## model npar AICc DeltaAICc weight Deviance
## 7 p(~time + c)c()f0(~1) 8 -960.3660 0.000000 9.869767e-01 139.3618
## 6 p(~time)c()f0(~1) 7 -951.7102 8.655815 1.302329e-02 150.0329
## 1 p(~t2)c(~1)f0(~1) 4 -904.4809 55.885099 7.227895e-13 203.2965
## 2 p(~t2 + t3)c(~1)f0(~1) 5 -902.6037 57.762262 2.827425e-13 203.1642
## 4 p(~1)c(~1)f0(~1) 3 -901.7478 58.618258 1.842951e-13 208.0373
## 3 p(~1)c()f0(~1) 2 -889.2287 71.137277 0.000000e+00 222.5621
## 5 p(~temp)c()f0(~1) 3 -888.6865 71.679478 0.000000e+00 221.0985# examine model names and find the name of the top model
names(caps.results)## [1] "p.bh.2p" "p.bh.3p" "p.dot"
## [4] "p.dot.behav" "p.temp" "p.time"
## [7] "p.time.behav.add" "model.table"# examine the output from top-ranked model (#7)
caps.results$p.time.behav.add$results## $lnl
## [1] -976.4347
##
## $deviance
## [1] 139.3618
##
## $deviance.df
## [1] 38
##
## $npar
## [1] 8
##
## $n
## [1] 2106
##
## $AICc
## [1] -960.366
##
## $beta
## estimate se lcl ucl
## p:(Intercept) -1.7195022 0.4628964 -2.6267791 -0.8122254
## p:time2 -0.3843736 0.1860089 -0.7489510 -0.0197962
## p:time3 -0.8840206 0.2245540 -1.3241464 -0.4438947
## p:time4 0.1847690 0.2530719 -0.3112518 0.6807898
## p:time5 -0.7517652 0.2840413 -1.3084862 -0.1950442
## p:time6 -0.5226170 0.2975576 -1.1058299 0.0605960
## p:c 1.2839083 0.6630643 -0.0156978 2.5835145
## f0:(Intercept) 5.7897150 0.7920030 4.2373891 7.3420408
##
## $real
## estimate se lcl ucl fixed note
## p g1 t1 0.1519353 0.0596446 0.0674347 0.3074165
## p g1 t2 0.1087207 0.0542606 0.0391130 0.2676947
## p g1 t3 0.0689120 0.0390366 0.0219675 0.1960657
## p g1 t4 0.1773022 0.0968689 0.0553907 0.4419856
## p g1 t5 0.0778971 0.0496739 0.0213167 0.2467877
## p g1 t6 0.0960314 0.0616423 0.0257335 0.2993600
## c g1 t2 0.3057706 0.0366059 0.2390380 0.3817861
## c g1 t3 0.2108824 0.0255312 0.1651604 0.2652404
## c g1 t4 0.4376205 0.0268714 0.3858411 0.4907956
## c g1 t5 0.2337316 0.0214465 0.1943388 0.2783506
## c g1 t6 0.2772366 0.0234123 0.2337564 0.3253701
## f0 g1 a0 t1 326.9198300 258.9214700 83.2674450 1283.5337000
##
## $beta.vcv
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.21427304 0.03234258 0.05272364 0.08135166 0.09164817
## [2,] 0.03234258 0.03459931 0.02659765 0.03266217 0.03556763
## [3,] 0.05272364 0.02659765 0.05042451 0.04270207 0.04705619
## [4,] 0.08135166 0.03266217 0.04270207 0.06404536 0.06133912
## [5,] 0.09164817 0.03556763 0.04705619 0.06133912 0.08067946
## [6,] 0.10070548 0.03756570 0.04999022 0.06554541 0.07285078
## [7,] -0.29515357 -0.06658358 -0.09797679 -0.13922375 -0.15662851
## [8,] -0.35668715 -0.07701017 -0.11285224 -0.16319738 -0.18130478
## [,6] [,7] [,8]
## [1,] 0.10070548 -0.29515357 -0.35668715
## [2,] 0.03756570 -0.06658358 -0.07701017
## [3,] 0.04999022 -0.09797679 -0.11285224
## [4,] 0.06554541 -0.13922375 -0.16319738
## [5,] 0.07285078 -0.15662851 -0.18130478
## [6,] 0.08854054 -0.16995986 -0.19723293
## [7,] -0.16995986 0.43965431 0.51905599
## [8,] -0.19723293 0.51905599 0.62726870
##
## $derived
## $derived$`N Population Size`
## estimate lcl ucl
## 1 677.9198 434.2674 1634.534
##
##
## $derived.vcv
## $derived.vcv$`N Population Size`
## [,1]
## [1,] 67040.33
##
##
## $covariate.values
## NULL
##
## $singular
## NULL
##
## $real.vcv
## NULL# 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)