ASM Chapter 4

code
julia
statistics
asm
ecology
Author

Alun ap Rhisiart

Published

October 12, 2026

We turn to Chapter 4 of Applied Statistical Modelling for Ecologists, by Marc Kéry and Kenneth Kellner. In this chapter we turn to a comparison of different statistical modelling approaches for ecological data. The book considers jagsUI, NIMBLE, TMB, Stan, and Maximum Likelihood. I add Turing.jl to the mix, with a brief discussion of other Julia packages. This discussion of Maximum Likelihood Estimation in this chapter is the clearest I have encountered.

First we pull in any packages we will need.

using DataFramesMeta, Distributions, StatsBase, StatsPlots, RCall, MCMCChains

Next, let’s make sure we have the exact same random variables as the book, for direct comparison. To ensure this, we actually call the R code itself, rather than using any of Julia’s RNGs.

R"""
set.seed(4)
y10 <- rnorm(n = 10, mean = 600, sd = 30) # Sample of 10 birds
y1000 <- rnorm(n = 1000, mean = 600, sd = 30) 
"""
@rget y10 y1000

truth = (mean=600, sd=30)
first(y1000, 10)
10-element Vector{Float64}:
 616.9981349454094
 600.4715836201372
 611.4917201555145
 598.6458865226008
 601.0305572219098
 605.0708032265492
 634.9508051708302
 598.6738800824373
 596.9889467224228
 591.4966629337923

Next, eyeball the data, which will check that it looks the same as in the book.

p1 = histogram(y10,
    color=:grey,
    xlims=(500, 700),
    bins=5,
    xlabel="y10",
    ylabel="Frequency",
    legend=false
)

p2 = histogram(y1000,
    color=:grey,
    bins=10,
    xlabel="y1000",
    legend=false
)

# Combine in a 1x2 layout
plot(p1, p2, layout=(1, 2), size=(800, 400))

For the ‘canned’ example, the book uses lm(). We can use GLM.jl for this. It takes an object that implements the Tables.jl inteface, not a vector, so we just wrap it in a DataFrame.

using GLM

y10lm = lm(@formula(y10 ~ 1), DataFrame(y10=y10))
y1000lm = lm(@formula(y1000 ~ 1), DataFrame(y1000=y1000))
y1000lm
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}

y1000 ~ 1

Coefficients:
────────────────────────────────────────────────────────────────────────
               Coef.  Std. Error       t  Pr(>|t|)  Lower 95%  Upper 95%
────────────────────────────────────────────────────────────────────────
(Intercept)  598.846    0.917595  652.63    <1e-99    597.045    600.647
────────────────────────────────────────────────────────────────────────

Maximum Likelihood Estimation

To keep in line with the book we will also do this manually, using idiomatic Julia where appropriate.

using Random, Optim, ForwardDiff, ADTypes

# Optim requires arrays of floats, not integers
starts = [550.0, 25.0] 

# Objective Function (Negative Log-Likelihood)
function NLL(param, y)
    pop_mean = param[1]
    pop_sd = param[2]
    
    # Catch negative standard deviations so the optimizer doesn't crash
    if pop_sd <= 0.0
        return Inf
    end
    
    # Equivalent to -sum(dnorm(y, mean=pop_mean, sd=pop_sd, log=TRUE))
    return -loglikelihood(Normal(pop_mean, pop_sd), y)
end

# Optimize
# We use an anonymous function `p -> NLL(p, y1000)` to pass the data
# autodiff = AutoForwardDiff() calculates the gradients automatically for faster convergence
fm = optimize(p -> NLL(p, y1000), starts, LBFGS(), autodiff = AutoForwardDiff())

# Extract the maximum likelihood estimates
mle_estimates = Optim.minimizer(fm)
println("MLE Estimates (Mean, SD): ", mle_estimates)

# Calculate the Hessian at the minimum (equivalent to R's hessian=TRUE)
H = ForwardDiff.hessian(p -> NLL(p, y1000), mle_estimates)
println("Hessian Matrix:\n", H)
MLE Estimates (Mean, SD): [598.8460950752325, 29.002390678087444]
Hessian Matrix:
[1.188864620730756 6.827315404646025e-10; 6.827315179674075e-10 2.377729241796953]

The book now checks convergence, inverts the ‘Hessian’, and gets the standard errors for the parameter estimates. We can do the same in Julia as follows:

using LinearAlgebra

# Check convergence
converged = Optim.converged(fm)
println("Converged: ", converged)

# Invert the Hessian to get the covariance matrix
cov_matrix = inv(H)
println("Covariance Matrix:\n", cov_matrix)

# Standard errors for the parameter estimates
std_errors = sqrt.(diag(cov_matrix))
println("Standard Errors (Mean, SD): ", std_errors)
Converged: true
Covariance Matrix:
[0.8411386650444126 -2.415211481758575e-10; -2.415211402173144e-10 0.420569332462874]
Standard Errors (Mean, SD): [0.9171361213279153, 0.6485131706163523]

The Turing.jl model

This is the main point about using Julia for Bayesian statistics. Let’s build the same model in Turing.jl. Also, with the latest version of Turing, the default chain type has changed, so you need explicitly set it to Chains in the sampling call to maintain compatibility with the MCMCDiagnosticsTools functions, which are not yet updated to FlexiChains.VNChain format. Here I will just use the diagnostics available in MCMCChains itself.

using Turing

@model function model_4_4(mass)
    # Priors
    # Note: BUGS dnorm(0, 0.0001) uses precision τ = 0.0001, 
    # which corresponds to variance 10,000 or SD = 100 in Turing.
    pop_mean ~ Uniform(0, 5000)
    pop_sd ~ Uniform(1, 100)

    # Add derived quantity for variance if required
    pop_var := pop_sd^2

    # Likelihood (Vectorized)
    mass .~ Normal(pop_mean, pop_sd)
end

# Sampling
model = model_4_4(y1000)
chains = sample(model, NUTS(), MCMCThreads(), 10000, 4; chain_type=Chains)

Keep in mind that this is a HMC/NUTS model, so equivalent to the Stan example. To me, at least, the model code is far more readable, and the steps needed to compile and run the model also short, simple, and easy to understand. More so than NIMBLE: no configure/build/compile/compile/run sequence, just define the model and sample. Yes, there is a single function to run NIMBLE, too, but then you can’t run in parallel. In Turing, the sample() function can take MCMCThreads() or MCMCDistributed(). And the model itself has no declaring of types, no defined sections for params and data, and the code is Julia, so it has the full range of Julia functions available. Also, there is no likelihood loop, because it can do broadcasting as normal, so it is even shorter than the JAGS model code.

Note

One thing to note is that, as it is Julia code, the order matters, unlike JAGS. You have to define variables before you can refer to them.

Now we can examine the posterior distribution of the parameters.

describe(chains)
Chains MCMC chain (10000×17×4 Array{Float64, 3}):

Iterations        = 1001:1:11000
Number of chains  = 4
Samples per chain = 10000
Wall duration     = 3.85 seconds
Compute duration  = 2.21 seconds
parameters        = pop_mean, pop_sd, pop_var
internals         = n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size, logprior, loglikelihood, logjoint

Summary Statistics

  parameters       mean       std      mcse     ess_bulk     ess_tail      rha ⋯
      Symbol    Float64   Float64   Float64      Float64      Float64   Float6 ⋯

    pop_mean   598.8509    0.9156    0.0051   32136.6764   29319.9822    1.000 ⋯
      pop_sd    29.0501    0.6473    0.0049   17147.4235   19526.1873    1.000 ⋯
     pop_var   844.3279   37.6630    0.2880   17147.4235   19526.1873    1.000 ⋯

                                                               2 columns omitted

Quantiles

  parameters       2.5%      25.0%      50.0%      75.0%      97.5% 
      Symbol    Float64    Float64    Float64    Float64    Float64 

    pop_mean   597.0497   598.2408   598.8497   599.4681   600.6542
      pop_sd    27.8258    28.6133    29.0415    29.4722    30.3635
     pop_var   774.2741   818.7192   843.4098   868.6116   921.9396

Trace plots

Getting trace and density plots from the chains is very straightforward.

plot(chains; size=(800,600))

There are lots of other plotting options, shown here: https://turinglang.org/MCMCChains.jl/stable/statsplots/

Diagnostics any of the following (and more):

  • rstar(chains) # updated from gelmandiag
  • ess(chains) # effective sample size for each parameter
  • rhat(chains) # potential scale reduction factor for each parameter
  • mcse(chains) # Monte Carlo standard error for each parameter

and more, or the Highest Posterior Density (HPD) intervals.

hpd(chains)
HPD

  parameters      lower      upper 
      Symbol    Float64    Float64 

    pop_mean   597.1196   600.7176
      pop_sd    27.8173    30.3535
     pop_var   771.2632   918.7832

Finally, we can do Prior Predictive Checks and Posterior Predictive Checks. But as the chapter doesn’t mention them we will leave to a later model.

Logical tests

For jagsUI (section 4.4.2) the book shows how to do straightforward logical tests, such as mean(sims[,"pop.mean"] < 597). We can do a similar thing with Turing.

The direct vector extraction way

# Extract flattened vector of all posterior draws (e.g. 10,000 x 4 = 40,000 draws)
pop_mean_draws = vec(Array(chains[:pop_mean]))

# Calculate posterior probability P(pop_mean < 597)
prob = mean(pop_mean_draws .< 597)
0.021975

As a one-liner

# Direct one-liner
prob = mean(chains[:pop_mean] .< 597)
0.021975

Other options besides Turing

The Julia ecosystem is a lot less mature than R. Turing itself is mature and stable, but other packages for Bayesian modeling may not be as fully developed or as user-friendly as their R counterparts.

TuringGLM

I tried JuliaBUGS and TuringGLM as well, and could get neither to work. TuringGLM has just been archived, end of August 2026, in favour of newer packages. It would be the equivalent of brms. It is being replaced by TuringRegressions.jl and/or BayesianRegressionModels.jl, which I have not yet tried. There is also Soss.jl which seems to be a Julia interface to the Python library, so will not be considered here. I will try the other two packages and if successful, I will report my findings here.

JuliaBUGS

If TuringGLM is the equivalent of brms, then JuliaBUGS is the equivalent of rjags or jagsui, providing an interface to the JAGS modelling language from within Julia. It does not call the R programs, but runs on top of Turing.jl. This sounds great for someone coming from R. Even better, it refers to a subpackage in R which would allow R users to run their model in Julia. However, the documentation is messy and I could not get either this simple model to work nor the examples in the documentation. Therefore, I will not be considering JuliaBUGS further in this series.

Gibbs Sampler

It is a shame that there is no direct equivalent to a Gibbs sampler in Julia, as they do make handling of discrete variables such as N-Mixture or occupancy models much easier. There is a Gibbs sampler, but it seems to be for combining different samplers within a single framework. However, there is an argument to be made that marginalisation is the best approach.

I note that Turing does support Particle Gibbs, but I don’t think that is quite the same use case. They provide AbstractMCMC as a common interface for MCMC samplers, which allows for flexibility in choosing different sampling algorithms, and it would seem not too difficult to have a Gibbs sampler for common ecological use cases to be built using it. There is, after all, a Metropolis-Hastings sampler. I’ll look into that possibility.