1. Model Fit, Selection, and Diagnostics

Deviance

For model with likelihood \(f(y|\theta)\), deviance is given by

\[ D(\theta) = -2\log(f(y|\theta)) \]

When comparing 2 models with deviance, the model with smaller deviance is better, in that it better fits the data. An alternate method for calculating deviance is

\[ D_s(\theta) = -2(\log(f(y|\theta)) - \log(f(y|\theta_s))) \]

where \(f(y|\theta_s)\) is the saturated model, that is a model that perfectly fits the data. A few examples of saturated deviance are

\[\begin{split} \begin{align*} y_i \sim Bin(n_i,\theta_i): & \space \space D_s(\underset{\sim}{\theta}) = 2 \sum_i \left ( y_i \log \frac{y_i/n_i}{\theta_i} + (n_i - y_i)\log \frac{1-y_i/n_i}{1-\theta_i} \right ) \\ y_i \sim Poi(\theta_i): & \space \space D_s(\underset{\sim}{\theta}) = 2 \sum_i \left ( y_i \log \frac{y_i}{\theta_i} - (y_i - \theta_i) \right ) \\ y_i \sim N(\theta_i,\delta_i^2): & \space \space D_s(\underset{\sim}{\theta}) = \sum_i \left ( \frac{y_i - \theta_i}{\delta_i} \right ) ^2 \\ \end{align*} \end{split}\]

Deviance example

Consider 2 models, Weibull and Exponential, for the data

\[ y = (1,1,2,2,3,4,4,5,5,8) \]

with non-informative priors. Which model has the smaller deviance?

import pymc as pm
import arviz as az
import numpy as np

rng = np.random.default_rng(6420)

y = np.array([1, 1, 2, 2, 3, 4, 4, 5, 5, 8])

# Weibull Model
with pm.Model() as m1:
    # priors
    alpha = pm.Exponential("alpha", 1)
    gamma = pm.Gamma("gamma", 0.001, 0.001)

    # convert Weibull parameterization
    beta = gamma ** (-1 / alpha)

    # likelihood
    like = pm.Weibull("like", alpha, beta, observed=y)

    # log_like = pm.math.log(
    #    alpha * y ** (alpha - 1) * pm.math.exp(-((y / beta) ** alpha)) / (beta**alpha)
    # )
    log_like = pm.logp(pm.Weibull.dist(alpha=alpha, beta=beta), y)
    deviance = pm.Deterministic("deviance", -2 * log_like.sum())

    trace1 = pm.sample(target_accept=0.9, random_seed=rng)

az.summary(trace1, var_names=["deviance"])

NUTS[nutpie]: [alpha, gamma]

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
deviance	42.8	1.8	41	47	1437	1626	1.00	0.051	0.062

# Exponential Model
with pm.Model() as m2:
    # prior
    lambda_ = pm.Gamma("lambda_", 0.001, 0.001)

    # likelihood
    like = pm.Exponential("like", lambda_, observed=y)

    # log_like = pm.math.log(lambda_ * pm.math.exp(-lambda_ * y))
    log_like = pm.logp(pm.Exponential.dist(lam=lambda_), y)
    deviance = pm.Deterministic("deviance", -2 * log_like.sum())

    trace2 = pm.sample(random_seed=rng)

az.summary(trace2, var_names="deviance")

NUTS[nutpie]: [lambda_]

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
deviance	46.08	1.4	45	49	1816	2031	1.00	0.034	0.057

We see the Weibull model has a lower deviance, and therefore fits the data better than the exponential model. This is not too surprising, given the extra parameter.

Extra learning for PyMC

Deviance measures model fit, but doesn’t tell us how well a model predicts new data. A more modern approach is Pareto-smoothed importance-sampling leave-one-out cross-validation (PSIS-LOO) (Vehtari, Gelman, and Gabry [2017]), which estimates predictive accuracy using posterior samples. ArviZ provides az.loo() and az.compare(). Higher ELPD is better, but the metric gets complicated by some diagnostics that we’ll talk about later.

See the PyMC model comparison guide and Vehtari, Simpson, Gelman, Yao, and Gabry [2024]. The CPO discussed later in lectures is also based on leave-one-out prediction.

# need log-likelihoods for PSIS-LOO
with m1:
    pm.compute_log_likelihood(trace1)

with m2:
    pm.compute_log_likelihood(trace2)

az.compare(
    {
        "Weibull": trace1,
        "Exponential": trace2,
    }
)

/Users/aaron/miniforge3/envs/pymc_f26/lib/python3.14/site-packages/arviz_stats/loo/loo_helper.py:1146: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.70 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

	rank	elpd_diff	dse	p_worse	diag_diff	diag_elpd	p	elpd	se	weight
Weibull	0	0.0	0.0	NaN		1 k̂ > 0.70	1.7	-22.0	2.3	0.91
Exponential	1	-1.0	1.4	0.73	N < 100		0.4	-23.0	1.9	0.09