using RCall, DataFrames, StatsPlots
n1 = 60
n2 = 40
n = n1 + n2
mu1 = 105.0
mu2 = 77.5
sigma = 2.75
@rput mu1 mu2 sigma n1 n2
R"""
set.seed(61)
y1 <- rnorm(n1, mu1, sigma)
y2 <- rnorm(n2, mu2, sigma)
y <- c(y1, y2)
"""
x = vcat(zeros(Int, n1), ones(Int, n2))
@rget yChapter 6 of Applied Statistical Modelling for Ecologists, by Marc Kéry and Kenneth Kellner, considers the case of two predictors in a linear regression model. Traditionally this would be a t-test if variance is equal, or the Welch test if they are unequal variances.
6.2 Equal Variance
The simulated data in this chapter is wing lengths of Peregrine falcons and tiercels. In the first section we assume (and thus simulate) that the data has equal variance. In the previous chapter x was continuous, but here it is a categorical dummy variable (‘maleness’).
The base model is the same:
\(y_i = \alpha + \beta x_i + \epsilon_i\)
\(\epsilon_i \sim \mathcal{N}(0, \sigma^2)\)
Since x is categorical, we can split the linear part into group-specific means:
\(y_i = \begin{cases} \alpha + \beta\cdot 0 = \alpha & \text{if } x_i = 0 \\ \alpha + \beta\cdot 1 = \alpha + \beta & \text{if } x_i = 1 \end{cases}\)
‘Slope’ is thus the difference between the group means. Therefore this is effects parameterization.
6.2.1 Data Generation
As with previous chapters, I drop into R to build the simulation dataset to exactly match the book. Otherwise I would generate it directly in Julia.
As Kéry and Kellner describe, the way we simulate the data here corresponds to means parameterization. The variables y1 and y2 are built with their own means, rather than the mean for y1 and the difference for y2.
An aside on rep()
In R, the rep() function is used to replicate the values in a vector. For example, rep(c(0,1), c(n1, n2)) creates a vector of n1 zeros followed by n2 ones. This is equivalent to the vcat(zeros(Int, n1), ones(Int, n2)) approach we used in Julia. There are different ways of reproducing rep()’s functionality in Julia, depending on the use case.
In the case we have here, where we are replicating 0s and 1s, the approach using vcat(zeros(Int, n1), ones(Int, n2)) is straightforward and efficient. For more general cases, we can use fill, like this:
```{julia}
n1 = 3
n2 = 2
x = vcat(fill(0, n1), fill(1, n2))
# Output: [0, 0, 0, 1, 1]
x_group = vcat(fill("Control", n1), fill("Treatment", n2))
# Output: ["Control", "Control", "Control", "Treatment", "Treatment"]
```To replicate rep’s ability to take a vector of values and repeat each according to a vector of counts, we can do the following:
```{julia}
# General R-style rep function in 1 line
rep(vals, counts) = reduce(vcat, fill.(vals, counts))
# Usage:
x = rep([0, 1], [n1, n2])
# Output: [0, 0, 0, 1, 1]
x_multi = rep([10, 20, 30], [2, 1, 3])
# Output: [10, 10, 20, 30, 30, 30]
```Here are some other cases:
| R Syntax | Julia |
|---|---|
rep(x, each = n) |
repeat(x, inner = n) |
rep(x, times = n) |
repeat(x, n) |
rep(x, each = a, times = b) |
repeat(x, inner = a, outer = b) |
rep(c(a,b), c(n1, n2)) |
vcat(fill(a, n1), fill(b, n2)) |
rep(x, length.out = n) |
repeat(x, ceil(Int, n / length(x)))[1:n] |
6.2.2 Likelihood analysis with canned functions
First, we use GLM.jl to fit a linear model to the simulated data.
using GLM
y_obs = DataFrame(y = y, x = x)
model = lm(@formula(y ~ x), y_obs)
modelStatsModels.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) 104.771 0.357902 292.74 <1e-99 104.06 105.481
x -27.6184 0.565893 -48.80 <1e-69 -28.7414 -26.4954
─────────────────────────────────────────────────────────────────────────The book then uses the canned t.test function. In Julia, HypothesisTests.jl provides similar functionality with the EqualVarianceTTest and UnequalVarianceTTest functions. However, the output is different in many ways; the value of t and the degrees of freedom, etc are the same, but it does not show the means of each group directly. If I really wanted it to show exactly the same I would drop into R and use t.test.
6.2.5 Bayesian analysis
Since Turing is an HMC sampler we will refer to the Stan section for model code comparison, but the output is shown for JAGS in the book. As before, we use the Stan code’s prior of half-Cauchy for the standard deviation. Turing uses 1000 burn-in, so 2000 saved with 4 chains should match the number saved in JAGS.
using Turing, FlexiChains, MCMCDiagnosticTools
@model function model6_2_5(y, x)
α ~ Normal(0, 100)
β ~ Normal(0, 100)
σ ~ truncated(Cauchy(0, 10); lower=0)
μ = @. α + β * x
y ~ arraydist(Normal.(μ, σ))
# Return what you want to track, in this case the residuals. This is for method 2, below
return y .- μ
end
model_turing25 = model6_2_5(y_obs.y, y_obs.x)
chains = sample(model_turing25, NUTS(), MCMCThreads(), 2000, 4; 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 … q5 q50 q95 │ │ α 104.7683 0.3600 0.0054 104.1758 104.7717 105.3558 │ │ β -27.6198 0.5768 0.0089 -28.5485 -27.6240 -26.6754 │ │ σ 2.8108 0.2039 0.0024 2.4981 2.7986 3.1750 │ ╰──────────────────────────────────────────────────────────────────────────────╯
An almost exact match with the JAGS output, but again in Turing we don’t put derived quantities into the model.
Means and Residuals
The book shows means and residuals. First the means, which are straightforward. These are extracted as column vectors (S x 1).
using Statistics
α_post = vec(chains[:α])
β_post = vec(chains[:β])
println("μ1: ", mean(α_post))
println("μ2: ", mean(α_post .+ β_post))μ1: 104.76828999825054
μ2: 77.1484998431667Then we calculate the residuals. There are two ways to do this. The first is quick, straightforward, and appropriate here, and uses matrix broadcasting; no need to add any tracked quantities to the model (no ‘return’). While we are about it, we also get the predictions.
# Calculate μ for all samples and all data points
# Data vectors are transposed to row vectors (1 x N)
# The result is an S x N matrix of predictions
μ_post = α_post .+ β_post .* y_obs.x'
# Calculate the residuals (S x N matrix)
# Each row is a posterior sample, each column is a data point
residuals_matrix = y_obs.y' .- μ_post
# If you want the mean residual for each observation:
mean_residuals = vec(mean(residuals_matrix, dims=1))
# Or the 95% credible intervals for the residuals:
using StatsBase
lower_ci = vec(mapslices(x -> percentile(x, 2.5), residuals_matrix, dims=1))
upper_ci = vec(mapslices(x -> percentile(x, 97.5), residuals_matrix, dims=1))
print("Mean residuals: ", first(mean_residuals, 5))Mean residuals: [-0.8150204675676441, -0.799344617312807, -4.501247947095701, 1.1955242844650562, -3.5744552334882465]The second way, the more ‘Turing’ way, is using generated_quantities. This is less transparent, but much less complicated if we are dealing with hierarchical models. The return calculation we added to the model allows us to track, but they are not calculated during the run (so no slowdown), they are calculated when we call generated_quantities afterwards.
# This returns a matrix of arrays containing the residuals
gen_residuals = generated_quantities(model_turing25, chains)
# Because gen_residuals is an array of arrays (Chain length × Chains),
# you will usually want to flatten/format it.
# Here is a quick way to stack it into a standard S x N matrix:
residuals_matrix = reduce(vcat, transpose.(vec(gen_residuals)))
mean_residuals = vec(mean(residuals_matrix, dims=1))
print("Mean residuals: ", first(mean_residuals, 5 ))Mean residuals: [-0.8150204675676441, -0.799344617312807, -4.501247947095701, 1.1955242844650562, -3.5744552334882465]Trace plots
Then we can check for convergence.
plot(chains; size = (800, 600))
Then the residuals plots (not shown in the book). These look fine, no obvious patterns or issues are visible.
using StatsPlots
# Create the scatter plot
p1 = scatter(1:length(mean_residuals), mean_residuals,
color=:black,
legend=false,
xlabel="Index",
ylabel="Mean Residual"
)
# Add the h=0 line (the '!' adds it to the current plot)
hline!(p1, [0], color=:black, linestyle=:dash)
# Create the boxplot
p2 = boxplot(string.(y_obs.x), mean_residuals,
color=:grey,
legend=false,
xlabel="Male",
ylabel="Wingspan residuals (cm)"
)
# Add the h=0 line
hline!(p2, [0], color=:black, linestyle=:dash)
# Combine p1 and p2 side-by-side
# The layout=(1, 2) creates a grid with 1 row and 2 columns
# We also increase the size so the plots don't look squished
combined_plot = plot(p1, p2, layout=(1, 2), size=(800, 400))
# Display the final combined figure
display(combined_plot)
6.2.3 Do-it-yourself maximum likelihood estimates
This is the same as Chapter 5, which I have already covered in detail.
6.3 Comparing two groups with unequal variance
Here we consider that the variance in wingspan may differ between the two groups. This is important because assuming equal variance when it is not true can lead to incorrect inferences. The new model allows each group to have its own variance, typically called \(\sigma_{\text{group}}^2\).
\(y_i = \alpha + \beta x_i + \epsilon_i\)
\(y_i = \begin{cases} \epsilon_i \sim{Normal(0, \sigma_1^2)} & \text{if } x_i = 0 \text{ [females]} \\ \epsilon_i \sim{Normal(0, \sigma_2^2)} & \text{if } x_i = 1 \text{ [males]} \end{cases}\)
So now we have to estimate an additioal parameter for the variance of each group, namely \(\sigma_1^2\) for females and \(\sigma_2^2\) for males.
6.3.2 Data generation
Use exact R code as usual, except I get tired trying to keep in mind what 1 and 2 mean, so I will use descriptive variable names for the groups.
# Number of observations per group
n_females = 60
n_males = 40
@rput n_females n_males
R"""
# True parameters
mu_females <- 105
mu_males <- 77.5
sigma_females <- 3
sigma_males <- 2.5
n <- n_females + n_males
y1 <- rnorm(n_females, mu_females, sigma_females)
y2 <- rnorm(n_males, mu_males, sigma_males)
y <- c(y1, y2)
"""
@rget y mu_females mu_males sigma_females sigma_males
x = vcat(zeros(Int, n_females), ones(Int, n_males))
# Create a DataFrame
y_obs = DataFrame(x=x, y=y)
y1_obs = y_obs.y[y_obs.x .== 0]
y2_obs = y_obs.y[y_obs.x .== 1]Lets see a boxplot
boxplot(string.(y_obs.x), y_obs.y,
color=:grey,
legend=false,
xlabel="Male",
ylabel="Wingspan (cm)"
)
6.3.5 Bayesian analysis
As mentioned, the HypothesisTests.jl functions don’t have quite the same output. I’ll move straight to Bayesian analysis.
@model function model6_3_5(y1, y2)
# priors
μ1 ~ Normal(0, 100)
μ2 ~ Normal(0, 100)
σ1 ~ truncated(Cauchy(0, 10); lower=0)
σ2 ~ truncated(Cauchy(0, 10); lower=0)
# Likelihood
y1 .~ Normal(μ1, σ1) # females
y2 .~ Normal(μ2, σ2) # males
delta := μ2 - μ1
end
model_turing35 = model6_3_5(y1_obs, y2_obs)
chains = sample(model_turing35, NUTS(), MCMCThreads(), 3000, 4; progress=false)summarystats(chains)╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮ │ iter collapsed │ │ chain collapsed │ │ ↓ stat = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95] │ │ │ │ Parameters (5) ── AbstractPPL.VarName │ │ Float64 μ1, μ2, σ1, σ2, delta │ │ │ │ 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 … q5 q50 q95 │ │ μ1 105.4417 0.4189 0.0032 104.7564 105.4395 106.1270 │ │ μ2 77.8711 0.4181 0.0036 77.1862 77.8679 78.5648 │ │ σ1 3.2498 0.3123 0.0025 2.7786 3.2260 3.7963 │ │ σ2 2.6159 0.3074 0.0027 2.1686 2.5860 3.1625 │ │ delta -27.5706 0.5938 0.0049 … -28.5365 -27.5639 -26.5951 │ ╰──────────────────────────────────────────────────────────────────────────────╯
Which matches very well. Note that there is a difference in the way the likelihood is written here compared with previous examples. Previously, we had something like
μ ~ @. α + β * x
y ~ arraydist(Normal.(μ, σ))
Here, μ is a vector of means, and arraydist then constructs the joint likelihood as the product of the individual Normal densities.
In the two groups case, however, we have two separate groups with their own means and standard deviations. Therefore, we can directly write the likelihood for each group using the element-wise broadcasting syntax .~ without needing to construct a vector of means. This simplifies the model and makes it more interpretable for the two-group comparison.
6.3.6 Do-it-yourself maximum likelihood estimation
using Distributions
using Optim
# Define the NLL function
# params is a vector: [mu1, mu2, sigma1, sigma2]
function nll(params, y1, y2)
mu1 = params[1]
mu2 = params[2]
sigma1 = params[3]
sigma2 = params[4]
# Optimizers explore wildly; we must return Infinity if they guess
# a negative or zero standard deviation to reject that path.
if sigma1 <= 0 || sigma2 <= 0
return Inf
end
# logpdf.(distribution, data) calculates the log-likelihood for each data point
# sum() adds them up for the whole vector
ll_y1 = sum(logpdf.(Normal(mu1, sigma1), y1))
ll_y2 = sum(logpdf.(Normal(mu2, sigma2), y2))
# Return the negative log-likelihood
return -(ll_y1 + ll_y2)
endnll (generic function with 1 method)Now we can optimize this negative log-likelihood function to find the maximum likelihood estimates for the parameters. As in previous chapters we need to wrap the parameters in an anonymous function for optim(). Same initial guesses as the book.
# 1. Set your initial guesses for [mu1, mu2, sigma1, sigma2]
initial_guesses = [50.0, 50.0, 1.0, 1.0]
# 2. Create a closure so the optimizer only sees a function of `p` (the parameters)
objective_function = p -> nll(p, y1_obs, y2_obs)
# 3. Run the optimization
# You can explicitly ask for NelderMead(), or let it default
result = optimize(objective_function, initial_guesses, NelderMead())
# 4. Extract your results
println(result)
# Get the exact maximum likelihood estimates
mle_estimates = Optim.minimizer(result)
println("MLE for mu1: ", mle_estimates[1])
println("MLE for mu2: ", mle_estimates[2])
println("MLE for sigma1: ", mle_estimates[3])
println("MLE for sigma2: ", mle_estimates[4]) * Status: success
* Candidate solution
Final objective value: 2.475664e+02
* Found with
Algorithm: Nelder-Mead
* Convergence measures
√(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08
* Work counters
Seconds run: 0 (vs limit Inf)
Iterations: 189
f(x) calls: 334
MLE for mu1: 105.43962547971826
MLE for mu2: 77.87025749126771
MLE for sigma1: 3.1584887401367543
MLE for sigma2: 2.5009815159024686Efficiency 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.
There are two models in this chapter; I just compare the second one. 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 6.3.3 section:
```{r}
# Output metrics for comparison with Turing
compare.params = c("mu1", "mu2", "sigma1", "sigma2")
# 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 = "model6.3.3.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_mcmc6_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 3000 iterations on 4 chains to make it comparable. JAGS saved 12,000 draws (3,000 per chain for 4 chains). By setting Turing to 3000 iterations it also saved 12,000 draws. Adaptation is 1,000 draws in both. This effects overall time, but not ess per second.
include("../../files/model_comparison.jl")
using .ModelComparison
using CSV
params_to_track = ["μ1", "μ2", "σ1", "σ2"]
n_chains = 4
n_iter = 3000
jags_metrics = CSV.read("../../files/jags_mcmc6_metrics.csv", DataFrame)
mcmcresult = compare_model_efficiency(model_turing35, n_chains, params_to_track, jags_metrics, n_iter)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 14931.6 16114.3 0.169898 87885.8 94846.8 2 │ JAGS 4625.52 6255.88 3.496 1323.09 1789.44
The Turing model is much faster, which is as expected. As the model becomes more complex, and more samples are saved, the discrepancy becomes greater, it seems. This isn’t the fastest it could be, because I use the := syntax to store delta; this calculates it during the sampling run and stores the results immediately in the chain, which is convenient. This is a simple calculation so not costly, but better is to use return delta = μ2 - μ1. This is not run during sampling, but only afterwards when you run generate_quantities(chains).
Note
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 6.