ASM Chapter 5

code
julia
statistics
asm
ecology
Author

Dr Alun ap Rhisiart

Published

October 21, 2026

Chapter 5 of Applied Statistical Modelling for Ecologists, by Marc Kéry and Kenneth Kellner, focuses on Ordinary Linear Regression. This chapter also gets in the interesting subject of Posterior Predictive Checks and Goodness Of Fit.

In this chapter, we use a simulated dataset about wallcreepers (which we don’t have in the UK, grrrh). I guess this is very familiar to most people. If the value of our dependent variable \(y_i\) varies according to independent variable \(x_i\), then \[\mu = \alpha + \beta\cdot x\] We take this deterministic, linear part to be the mean, or expected value. \[y_i\sim{Normal(\alpha + \beta\cdot x_i, \sigma^2)}\] We can also express this as a linear equation for the mean plus an error term. \[y_i = \alpha + \beta\cdot x_i + \epsilon_i\] \[\epsilon_i\sim{Normal(0, \sigma^2)}\] So, we have an overall intercept, a coefficient for the effect of x, which will be a regression line, and some scatter around that line which is normally distributed but centered on the line itself.

5.2 Data generation

As usual, the data is simulated, so we know the real value of the parameters, and can check that our model can recover them.

n = 16     # number of years 
var = 25 # residual variance

truth = (α = 40, β = -0.5, σ = sqrt(var))
(α = 40, β = -0.5, σ = 5.0)

Next, we plot the generated data. As in previous chapters, although we can get random numbers from the Normal distribution, they won’t exactly match the book output, because it is a different RNG. So we will just drop into R to get the value of epsilon. I fudged putting the sd into R because R doesn’t know Unicode; I could just use roman letters (eps) for the variable, but it isn’t important here, since I just want to get the same numbers. In reality, I would use the Julia RNG anyway.

using RCall

@rput n
R"set.seed(5);eps <- rnorm(n, 0, sd = 5)"
@rget eps
using StatsPlots, Distributions

ϵ = eps
x = 1:n
y = @. truth.α + truth.β * x + ϵ
scatter(1989 .+ x, y, xlabel="Year", ylabel="Wallcreeper Count", title="Simulated Wallcreeper Data")

5.3 Analysis with Canned Functions in Julia

This section is going to fit a linear regression model to the simulated wallcreeper data using Julia’s built-in functions.

using GLM, StatsBase

data = (; x = x, y = y)
model = lm(@formula(y ~ x), data)
model
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}, Matrix{Float64}}

y ~ 1 + x

Coefficients:
─────────────────────────────────────────────────────────────────────────
                 Coef.  Std. Error      t  Pr(>|t|)  Lower 95%  Upper 95%
─────────────────────────────────────────────────────────────────────────
(Intercept)  40.7426      2.40168   16.96    <1e-10   35.5915   45.8937
x            -0.690722    0.248376  -2.78    0.0147   -1.22344  -0.158009
─────────────────────────────────────────────────────────────────────────

5.3.2 Goodness-of-fit assessment

The book prints out the residuals, which we can do also, but it also gives the code, but not the plot, of a histogram. I’ll display it here, although it isn’t that useful for such a small dataset.

histogram(residuals(model),
    bins = 20,
    xlabel = "Residuals",
    ylabel = "Frequency",
    title = "Histogram of Residuals"
    )

5.3.3 Computing Predictions

Next we get the predictions from the model. In this case there is an identity link function, so the predictions come directly from the linear predictor without any transformation. That is, \(\hat{y} = \alpha + \beta x_i\).

Note

Kéry & Kellner’s R code manually calculates the interval using 1.96 (pred$fit - 1.96 * pred$se.fit). That is an asymptotic approximation (assuming infinite degrees of freedom). Julia’s predict(…, interval=:confidence) uses the exact Student’s t-distribution critical value based on the model’s actual residual degrees of freedom. This means Julia’s frequentist ribbon will be slightly more statistically strictly correct than the manual 1.96 calculation in the ASM book (though practically, with enough data points, the visual difference is negligible).

using DataFrames

# Create a dense grid of x-values for prediction
xpred = range(1, 16, length=100)
newdata = DataFrame(x = xpred)

# Get predictions and 95% Confidence Intervals directly
# GLM returns a DataFrame with columns: :prediction, :lower, and :upper
preds = predict(model, newdata, interval=:confidence, level=0.95)

# Assemble the equivalent of pred.table
pred_table = DataFrame(
    x = newdata.x,
    point_est = preds.prediction,
    LCL = preds.lower,
    UCL = preds.upper
)

first(pred_table, 5)
5×4 DataFrame
Row x point_est LCL UCL
Float64 Float64? Float64? Float64?
1 1.0 40.0519 35.3622 44.7416
2 1.15152 39.9472 35.3261 44.5683
3 1.30303 39.8426 35.2896 44.3955
4 1.45455 39.7379 35.2527 44.2231
5 1.60606 39.6333 35.2154 44.0511

5.4 Bayesian Analysis

We will translate the JAGS version of the linear regression model into a Bayesian analysis using Julia’s Turing package. Now we start to see some major differences between the way you would approach this in BUGS/JAGS and how you do it in Turing.

In the JAGS code, you put everything into the model block, including priors, likelihood, posterior predictive checks, derived residuals, etc. In Turing, this slows down or may even break the NUTS sampler. In JAGS, the model string is a scripting language telling JAGS what nodes to track during its Gibbs sampling loops. In Turing, the @model macro defines a strict probabilistic log-density target for the sampler. Everything else, residuals, fit statistics, predictive datasets, etc, is just standard Julia code applied to the resulting DataFrame or Chain (see 5.4.2 below). Since Julia is fast, we don’t have to convert this into a different language such as C++.

NB the @. macro at the start of a statement just means to broadcast each of the operations element-wise over arrays. You can also, (and you see I do) just add the . broadcast operator for each one.

Note

The Turing philosophy is: fit the parameters in the model, generate predictions outside the model.

using Turing, FlexiChains, MCMCDiagnosticTools

@model function linear_regression(x, y)
  # Priors. tau = 0.0001 is sd = 100
    α ~ Normal(0, 100)
    β ~ Normal(0, 100)
    σ ~ truncated(Cauchy(0, 10); lower=0)

  # Likelihood
    μ = @. α + β * x
    y ~ arraydist(Normal.(μ, σ))
end

n_chains = 4
model = linear_regression(x, y)
chains = sample(model, NUTS(), MCMCThreads(), 10000, n_chains; progress=false)
summarystats(chains)
╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95]          │
│                                                                              │
│ Parameters (3) ── AbstractPPL.VarName                                        │
│  Float64  α, β, σ                                                            │
│                                                                              │
│ Extras (14)                                                                  │
│  Float64  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                                                                      │
│  param  mean     std     mcse    …  rhat    q5       q50      q95            │
│  α      40.6954  2.6411  0.0223     1.0000  36.3268  40.7183  44.9378        │
│  β      -0.6852  0.2728  0.0023     1.0000  -1.1272  -0.6863  -0.2351        │
│  σ      4.9408   1.0223  0.0079     1.0001  3.5819   4.7914   6.8246         │
╰──────────────────────────────────────────────────────────────────────────────╯

These figures closely match the output from JAGS.

Now we can check for convergence.

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

5.4.2 GOF and Posterior Predictive Checks

To replicate the fit and fit.new discrepancy measures from JAGS, we extract the posterior samples and calculate the sums-of-squares in standard Julia. This approach is lightning-fast and mathematically identical to the JAGS output.


# Extract parameter vectors from the chain
α_post = vec(Array(chains[:α]))
β_post  = vec(Array(chains[:β]))
σ_post = vec(Array(chains[:σ]))

S = length(α_post) # Number of posterior draws
fit = zeros(S)         # Sum of squares for actual data
fit_new = zeros(S)     # Sum of squares for new data

# Generate y.new and compute discrepancies for every posterior draw
for s in 1:S
    # Expected values for this iteration
    μ = @. α_post[s] + β_post[s] * x
    
    # Generate replicate dataset
    y_new = rand.(Normal.(μ, σ_post[s]))
    
    # Calculate squared residuals
    fit[s]     = sum((y .- μ).^2)
    fit_new[s] = sum((y_new .- μ).^2)
end

# Calculate the Bayesian p-value (proportion of times fit_new > fit)
bp_val = mean(fit_new .> fit)
println("Bayesian p-value: ", round(bp_val, digits=3))
Bayesian p-value: 0.526

As with the book, we have a value close to 0.5, so a decent fit.

# Replicate the JAGS goodness-of-fit scatter plot
scatter(fit, fit_new, 
    color = (:black), 
    markershape = :circle, 
    markerstrokewidth = 0,
    xlabel = "Discrepancy actual data (fit)",
    ylabel = "Discrepancy replicate data (fit.new)",
    title = "Bayesian p-value = $(round(bp_val, digits=3))",
    legend = false,
    framestyle = :box,
    aspect_ratio = :equal # Keeps the 1:1 line perfectly diagonal
)

# Add 1:1 reference line
min_val = min(minimum(fit), minimum(fit_new))
max_val = max(maximum(fit), maximum(fit_new))
plot!([min_val, max_val], [min_val, max_val], color=:red, linewidth=2)

Finally, we will plot the Credible Intervals (CRI) like figure 5.4 in the book. First obtain predictions and calculate offset distances for the plot ribbon. Note that in R the plots are built using polygon and absolute values. In Julia we use the ribbon keyword and provide the lower and upper offsets relative to the mean prediction.

using Statistics

# Matrix of predictions: rows = length(xpred), cols = number of MCMC draws
predictions = [α_post[s] + β_post[s] * xp for xp in xpred, s in 1:length(α_post)]

# Calculate Bayesian point estimates (post. means) and 95% CRI
BPE = vec(mean(predictions, dims=2))
LPB = [quantile(row, 0.025) for row in eachrow(predictions)]
UPB = [quantile(row, 0.975) for row in eachrow(predictions)]

# Calculate offset distances required by Plots.jl 'ribbon'
BPE_rib_lower = BPE .- LPB
BPE_rib_upper = UPB .- BPE

Now we can recreate figure 5.4

# Base plot: Raw data
p = scatter(x .+ 1989, y,
    xlabel = "Year",
    ylabel = "Prop. occupied (%)",
    ylims = (25, 50),
    color = :black, 
    fillalpha = 0.5,
    markershape = :circle, # pch = 16
    markersize = 6,        # cex = 1.5
    markerstrokewidth = 0,
    framestyle = :axes,    # frame = FALSE
    legend = :bottomleft   # Places legend automatically
)

# CIs from frequentist inference with GLM/lm() (Section 5.3.3)
# Requires calculating offsets: freq_lower_dist = freq_fit .- freq_lower, etc.
freq_fit = pred_table.point_est
freq_lower_dist = freq_fit .- pred_table.LCL
freq_upper_dist = pred_table.UCL .- freq_fit

plot!(p, xpred .+ 1989, freq_fit,
    ribbon = (freq_lower_dist, freq_upper_dist),
    color = :blue,
    fillalpha = 0.1,       # rgb(0, 0, 1, 0.1)
    linewidth = 2,
    label = "Maximum likelihood"
)

# CRIs from Bayesian inference
plot!(p, xpred .+ 1989, BPE,
    ribbon = (BPE_rib_lower, BPE_rib_upper),
    color = :red,
    fillalpha = 0.2,       # rgb(1, 0, 0, 0.2)
    linewidth = 2,
    label = "Posterior inference"
)

display(p)

5.7 Maximum Likelihood Estimation

We will do a more-or-less direct translation of the R code for maximum log-likelihood estimation into Julia.

using Optim

# Define the NLL function
function NLL(params, y, x)
    α = params[1]
    β  = params[2]
    σ = exp(params[3]) # convert from log scale
    n = length(y)

    # Vectorized calculation using Julia's broadcasting dot (.) syntax
    μ = @. α + β * x
    LL = logpdf.(Normal.(μ, σ), y)
    
    return -sum(LL)
end

# Minimize the NLL to find MLEs
inits = [0.0, 0.0, 0.0]

# Wrap NLL in a function for Optim
objective(params) = NLL(params, y, x)

# Run optimization using BFGS
# Autodiff is optional but highly recommended in Julia for exact gradients
out5_7 = optimize(objective, inits, BFGS(), autodiff = AutoForwardDiff())

# Extract parameters
# Optim.minimizer returns the optimal parameters
diy_params = Optim.minimizer(out5_7)

diy_est = [
    diy_params[1],        # alpha
    diy_params[2],        # beta
    diy_params[3]         # log.sigma
]

# To get the Standard Errors (SEs), you can extract the inverse Hessian:
# hessian = ForwardDiff.hessian(objective, diy_params)
# ses = sqrt.(diag(inv(hessian)))
3-element Vector{Float64}:
 40.74259614225698
 -0.6907222556084501
  1.4548947065085656

Once again we have a close match with the R DIY-MLE output

Efficiency Comparison

In the last section of each chapter I compare the efficiency of the HMC/NUTS model in Julia with the JAGs model in R. All run on a MacBook Pro with M5 Pro chip, 24GB RAM.

I ran the jagsUI model in RStudio and saved the MCMC metrics to a CSV file, which we will read into Julia for comparison. Here is the R code, added at the end of the 5.4.1 section:

```{r}
# Output metrics for comparison with Turing
compare.params = c("alpha", "beta", "sigma")

# Time the whole run (wall clock, includes adaptation + burnin)
elapsed <- system.time({
    fit <- jags(data     = dataList,
            inits        = inits,
            parameters.to.save = compare.params,  # keep this list IDENTICAL to Julia
            model.file   = "model5.4.txt",
            n.chains     = nc,
            n.iter       = ni,
            n.burnin     = nb,
            n.thin       = nt,
            parallel     = TRUE)
})[["elapsed"]]

# Per-parameter effective sample size, dropping derived quantities
neff <- fit$summary[, "n.eff"]
neff <- neff[!names(neff) %in% "deviance"]

jags_metrics <- data.frame(sampler  = "JAGS",
                           min_ess  = min(neff),
                           mean_ess = mean(neff),
                           seconds  = elapsed)
jags_metrics$min_ess_s  <- jags_metrics$min_ess  / elapsed
jags_metrics$mean_ess_s <- jags_metrics$mean_ess / elapsed

write.csv(jags_metrics, "jags_mcmc5_metrics.csv", row.names = FALSE)
```

Back in Julia, we read the R output file and compare it with the Turing model metrics. Here, Turing ran 5000 iterations on 4 chains to make it comparable. JAGS saved 20,000 draws (5,000 per chain for 4 chains, plus 1,000 adaptation). By setting Turing to 5000 iterations it also saved 20,000 draws. This effects overall time, but not ess per second.

include("../../files/model_comparison.jl")
using .ModelComparison
using CSV

params_to_track = ["α", "β", "σ"]
jags_metrics = CSV.read("../../files/jags_mcmc5_metrics.csv", DataFrame)
mcmcresult = compare_model_efficiency(model, n_chains, params_to_track, jags_metrics)
println(mcmcresult)
2×6 DataFrame
 Row │ sampler  min_ess   mean_ess  seconds   min_ess_s  mean_ess_s 
     │ String   Float64   Float64   Float64   Float64    Float64    
─────┼──────────────────────────────────────────────────────────────
   1 │ Turing    7964.79   8259.91  0.250552   31789.0     32966.9
   2 │ JAGS     12230.0   17410.0   3.519       3475.42     4947.43

The Turing model is significantly faster, which is as expected. In jagsUI, parallel=TRUE starts up multiple R processes to run the chains in parallel, which adds some overhead. Turing uses threads within a single Julia process, which is more efficient. This effects the wallclock figure. The key metric to focus on is ESS per second, which accounts for both sampling efficiency and computational time. A mean of around 8,000 effective sample size should be enough for estimating fine tail behaviour. ~1000 ESS is probably good enough for general parameter estimation.

In the Julia code I do not use GPU acceleration, but Turing.jl does support it. Since I’m on Apple Silicon, I could use Metal acceleration for certain computations (CUDA is also supported), which might further improve performance. I wouldn’t bother in this case, but if the dataset is > 100,000 or if it is a high-dimensional matrix model (spatial models, GLMMs, etc), I would experiment with hardware acceleration. Matrix multiplications scale significantly faster on GPU cores.

Let’s see a plot.


metrics = ["min ESS/s", "mean ESS/s"]
ymat    = permutedims(hcat(mcmcresult.min_ess_s, mcmcresult.mean_ess_s)) 

groupedbar(metrics, ymat,
           labels     = ["Turing" "JAGS"],  
           ylabel     = "Effective sample size per second",
           title      = "MCMC efficiency: jagsUI vs Turing.jl",
           framestyle = :axes,
           legend     = :topleft)

That concludes Chapter 5. I haven’t tried to reproduce all the DHARMa output nor all of the plots, but we have covered most of it.