Kaslo Bull Trout Productivity 2020

The suggested citation for this analytic appendix is:

Thorley, J.L., Amies-Galonski, E. and Dalgarno, S.I.J. (2020) Kaslo Bull Trout Productivity 2020. A Poisson Consulting Analysis Appendix. URL: https://www.poissonconsulting.ca/f/1283025452.

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 2020, field crews have night-snorkeled these systems in the fall and recorded all juvenile Bull Trout less than 350 mm in length. Keen Creek was not snorkelled in 2020 due to visibility issues. 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, with the exception of 2020, in which Keen Creek was not surveyed. 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 4.0.3 (R Core Team 2015).

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, lower and upper 95% credible limits (CLs) and the surprisal s-value (Greenland 2019). The estimate is the median (50th percentile) of the MCMC samples while the 95% CLs are the 2.5th and 97.5th percentiles. The s-value can be considered a test of directionality. More specifically it indicates how surprising (in bits) it would be to discover that the true value of the parameter is in the opposite direction to the estimate. An s-value of 4.3 bits, which is equivalent to a p-value (Kery and Schaub 2011; Greenland and Poole 2013) of 0.05, indicates that the surprise would be equivalent to throwing 4.3 heads in a row.

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\)) (Burnham and Anderson 2002). The best supported model was then refitted as its Bayesian equivalent.

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 4.0.3 (R Core Team 2015) and the mbr family of packages.

Model Descriptions

Length Correction

The annual bias (inaccuracy) and error (imprecision) in observer’s fish length estimates when spotlighting (standing) and snorkeling were quantified from the divergence of their length distribution from the length distribution for JLT and SH (the two most experience snorkelers) in that year. More specifically, the length correction that minimised the Jensen-Shannon divergence (Lin 1991) provided a measure of the inaccuracy while the minimum divergence (the Jensen-Shannon divergence was calculated with log to base 2 which means it lies between 0 and 1) provided a measure of the imprecision.

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:

  • The observer efficiency varies by the fork length as a second order polynomial.
  • The observer efficiency varies by observer.
  • The observer efficiency varies between systems.
  • The observer efficiency varies by the gradient, sinuosity and river kilometre.

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

The final Bayesian logistic regression assumed that:

  • The observer efficiency varies by the fork length.

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:

  • The lineal density varies between systems and years.
  • The lineal density varies randomly by site.
  • The lineal density varies by site sinuosity, gradient and river kilometre.
  • The observer efficiency for each system is as estimated by the observer efficiency model.

Preliminary analysis indicated that river kilometre was better supported as a predictor of lineal density than watershed area. Preliminary analysis also indicated that autocorrelation (modeled as an AR1 process) was not significant.

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:

  • The strength of density-dependence varies between systems.

Model Templates

Observer Efficiency

Maximum Likelihood
.#include <TMB.hpp>

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

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

PARAMETER_VECTOR(bObserver);

vector<Type> eObserved = Observed;

int n = Observed.size();

Type nll = 0.0;

for(int i = 0; i < n; i++){
eObserved(i) = invlogit(bObserver(Observer(i)) + 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;

Block 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])
  }

Block 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;

Block 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])
  }

Block 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

Block 5. Final model.

Results

Tables

Coverage

Table 1. 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]
Kaslo River 2018 8.44 [km]
Kaslo River 2019 7.19 [km]
Kaslo River 2020 10.42 [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]
Keen Creek 2018 3.37 [km]
Keen Creek 2019 1.48 [km]

Observer Efficiency

Table 2. 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 3. Model averaged parameter estimates and Akaike weights.

term estimate lower upper svalue
bGradient 0.0918623 -0.2039461 0.3876706 0.8816409
bIntercept -1.0808780 -2.3295953 0.1678392 3.4773626
bLength 0.6677614 0.2681462 1.0673766 9.8869492
bLength2 -0.0425534 -0.2839225 0.1988157 0.4546499
bObserver[1] -2.3082666 -55.9343683 51.3178350 0.1004115
bObserver[2] -0.3257582 -1.5238149 0.8722984 0.7512576
bObserver[3] -0.3072008 -1.4105807 0.7961792 0.7728019
bObserver[4] -2.2142360 -51.7735269 47.3450548 0.1043564
bObserver[5] -0.3711566 -2.0796860 1.3373728 0.5771860
bRkm -0.0126214 -0.2351652 0.2099224 0.1336998
bSinuosity -0.0094249 -0.1826084 0.1637587 0.1280685
bSystem -0.0278483 -0.5354908 0.4797943 0.1291396
Bayesian

Table 4. Final parameter estimates.

term estimate lower upper svalue
bIntercept -1.482489 -1.8471861 -1.130028 10.55171
bLength 0.695636 0.3598602 1.076327 10.55171

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

System Efficiency EfficiencySD EfficiencyLower EfficiencyUpper
Kaslo River 0.146685 0.0280345 0.0980146 0.2080916

Table 6. Final model summary.

n K nchains niters nthin ess rhat converged
220 2 3 500 1 730 1.003 TRUE

Lineal Density

Table 7. 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 8. Model averaged parameter estimates and Akaike weights.

term estimate lower upper svalue
bGradient 0.0481677 -0.0070310 0.1033664 3.519391
bRkm 0.4821795 0.4803001 0.4840588 Inf
bSinuosity 0.0522280 -0.0014249 0.1058808 4.148099
bSystemYear[1,1] -2.6189795 -2.6319910 -2.6059679 Inf
bSystemYear[2,1] -1.1819174 -1.2211456 -1.1426892 Inf
bSystemYear[1,2] -2.8299663 -2.8359006 -2.8240320 Inf
bSystemYear[2,2] -1.5532761 -1.6287772 -1.4777751 Inf
bSystemYear[1,3] -2.9652543 -2.9733626 -2.9571459 Inf
bSystemYear[2,3] -1.1165560 -1.2021444 -1.0309676 476.602014
bSystemYear[1,4] -2.3892524 -2.4011550 -2.3773498 Inf
bSystemYear[2,4] -0.5613352 -0.6416158 -0.4810546 139.586467
bSystemYear[1,5] -2.8324690 -2.8413748 -2.8235633 Inf
bSystemYear[2,5] -2.6859777 -2.7675664 -2.6043891 Inf
bSystemYear[1,6] -2.1253187 -2.1336582 -2.1169792 Inf
bSystemYear[2,6] -1.3239973 -1.3974273 -1.2505673 906.352094
bSystemYear[1,7] -2.5239050 -2.5336409 -2.5141691 Inf
bSystemYear[2,7] -1.6582903 -1.7308275 -1.5857532 Inf
bSystemYear[1,8] -2.6064297 -2.6154131 -2.5974463 Inf
bSystemYear[2,8] -1.7108347 -1.7730932 -1.6485762 Inf
bSystemYear[1,9] -2.5064879 -2.5150875 -2.4978884 Inf
bSystemYear[2,9] -2.5000000 -2.5000000 -2.5000000 Inf
log_sDispersion -0.6680712 -0.6688614 -0.6672810 Inf
log_sSite -0.8655976 -0.8783731 -0.8528222 Inf
Bayesian

Table 9. Parameter estimates.

term estimate lower upper svalue
bEfficiency 0.1433574 0.0776347 0.1994434 10.5517083
bRkm 0.4854707 0.3825696 0.5918251 10.5517083
bSystemYear[1,1] -2.5895614 -3.0108882 -1.9734800 10.5517083
bSystemYear[2,1] -1.0899258 -1.7006902 -0.3491926 7.7443533
bSystemYear[1,2] -2.8079909 -3.2031556 -2.1574889 10.5517083
bSystemYear[2,2] -1.4344209 -2.2396465 -0.5601481 8.2297802
bSystemYear[1,3] -2.9489954 -3.3772655 -2.3043877 10.5517083
bSystemYear[2,3] -0.9883088 -1.8059546 -0.0885136 5.0598552
bSystemYear[1,4] -2.3473965 -2.7798670 -1.6936240 10.5517083
bSystemYear[2,4] -0.4377502 -1.1189441 0.3850282 1.9114633
bSystemYear[1,5] -2.8119805 -3.2278659 -2.1682952 10.5517083
bSystemYear[2,5] -2.6248922 -3.8121141 -1.5371765 10.5517083
bSystemYear[1,6] -2.1045876 -2.4989584 -1.4775135 10.5517083
bSystemYear[2,6] -1.2136227 -1.8452623 -0.4608058 10.5517083
bSystemYear[1,7] -2.4784520 -2.9054474 -1.8367602 10.5517083
bSystemYear[2,7] -1.5623849 -2.1107687 -0.8570326 10.5517083
bSystemYear[1,8] -2.5725677 -3.0121570 -1.9277275 10.5517083
bSystemYear[2,8] -1.6182459 -2.2834911 -0.8053585 10.5517083
bSystemYear[1,9] -2.4768057 -2.8845415 -1.8388689 10.5517083
bSystemYear[2,9] -0.0126652 -9.8285520 9.2786919 0.0057785
log_sDispersion -0.5879317 -0.7658651 -0.4217460 10.5517083
log_sSite -0.8565095 -1.1981595 -0.6182791 10.5517083

Table 10. Model summary.

n K nchains niters nthin ess rhat converged
1101 22 3 500 20 110 1.008 FALSE

Stock-Recruiment

Table 11. 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 12. Final parameter estimates.

term estimate lower upper svalue
bAlpha 100.4818607 32.3831554 15569.5777983 10.5517083
bBeta[1] 0.4610103 0.0718305 91.5513575 10.5517083
bBeta[2] 0.2240096 0.0027319 53.6633599 10.5517083
bCarryingCapacity[1] 218.8600050 143.6379897 459.8148416 10.5517083
bCarryingCapacity[2] 441.8676218 229.9165126 10995.6324474 10.5517083
log_bAlpha 4.6099745 3.4776377 9.6530696 10.5517083
log_bBeta[1] -0.7743362 -2.6336163 4.5168063 0.9611212
log_bBeta[2] -1.4960666 -5.9030147 3.9825656 1.5830415
log_sRecruits -0.8356215 -1.1613737 -0.3544944 8.2297802
sRecruits 0.4336049 0.3130558 0.7015290 10.5517083

Table 13. Final model summary.

n K nchains niters nthin ess rhat converged
17 10 3 500 100 162 1.022 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.

Length Correction

figures/length/corrected.png
Figure 7. Corrected length-frequency histogram by year and observation type.

Fish

figures/fish/capture.png
Figure 8. Length-frequency plot of marked Bull Trout by year and system.
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 and tag color.
figures/observer/jmb/length.png
Figure 11. Estimated observer efficiency by fork length and system (with 95% CIs).

Lineal Density

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

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

  • Observer efficiency is approximately 17% for age-1 Bull Trout.
  • Age-1 Bull Trout are relatively evenly distributed with respect to mesohabitat.
  • Lineal densities of age-1 Bull Trout increase with river kilometre in both systems.
  • The age-1 carrying capacity is estimated to be around 220 fish per km in Kaslo River and around 440 fish per km in Keen Creek.

Acknowledgements

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

  • Habitat Conservation Trust Foundation
  • Ministry of Forest, Lands and Natural Resource Operations
    • Greg Andrusak
  • Ministry of Environment
    • Jen Sarchuk
  • BC Fish and Wildlife
    • Trina Radford
  • Stephan Himmer
  • Gillian Sanders
  • Jeff Berdusco
  • Vicky Lipinski
  • Jimmy Robbins
  • Jason Bowers

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.
Greenland, Sander. 2019. “Valid p -Values Behave Exactly as They Should: Some Misleading Criticisms of p -Values and Their Resolution With s -Values.” The American Statistician 73 (sup1): 106–14. https://doi.org/10.1080/00031305.2018.1529625.
Greenland, Sander, and Charles Poole. 2013. “Living with p Values: Resurrecting a Bayesian Perspective on Frequentist Statistics.” Epidemiology 24 (1): 62–68. https://doi.org/10.1097/EDE.0b013e3182785741.
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.
Lin, J. 1991. “Divergence Measures Based on the Shannon Entropy.” IEEE Transactions on Information Theory 37 (1): 145–51. https://doi.org/10.1109/18.61115.
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/.