This markdown file runs 3 robust-design models of interest in RMark. The models are described in section 15.6.1 of Cooch & White’ “Gentle Introduction to MARK”. They represent models with either (i) no temporary emigration out of the study population, (ii) random temporary emigration, or (iii) Markovian temporary emigration.
library(RMark)## This is RMark 2.2.4
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.ziprd.inp <- convert.inp("http://www.montana.edu/rotella/documents/502/rd_simple1.inp")
# set up time intervals for 5 primary occasions each with 3 secondary occasions
# these values are the interval lengths between occasions
# e.g, 0 time between occ 1 & 2, 0 between 2 & 3, 1 year between 3 & 4
time.intervals <- c(0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0)
rd <- process.data(data = rd.inp,
model = "Robust",
time.intervals = time.intervals)
names(rd)## [1] "data" "model" "mixtures"
## [4] "freq" "nocc" "nocc.secondary"
## [7] "time.intervals" "begin.time" "age.unit"
## [10] "initial.ages" "group.covariates" "nstrata"
## [13] "strata.labels" "counts" "reverse"
## [16] "areas"# take a look at the data
head(rd$data)## ch freq
## 1 111000000000000 146
## 2 010101000000000 11
## 3 110110000000000 8
## 4 101001110001001 1
## 5 100010000000000 10
## 6 110101000000000 8# check the number of primary occasions
rd$nocc## [1] 5# check the number of secondary occasions
rd$nocc.secondary## [1] 3 3 3 3 3For this week’s lab, we’ll set up specific models in the RMark function. This is done because here we don’t just want all possible combinations of the parameterizations for S, p, and gamma. Rather, we want only specific combinations and so need to dictate which ones we want one model at a time.
run.robust <- function() {
# In all 3 models, S varies by year, p varies by primary occasion but not secondary,
# So we set up structures for S and p before building competing models of TE
# Apparent survival varies by year
S.time = list(formula = ~ time)
# p varies by primary session but not among secondaries within primary
# session = varies by primary, time=varies by secondary
# session:time = varies by both
# p=c due to use of "share=TRUE"
p.session = list(formula = ~ session, share = TRUE)
# set up Gamma structure for Markovian time-varying temporary emigration
# With time-varying S,p, and TE, we need to constrain some gamma's
# A conventional way to do it is to constrain
# gamma"_(k) = gamm"_(k-1) & gamma'_(k) = gamm'_(k-1)
# so that all S and remaining gamma's are identifiable.
GammaDoublePrime.markov = list(formula = ~ time)
GammaPrime.time = list(formula = ~ time)
# Model 1: TE is Markovian & time-varying, S(t), p(varies by primary)
# use design parameters to constrain some of the gammas
# specifically, this uses time-binning to constrain gamma" and gamma'
# bins need to be same for each type of gamma to make this happen
mod.markov = mark(data = rd,
design.parameters = list(
GammaDoublePrime = list(time.bins = c(1, 2, 3, 5)),
GammaPrime = list(time.bins = c(2, 3, 5))),
right = FALSE,
model.parameters = list(
S = S.time,
p = p.session,
GammaDoublePrime = GammaDoublePrime.markov,
GammaPrime = GammaPrime.time))
# Model 2: random, time-varying TE, S(t), p(varies by primary)
# set up Gamma structure of random temporary emigration
# by using "share=TRUE"
# TE rate varies by year
# Even though Gamma' is shared with Gamma" such that don't
# provide model.parameters explicitly for GammaPrime,
# still need to provide time bins for GammaPrime and need to
# make them the same as those used for GammaDoublePrime
GammaDoublePrime.random = list(formula = ~ time, share = TRUE)
mod.random = mark(data = rd,
design.parameters = list(
GammaDoublePrime = list(time.bins = c(1, 2, 3, 5)),
GammaPrime = list(time.bins = c(2, 3, 5))),
right = FALSE,
model.parameters = list(
S = S.time,
p = p.session,
GammaDoublePrime = GammaDoublePrime.random))
# Model 3 = no TE
# fix Gamma's = 0
GammaDoublePrime.zeroTE = list(formula = ~ 1,
fixed = 0)
GammaPrime.zero = list(formula = ~ 1,
fixed = 0)
mod.ZeroTE = mark(data = rd,
model.parameters = list(
S = S.time,
p = p.session,
GammaDoublePrime = GammaDoublePrime.zeroTE,
GammaPrime = GammaPrime.zero))
# Return model table and list of models
return(collect.models())
}It’s very simple to then use the function to run 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.
robust.results=run.robust()##
## Output summary for Robust model
## Name : S(~time)Gamma''(~time)Gamma'(~time)p(~session)c()f0(~session)
##
## Npar : 19
## -2lnL: -50048.47
## AICc : -50010.41
##
## Beta
## estimate se lcl ucl
## S:(Intercept) 0.8048548 0.0503436 0.7061813 0.9035284
## S:time2 0.5010670 0.1158688 0.2739641 0.7281700
## S:time3 1.3946568 0.2809366 0.8440211 1.9452925
## S:time4 1.6863700 0.4347458 0.8342683 2.5384718
## GammaDoublePrime:(Intercept) -1.4190880 0.0839355 -1.5836016 -1.2545743
## GammaDoublePrime:time[2,3) 0.6462300 0.1232970 0.4045679 0.8878921
## GammaDoublePrime:time[3,5] 0.5442858 0.1356674 0.2783777 0.8101940
## GammaPrime:(Intercept) -1.3978626 0.2365017 -1.8614060 -0.9343192
## GammaPrime:time[3,5] 0.8668237 0.3131704 0.2530097 1.4806378
## p:(Intercept) 0.0132579 0.0278658 -0.0413591 0.0678749
## p:session2 0.4106212 0.0444331 0.3235323 0.4977102
## p:session3 0.3716419 0.0503407 0.2729742 0.4703097
## p:session4 -0.0393964 0.0560363 -0.1492276 0.0704347
## p:session5 0.8660128 0.0540865 0.7600032 0.9720223
## f0:(Intercept) 5.9004951 0.0736910 5.7560607 6.0449294
## f0:session2 -1.2673029 0.1427056 -1.5470058 -0.9875999
## f0:session3 -1.5740309 0.1618332 -1.8912240 -1.2568377
## f0:session4 -1.0348216 0.1456652 -1.3203254 -0.7493177
## f0:session5 -2.7759624 0.2492861 -3.2645632 -2.2873615
##
##
## Real Parameter S
##
## 1 2 3 4
## 1 0.691012 0.7868299 0.9002056 0.9235244
## 2 0.7868299 0.9002056 0.9235244
## 3 0.9002056 0.9235244
## 4 0.9235244
##
##
## Real Parameter GammaDoublePrime
##
## 1 2 3 4
## 1 0.1948046 0.3158612 0.2942561 0.2942561
## 2 0.3158612 0.2942561 0.2942561
## 3 0.2942561 0.2942561
## 4 0.2942561
##
##
## Real Parameter GammaPrime
##
## 2 3 4
## 1 0.1981555 0.3702746 0.3702746
## 2 0.3702746 0.3702746
## 3 0.3702746
##
##
## Real Parameter p
## Session:1
## 1 2 3
## 0.5033144 0.5033144 0.5033144
##
## Session:2
## 1 2 3
## 0.6044111 0.6044111 0.6044111
##
## Session:3
## 1 2 3
## 0.5950543 0.5950543 0.5950543
##
## Session:4
## 1 2 3
## 0.4934657 0.4934657 0.4934657
##
## Session:5
## 1 2 3
## 0.7066711 0.7066711 0.7066711
##
##
## Real Parameter c
## Session:1
## 2 3
## 0.5033144 0.5033144
##
## Session:2
## 2 3
## 0.6044111 0.6044111
##
## Session:3
## 2 3
## 0.5950543 0.5950543
##
## Session:4
## 2 3
## 0.4934657 0.4934657
##
## Session:5
## 2 3
## 0.7066711 0.7066711
##
##
## Real Parameter f0
## Session:1
## 0
## 365.2182
##
## Session:2
## 0
## 102.8418
##
## Session:3
## 0
## 75.67623
##
## Session:4
## 0
## 129.7583
##
## Session:5
## 0
## 22.74926
##
## Output summary for Robust model
## Name : S(~time)Gamma''(~time)Gamma'()p(~session)c()f0(~session)
##
## Npar : 17
## -2lnL: -50038.02
## AICc : -50003.97
##
## Beta
## estimate se lcl ucl
## S:(Intercept) 0.8258660 0.0507790 0.7263392 0.9253928
## S:time2 0.4063298 0.1023203 0.2057820 0.6068775
## S:time3 1.1517761 0.1799759 0.7990233 1.5045289
## S:time4 1.6270581 0.3972188 0.8485093 2.4056070
## GammaDoublePrime:(Intercept) -1.3864256 0.0824090 -1.5479473 -1.2249038
## GammaDoublePrime:time[2,3) 0.4805955 0.1107678 0.2634906 0.6977005
## GammaDoublePrime:time[3,5] 0.5014666 0.1278346 0.2509107 0.7520224
## p:(Intercept) 0.0132591 0.0278658 -0.0413579 0.0678760
## p:session2 0.4106220 0.0444331 0.3235331 0.4977110
## p:session3 0.3716336 0.0503417 0.2729639 0.4703033
## p:session4 -0.0394020 0.0560366 -0.1492338 0.0704297
## p:session5 0.8660091 0.0540865 0.7599994 0.9720187
## f0:(Intercept) 5.9004971 0.0736907 5.7560632 6.0449309
## f0:session2 -1.2673127 0.1427063 -1.5470170 -0.9876085
## f0:session3 -1.5740297 0.1618348 -1.8912259 -1.2568334
## f0:session4 -1.0348265 0.1456659 -1.3203317 -0.7493214
## f0:session5 -2.7759688 0.2492881 -3.2645735 -2.2873641
##
##
## Real Parameter S
##
## 1 2 3 4
## 1 0.6954801 0.7742027 0.8784296 0.920775
## 2 0.7742027 0.8784296 0.920775
## 3 0.8784296 0.920775
## 4 0.920775
##
##
## Real Parameter GammaDoublePrime
##
## 1 2 3 4
## 1 0.199979 0.2878539 0.2921512 0.2921512
## 2 0.2878539 0.2921512 0.2921512
## 3 0.2921512 0.2921512
## 4 0.2921512
##
##
## Real Parameter GammaPrime
##
## 2 3 4
## 1 0.2878539 0.2921512 0.2921512
## 2 0.2921512 0.2921512
## 3 0.2921512
##
##
## Real Parameter p
## Session:1
## 1 2 3
## 0.5033147 0.5033147 0.5033147
##
## Session:2
## 1 2 3
## 0.6044116 0.6044116 0.6044116
##
## Session:3
## 1 2 3
## 0.5950526 0.5950526 0.5950526
##
## Session:4
## 1 2 3
## 0.4934646 0.4934646 0.4934646
##
## Session:5
## 1 2 3
## 0.7066705 0.7066705 0.7066705
##
##
## Real Parameter c
## Session:1
## 2 3
## 0.5033147 0.5033147
##
## Session:2
## 2 3
## 0.6044116 0.6044116
##
## Session:3
## 2 3
## 0.5950526 0.5950526
##
## Session:4
## 2 3
## 0.4934646 0.4934646
##
## Session:5
## 2 3
## 0.7066705 0.7066705
##
##
## Real Parameter f0
## Session:1
## 0
## 365.219
##
## Session:2
## 0
## 102.841
##
## Session:3
## 0
## 75.67648
##
## Session:4
## 0
## 129.7579
##
## Session:5
## 0
## 22.74916
##
## Output summary for Robust model
## Name : S(~time)Gamma''(~1)Gamma'(~1)p(~session)c()f0(~session)
##
## Npar : 14
## -2lnL: -49178.71
## AICc : -49150.68
##
## Beta
## estimate se lcl ucl
## S:(Intercept) 0.6870743 0.0436630 0.6014948 0.7726538
## S:time2 0.5007340 0.0785131 0.3468484 0.6546196
## S:time3 1.1313854 0.1177992 0.9004989 1.3622719
## S:time4 0.1535542 0.0833413 -0.0097947 0.3169031
## p:(Intercept) 0.0132540 0.0278661 -0.0413634 0.0678715
## p:session2 0.0379510 0.0395811 -0.0396279 0.1155299
## p:session3 -0.2238613 0.0420872 -0.3063522 -0.1413704
## p:session4 -0.5240347 0.0458981 -0.6139950 -0.4340745
## p:session5 0.8660159 0.0540868 0.7600057 0.9720261
## f0:(Intercept) 5.9005050 0.0736912 5.7560703 6.0449397
## f0:session2 -0.5810066 0.1156812 -0.8077418 -0.3542713
## f0:session3 -0.5219328 0.1166189 -0.7505059 -0.2933598
## f0:session4 -0.2616709 0.1144630 -0.4860184 -0.0373235
## f0:session5 -2.7759713 0.2492870 -3.2645738 -2.2873688
##
##
## Real Parameter S
##
## 1 2 3 4
## 1 0.6653158 0.7663488 0.8603812 0.6985976
## 2 0.7663488 0.8603812 0.6985976
## 3 0.8603812 0.6985976
## 4 0.6985976
##
##
## Real Parameter GammaDoublePrime
##
## 1 2 3 4
## 1
## 2
## 3
## 4
##
##
## Real Parameter GammaPrime
##
## 2 3 4
## 1
## 2
## 3
##
##
## Real Parameter p
## Session:1
## 1 2 3
## 0.5033135 0.5033135 0.5033135
##
## Session:2
## 1 2 3
## 0.5127985 0.5127985 0.5127985
##
## Session:3
## 1 2 3
## 0.4475419 0.4475419 0.4475419
##
## Session:4
## 1 2 3
## 0.3750105 0.3750105 0.3750105
##
## Session:5
## 1 2 3
## 0.7066709 0.7066709 0.7066709
##
##
## Real Parameter c
## Session:1
## 2 3
## 0.5033135 0.5033135
##
## Session:2
## 2 3
## 0.5127985 0.5127985
##
## Session:3
## 2 3
## 0.4475419 0.4475419
##
## Session:4
## 2 3
## 0.3750105 0.3750105
##
## Session:5
## 2 3
## 0.7066709 0.7066709
##
##
## Real Parameter f0
## Session:1
## 0
## 365.2219
##
## Session:2
## 0
## 204.2814
##
## Session:3
## 0
## 216.7126
##
## Session:4
## 0
## 281.1347
##
## Session:5
## 0
## 22.74928options(width = 90)
robust.results## model npar AICc DeltaAICc
## 1 S(~time)Gamma''(~time)Gamma'(~time)p(~session)c()f0(~session) 19 -50010.41 0.000000
## 2 S(~time)Gamma''(~time)Gamma'()p(~session)c()f0(~session) 17 -50003.97 6.441603
## 3 S(~time)Gamma''(~1)Gamma'(~1)p(~session)c()f0(~session) 14 -49150.68 859.730822
## weight Deviance
## 1 0.96160961 5205.979
## 2 0.03839039 5216.432
## 3 0.00000000 6075.736names(robust.results)## [1] "mod.markov" "mod.random" "mod.ZeroTE" "model.table"round(robust.results$mod.markov$results$real[, 1:4], 3)## estimate se lcl ucl
## S g1 c1 a0 t1 0.691 0.011 0.670 0.712
## S g1 c1 a1 t2 0.787 0.016 0.755 0.816
## S g1 c1 a2 t3 0.900 0.025 0.840 0.939
## S g1 c1 a3 t4 0.924 0.030 0.838 0.966
## Gamma'' g1 c1 a0 t1 0.195 0.013 0.170 0.222
## Gamma'' g1 c1 a1 t2 0.316 0.018 0.281 0.352
## Gamma'' g1 c1 a2 t3 0.294 0.022 0.253 0.339
## Gamma' g1 c1 a1 t2 0.198 0.038 0.135 0.282
## Gamma' g1 c1 a2 t3 0.370 0.053 0.273 0.479
## p g1 s1 t1 0.503 0.007 0.490 0.517
## p g1 s2 t1 0.604 0.008 0.588 0.621
## p g1 s3 t1 0.595 0.010 0.575 0.615
## p g1 s4 t1 0.493 0.012 0.470 0.517
## p g1 s5 t1 0.707 0.010 0.687 0.725
## f0 g1 a0 s1 t0 365.218 26.913 316.162 421.885
## f0 g1 a0 s2 t0 102.842 12.568 81.008 130.560
## f0 g1 a0 s3 t0 75.676 10.904 57.140 100.225
## f0 g1 a0 s4 t0 129.758 16.304 101.531 165.834
## f0 g1 a0 s5 t0 22.749 5.418 14.356 36.049round(robust.results$mod.markov$results$derived$`N Population Size`, 3)## estimate lcl ucl
## 1 2984.218 2935.162 3040.885
## 2 1668.842 1647.008 1696.560
## 3 1146.676 1128.140 1171.225
## 4 1001.758 973.531 1037.834
## 5 920.749 912.356 934.049round(robust.results$mod.random$results$real[, 1:4], 3)## estimate se lcl ucl
## S g1 c1 a0 t1 0.695 0.011 0.674 0.716
## S g1 c1 a1 t2 0.774 0.014 0.747 0.800
## S g1 c1 a2 t3 0.878 0.018 0.837 0.910
## S g1 c1 a3 t4 0.921 0.029 0.843 0.962
## Gamma'' g1 c1 a0 t1 0.200 0.013 0.175 0.227
## Gamma'' g1 c1 a1 t2 0.288 0.015 0.259 0.318
## Gamma'' g1 c1 a2 t3 0.292 0.020 0.254 0.333
## p g1 s1 t1 0.503 0.007 0.490 0.517
## p g1 s2 t1 0.604 0.008 0.588 0.621
## p g1 s3 t1 0.595 0.010 0.575 0.615
## p g1 s4 t1 0.493 0.012 0.470 0.517
## p g1 s5 t1 0.707 0.010 0.687 0.725
## f0 g1 a0 s1 t0 365.219 26.913 316.163 421.886
## f0 g1 a0 s2 t0 102.841 12.568 81.008 130.559
## f0 g1 a0 s3 t0 75.676 10.904 57.140 100.225
## f0 g1 a0 s4 t0 129.758 16.304 101.530 165.833
## f0 g1 a0 s5 t0 22.749 5.418 14.356 36.049round(robust.results$mod.random$results$derived$`N Population Size`, 3)## estimate lcl ucl
## 1 2984.219 2935.163 3040.886
## 2 1668.841 1647.008 1696.559
## 3 1146.676 1128.140 1171.225
## 4 1001.758 973.530 1037.833
## 5 920.749 912.356 934.049round(robust.results$mod.ZeroTE$results$real[, 1:4], 3)## estimate se lcl ucl
## S g1 c1 a0 t1 0.665 0.010 0.646 0.684
## S g1 c1 a1 t2 0.766 0.011 0.744 0.787
## S g1 c1 a2 t3 0.860 0.013 0.833 0.884
## S g1 c1 a3 t4 0.699 0.015 0.669 0.727
## p g1 s1 t1 0.503 0.007 0.490 0.517
## p g1 s2 t1 0.513 0.007 0.499 0.527
## p g1 s3 t1 0.448 0.008 0.432 0.463
## p g1 s4 t1 0.375 0.009 0.358 0.392
## p g1 s5 t1 0.707 0.010 0.687 0.725
## f0 g1 a0 s1 t0 365.222 26.914 316.165 421.890
## f0 g1 a0 s2 t0 204.281 18.216 171.583 243.211
## f0 g1 a0 s3 t0 216.713 19.588 181.595 258.622
## f0 g1 a0 s4 t0 281.135 24.624 236.866 333.678
## f0 g1 a0 s5 t0 22.749 5.418 14.356 36.049
## Gamma'' g1 c1 a0 t1 0.000 0.000 0.000 0.000round(robust.results$mod.ZeroTE$results$derived$`N Population Size`, 3)## estimate lcl ucl
## 1 2984.222 2935.165 3040.890
## 2 1770.281 1737.583 1809.211
## 3 1287.713 1252.595 1329.622
## 4 1153.135 1108.866 1205.678
## 5 920.749 912.356 934.049