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.

Bring in the data

library(RMark)
## This is RMark 2.2.4
##  Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
rd.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 3

Set up models to run

For 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())
}

Run the models and examine the output

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

Examine model-selection table

options(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.736
names(robust.results)
## [1] "mod.markov"  "mod.random"  "mod.ZeroTE"  "model.table"

Examine output from each model

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.049
round(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.049
round(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.049
round(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.049
round(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.000
round(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