Kaslo Bull Trout Productivity 2017

The suggested citation for this analytic report is:

Thorley, J.L. and Dalgarno, S.I.J. (2018) Kaslo Bull Trout Productivity 2017. A Poisson Consulting Analysis Report. URL: http://www.poissonconsulting.ca/f/443864882.

Background

The Kaslo River and Keen Creek, which is a tributary of the Kaslo River, are important Bull Trout spawning and rearing tributaries on Kootenay Lake. From 2012 to 2017, field crews have night-snorkeled both systems in the fall and recorded all juvenile Bull Trout between 30 mm and 350 mm in length. Snorkel and electrofishing marking crews have also captured and tagged juvenile Bull Trout for the snorkel crews to resight. Redd counts have been conducted in both systems since 2006.

The primary goal of the current analyses is to answer the following questions:

What is the observer efficiency when night-snorkeling for juvenile Bull Trout in the Kaslo River and Keen Creek?

What are the densities of age-1 Bull Trout in the Kaslo River and Keen Creek?

What is the relationship between the number of redds and the resultant densities of age-1 Bull Trout two years later?

Data Preparation

The data were cleaned, tidied and databased using R version 3.4.4 (R Core Team 2015).

Table 1. Bull Trout length cutoffs by age-class.

Age Class Length (mm)
Age-0 30 - 90
Age-1 91 - 150
Age-2+ 151 - 350

Statistical Analysis

Model parameters were estimated using Maximum-Likelihood (ML) and Bayesian methods. The Maximum-Likelihood Estimates (MLE) were obtained using TMB (Kristensen et al. 2016) while the Bayesian estimates were produced using JAGS (Plummer 2015). For additional information on ML and Bayesian estimation the reader is referred to Millar (2011) and McElreath (2016), respectively.

Unless indicated otherwise, the Bayesian analyses used uninformative normal prior distributions (Kery and Schaub 2011, 36). The posterior distributions were estimated from 1500 Markov Chain Monte Carlo (MCMC) samples thinned from the second halves of 3 chains (Kery and Schaub 2011, 38–40). Model convergence was confirmed by ensuring that \(\hat{R} \leq 1.1\) (Kery and Schaub 2011, 40) and \(\textrm{ESS} \geq 150\) for each of the monitored parameters (Kery and Schaub 2011, 61). Where \(\hat{R}\) is the potential scale reduction factor and \(\textrm{ESS}\) is the effective sample size.

The parameters are summarised in terms of the point estimate, standard deviation (sd), the z-score, lower and upper 95% confidence/credible limits (CLs) and the p-value (Kery and Schaub 2011, 37, 42). For ML models, the point estimate is the MLE, the standard deviation is the standard error, the z-score is \(\mathrm{sd}/\mathrm{MLE}\) and the 95% CLs are the \(\mathrm{MLE} \pm 1.96 \cdot \mathrm{sd}\). For Bayesian models, the estimate is the median (50th percentile) of the MCMC samples, the z-score is \(\mathrm{sd}/\mathrm{mean}\) and the 95% CLs are the 2.5th and 97.5th percentiles. A p-value of 0.05 indicates that the lower or upper 95% CL is 0.

Model averaging and selection was implemented using an information theoretic approach (Burnham and Anderson 2002). ML models were compared using the marginal Akaike’s Information Criterion corrected for small sample size (AICc, Burnham and Anderson 2002). The support for fixed terms in models with identical random (Kery and Schaub 2011, 75) effects was assessed using the Akaike’s weights (\(w_i\)) which were also used for model averaging (Burnham and Anderson 2002).

Where relevant, model adequacy was confirmed by examination of residual plots for the full model(s).

The results are displayed graphically by plotting the modeled relationships between particular variables and the response(s) with the remaining variables held constant. In general, continuous and discrete fixed variables are held constant at their mean and first level values, respectively, while random variables are held constant at their typical values (expected values of the underlying hyperdistributions) (Kery and Schaub 2011, 77–82). When informative the influence of particular variables is expressed in terms of the effect size (i.e., percent change in the response variable) with 95% confidence/credible intervals (CIs, Bradford, Korman, and Higgins 2005).

The analyses were implemented using R version 3.4.4 (R Core Team 2015) and the mbr family of packages.

Model Descriptions

Observer Efficiency

All resighted fish with a tag were allocated to the closest unallocated marked fish (with the same colour tag) by fork length and distance. The marked fish were analysed using a logistic regression model. Key assumptions of the full Maximum Likelihood logistic regression include:

Preliminary analysis indicated that river kilometre received similar support to watershed area.

The final Bayesian logistic regression assumed that:

Lineal Density

Both systems were broken into 100 m sites by bank. The lineal density at each site was estimated using an over-dispersed Poisson Generalized Linear Mixed Model. Key assumptions of the full Maximum Likelihood GLMM include:

Preliminary analysis indicated that river kilometre was better supported as a predictor of lineal density than watershed area.

The final Bayesian GLMM dropped sinuosity and gradient and took into account the uncertainty in the observer efficiencies as estimated by the observer efficiency model.

Stock-Recruitment

The stock-recruitment relationship was estimated using a Bayesian Beverton-Holt stock-recruitment curve. Key assumptions of the final BH SR model include:

Model Templates

Observer Efficiency

Maximum Likelihood

#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Observed);
DATA_FACTOR(System);
DATA_VECTOR(Tagged);
DATA_VECTOR(Length);
DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Rkm);

PARAMETER(bIntercept);
PARAMETER(bSystem);
PARAMETER(bLength);
PARAMETER(bLength2);
PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bRkm);

vector<Type> eObserved = Observed;

int n = Observed.size();

Type nll = 0.0;

for(int i = 0; i < n; i++){
eObserved(i) = invlogit(bIntercept + bSystem * System(i) + bLength * Length(i) + bLength2 * pow(Length(i), 2) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) + bRkm * Rkm(i));
nll -= dbinom(Observed(i), Tagged(i), eObserved(i), true);
return nll;
..

Template 1. Full model.

Bayesian

model{
  bIntercept ~ dnorm(0, 5^-2)
  bLength ~ dnorm(0, 5^-2)

  for(i in 1:length(Observed)){
    logit(eObserved[i]) <-  bIntercept + bLength * Length[i]
    Observed[i] ~ dbern(eObserved[i])
  }
..

Template 2. Final model.

Lineal Density

Maximum Likelihood

#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() () {
DATA_VECTOR(Count);
DATA_VECTOR(Length);
DATA_VECTOR(Coverage);
DATA_VECTOR(Efficiency);
DATA_VECTOR(Dispersion);

DATA_VECTOR(Gradient);
DATA_VECTOR(Sinuosity);
DATA_VECTOR(Rkm);

DATA_FACTOR(System);

DATA_FACTOR(Site);
DATA_INTEGER(nSite);
DATA_FACTOR(Year);

PARAMETER_MATRIX(bSystemYear);
PARAMETER_VECTOR(bSite);

PARAMETER(bGradient);
PARAMETER(bSinuosity);
PARAMETER(bRkm);

PARAMETER_VECTOR(bDispersion);
DATA_INTEGER(nDispersion);
PARAMETER(log_sDispersion);
PARAMETER(log_sSite);

Type sSite = exp(log_sSite);
Type sDispersion = exp(log_sDispersion);

ADREPORT(sSite);
ADREPORT(sDispersion);

vector<Type> eDensity = Count;
vector<Type> eCount = Count;
vector<Type> eNorm = pnorm(bDispersion, Type(0), Type(1));
vector<Type> eDispersion = qgamma(eNorm, pow(sDispersion, -2), pow(sDispersion, 2));

Type nll = 0.0;

for(int i = 0; i < nSite; i++){
nll -= dnorm(bSite(i), Type(0), sSite, true);

for(int i = 0; i < nDispersion; i++){
nll -= dnorm(bDispersion(i), Type(0), Type(1), true);

eDensity(i) = exp(bSystemYear(System(i), Year(i)) + bGradient * Gradient(i) + bSinuosity * Sinuosity(i) +  bRkm * Rkm(i) + bSite(Site(i)));
eCount(i) = eDensity(i) * Length(i) * Coverage(i);

nll -= dpois(Count(i), eCount(i) * Efficiency(i) * eDispersion(i), true);

return nll;
..

Template 3. The full model.

Bayesian

model{

  bEfficiency ~ dnorm(Efficiency[1], EfficiencySD[1]^-2) T(0, )

  for(i in 1:nSystem) {
    for(j in 1:nYear) {
      bSystemYear[i,j] ~ dnorm(0, 5^-2)
    }
  }

  log_sSite ~ dnorm(0, 5^-2)
  log(sSite) <- log_sSite
  for(i in 1:nSite) {
    bSite[i] ~ dnorm(0, sSite^-2)
  }

  bRkm ~ dnorm(0, 5^-2)

  log_sDispersion ~ dnorm(0, 5^-2)
  log(sDispersion) <- log_sDispersion


  for(i in 1:length(Count)) {
    log(eDensity[i]) <- bSystemYear[System[i],Year[i]] + bRkm * Rkm[i] + bSite[Site[i]]
    eDispersion[i] ~ dgamma(sDispersion^-2, sDispersion^-2)
    Count[i] ~ dpois(eDensity[i] * Length[i] * Coverage[i] * bEfficiency * eDispersion[i])
  }
..

Template 4. The final model.

Stock-Recruiment

model {
  log_bAlpha ~ dnorm(5, 5^-2)
  log(bAlpha) <- log_bAlpha

  for(i in 1:nSystem) {
    log_bBeta[i] ~ dnorm(0, 5^-2)
    log(bBeta[i]) <- log_bBeta[i]
  }

  log_sRecruits ~ dnorm(0, 5^-2)
  log(sRecruits) <- log_sRecruits

  for(i in 1:length(Recruits)){
    eRecruits[i] <- bAlpha * Stock[i] / (1 + bBeta[System[i]] * Stock[i])
    Recruits[i] ~ dlnorm(log(eRecruits[i]), sRecruits^-2)
  }
  bCarryingCapacity <- bAlpha / bBeta
..

Template 5. Final model.

Results

Tables

Coverage

Table 2. Total length of river bank counted (including replicates) by system and year.

System Year Length (km)
Kaslo River 2012 10.57 km
Kaslo River 2013 13.43 km
Kaslo River 2014 11.45 km
Kaslo River 2015 7.32 km
Kaslo River 2016 11.59 km
Kaslo River 2017 9.79 km
Keen Creek 2012 1.44 km
Keen Creek 2013 0.95 km
Keen Creek 2014 0.67 km
Keen Creek 2015 0.72 km
Keen Creek 2016 0.85 km
Keen Creek 2017 1.68 km

Fish

Table 3. Number of fish observed by system, age-class and species. All species are assumed to have the Bull Trout age-class length cutoffs.

System AgeClass Year Bull Trout Brook Trout Rainbow Trout
Kaslo River Age-0 2012 104 0 1
Kaslo River Age-0 2013 104 1 41
Kaslo River Age-0 2014 253 1 49
Kaslo River Age-0 2015 89 1 33
Kaslo River Age-0 2016 187 4 138
Kaslo River Age-0 2017 172 9 143
Kaslo River Age-1 2012 151 2 86
Kaslo River Age-1 2013 161 12 63
Kaslo River Age-1 2014 113 10 66
Kaslo River Age-1 2015 166 3 81
Kaslo River Age-1 2016 172 20 277
Kaslo River Age-1 2017 245 25 233
Kaslo River Age-2+ 2012 173 29 205
Kaslo River Age-2+ 2013 196 14 174
Kaslo River Age-2+ 2014 207 21 140
Kaslo River Age-2+ 2015 69 3 61
Kaslo River Age-2+ 2016 114 4 188
Kaslo River Age-2+ 2017 169 36 361
Keen Creek Age-0 2012 8 0 0
Keen Creek Age-0 2013 4 0 0
Keen Creek Age-0 2014 8 0 0
Keen Creek Age-0 2015 2 1 4
Keen Creek Age-0 2016 1 0 0
Keen Creek Age-0 2017 4 0 5
Keen Creek Age-1 2012 40 0 3
Keen Creek Age-1 2013 20 0 0
Keen Creek Age-1 2014 21 0 1
Keen Creek Age-1 2015 43 0 4
Keen Creek Age-1 2016 4 0 0
Keen Creek Age-1 2017 47 3 19
Keen Creek Age-2+ 2012 88 1 24
Keen Creek Age-2+ 2013 61 3 10
Keen Creek Age-2+ 2014 42 0 7
Keen Creek Age-2+ 2015 40 2 9
Keen Creek Age-2+ 2016 19 1 10
Keen Creek Age-2+ 2017 53 8 36

Observer Efficiency

Table 4. Parameter descriptions.

Parameter Description
bGradient The effect of Gradient on bIntercept
bIntercept The intercept for logit(eObserved)
bLength The effect of Length on bIntercept
bLength2 The effect of Length on the effect of Length on bIntercept
bRkm The effect of Rkm on bIntercept
bSinuosity The effect of Sinuosity on bIntercept
eObserved The expected probability of observing an individual
Gradient The standardized site-level gradient
Length The standardized fork length
Observed The number of individuals observed (0 or 1)
Rkm The standardized Rkm
Sinuosity The standardized site-level sinuosity
Tagged The number of tagged individuals (1)

Maximum Likelihood

Table 5. Model averaged parameter estimates and Akaike weights.

term estimate sd zscore lower upper pvalue nmodels proportion ICWt
bGradient 0.1073715 0.1639901 0.6547440 -0.2140432 0.4287862 0.5126326 48 0.5000000 0.49
bIntercept -1.3314248 0.2244021 -5.9332104 -1.7712448 -0.8916048 0.0000000 48 1.0000000 1.00
bLength 0.7757319 0.2170603 3.5738087 0.3503016 1.2011623 0.0003518 48 0.6666667 1.00
bLength2 -0.0195316 0.1099985 -0.1775623 -0.2351247 0.1960615 0.8590667 48 0.3333333 0.28
bRkm -0.0296348 0.1285408 -0.2305475 -0.2815702 0.2223006 0.8176663 48 0.5000000 0.31
bSinuosity -0.0190170 0.1021890 -0.1860967 -0.2193037 0.1812697 0.8523689 48 0.5000000 0.29
bSystem 0.0602404 0.2896767 0.2079574 -0.5075155 0.6279963 0.8352622 48 0.5000000 0.30

Bayesian

Table 6. Final parameter estimates.

term estimate sd zscore lower upper pvalue
bIntercept -1.3040615 0.1946852 -6.746730 -1.7185581 -0.9392931 7e-04
bLength 0.8002934 0.2041332 3.933018 0.4097444 1.2043828 7e-04

Table 7. Observer Efficiency estimates for a 123 mm Bull Trout.

System Efficiency EfficiencySD EfficiencyLower EfficiencyUpper
Kaslo River 0.1644315 0.0322458 0.1060672 0.2318914

Table 8. Final model summary.

n K nchains niters nthin ess rhat converged
190 2 3 500 1 540 1 TRUE

Lineal Density

Table 9. Parameter descriptions.

Parameter Description
bDispersion The random effect of overdispersion
bGradient The effect of Gradient on bSystemYear
bRkm The effect of Rkm on bSystemYear
bSinuosity The effect of Sinuosity on bSystemYear
bSite The random effect of Site on bSystemYear
bSystemYear The effect of System and Year on log(eDensity)
Count Number of fish counted
Coverage Proportion of site surveyed
Dispersion Factor for random effect of overdispersion
eCount The expected Count
eDensity The expected lineal density
Efficiency The observer efficiency from the observer efficiency model
Gradient Standardized gradient
Length Length of site (m)
log_sDispersion log(sDispersion)
log_sSite log(sSite)
Rkm Standardized Rkm
sDispersion The SD of bDispersion
Sinuosity Standardized sinuosity
Site The site
sSite The SD of bSite
System The system
Year The year

Maximum Likelihood

Table 10. Model averaged parameter estimates and Akaike weights.

term estimate sd zscore lower upper pvalue nmodels proportion ICWt
bGradient 0.0219370 0.0424145 0.5172052 -0.0611939 0.1050678 0.6050129 8 0.5 0.39
bRkm 0.3753258 0.0581847 6.4505957 0.2612859 0.4893656 0.0000000 8 0.5 1.00
bSinuosity 0.0113107 0.0318410 0.3552244 -0.0510965 0.0737179 0.7224215 8 0.5 0.32
bSystemYear[1,1] -2.7708561 0.1217848 -22.7520757 -3.0095499 -2.5321624 0.0000000 8 1.0 1.00
bSystemYear[2,1] -1.4433005 0.2625496 -5.4972499 -1.9578882 -0.9287129 0.0000000 8 1.0 1.00
bSystemYear[1,2] -2.8854328 0.1104037 -26.1352859 -3.1018201 -2.6690455 0.0000000 8 1.0 1.00
bSystemYear[2,2] -1.5938008 0.3439549 -4.6337491 -2.2679401 -0.9196616 0.0000036 8 1.0 1.00
bSystemYear[1,3] -3.1408441 0.1260153 -24.9243169 -3.3878294 -2.8938587 0.0000000 8 1.0 1.00
bSystemYear[2,3] -1.2106845 0.3641526 -3.3246626 -1.9244105 -0.4969586 0.0008853 8 1.0 1.00
bSystemYear[1,4] -2.2830433 0.1292692 -17.6611514 -2.5364063 -2.0296803 0.0000000 8 1.0 1.00
bSystemYear[2,4] -0.4973132 0.3259981 -1.5255094 -1.1362577 0.1416314 0.1271321 8 1.0 1.00
bSystemYear[1,5] -2.6779658 0.1144346 -23.4017121 -2.9022536 -2.4536781 0.0000000 8 1.0 1.00
bSystemYear[2,5] -3.1198760 0.5727878 -5.4468272 -4.2425195 -1.9972326 0.0000001 8 1.0 1.00
bSystemYear[1,6] -2.1695035 0.1041407 -20.8324317 -2.3736155 -1.9653916 0.0000000 8 1.0 1.00
bSystemYear[2,6] -1.4367724 0.2572289 -5.5855792 -1.9409318 -0.9326130 0.0000000 8 1.0 1.00
log_sDispersion -0.5637376 0.1072507 -5.2562592 -0.7739451 -0.3535300 0.0000001 8 1.0 1.00
log_sSite -0.8380570 0.1586402 -5.2827533 -1.1489860 -0.5271279 0.0000001 8 1.0 1.00

Bayesian

Table 11. Parameter estimates.

term estimate sd zscore lower upper pvalue
bEfficiency 0.1583234 0.0320791 4.9908766 0.1023206 0.2247109 0.0007
bRkm 0.3787761 0.0579048 6.5537086 0.2642380 0.4962428 0.0007
bSystemYear[1,1] -2.6910260 0.2361624 -11.3961262 -3.1263956 -2.2159095 0.0007
bSystemYear[2,1] -1.3574257 0.3420859 -3.9344371 -1.9901016 -0.6522819 0.0007
bSystemYear[1,2] -2.8065407 0.2320346 -12.1051927 -3.2538568 -2.3341389 0.0007
bSystemYear[2,2] -1.4690385 0.4024343 -3.6494544 -2.2516102 -0.7195669 0.0027
bSystemYear[1,3] -3.0661694 0.2391919 -12.8086954 -3.5342248 -2.5983214 0.0007
bSystemYear[2,3] -1.0931471 0.4182975 -2.5855224 -1.9154105 -0.2378714 0.0107
bSystemYear[1,4] -2.2086802 0.2362729 -9.3479884 -2.6667646 -1.7470974 0.0007
bSystemYear[2,4] -0.3694095 0.3797515 -0.9805993 -1.1532209 0.3684330 0.3200
bSystemYear[1,5] -2.6082042 0.2307494 -11.2686788 -3.0291092 -2.1198527 0.0007
bSystemYear[2,5] -3.0325974 0.5996325 -5.0925603 -4.3506536 -1.9780550 0.0007
bSystemYear[1,6] -2.0985962 0.2337713 -8.9589431 -2.5122642 -1.6089384 0.0007
bSystemYear[2,6] -1.3305624 0.3253582 -4.0528784 -1.9297537 -0.6908107 0.0007
log_sDispersion -0.4928908 0.1056907 -4.7088345 -0.7210984 -0.3098148 0.0007
log_sSite -0.8466801 0.1656018 -5.1588087 -1.2155074 -0.5537102 0.0007

Table 12. Model summary.

n K nchains niters nthin ess rhat converged
758 16 3 500 20 225 1.02 TRUE

Stock-Recruiment

Table 13. Parameter descriptions.

Parameter Description
bBeta Density-dependence
eRecruits Expected value of Recruits
Recruits Number of age-1 Bull Trout
Stock Number of Bull Trout redds
bAlpha Recruits per stock at low stock density
sRecruits SD of residual variation in Recruits

Bayesian

Table 14. Final parameter estimates.

term estimate sd zscore lower upper pvalue
bAlpha 60.6428203 2.015650e+04 0.1301767 24.6287527 1.510198e+04 0.0007
bBeta[1] 0.2577378 1.171675e+02 0.1269806 0.0154937 8.720062e+01 0.0007
bBeta[2] 0.0653216 7.132601e+01 0.1221502 0.0000299 5.018846e+01 0.0007
bCarryingCapacity[1] 238.8018290 7.010678e+04 0.0493589 111.2201443 1.605357e+03 0.0007
bCarryingCapacity[2] 905.2124317 5.486044e+08 0.0330460 234.0531743 1.327449e+06 0.0007
log_bAlpha 4.1050012 1.615337e+00 2.8894584 3.2039118 9.622578e+00 0.0007
log_bBeta[1] -1.3558127 2.024859e+00 -0.4648962 -4.1673276 4.468051e+00 0.4173
log_bBeta[2] -2.7284332 3.422965e+00 -0.8745201 -10.4179713 3.915781e+00 0.2880
log_sRecruits -0.4741538 2.450922e-01 -1.8411415 -0.8729008 1.173366e-01 0.1000
sRecruits 0.6224115 1.772484e-01 3.7079487 0.4177382 1.124500e+00 0.0007

Table 15. Final model summary.

n K nchains niters nthin ess rhat converged
12 10 3 500 100 394 1.012 TRUE

Figures

Systems

figures/rkm/map.png
Figure 1. Spatial distribution of fish-bearing channel.
figures/rkm/elevation.png
Figure 2. Channel elevation by river kilometre and system.
figures/rkm/area.png
Figure 3. Catchment area by river kilometre and system.

Sites

figures/site/gradient.png
Figure 4. Site gradient by river kilometre and system.
figures/site/sinuosity.png
Figure 5. Site sinuosity by river kilometre and system.

Coverage

figures/visit/count.png
Figure 6. Snorkel count coverage by year and bank.

Fish

figures/fish/capture.png
Figure 7. Length-frequency plot of marked Bull Trout by year and system.
figures/fish/count.png
Figure 8. Length-frequency plot of counted Bull Trout by observer.
figures/fish/freq.png
Figure 9. Length-frequency plot of observed Bull Trout by year.

Observer Efficiency

Bayesian

figures/observer/jmb/capture.png
Figure 10. Distribution of marked juvenile Bull Trout by year.
figures/observer/jmb/length.png
Figure 11. Estimated observer efficiency by fork length and system (with 95% CIs).

Lineal Density

Bayesian

figures/density/jmb/coverage.png
Figure 12. Percent coverage by year and system.
figures/density/jmb/year.png
Figure 13. Estimated density by year and system (with 95% CIs).
figures/density/jmb/abundance.png
Figure 14. The estimated abundance by year and system (with 95% CIs).
figures/density/jmb/site.png
Figure 15. The estimated density by site and system (with 95% CIs).

Stock-Recruiment

Bayesian

figures/redds/jmb/redds.png
Figure 16. Complete redds by system and spawn year.
figures/redds/jmb/stock.png
Figure 17. Estimated stock-recruitment relationship (with 95% CIs). The points are labelled by spawn year.

Conclusions

Recommendations

Acknowledgements

The organisations and individuals whose contributions have made this analysis report possible include:

References

Bradford, Michael J, Josh Korman, and Paul S Higgins. 2005. “Using Confidence Intervals to Estimate the Response of Salmon Populations (Oncorhynchus Spp.) to Experimental Habitat Alterations.” Canadian Journal of Fisheries and Aquatic Sciences 62 (12): 2716–26. https://doi.org/10.1139/f05-179.

Burnham, Kenneth P., and David R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Vol. Second Edition. New York: Springer.

Kery, Marc, and Michael Schaub. 2011. Bayesian Population Analysis Using WinBUGS : A Hierarchical Perspective. Boston: Academic Press. http://www.vogelwarte.ch/bpa.html.

Kristensen, Kasper, Anders Nielsen, Casper W. Berg, Hans Skaug, and Bradley M. Bell. 2016. “TMB : Automatic Differentiation and Laplace Approximation.” Journal of Statistical Software 70 (5). https://doi.org/10.18637/jss.v070.i05.

McElreath, Richard. 2016. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. Chapman & Hall/CRC Texts in Statistical Science Series 122. Boca Raton: CRC Press/Taylor & Francis Group.

Millar, R. B. 2011. Maximum Likelihood Estimation and Inference: With Examples in R, SAS, and ADMB. Statistics in Practice. Chichester, West Sussex: Wiley.

Plummer, Martyn. 2015. “JAGS Version 4.0.1 User Manual.” http://sourceforge.net/projects/mcmc-jags/files/Manuals/4.x/.

R Core Team. 2015. “R: A Language and Environment for Statistical Computing.” Vienna, Austria: R Foundation for Statistical Computing. http://www.R-project.org/.