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’.
library(RMark)## This is RMark 2.2.4
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zipcoho <- 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 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"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())
}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.39506Once 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# 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.3472To 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.9999For 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 NAChapter 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.