Getting familiar with INLA
The Getting started tutorial introduced the surgery model and ran it through all three engines in a few lines each. This tutorial stays with that same model and digs into the INLA result object: how to read the marginals back through the accessor functions, how to use the model-comparison criteria INLA computes along the way, and how to draw joint posterior samples.
The surgery model
We rebuild the model from Getting started unchanged: a Binomial likelihood for the per-hospital death counts, an overall logit-scale intercept β, and a sum-to-zero IID hospital effect u with precision τ_u.
using DataFrames
surg_data = DataFrame(
hospital = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"],
n = [47, 148, 119, 810, 211, 196, 148, 215, 207, 97, 256, 360],
r = [0, 18, 8, 46, 8, 13, 9, 31, 14, 8, 29, 24]
)
using Latte
using Distributions
using GaussianMarkovRandomFields: IIDModel
using StatsFuns: logistic
using LinearAlgebra
@latte function surg_mortality(r, n_trials, hospital_idx, H)
τ_u ~ Gamma(2.0, 1.0)
β ~ MvNormal(zeros(1), 100.0 * I(1))
u ~ IIDModel(H, constraint = :sumtozero)(τ = τ_u)
for i in eachindex(r)
r[i] ~ Binomial(
n_trials[i], logistic(β[1] + u[hospital_idx[i]])
)
end
end
H = length(surg_data.hospital)
hospital_idx = 1:H |> collect
lgm = surg_mortality(surg_data.r, surg_data.n, hospital_idx, H)
inla_result = inla(lgm, surg_data.r)INLAResult:
Model: LatentGaussianModel{HyperparameterSpec{@NamedTuple{τ_u::Hyperparameter{Base.Fix1{typeof(broadcast), typeof(log)}, :natural}}, @NamedTuple{}}, Latte._PatternAugmentedLatentModel{Latte.RoutedLatentModel{CombinedModel{LinearSolve.CHOLMODFactorization{Nothing}}, @NamedTuple{τ_iid::Symbol}}, SparseArrays.SparseMatrixCSC{Int64, Int64}}, LinearlyTransformedObservationModel{Latte.BinomialTrialsObservationModel{ExponentialFamily{Distributions.Binomial, LogitLink, Nothing, Nothing}, Vector{Int64}}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing}}
Hyperparameters: 1
Latent variables: 13
Mode: (τ_u=3.4264)
Convergence: ✓
Total time: 0.01 seconds
Exploration: 4 points (4 integration)
Model comparison metrics:
Deviance Information Criterion (DIC):
DIC: 54.26
Effective parameters (p_D): -0.1
Mean deviance (D̄): 54.36
Deviance at mode: 54.47
Marginal Log-Likelihood:
log p(y): -44.42
Watanabe-Akaike Information Criterion (WAIC):
WAIC: 67.2
Effective parameters (p_WAIC): 3.54
Log pointwise predictive density (lppd): -30.06
Conditional Predictive Ordinates (CPO):
LPML: -101.09
Mean CPO: 0.0027
Min CPO: 0.0
Unreliable observations: 12 / 12
Max failure score: 4.52
PIT computed: 12 values
PIT mean: 0.4976 (ideal: 0.5)
Approximation quality (KLD):
Max SKLD: 0.0493 (variable 1)
Mean SKLD: 0.0057
Use .hyperparameter_marginals, .latent_marginals, .accumulators for analysisMarginals in depth
Every marginal Latte computes implements the Distributions.jl interface, so you can call Distributions.jl methods on these objects directly. The rule of thumb is to read marginals back out through the accessor functions rather than the result's fields: the accessors work the same way regardless of how the result was stored internally.
Each accessor returns a block of marginals. hyperparameter_marginals takes a hyperparameter name and returns the marginals for that named block, here a single distribution for τ_u:
τ_marginal = hyperparameter_marginals(inla_result, :τ_u)[1]SplineMarginalDistribution(mean=3.186, sd=0.9639, support=[1.449, 5.269])Because it is an ordinary Distributions.jl object, the usual methods apply. Read off the median and a 95% credible interval:
median(τ_marginal), quantile(τ_marginal, 0.025), quantile(τ_marginal, 0.975)(3.1073564818454895, 1.5976887018701216, 5.07406012330638)Sample from it:
rand(τ_marginal, 3)3-element Vector{Float64}:
3.5442355470638476
2.0265855920537708
2.1357538599313286Or evaluate its density and CDF at a point:
pdf(τ_marginal, 1.0), cdf(τ_marginal, 1.0)(0.07143329576647237, 0.0)Calling an accessor without a name returns every block at once. summary_df turns a block into a tidy DataFrame of key statistics — the quickest way to scan a whole posterior:
summary_df(hyperparameter_marginals(inla_result))1×6 DataFrame
Row │ mode median q2_5 q97_5 mean std
│ Float64 Float64 Float64 Float64 Float64 Float64
─────┼──────────────────────────────────────────────────────
1 │ 2.81048 3.10736 1.59769 5.07406 3.186 0.96391The same works for the latent field. These are exactly the variables in the model: the intercept β plus the 12 hospital effects, 13 in total. Pass a name to slice out a single latent block:
summary_df(latent_marginals(inla_result, :u))12×6 DataFrame
Row │ mode median q2_5 q97_5 mean std
│ Float64 Float64 Float64 Float64 Float64 Float64
─────┼─────────────────────────────────────────────────────────────────────
1 │ -0.599959 -0.643315 -1.64449 0.178207 -0.666569 0.46145
2 │ 0.478828 0.481841 0.0283837 0.946088 0.483255 0.234046
3 │ -0.0582782 -0.0576395 -0.622114 0.509027 -0.0573584 0.288136
4 │ -0.21602 -0.215365 -0.522049 0.0937247 -0.215053 0.157036
5 │ -0.483302 -0.488005 -1.0342 0.0404137 -0.49033 0.273689
6 │ -0.0659345 -0.0652301 -0.542318 0.414715 -0.0648606 0.243937
7 │ -0.133422 -0.133621 -0.671866 0.403988 -0.133708 0.274092
8 │ 0.698258 0.700337 0.321449 1.08697 0.701353 0.195276
9 │ -0.0491543 -0.0483063 -0.513808 0.420506 -0.0478774 0.238164
10 │ 0.0901248 0.0926068 -0.48309 0.67785 0.0938474 0.295757
11 │ 0.4526 0.454407 0.0776654 0.837972 0.455298 0.193916
12 │ -0.0592475 -0.05841 -0.442781 0.329136 -0.0579978 0.196831We can also ask for the per-hospital linear predictors η = β + u[hospital]. These aren't latent variables we sampled directly, so Latte derives them from the latent posterior on demand:
summary_df(linear_predictor_marginals(inla_result))12×6 DataFrame
Row │ mode median q2_5 q97_5 mean std
│ Float64 Float64 Float64 Float64 Float64 Float64
─────┼────────────────────────────────────────────────────────────
1 │ -3.17803 -3.22566 -4.32341 -2.32433 -3.2511 0.506264
2 │ -2.10304 -2.10186 -2.57024 -1.62903 -2.10128 0.240073
3 │ -2.63967 -2.64116 -3.24966 -2.03822 -2.64189 0.308679
4 │ -2.79927 -2.79948 -3.08592 -2.51383 -2.79958 0.145932
5 │ -3.06496 -3.07159 -3.6605 -2.50766 -3.07486 0.293698
6 │ -2.64828 -2.64901 -3.15262 -2.14831 -2.64939 0.256056
7 │ -2.71522 -2.71724 -3.29514 -2.14702 -2.71824 0.292585
8 │ -1.88395 -1.88343 -2.25751 -1.50739 -1.88318 0.19136
9 │ -2.63152 -2.63211 -3.12144 -2.14504 -2.63241 0.248957
10 │ -2.4912 -2.49086 -3.11103 -1.86933 -2.49068 0.316418
11 │ -2.12991 -2.12945 -2.50337 -1.75384 -2.12923 0.191198
12 │ -2.64213 -2.6424 -3.03068 -2.25512 -2.64253 0.197807And the posterior predictive marginals, which map the linear predictors through the inverse link function onto the mortality-rate scale:
summary_df(observation_marginals(inla_result))12×6 DataFrame
Row │ mode median q2_5 q97_5 mean std
│ Float64 Float64 Float64 Float64 Float64 Float64
─────┼───────────────────────────────────────────────────────────────────
1 │ 0.0322651 0.0382114 0.0130812 0.0891276 0.041518 0.0197637
2 │ 0.104474 0.108916 0.0710787 0.163964 0.111157 0.0238167
3 │ 0.0616866 0.0665357 0.0373392 0.115249 0.0690664 0.0200402
4 │ 0.0563497 0.0573522 0.043692 0.0748944 0.0578572 0.00797633
5 │ 0.0413621 0.0442946 0.0250747 0.0753228 0.0458246 0.0129135
6 │ 0.0626687 0.06605 0.0409881 0.104489 0.0677874 0.0162952
7 │ 0.0578649 0.061964 0.0357384 0.10461 0.0640939 0.0177017
8 │ 0.128858 0.131995 0.0947033 0.181325 0.133563 0.0221697
9 │ 0.0638515 0.0671 0.0422316 0.104796 0.0687663 0.01605
10 │ 0.0707098 0.0765015 0.0426546 0.133619 0.0795203 0.0234072
11 │ 0.103508 0.106267 0.0756224 0.147564 0.107654 0.0184138
12 │ 0.0643967 0.0664593 0.0460589 0.0949084 0.0675058 0.0125073Model-comparison criteria
INLA computes several model-comparison quantities along the way, evaluated at each integration point of the hyperparameter posterior. By default Latte attaches four accumulators: the Deviance Information Criterion (DIC), the marginal likelihood, WAIC, and CPO. Their values appear in the INLAResult summary, and you can also reach them directly through the accumulators field.
The first accumulator carries the DIC and its effective number of parameters:
inla_result.accumulators[1].DIC, inla_result.accumulators[1].p_D(54.257451410310246, -0.10489928634216028)The marginal likelihood has its own accessor:
log_marginal_likelihood(inla_result)-44.423425943588185Quantities like these are the basis of Bayesian model selection: fit competing models to the same data and compare them on DIC or WAIC (lower is better), or on the marginal likelihood (higher is better). The marginal likelihood is also what turns several candidate models into a weighted average. The Bayesian model averaging tutorial works through a comparison of two competing trend models on a shared dataset, built from exactly these quantities.
Posterior sampling
The marginals above are one-dimensional summaries. When you need joint draws — for instance to propagate posterior uncertainty through a downstream function of several latent variables — rand draws from the full approximate posterior. Each draw picks a hyperparameter configuration from the integration grid, then a joint latent field from the inner Gaussian at that configuration.
using Random
Random.seed!(0)
samples = rand(inla_result, 1000);rand(result, n) returns a PosteriorSamples object holding row-aligned matrices: one row per draw, with θ the hyperparameters and x the latent field. Their shapes:
size(samples.θ), size(samples.x)((1000, 1), (1000, 13))From the joint draws we can compute any functional of the posterior. The columns of samples.x follow the model's latent layout, which latent_groups reports as a name-to-column-range map. We use it to locate the hospital block rather than hard-coding offsets:
u_cols = latent_groups(inla_result)[:u]2:13With that, the posterior probability that hospital H (the eighth) has a higher logit-scale effect than hospital A (the first):
mean(samples.x[:, u_cols[8]] .> samples.x[:, u_cols[1]])0.997Passing include_y = true additionally draws posterior-predictive observations, the basis of posterior predictive checks. The posterior predictive checks tutorial walks through that workflow end to end.
Conclusion
That covers the INLA result object in depth: the marginal accessors and their Distributions.jl interface, the DIC / WAIC / marginal-likelihood accumulators for model comparison, and joint posterior sampling. The same accessors work on tmb and hmc_laplace results too, so once you are comfortable here the other engines read the same way.
References
The original INLA paper: deterministic approximate Bayesian inference for latent Gaussian models via nested Laplace approximations and numerical integration over the hyperparameters.
A modern review of the INLA methodology, the SPDE approach, and the R-INLA ecosystem.
This page was generated using Literate.jl.