This markdown document runs the POPAN models for data on Coho salmon that you analyzed directly in MARK in lab 8. Here, we’re dealing with 4 types of parameters: ‘Phi’, ‘p’, ‘pent’, and ‘N’.

Bring in the Data

library(RMark)
## This is RMark 2.2.4
##  Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
coho <- convert.inp("http://www.montana.edu/rotella/documents/502/chaseAD.inp",
                 group.df = NULL)
head(coho)
##         ch freq
## 1 11110000   -1
## 2 11010000   -1
## 3 11010000    1
## 4 11001001   -1
## 5 11000000    1
## 6 11000000   -1

Process the data

# Process data
# NOTE: for this data set, we need to provide the time
# intervals because they are not all the same length
coho.pr <- process.data(coho, begin.time = 1, model = "POPAN",
                        time.intervals = c(1.5, 1, 1, 1, 1, 1, 1.5))
#  
# Create default design data
coho.ddl = make.design.data(coho.pr)

# examine the design data 
coho.ddl
## $Phi
##   par.index model.index group age time Age Time
## 1         1           1     1   0    1 0.0  0.0
## 2         2           2     1 1.5  2.5 1.5  1.5
## 3         3           3     1 2.5  3.5 2.5  2.5
## 4         4           4     1 3.5  4.5 3.5  3.5
## 5         5           5     1 4.5  5.5 4.5  4.5
## 6         6           6     1 5.5  6.5 5.5  5.5
## 7         7           7     1 6.5  7.5 6.5  6.5
## 
## $p
##   par.index model.index group age time Age Time
## 1         1           8     1   0    1 0.0  0.0
## 2         2           9     1 1.5  2.5 1.5  1.5
## 3         3          10     1 2.5  3.5 2.5  2.5
## 4         4          11     1 3.5  4.5 3.5  3.5
## 5         5          12     1 4.5  5.5 4.5  4.5
## 6         6          13     1 5.5  6.5 5.5  5.5
## 7         7          14     1 6.5  7.5 6.5  6.5
## 8         8          15     1   8    9 8.0  8.0
## 
## $pent
##   par.index model.index group age time Age Time
## 1         1          16     1 1.5  2.5 1.5  0.0
## 2         2          17     1 2.5  3.5 2.5  1.0
## 3         3          18     1 3.5  4.5 3.5  2.0
## 4         4          19     1 4.5  5.5 4.5  3.0
## 5         5          20     1 5.5  6.5 5.5  4.0
## 6         6          21     1 6.5  7.5 6.5  5.0
## 7         7          22     1   8    9 8.0  6.5
## 
## $N
##   par.index model.index group age time Age Time
## 1         1          23     1   0    1   0    0
## 
## $pimtypes
## $pimtypes$Phi
## $pimtypes$Phi$pim.type
## [1] "all"
## 
## 
## $pimtypes$p
## $pimtypes$p$pim.type
## [1] "all"
## 
## 
## $pimtypes$pent
## $pimtypes$pent$pim.type
## [1] "all"
## 
## 
## $pimtypes$N
## $pimtypes$N$pim.type
## [1] "all"

Build function for creating models

Here, we set up a function that contains the structures for ‘Phi’, ‘p’, ‘pent’ and ‘N’. It will run all combinations of the structures for each parameter type when executed.

run.coho <- function() {
  #
  # Define range of models for each parameter type. Here,
  # we only have a few models as the labs focused more on how
  # how work with the output than on setting up a wide variety
  # of competing models
  
  # Define model for Phi
  Phi.time = list(formula =  ~ time)
  
  #  Define range of models for p
  p.dot = list(formula =  ~ 1)
  p.time = list(formula =  ~ time)
  
  #  Define model for pent (probability of entry)
  pent.time = list(formula =  ~ time)
  
  # Define model for N
  N.dot = list(formula =  ~ 1)
  
  # Run assortment of models
  Phi.time.p.dot.pent.time.N = mark(
    coho.pr,
    coho.ddl,
    model.parameters = list(
      Phi = Phi.time,
      p = p.dot,
      pent = pent.time,
      N = N.dot))
  
  # Note this initial run will penalize for too many parameters because
  #  we are not adjusting the parameter count for confounded parameters
  Phi.time.p.time.pent.time.N = mark(
    coho.pr,
    coho.ddl,
    model.parameters = list(
      Phi = Phi.time,
      p = p.time,
      pent = pent.time,
      N = N.dot))
    
    # Here, we adjust the parameter count to the correct value
    Phi.time.p.time.pent.time.N =
    adjust.parameter.count(Phi.time.p.time.pent.time.N, 21)
  
  # Return model table and list of models
  return(collect.models())
}

Run the Models and Examine the Output

It’s very simple to then run the function and each of the models. As implemented here, the output from each of the models appears in the console where it can be reviewed when the function written above is called. If you don’t want the output from every model to appear automatically, see the note above the “mark.wrapper” command for how to set the option of suppressing the output. One can also examine model-specific output in other ways as shown below.

coho.results <- run.coho()
## 
## Note: only 14 parameters counted of 16 specified parameters
## AICc and parameter count have been adjusted upward
## 
## Output summary for POPAN model
## Name : Phi(~time)p(~1)pent(~time)N(~1) 
## 
## Npar :  16  (unadjusted=14)
## -2lnL:  489.9401
## AICc :  524.5306  (unadjusted=519.92125)
## 
## Beta
##                     estimate          se           lcl          ucl
## Phi:(Intercept)    0.3405772   0.4146701    -0.4721763    1.1533307
## Phi:time2.5        2.1925248   2.7452806    -3.1882253    7.5732749
## Phi:time3.5        1.0292293   0.9068655    -0.7482271    2.8066858
## Phi:time4.5        2.1249729   1.8553331    -1.5114802    5.7614259
## Phi:time5.5       -0.0208952   0.6894135    -1.3721456    1.3303552
## Phi:time6.5       -1.0457876   0.7023157    -2.4223265    0.3307512
## Phi:time7.5        0.1407362   0.7245702    -1.2794215    1.5608938
## p:(Intercept)     -0.7751083   0.1958157    -1.1589070   -0.3913096
## pent:(Intercept)  -2.0593699   1.2107900    -4.4325184    0.3137786
## pent:time3.5       2.0009395   1.2211165    -0.3924489    4.3943279
## pent:time4.5       0.9016539   1.2798597    -1.6068711    3.4101789
## pent:time5.5     -12.5341290 624.3467100 -1236.2537000 1211.1854000
## pent:time6.5      -0.4771468   1.8307879    -4.0654912    3.1111976
## pent:time7.5       1.0721899   1.1726052    -1.2261162    3.3704961
## pent:time9       -13.6127300 584.3804600 -1158.9985000 1131.7730000
## N:(Intercept)      4.9487594   0.2062252     4.5445581    5.3529607
## 
## 
## Real Parameter Phi
##          1       2.5       3.5       4.5       5.5      6.5      7.5
##  0.5843307 0.9264301 0.7973489 0.9216912 0.5792467 0.330658 0.618058
## 
## 
## Real Parameter p
##          1       2.5       3.5       4.5       5.5       6.5       7.5
##  0.3153751 0.3153751 0.3153751 0.3153751 0.3153751 0.3153751 0.3153751
##          9
##  0.3153751
## 
## 
## Real Parameter pent
##        2.5       3.5       4.5          5.5       6.5       7.5
##  0.0449579 0.3325087 0.1107616 1.619209e-07 0.0278987 0.1313567
##             9
##  5.506461e-08
## 
## 
## Real Parameter N
##         1
##  331.9999
## 
## Note: only 16 parameters counted of 23 specified parameters
## 
## AICc and parameter count have been adjusted upward
## 
## Output summary for POPAN model
## Name : Phi(~time)p(~time)pent(~time)N(~1) 
## 
## Npar :  23  (unadjusted=16)
## -2lnL:  486.8046
## AICc :  538.243  (unadjusted=521.39506)
## 
## Beta
##                     estimate          se          lcl        ucl
## Phi:(Intercept)    0.2889665   0.3991072   -0.4932837   1.071217
## Phi:time2.5       69.7831780   0.0000000   69.7831780  69.783178
## Phi:time3.5        0.6512625   0.7679799   -0.8539780   2.156503
## Phi:time4.5       41.7986950   0.0000000   41.7986950  41.798695
## Phi:time5.5        0.9926691   2.8083544   -4.5117055   6.497044
## Phi:time6.5       -1.2225501   1.2208442   -3.6154048   1.170305
## Phi:time7.5       20.2824040   0.0000000   20.2824040  20.282404
## p:(Intercept)     10.2386410 125.7270800 -236.1864400 256.663730
## p:time2.5        -10.7490360 125.7283400 -257.1766000 235.678520
## p:time3.5        -11.4209800 125.7272400 -257.8463700 235.004410
## p:time4.5        -10.7720870 125.7272900 -257.1975900 235.653410
## p:time5.5        -11.0492700 125.7272800 -257.4747500 235.376210
## p:time6.5        -11.5355790 125.7294100 -257.9652300 234.894070
## p:time7.5        -11.2582550 125.7288900 -257.6868900 235.170380
## p:time9          -12.1405600 125.7288300 -258.5690600 234.287940
## pent:(Intercept)   0.1421701   0.4732481   -0.7853961   1.069736
## pent:time3.5       1.3465751   0.5541140    0.2605116   2.432639
## pent:time4.5     -24.0561900   0.0000000  -24.0561900 -24.056190
## pent:time5.5     -36.1883660   0.0000000  -36.1883660 -36.188366
## pent:time6.5      -0.7727439   1.3487229   -3.4162408   1.870753
## pent:time7.5       0.1271767   0.7335699   -1.3106204   1.564974
## pent:time9       -61.0376280   0.0000000  -61.0376280 -61.037628
## N:(Intercept)      4.7938177   0.3177971    4.1709353   5.416700
## 
## 
## Real Parameter Phi
##          1 2.5       3.5 4.5       5.5       6.5 7.5
##  0.5717431   1 0.7191459   1 0.7827281 0.2821982   1
## 
## 
## Real Parameter p
##          1       2.5       3.5       4.5       5.5       6.5       7.5
##  0.9999642 0.3751008 0.2346318 0.3697136 0.3077565 0.2146807 0.2651026
##          9
##  0.1298914
## 
## 
## Real Parameter pent
##        2.5       3.5          4.5          5.5       6.5      7.5
##  0.1368163 0.5259541 4.882778e-12 2.628633e-17 0.0631741 0.155371
##             9
##  4.244559e-28
## 
## 
## Real Parameter N
##         1
##  311.7615
## 
## Number of parameters adjusted from 23 to 21
## Adjusted AICc = 533.311897073171
## Unadjusted AICc = 521.39506

Examine Model-Selection Table

Once the models are run, we can examine the model-selection table. We can also examine model-specific output.

options(width = 150)
coho.results
##                                model npar     AICc DeltaAICc     weight  Deviance
## 1    Phi(~time)p(~1)pent(~time)N(~1)   16 524.5306  0.000000 0.98775903 -517.7350
## 2 Phi(~time)p(~time)pent(~time)N(~1)   21 533.3119  8.781301 0.01224097 -520.8706

Look at output from top model

# View the model names
names(coho.results)
## [1] "Phi.time.p.dot.pent.time.N"  "Phi.time.p.time.pent.time.N" "model.table"
# examine the output from top-ranked model (#1) and
# store top model in object with short name 
top <- coho.results$Phi.time.p.dot.pent.time.N

# Check which link functions were used for the parameters
# Notice that RMark used parameter-specific link functions that 
# make sense here, i.e., the logit link for Phi & p, the mlogit for pent,
# and a log link for N.
top$links
##  [1] "Logit"     "Logit"     "Logit"     "Logit"     "Logit"     "Logit"     "Logit"     "Logit"     "Log"       "mlogit(1)" "mlogit(1)" "mlogit(1)"
## [13] "mlogit(1)" "mlogit(1)" "mlogit(1)" "mlogit(1)"
# look at summary of top model's output
summary(top, showall = FALSE)
## Output summary for POPAN model
## Name : Phi(~time)p(~1)pent(~time)N(~1) 
## 
## Npar :  16  (unadjusted=14)
## -2lnL:  489.9401
## AICc :  524.5306  (unadjusted=519.92125)
## 
## Beta
##                     estimate          se           lcl          ucl
## Phi:(Intercept)    0.3405772   0.4146701    -0.4721763    1.1533307
## Phi:time2.5        2.1925248   2.7452806    -3.1882253    7.5732749
## Phi:time3.5        1.0292293   0.9068655    -0.7482271    2.8066858
## Phi:time4.5        2.1249729   1.8553331    -1.5114802    5.7614259
## Phi:time5.5       -0.0208952   0.6894135    -1.3721456    1.3303552
## Phi:time6.5       -1.0457876   0.7023157    -2.4223265    0.3307512
## Phi:time7.5        0.1407362   0.7245702    -1.2794215    1.5608938
## p:(Intercept)     -0.7751083   0.1958157    -1.1589070   -0.3913096
## pent:(Intercept)  -2.0593699   1.2107900    -4.4325184    0.3137786
## pent:time3.5       2.0009395   1.2211165    -0.3924489    4.3943279
## pent:time4.5       0.9016539   1.2798597    -1.6068711    3.4101789
## pent:time5.5     -12.5341290 624.3467100 -1236.2537000 1211.1854000
## pent:time6.5      -0.4771468   1.8307879    -4.0654912    3.1111976
## pent:time7.5       1.0721899   1.1726052    -1.2261162    3.3704961
## pent:time9       -13.6127300 584.3804600 -1158.9985000 1131.7730000
## N:(Intercept)      4.9487594   0.2062252     4.5445581    5.3529607
## 
## 
## Real Parameter Phi
##          1       2.5       3.5       4.5       5.5      6.5      7.5
##  0.5843307 0.9264301 0.7973489 0.9216912 0.5792467 0.330658 0.618058
## 
## 
## Real Parameter p
##          1       2.5       3.5       4.5       5.5       6.5       7.5         9
##  0.3153751 0.3153751 0.3153751 0.3153751 0.3153751 0.3153751 0.3153751 0.3153751
## 
## 
## Real Parameter pent
##        2.5       3.5       4.5          5.5       6.5       7.5            9
##  0.0449579 0.3325087 0.1107616 1.619209e-07 0.0278987 0.1313567 5.506461e-08
## 
## 
## Real Parameter N
##         1
##  331.9999
# store and examine estimates of Phi, p, pent, and N

# First examine estimates of the real parameters
# to see how results are stored. We drop 'fixed' & 'note' columns, so
# we can use the 'round' command
round(top$results$real[ , 1:4], 5)
##                    estimate       se       lcl       ucl
## Phi g1 a0 t1        0.58433  0.10072   0.38410   0.76012
## Phi g1 a1.5 t2.5    0.92643  0.17596   0.07400   0.99950
## Phi g1 a2.5 t3.5    0.79735  0.13284   0.43989   0.95172
## Phi g1 a3.5 t4.5    0.92169  0.13129   0.24981   0.99760
## Phi g1 a4.5 t5.5    0.57925  0.13726   0.31343   0.80589
## Phi g1 a5.5 t6.5    0.33066  0.12581   0.13952   0.60082
## Phi g1 a6.5 t7.5    0.61806  0.14054   0.33501   0.83865
## p g1 a0 t1          0.31538  0.04228   0.23887   0.40340
## N g1 a0 t1        331.99993 29.07773 285.51508 401.34718
## pent g1 a1.5 t2.5   0.04496  0.05065   0.00464   0.32216
## pent g1 a2.5 t3.5   0.33251  0.07181   0.20899   0.48434
## pent g1 a3.5 t4.5   0.11076  0.07183   0.02896   0.34217
## pent g1 a4.5 t5.5   0.00000  0.00010   0.00000   1.00000
## pent g1 a5.5 t6.5   0.02790  0.04113   0.00147   0.35932
## pent g1 a6.5 t7.5   0.13136  0.03557   0.07585   0.21790
## pent g1 a8 t9       0.00000  0.00003   0.00000   1.00000
# for Phi in top model there are 7 estimates
top.Phi <- top$results$real[1:7, ]
top.Phi
##                   estimate        se       lcl       ucl fixed    note
## Phi g1 a0 t1     0.5843307 0.1007185 0.3841013 0.7601188              
## Phi g1 a1.5 t2.5 0.9264301 0.1759566 0.0739955 0.9994963              
## Phi g1 a2.5 t3.5 0.7973489 0.1328446 0.4398944 0.9517177              
## Phi g1 a3.5 t4.5 0.9216912 0.1312869 0.2498095 0.9976020              
## Phi g1 a4.5 t5.5 0.5792467 0.1372581 0.3134262 0.8058886              
## Phi g1 a5.5 t6.5 0.3306580 0.1258054 0.1395154 0.6008240              
## Phi g1 a6.5 t7.5 0.6180580 0.1405429 0.3350128 0.8386511
# for p in top model there is 1 estimate
# and it's in the 8th row of output
top.p <- top$results$real[8, ]
top.p
##             estimate        se       lcl       ucl fixed    note
## p g1 a0 t1 0.3153751 0.0422793 0.2388659 0.4034021
# Store and examine the estimates of pent
top.pent <- top$results$real[10:16, ]
round(top.pent[ ,1:4], 5)
##                   estimate      se     lcl     ucl
## pent g1 a1.5 t2.5  0.04496 0.05065 0.00464 0.32216
## pent g1 a2.5 t3.5  0.33251 0.07181 0.20899 0.48434
## pent g1 a3.5 t4.5  0.11076 0.07183 0.02896 0.34217
## pent g1 a4.5 t5.5  0.00000 0.00010 0.00000 1.00000
## pent g1 a5.5 t6.5  0.02790 0.04113 0.00147 0.35932
## pent g1 a6.5 t7.5  0.13136 0.03557 0.07585 0.21790
## pent g1 a8 t9      0.00000 0.00003 0.00000 1.00000
# Store estimate for N
top.N <- top$results$real[9, ]
top.N
##            estimate       se      lcl      ucl fixed    note
## N g1 a0 t1 331.9999 29.07773 285.5151 401.3472

Calculate the estimated number of births

To calculate the estimates of B[i], you need to multiply N by b[i]. Remember, to get b[0], you subtract the sum of all the b[i] for i = 1, 2, …, 7 from 1.

# calculate the proportion of the population present at the start. 
# To do this we need the b[i] for each occasion
b <- top.pent$estimate
b0 <- 1 - sum(top.pent$estimate)
#
# Examine the b[i] including b[0]
round(c(b0, b), 4)
## [1] 0.3525 0.0450 0.3325 0.1108 0.0000 0.0279 0.1314 0.0000
# Next, calculate the net births by occasion
B.occ <- b*top.N$estimate 
round(B.occ, 4)
## [1]  14.9260 110.3929  36.7728   0.0001   9.2624  43.6104   0.0000
# B[0] is simply the proportion of the superpopulation that was 
# present at the start of the study (b[0]*N-hat)
(B0 <- b0*top.N$estimate)
## [1] 117.0353
# Sum up all the B[i] including B0
sum(B0, B.occ)
## [1] 331.9999

Work out N[i] values

For this, we have to deal with the losses on capture and the fact that not all time intervals are the same between the various occasions. The details are provided on page 12-20 of the C&W chapter.

# review values for B[i]
round(c(B0, B.occ), 4)
## [1] 117.0353  14.9260 110.3929  36.7728   0.0001   9.2624  43.6104   0.0000
# calculate and store values of N[i] for each occasion iteratively 
# set up vector for containing N by occ
N.occ <- vector("numeric", 8)
# Store starting number
N.occ[1] <- B0

# Store time intervals for use in calculations
time.intervals <- c(1.5, 1, 1, 1, 1, 1, 1.5)
# Store the losses on capture (see Table 12.6 on page 12-15 of C&W)
losses <- c(0, 1, 11, 2, 8, 8, 6, 10)

# Calculate the N.occ values
for (i in 2:8) {
  N.occ[i] = (N.occ[i - 1] - losses[i - 1]) * 
    (top.Phi$estimate[i - 1]^time.intervals[i - 1]) +
    B.occ[i - 1]
}

# examine estimated numbers 
round(cbind(N.occ = N.occ, B.occ = c(B.occ, NA)), 4)
##         N.occ    B.occ
## [1,] 117.0353  14.9260
## [2,]  67.2024 110.3929
## [3,] 171.7247  36.7728
## [4,] 164.9265   0.0001
## [5,] 150.1680   9.2624
## [6,]  91.6127  43.6104
## [7,]  71.2576   0.0000
## [8,]  31.7085       NA

Further work that can be done

Chapter 12 goes through more modeling that can be done including additional models of the adults and combined modeling for adults and jacks. Also, you could work with the delta method to obtain estimates of the SE’s for the quantities we derived above for B.occ and N.occ.