7 Mortality Components
Cobb 1917.
7.1 Exploitation rate
A tagging study can provide a real-time estimate of the exploitation rate, or the fraction of the population being harvested. Similar to the tag-return method of estimating survival (Section 6.3), the only field work required is the initial tagging effort. Here there is a single release of fish with external tags at the start of the study, and the number of returned tags is recorded over a specified interval. The study duration could be an entire year, or a fishing season for a population managed by seasonal harvest. A high exploitation rate might indicate that harvest regulations such as a size limit or creel limit could be worthwhile. Alternatively a low estimate might indicate that other factors such as natural mortality or recruitment may be more important in regulating population size. With only a single release, we cannot estimate survival (or its complement, mortality) but it is nevertheless an important first step in judging the relative impact of fishing.
There are several practical issues to consider in planning an exploitation rate study. We want to ensure that essentially all tags of harvested fish are reported, so high-reward tags (e.g., $100) could be used (Pollock et al. 2001). Tag loss is another potential source of bias. It can be estimated by double-tagging (giving some or all fish two tags, and recording returns of two versus only one tag), but it greatly complicates the analysis. For this example, we will keep things simple and assume that tag loss is negligible. We also assume that there is no immediate tagging mortality (i.e., due to capture, handling and tagging). This can be challenging to investigate. It is sometimes estimated by holding a sample of fish in a cage after tagging but this can itself be a source of mortality (Pollock and Pine 2007). For now, we assume that there is no short- or long-term mortality associated with tagging. We assume that all caught fish are harvested. If catch-and-release was occurring, our model could be extended to estimate rates of harvest and catch-and-release. The decision of how many fish to tag will depend on the required precision of the study, and as usual can be examined through simulation. Finally, it is important that tagged fish be representative of the entire population. This refers not only to size or age, but also the spatial pattern. Tagging could be done at randomly selected locations, or when the entire study population is concentrated within a single area (e.g., winter holding area or spawning location). The key is for all fish, including the tagged subset, to have the same probability of being harvested.
Our simulation code uses an assumed exploitation rate (u) of 0.4. The value could be based on a literature review or prior work from the study area or similar areas regionally. The number of returned tags has a binomial distribution:
rm(list=ls()) # Clear Environment
# Tagging study to estimate exploitation rate
u.true <- 0.4
n.tagged <- 100
n.returned <- rbinom(n=1, prob=u.true, size=n.tagged)We need only a single observation for the returned number of tags. What two lines of code could you add to look at the distribution of returned tags (using a new variable for the vector of tag returns)?
The JAGS code is also extremely simple. We have one line of code for the uninformative prior distribution for the exploitation rate, and one for the likelihood.
# Load necessary library
library(rjags)
library(R2jags)
# JAGS code
sink("ExpRate.txt")
cat("
model{
# Priors
u.est ~ dunif(0,1) # Exploitation rate
# Likelihood
n.returned ~ dbin(u.est, n.tagged)
}
",fill=TRUE)
sink()
# Bundle data
jags.data <- list("n.returned", "n.tagged")
# Initial values
jags.inits <- function(){ list(u.est=runif(n=1, min=0, max=1))}
model.file <- 'ExpRate.txt'
# Parameters monitored
jags.params <- c("u.est")
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=jags.inits,
n.chains = 3, n.thin=1, n.iter = 4000, n.burnin=2000,
model.file)
print(jagsfit)
plot(jagsfit)A single run with Ntags=100 provided a credible interval of 0.27-0.46 (from a random return of 36 tags). Experiment with different sample sizes and consider what level of precision would be adequate for management purposes. As always, this is an optimistic scenario, with the number of returned tags perfectly conforming to a binomial distribution.
7.2 Completed tagging study
A single release of tags (Section 7.1) can also be used to estimate the rate of natural death. Hearn et al. (1987) provide a method that can be used in completed tagging studies; that is, a study of sufficient duration that no live tagged fish remain (and tag returns have ceased). Their method requires only that fishing mortality in each period be greater than zero, in order to establish the study’s endpoint. It is assumed that tags of harvested fish are always reported, there is no mortality associated with tagging, and there is no tag loss. The model assumes a constant rate of natural mortality, so the estimated rate would in practice represent the average over the study period.
The following Bayesian version is a modification of the Hearn et al. (1987) method. We allow for zero harvest but assume that tag returns are monitored for a sufficiently long time to confirm the true study length (i.e., to allow for the possibility of an interspersed 0). The Hearn et al. (1987) model uses instantaneous rates, which we have to this point avoided. The relationship between probabilities (survival, death, exploitation, natural death) and instantaneous rates is as follows. The instantaneous rate of population decline is \(dN/dt = -Zt\), where N is population size and Z is the instantaneous total mortality rate (rate of decline due to fishing and natural sources). An advantage of instantaneous rates is that they are additive. We can partition Z into whatever categories are useful and estimable; for example, fishing versus natural mortality, commmercial versus recreational fishing, longline versus trawl, etc. The most common partitioning is for fishing (F) versus natural (M) sources: \(dN/dt = -(M+F)t\). Integrating this equation provides us with an expression for population size at any time t: \(N_t = N_0 * exp(-(M+F)*t)\) where N0 is initial abundance and \(S=exp(-(M+F)*t)\) is the survival rate. The probability of dying from all causes is A=1-S. The exploitation rate is obtained as u=FA/Z, which makes sense as the fraction of total deaths due to fishing (F/Z). Similarly the probability of natural death is v=MA/Z. These expressions make it possible to go back and forth between probabilities and instantaneous rates, depending on which is more useful or traditional for a particular model. These expressions apply for the most common situation where fishing and natural mortality occur simultaneously (referred to by Ricker (1975) as a Type 2 fishery). Ricker (1975) provides alternate expressions for the less common case (a Type 1 fishery) where natural mortality occurs after a short fishing season.
Unlike the probabilities for survival, harvest, etc., instantaneous rates are unbounded. In a Bayesian context, this makes them more challenging to estimate compared to probabilities, which are conveniently bounded by 0 and 1. Instantaneous rates are also less intuitive in terms of their magnitude compared to probabilities. There is typically more uncertainty about M than F because natural deaths are almost never observed. The joke about M is that it is often assumed to be 0.2, because .2 looks like a distored version of a question mark. In a review of actual published natural mortality estimates, Vetter (1988) found that a majority were less than 0.5 but values exceeding 2.0 were occasionally reported.
We begin with simulation code:
rm(list=ls()) # Clear Environment
n.tagged <- 100
Periods <- 20 # Set to high value to ensure that end of study is observed
F.true <- 0.4
M.true <- 0.2
S.true <- exp(-F.true - M.true)
u.true <- F.true * (1-S.true) / (F.true + M.true)
TagsAtRisk <- array(data=NA, dim=Periods)
TagsAtRisk[1] <- n.tagged
TagFates <- array(data=NA, dim=Periods) # Returns by period
for (j in 2:Periods){
TagsAtRisk[j] <- rbinom(n=1, size=TagsAtRisk[j-1], prob=S.true)
TagFates[j-1] <- rbinom(n=1, size=(TagsAtRisk[j-1]
-TagsAtRisk[j]),(F.true/(F.true+M.true))) # Random realization: fishing deaths
} #j
TagFates[Periods] <- rbinom(n=1, size=(TagsAtRisk[Periods-1]
-TagsAtRisk[Periods]),(F.true/(F.true+M.true)))
t.Completed <- max(which(TagFates[]>0)) # Locate period when study completed (last return)
NatDeaths <- n.tagged - sum(TagFates)
TagFates <- c(TagFates[1:t.Completed],NatDeaths) # Returns + all natural deathsThe vectors for tags-at-risk and fates are binomially distributed random realizations. Tags at risk are obtained using the number at risk in the previous period and the survival rate. Tag returns (stored in tag fates vector) are obtained from the number of deaths between periods and the fraction of deaths due to fishing (F/Z). The which() function provides a vector of periods when tag returns are greater than 0, and the max function chooses the final period with at least one tag return. Try which(TagFates[]>0) in the Console to make clear how this expression works. The final version of the tag fates vector contains tag returns for the completed study length plus a final column for natural deaths (any tagged fish not seen in the completed tagging study).
The JAGS code for fitting the model is similar to that for the tag-based method of estimating survival (Section 6.3), except that our parameters now are instantaneous rates. We use an arbitrary uniform prior distribution, with an upper bound (2) set high enough to ensure that the interval contains the true value. The likelihood uses a multinomial distribution, where cell probabilities represent the fraction returned in each period plus the final column for fish not seen again (natural deaths in our completed study).
# Load necessary library
library(rjags)
library(R2jags)
######### JAGS code for estimating model parameters
sink("Hearn.txt")
cat("
model {
# Priors
M.est ~ dunif(0, 2) # Instantaneous natural mortality rate
for (i in 1:t.Completed) {
F.est[i] ~ dunif(0,2) # Instantaneous fishing mortality rate
S.est[i] <- exp(-M.est - F.est[i])
u.est[i] <- F.est[i] * (1-S.est[i])/(F.est[i]+M.est) # FA/Z
} #i
# Cell probabilities
p.est[1] <- u.est[1]
for (i in 2:t.Completed) {
p.est[i] <- prod(S.est[1:(i-1)])*u.est[i]
} #i
p.est[t.Completed+1] <- 1 - sum(p.est[1:t.Completed]) # Prob of not being seen again
TagFates[1:(t.Completed+1)] ~ dmulti(p.est[1:(t.Completed+1)], n.tagged)
}
",fill=TRUE)
sink()
# Bundle data
jags.data <- list("TagFates", "t.Completed", "n.tagged")
# Initial values
jags.inits <- function(){list(M.est=runif(1, 0, 2), F.est=runif(t.Completed, 0, 2))}
model.file <- 'Hearn.txt'
# Parameters monitored
jags.params <- c("M.est", "F.est")
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=jags.inits,
n.chains = 3, n.iter = 4000, n.burnin=2000,
model.file)
print(jagsfit)
plot(jagsfit)This approach works because there are only two fates for tagged fish: caught and reported or natural death. Using a completed tagging study means that the number of natural deaths is known exactly (number of tag releases - total returns). The model adjusts the estimate of M in order to account for (get rid of) all the natural deaths by the end of the completed study. The length of the completed study depends on the assumed values for F and M. The estimate of M will be higher when the study length is short and when many tagged fish are not seen again.
For the default settings, estimates of M tend to be reliable, which is unsurprising given that it is modeled as a constant rate and estimated from multiple years of data. How does uncertainty about M change as you increase the number of tags released (e.g., 200 or 1000)?
Estimates of F are reliable in the first few years but uncertainty increases dramatically towards the end of the completed study. This makes sense because the estimates of tags at risk are further and further from the known initial tagged population. There would likely be some improvement from simplifying the model to estimate only the average F rather than period-specific estimates, given that the true simulated F values vary stochastically around a fixed value.
7.3 Multi-year tagging study
We return to the multi-year tagging study (Section 6.3) but now assume that we have information about the tag-reporting rate. This information allows us to partition mortality into fishing and natural sources. Tag-reporting rate is fixed at 1.0 in the following code, as could be appropriate when using high-reward tags. In a later section of code, we estimate it internally, using planted tags (Hearn et al. 2003) and an additional likelihood component. Information about the tag-reporting rate is very powerful. For example, a 100% tag-reporting rate means that (assuming no other assumptions are violated) you ultimately know the number of natural deaths because all other tags will be reported (as in Section 7.2). There is uncertainty in the short term because tags not returned could belong to fish still alive and at risk, so information about natural deaths is clear only after sufficient time has passed.
Our simulation code is similar to that of Section 6.3 except that we now specify the probabilities of harvest (u.true) and natural death (v.true). We use vectors so that the rates can vary by period. We have arbitrarily set the probabilities for harvest higher than for natural death, so it is useful to investigate whether we can reliably detect the predominant threat in our simulated study design.
rm(list=ls()) # Clear Environment
# Parameter values for three-yr experiment
Periods <- 3 # Equal number of release and return periods assumed
NTags <- 1000 # Number released each period
NumTagged <- rep(NTags, Periods) # Release NTags fish each period
u.true <- c(0.4, 0.25, 0.3) # Probability of fishing death
v.true <- c(0.1, 0.15, 0.2) # Probability of natural death
lam.true <- 1 # Reporting rate
# Calculated value(s), array creation
S.true <- 1 - u.true - v.true
TagsAtRisk <- array(data=NA, dim=c(Periods, Periods))
TagFate <- array(data=NA, dim=c(Periods, Periods+1)) # Extra column for tags not seen again
for (i in 1:Periods){ # Expected values for tags-at-risk, returned/not-seen-again
TagsAtRisk[i,i] <- NumTagged[i]
TagFate[i,i] <- TagsAtRisk[i,i] * lam.true * u.true[i] # Returns in release period
j <- i+1
while(j <= Periods) {
TagsAtRisk[i,j] <- TagsAtRisk[i,(j-1)]*S.true[j-1]
TagFate[i,j] <- TagsAtRisk[i,j] * lam.true * u.true[j]
j <- j+1
} #while
TagFate[i, Periods+1] <- NTags-sum(TagFate[i,i:Periods]) # Tags not seen again
} #i
RandRet <- array(data=NA, dim=c(Periods, Periods+1)) # Extra column for fish not seen again
#Random trial
for (i in 1:Periods) # Create upper diagonal matrix of random tag returns
{
RandRet[i,i:(Periods+1)] <- t(rmultinom(1, NumTagged[i], TagFate[i,i:(Periods+1)]))
} #iThe analysis is again similar to Section 6.3 except that now we now have a multinomial distribution with three possible fates (survival, harvest, natural death). Following Kéry and Schaub (2011) (their Section 9.6), we use three gamma-distributed parameters that are scaled by their sum, thus ensuring that they sum to 1. The parameters of interest (u, v, S) are calculated functions of the “a” parameters, which are used internally but not themselves of interest.
# Load necessary library
library(rjags)
library(R2jags)
######### JAGS code for estimating model parameters
sink("TagReturn.txt")
cat("
model {
# Priors
lam.est <- 1 # Reporting rate assumed known (e.g. high-reward tags only)
for (i in 1:Periods) {
for (j in 1:3) # Rows are return periods, columns are fates
{
a[i,j] ~ dgamma(1,1)
} #j
S.est[i] <- a[i,1]/sum(a[i,]) # Probability of survival
u.est[i] <- a[i,2]/sum(a[i,]) # Probability of harvest (exploitation rate)
v.est[i] <- a[i,3]/sum(a[i,]) # Probability of natural death
# v.est[3] could also be obtained as 1-S.est[3]-u.est[3]
} #i
# Cell probabilities
for (i in 1:Periods) {
p.est[i,i] <- lam.est * u.est[i]
for (j in (i+1):Periods) {
p.est[i,j] <- prod(S.est[i:(j-1)])*lam.est*u.est[j]
} #j
p.est[i,Periods+1] <- 1 - sum(p.est[i, i:Periods]) # Prob of not being seen again
} #i
for (i in 1:Periods) {
RandRet[i,i:(Periods+1)] ~ dmulti(p.est[i,i:(Periods+1)], NumTagged[i])
}#i
}
",fill=TRUE)
sink()
# Bundle data
jags.data <- list("RandRet", "Periods", "NumTagged")
# Initial values
a.init <- array(data=NA, dim=c(Periods, 3))
for (i in 1:3){
a.init[i,] <- runif(n=3, min=0, max=1)
a.init[i,] <- a.init[i,]/sum(a.init[i,]) # Rescale to sum to 1
}
jags.inits <- function(){ list(a=a.init)}
model.file <- 'TagReturn.txt'
# Parameters monitored
jags.params <- c(#"lam.est",
"u.est", "v.est")
#, "S.est")
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=jags.inits,
n.chains = 3, n.thin = 1, n.iter = 4000, n.burnin=2000,
model.file)
print(jagsfit)
plot(jagsfit)We provide arbitrary starting values for the “a” parameters by rescaling uniform (0,1) random variates. Using a large release on each occasion (1000), our estimates are generally close to the true values except that v.est[3] is imprecise and typically overestimated. It is the most difficult time-dependent parameter to estimate. We get direct information (from returned tags) about the final exploitation rate but there is much uncertainty about whether fish not seen were still at risk or were natural deaths. This is especially true for the final release because most fish are still at large (examine the TagFate and RandRet matrices). Note that this differs from a completed tagging study (Section 7.2), where the study continues until no live fish remain.
The analysis can be extended to include estimation of the tag-reporting rate. We use the simplest approach of a binomial experiment using planted tags. The fraction of planted tags that are returned is a direct estimate of the reporting rate. The only change required in the simulation code is to replace the assignment lam.true <- 1.0 with the following three lines:
lam.true <- 0.9 # Reporting rate
PlantedTags <- 30
PlantReturns <- rbinom(n=1, size=PlantedTags, prob=lam.true)The assumed reporting rate and number of planted tags are arbitrary and can be varied to see the effect on parameter uncertainty.
# Load necessary library
library(rjags)
library(R2jags)
######### JAGS code for estimating model parameters
sink("TagReturn.txt")
cat("
model {
# Priors
lam.est ~ dunif(0,1) # Reporting rate
for (i in 1:Periods) {
for (j in 1:3) # Rows are return periods, columns are fates
{
a[i,j] ~ dgamma(1,1)
} #j
S.est[i] <- a[i,1]/sum(a[i,]) # Probability of survival
u.est[i] <- a[i,2]/sum(a[i,]) # Probability of harvest (exploitation rate)
v.est[i] <- a[i,3]/sum(a[i,]) # Probability of natural death
# v.est[3] could also be obtained as 1-S.est[i]-u.est[i]
} #i
# Cell probabilities
for (i in 1:Periods) {
p.est[i,i] <- lam.est * u.est[i]
for (j in (i+1):Periods) {
p.est[i,j] <- prod(S.est[i:(j-1)])*lam.est*u.est[j]
} #j
p.est[i,Periods+1] <- 1 - sum(p.est[i, i:Periods]) # Prob of not being seen again
} #i
for (i in 1:Periods) {
RandRet[i,i:(Periods+1)] ~ dmulti(p.est[i,i:(Periods+1)], NumTagged[i])
}#i
PlantReturns ~ dbin(lam.est, PlantedTags) # Additional likelihood component
}
",fill=TRUE)
sink()
# Bundle data
jags.data <- list("RandRet", "Periods", "NumTagged", "PlantedTags", "PlantReturns")
# Initial values
a.init <- array(data=NA, dim=c(Periods, 3))
for (i in 1:3){
a.init[i,] <- runif(n=3, min=0, max=1)
a.init[i,] <- a.init[i,]/sum(a.init[i,]) # Rescale to sum to 1
}
jags.inits <- function(){ list(a=a.init)}
model.file <- 'TagReturn.txt'
# Parameters monitored
jags.params <- c("lam.est", "u.est", "v.est")
#, "S.est")
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=jags.inits,
n.chains = 3, n.thin = 1, n.iter = 10000, n.burnin=5000,
model.file)
print(jagsfit)
plot(jagsfit)The auxiliary planted tag experiment requires only a few changes in the JAGS code. We provide an uninformative prior distribution for the now estimated parameter lam.est, and pass in the number of planted tags and returns. We let JAGS generate an initial value for lam.est. The additional likelihood component describes the binomial experiment. Convergence may be slower for this version of the code, so the number of MCMC iterations was increased to 10,000.
Estimates have similar precision to the original example if the reporting rate is high. Precision is reduced at lower reporting rates because of the reduced sample size of returned tags and uncertainty about lambda.
7.4 Telemetry estimates of F and M
In Section 6.2, we used fixed receiving stations and telemetry tags to estimate the survival rate. Here we extend that approach to partition total mortality into fishing and natural components by giving telemetered fish an external high-reward tag (Design C of Hightower and Harris (2017)). Live fish “detect themselves” and provide information on survival by passing within range of a receiver. Harvest information comes from the reported high-reward tags. Thus we have observations on two of the three possible states, with natural mortality being the only unobserved true state. An important practical advantage of this approach is that field surveys are not needed after tagging is done (other than maintaining the receiver array and downloading detection data).
We begin with a detour to introduce the JAGS dcat distribution, used to estimate probabilities for observations from different categories.
rm(list=ls()) # Clear Environment
# Load necessary library
library(rjags)
library(R2jags)
# JAGS code for estimating model parameters
sink("dcat_example.txt")
cat("
model {
# Priors
for (j in 1:3) {
a[j] ~ dgamma(1,1)
} #j
p[1] <- a[1]/sum(a[])
p[2] <- a[2]/sum(a[])
p[3] <- a[3]/sum(a[]) # Could also be obtained as 1-p[1]-p[2]
# Likelihood
for (i in 1:N){
y[i] ~ dcat(p[])
} #i
}
",fill = TRUE)
sink()
# Bundle data
y <- c(1, 1, 3, 3, 2) # Each observation drawn from category 1, 2, or 3
N <- length(y)
jags.data <- list("N", "y")
model.file <- 'dcat_example.txt'
# Parameters monitored
jags.params <- c("p")
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=NULL,
n.chains = 3, n.thin=1, n.iter = 2000,
model.file)
print(jagsfit)
plot(jagsfit)Our dcat example has a y vector with five observations in the range of 1 to 3. The categories can be of any sort; for example, age or species code for individual fish. The model parameters of interest are calculated variables representing the probabilities of the three categories. The probabilities sum to 1, so as in Section 7.3, we calculate them using an “a” vector with an uninformative gamma prior distribution. (A Dirichlet distribution can also be used as an uninformative prior distribution for probabilities that must sum to 1 (McCarthy 2007).) Returning to the purpose of this Section, the p vector could represent probabilities of a live detection, harvest, or not being detected. The MCMC process adjusts the p estimates to achieve the best possible match between model predictions and the observations (y vector). For simplicity, we let JAGS provide the initial values for the a vector. The estimates have wide credible intervals because of the small sample size.
Returning to the telemetry example, we begin with simulation code, using arbitrary values for the number of occasions and sample size. Varying the patterns for rates of exploitation and natural death is helpful in judging the model’s ability to detect moderate time variation. This is a multistate model (Kéry and Schaub 2011), with a matrix (PSI.STATE) describing the various true states and the probabilities for transitioning between states over time. The three rows represent true states alive, harvest, and natural death. A fish that is alive (row 1) can remain alive (probability S.true), be harvested (probability u.true), or become a natural mortality (probability v.true). A harvested fish (row 2) remains in that state with probability 1, as does a natural death (row 3). Note that the third matrix dimension for PSI.STATE is time, as the rates are time-dependent.
The observation matrix (PSI.OBS) has a similar structure. Rows are true states, columns are observed states, and the third dimension is time. There are three possible observed states: detected alive, harvested, and not detected. A live fish (row 1) can either be detected or not, with probabilities p[t] and 1-p[t]. A harvested fish (row 2) is assumed to be always detected (100% reporting), and a natural death (row 3) is never detected.
The function for generating the simulated capture-history matrix (simul.ms) is complex, but key features are as follows. All fish are assumed to be tagged (transmitters and external tags) at time 1. Next we loop over occasions and use a multinomial trial to determine the true state transition between time t-1 and t (i.e., a live fish at time t-1 can remain alive, be harvested, or transition to natural death). To better understand the multinomial trial, try a simpler example in the Console. For example, rmultinom(n=1, size=1, prob=c(0.7, 0.25, 0.05)) is a single trial (n=1) for one individual (size=1), with three states at the specified probabilities. The function returns a matrix with a 1 in a randomly determined row. The which() function locates the “1” and returns the row index (= fish’s state). This example could represent probabilities that a live fish at time t-1 survives (0.7), is harvested (0.25), or is a natural death (0.05). Returning to the code, a second multinomial trial based on the fish’s true state draws a random observed outcome (e.g., a live fish is either detected or missed).
# Code modified from Kery and Schaub 2012, Section 9.5
rm(list=ls()) # Clear Environment
# Generation of simulated data
# Define probabilities as well as number of occasions, states, observations and released individuals
n.occasions <- 5
n.tagged <- 100 # Simplify from Kery's version, only an initial release, and no spatial states
u.true <- c(0.2, 0.1, 0.3, 0.2)
v.true <- c(0.1, 0.15, 0.15, 0.1)
S.true <- 1-u.true-v.true
p <- c(1, 0.8, 0.8, 0.6, 0.7)
# Detection probability fixed at 1 for first period, n.occasions-1 searches (e.g start of period 2, 3, ..)
n.states <- 3
n.obs <- 3
# Define matrices with survival, transition and detection probabilities
# 1. State process matrix. Three-dimensional matrix. 1=state, time t, 2=state, time t+1, 3=time
PSI.STATE <- array(NA, dim=c(n.states, n.states, n.occasions-1))
for (t in 1:(n.occasions-1)){
PSI.STATE[,,t] <- matrix(c(
S.true[t], u.true[t], v.true[t],
0, 1, 0,
0, 0, 1), nrow = n.states, byrow = TRUE)
} #t
# 2.Observation process matrix. 1=true state, 2=observed state, 3=time
PSI.OBS <- array(NA, dim=c(n.states, n.obs, n.occasions))
for (t in 1:n.occasions){
PSI.OBS[,,t] <- matrix(c(
p[t], 0, 1-p[t],
0, 1, 0, # Caught w/ high-reward tag, assume 100% reporting
0, 0, 1), nrow = n.states, byrow = TRUE) # Natural deaths never detected
} #t
# Define function to simulate multistate capture-recapture data
simul.ms <- function(PSI.STATE, PSI.OBS, n.tagged, n.occasions, unobservable = NA){
# Unobservable: number of state that is unobservable
CH <- CH.TRUE <- matrix(NA, ncol = n.occasions, nrow = n.tagged)
CH[,1] <- CH.TRUE[,1] <- 1 # All releases at t=1
for (i in 1:n.tagged){
for (t in 2:n.occasions){
# Multinomial trials for state transitions
CH.TRUE[i,t] <- which(rmultinom(n=1, size=1,
prob=PSI.STATE[CH.TRUE[i,t-1],,t-1])==1)
# which vector element=1; i.e., which state gets the random draw
# at time t given true state at t-1
# True state determines which row of PSI.STATE provides the probabilities
# Multinomial trials for observation process
CH[i,t] <- which(rmultinom(1, 1, PSI.OBS[CH.TRUE[i,t],,t])==1)
# which observation gets the 1 random draw, given true time t state.
} #t
} #i
return(list(CH=CH, CH.TRUE=CH.TRUE)) # True (CH.TRUE) and observed (CH)
} # simul.msThe JAGS code is also complex but does generally mirror the simulation code.
# Load necessary library
library(rjags)
library(R2jags)
# JAGS code for estimating model parameters
sink("Mort_Tel_HRTag.txt")
cat("
model {
# Priors
for (t in 1:(n.occasions-1)) {
for (j in 1:3) { # Rows are return periods, columns are fates
a[t,j] ~ dgamma(1,1)
} #j
S.est[t] <- a[t,1]/sum(a[t,]) # Probability of survival
u.est[t] <- a[t,2]/sum(a[t,]) # Probability of harvest (exploitation rate)
v.est[t] <- a[t,3]/sum(a[t,]) # Probability of natural death
# v.est[t] could also be obtained as 1-S.est[t]-u.est[t]
} #t
p[1] ~ dbern(1)
for (t in 2:n.occasions){
p[t] ~ dunif(0,1)
}#t
# Define state-transition and observation matrices
# Probabilities of true state at t+1 given state at time t
for (t in 1:(n.occasions-1)){
ps[1,t,1] <- S.est[t]
ps[1,t,2] <- u.est[t]
ps[1,t,3] <- v.est[t]
ps[2,t,1] <- 0
ps[2,t,2] <- 1
ps[2,t,3] <- 0
ps[3,t,1] <- 0
ps[3,t,2] <- 0
ps[3,t,3] <- 1
} #t
# Probabilities of observed states given true state
for (t in 2:n.occasions){
po[1,t,1] <- p[t] # Row 1 true state = alive
po[1,t,2] <- 0
po[1,t,3] <- 1-p[t]
po[2,t,1] <- 0 # Row 2 true state = harvested
po[2,t,2] <- 1
po[2,t,3] <- 0
po[3,t,1] <- 0 # Row 3 true state= natural death
po[3,t,2] <- 0
po[3,t,3] <- 1
} #t
# Likelihood
for (i in 1:n.tagged){
z[i,1] <- y[i,1] # Latent state known at time of tagging (t=1)
for (t in 2:n.occasions){
# State process: draw state at time t given state at time t-1
z[i,t] ~ dcat(ps[z[i,(t-1)], (t-1),])
# Observation process: draw observed state given true state at time t
y[i,t] ~ dcat(po[z[i,t], t,])
} #t
} #i
}
",fill = TRUE)
sink()
# Generate initial values for z. Modified from Kery code for JAGS inits, age-specific example 9.5.3
z.init <- function(ch) {
# State 1=Obs 1 (alive) State 2=Obs 2 (fishing death). Start w/ known states from obs capture-history,
# replace "3" with possible state(s)
ch[ch==3] <- -1 # Not observed so temporarily replace w/ neg value
ch[,1] <- NA # Initial value not needed for release period
for (i in 1:nrow(ch)){
if(max(ch[i,], na.rm=TRUE)<2){
ch[i, 2:ncol(ch)] <- 1 # Not detected dead so initialize as alive (after release period)
} else {
m <- min(which(ch[i,]==2)) # Period when fishing death first detected
if(m>2) ch[i, 2:(m-1)] <- 1 # Initialize as alive up to period prior to harvest
ch[i, m:ncol(ch)] <- 2 # Initialize as dead after detected fishing death
} # if/else
} # i
return(ch)
} # z.init
# Call function to get simulated true (CH.TRUE) and observed (CH) states
sim <- simul.ms(PSI.STATE, PSI.OBS, n.tagged, n.occasions)
y <- sim$CH
# Bundle data
jags.data <- list("n.occasions", "n.tagged", "y")
# Initial values.
jags.inits <- function(){ list(z = z.init(y))}
model.file <- 'Mort_Tel_HRTag.txt'
# Parameters monitored
jags.params <- c("u.est", "v.est", "p")
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=jags.inits,
n.chains = 3, n.thin=1, n.iter = 5000,
model.file)
print(jagsfit)
plot(jagsfit)The parameters for rates of survival, exploitation and natural death use the same approach as in the above dcat example.. The “a” matrix is estimated internally using a vague gamma prior distribution, then the three fractions estimate the probabilities of the three possible true states (survival, harvest, natural death), which sum to 1. The latent true states (z matrix) are estimated in the likelihood section of the code. The dcat (categorical distribution) is used to draw the true state at time t given state at time t-1, as well as the observation at each time given the true state.
One difficult part of fitting this model is obtaining initial values for the latent true states. JAGS will refuse to update model parameters if initial values are inconsistent with the observations (link). The function used here z.init() initializes all unknown true states to be “alive” unless a harvest occurs. This works because fish not detected could either be alive or a natural death, so “alive” is a valid initial value. Individuals that are ultimately harvested are known to be alive up to the period prior to harvest. It can be instructive to compare the matrices containing true and observed states with the initial values. For example, initial values can be saved by entering z.test <- z.init(y) in the Console. Entering sim$CH.TRUE[1:5,] and CH[1:5,] in the Console will display true and observed states for the first five individuals for comparison to z.test.
Results for the chosen simulation settings are moderately reliable. There are many parameters, as the model estimates true latent states as well as the probabilities for survival, harvest and natural death. It seems to consistently detect the increase in harvest rate in period 3. Uncertainty increases toward the end of the study, especially for probabilities of detection and natural death. This makes sense because fish not detected on the final occasion could either be alive or a natural death.
Perhaps the most obvious simulation settings to vary are the fairly robust sample size of tagged fish (100) and the detection probabilities (0.6 and above). Uncertainty increases markedly as detection probabilities decrease, so these simulations can be very helpful in planning the field study (e.g., number of receiver stations). The number of occasions can also be varied but it requires commensurate changes in the vectors for u, v, and p.
7.5 Exercises
For the exploitation rate study design (Section 7.1), modify the JAGS code to plot the posterior distribution for the exploitation rate. Include a vertical line showing the underlying true exploitation rate.
For the completed tagging study design (Section 7.2), modify the code to use a single F parameter (representing an average over periods).
For the multi-year tagging study (Section 7.3), run three replicate simulations of the version using planted tags, with tag releases of 1000 (current code), 100, and 30. How does uncertainty (e.g. credible interval width) vary? What sample size would you recommend if planning this study? Also, save credible intervals for v3 for the 1000 fish releases (for exercise 4 below).
Compare results for v3 (from exercise 3) with a modified design using five periods. Extend the u and v vectors by using the values from periods 1-2 for the final two periods. How do credible intervals for v3 change when the study is extended?
For the telemtry study (Section 7.4), compare results (particularly uncertainty) for parameter v4 when p5 takes on the following values: 0.2, 0.4, 0.7, 0.9. How might this affect planning for conducting the field work?