Compare and interpret models¶

Fitting one model is rarely the end of an analysis — the real question is usually which of several candidate models best explains the data. This page shows the HSSM pattern for that question:

fit each candidate with idata_kwargs=dict(log_likelihood=True) (comparison needs the pointwise log-likelihood),
rank the fits with az.compare, and
read the comparison table — including when it tells you the candidates are indistinguishable.

The worked example — two hierarchical regression models for drift rate — is adapted from the Winterbrain 2025 workshop (archived snapshot). The stored outputs come from a seeded, reproducible run of this notebook; the two hierarchical fits take a while if you re-execute them yourself.

Run this how-to¶

On Colab, uncomment and run the installation cell below once, then restart the runtime. For local setup, GPU extras, and troubleshooting see the Installation guide.

In [1]:

Copied!

# %pip install hssm
# %pip install hssm

Setup¶

In [2]:

Copied!





# Import modules
import arviz as az
import jax
import matplotlib as plt
import pytensor

import hssm

pytensor.config.floatX = "float32"
jax.config.update("jax_enable_x64", False)

plt.use("Agg")
%matplotlib inline

# hssm.set_floatX("float32")
# Import modules
import arviz as az
import jax
import matplotlib as plt
import pytensor

import hssm

pytensor.config.floatX = "float32"
jax.config.update("jax_enable_x64", False)

plt.use("Agg")
%matplotlib inline

# hssm.set_floatX("float32")

In [3]:

Copied!

import numpy as np
import pandas as pd

np.random.seed(20250801)  # reproducible covariates
import numpy as np
import pandas as pd

np.random.seed(20250801)  # reproducible covariates

The data: one outcome, two candidate explanations¶

We simulate a dataset in which each participant's drift rate depends on two (z-scored) brain-measure regressors (zSTN, zGPe) and a symptom-severity score (sevScore), with participants belonging to one of two diagnosis groups. The scientific question: does diagnosis moderate the brain–behavior relationship (a fixed interaction), or is it enough to let effects vary across participants nested within diagnosis groups?

In [4]:

Copied!





# Function to simulate data for one participant
def simulate_participant3(participant_id, sevScore, diagnosis, size=300):
    """Simulate DDM trial data for one participant with diagnosis-dependent drift."""
    # intercept = 0.5
    zSTN = np.random.normal(loc=1, scale=2, size=size)
    zGPe = np.random.normal(loc=1, scale=2, size=size)

    if diagnosis == "PD":
        intercept = 0.3
    else:
        intercept = 0.6

    ## Strength of neural modulation depends on PD severity score
    if sevScore <= 0.2:
        v = intercept + 0.8 * zSTN + 0.3 * zGPe
    elif sevScore > 0.2 and sevScore <= 0.4:
        v = intercept + 0.6 * zSTN + 0.2 * zGPe
    elif sevScore > 0.4 and sevScore < 0.6:
        v = intercept + 0.5 * zSTN + 0.1 * zGPe
    else:
        v = intercept + 0.4 * zSTN + 0.0 * zGPe

    # Assume `hssm.simulate_data` returns a DataFrame
    true_values = np.column_stack(
        [v, np.repeat([[1.5, 0.5, 0.5]], axis=0, repeats=size)]
    )
    dataset_reg_v = hssm.simulate_data(
        model="ddm",
        theta=true_values,
        size=1,
        # Generate one data point for each set of true values
        random_state=int(sevScore * 100) + (0 if diagnosis == "PD" else 1000),
    )

    # Adding additional variables to the dataset
    dataset_reg_v["zSTN"] = zSTN
    dataset_reg_v["zGPe"] = zGPe
    dataset_reg_v["participant_id"] = str(participant_id)
    dataset_reg_v["sevScore"] = sevScore
    dataset_reg_v["diagnosis"] = diagnosis

    return dataset_reg_v


# Simulate data for four participants
### note that we assume that STN & GPe

## PD patients:
# patients with low severity scores
subj_list = [1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8]
sevscore_list = [
    0.10,
    0.20,
    0.30,
    0.40,
    0.50,
    0.60,
    0.70,
    0.80,
    0.10,
    0.20,
    0.30,
    0.40,
    0.50,
    0.60,
    0.70,
    0.80,
]
diagnosis_list = [
    "PD",
    "PD",
    "PD",
    "PD",
    "PD",
    "PD",
    "PD",
    "PD",
    "DD",
    "DD",
    "DD",
    "DD",
    "DD",
    "DD",
    "DD",
    "DD",
]

# [simulate_participant3(subj, sevscore, diagnosis) for subj, sevscore, diagnosis
# in zip(subj_list, sevscore_list, diagnosis_list)]
# dataset2_participant1 = simulate_participant3(1,0.10,"PD")
# dataset2_participant2 = simulate_participant3(2,0.20,"PD")
# dataset2_participant3 = simulate_participant3(3,0.30,"PD")
# dataset2_participant4 = simulate_participant3(4,0.40,"PD")
# # patients with high severity scores
# dataset2_participant5 = simulate_participant3(5,0.50,"PD")
# dataset2_participant6 = simulate_participant3(6,0.60,"PD")
# dataset2_participant7 = simulate_participant3(7,0.70,"PD")
# dataset2_participant8 = simulate_participant3(8,0.80,"PD")

# ## Dystonia patients:
# # patients with low severity scores
# dataset2_participant9 = simulate_participant3(1,0.10,"DD")
# dataset2_participant10 = simulate_participant3(2,0.20,"DD")
# dataset2_participant11 = simulate_participant3(3,0.30,"DD")
# dataset2_participant12 = simulate_participant3(4,0.40,"DD")
# # patients with high severity scores
# dataset2_participant13 = simulate_participant3(5,0.50,"DD")
# dataset2_participant14 = simulate_participant3(6,0.60,"DD")
# dataset2_participant15 = simulate_participant3(7,0.70,"DD")
# dataset2_participant16 = simulate_participant3(8,0.80,"DD")

# Combine datasets into one DataFrame
combined_dataset3 = pd.concat(
    [
        simulate_participant3(subj, sevscore, diagnosis)
        for subj, sevscore, diagnosis in zip(subj_list, sevscore_list, diagnosis_list)
    ],
    ignore_index=True,
)
combined_dataset3
# Function to simulate data for one participant
def simulate_participant3(participant_id, sevScore, diagnosis, size=300):
    """Simulate DDM trial data for one participant with diagnosis-dependent drift."""
    # intercept = 0.5
    zSTN = np.random.normal(loc=1, scale=2, size=size)
    zGPe = np.random.normal(loc=1, scale=2, size=size)

    if diagnosis == "PD":
        intercept = 0.3
    else:
        intercept = 0.6

    ## Strength of neural modulation depends on PD severity score
    if sevScore <= 0.2:
        v = intercept + 0.8 * zSTN + 0.3 * zGPe
    elif sevScore > 0.2 and sevScore <= 0.4:
        v = intercept + 0.6 * zSTN + 0.2 * zGPe
    elif sevScore > 0.4 and sevScore < 0.6:
        v = intercept + 0.5 * zSTN + 0.1 * zGPe
    else:
        v = intercept + 0.4 * zSTN + 0.0 * zGPe

    # Assume `hssm.simulate_data` returns a DataFrame
    true_values = np.column_stack(
        [v, np.repeat([[1.5, 0.5, 0.5]], axis=0, repeats=size)]
    )
    dataset_reg_v = hssm.simulate_data(
        model="ddm",
        theta=true_values,
        size=1,
        # Generate one data point for each set of true values
        random_state=int(sevScore * 100) + (0 if diagnosis == "PD" else 1000),
    )

    # Adding additional variables to the dataset
    dataset_reg_v["zSTN"] = zSTN
    dataset_reg_v["zGPe"] = zGPe
    dataset_reg_v["participant_id"] = str(participant_id)
    dataset_reg_v["sevScore"] = sevScore
    dataset_reg_v["diagnosis"] = diagnosis

    return dataset_reg_v


# Simulate data for four participants
### note that we assume that STN & GPe

## PD patients:
# patients with low severity scores
subj_list = [1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8]
sevscore_list = [
    0.10,
    0.20,
    0.30,
    0.40,
    0.50,
    0.60,
    0.70,
    0.80,
    0.10,
    0.20,
    0.30,
    0.40,
    0.50,
    0.60,
    0.70,
    0.80,
]
diagnosis_list = [
    "PD",
    "PD",
    "PD",
    "PD",
    "PD",
    "PD",
    "PD",
    "PD",
    "DD",
    "DD",
    "DD",
    "DD",
    "DD",
    "DD",
    "DD",
    "DD",
]

# [simulate_participant3(subj, sevscore, diagnosis) for subj, sevscore, diagnosis
# in zip(subj_list, sevscore_list, diagnosis_list)]
# dataset2_participant1 = simulate_participant3(1,0.10,"PD")
# dataset2_participant2 = simulate_participant3(2,0.20,"PD")
# dataset2_participant3 = simulate_participant3(3,0.30,"PD")
# dataset2_participant4 = simulate_participant3(4,0.40,"PD")
# # patients with high severity scores
# dataset2_participant5 = simulate_participant3(5,0.50,"PD")
# dataset2_participant6 = simulate_participant3(6,0.60,"PD")
# dataset2_participant7 = simulate_participant3(7,0.70,"PD")
# dataset2_participant8 = simulate_participant3(8,0.80,"PD")

# ## Dystonia patients:
# # patients with low severity scores
# dataset2_participant9 = simulate_participant3(1,0.10,"DD")
# dataset2_participant10 = simulate_participant3(2,0.20,"DD")
# dataset2_participant11 = simulate_participant3(3,0.30,"DD")
# dataset2_participant12 = simulate_participant3(4,0.40,"DD")
# # patients with high severity scores
# dataset2_participant13 = simulate_participant3(5,0.50,"DD")
# dataset2_participant14 = simulate_participant3(6,0.60,"DD")
# dataset2_participant15 = simulate_participant3(7,0.70,"DD")
# dataset2_participant16 = simulate_participant3(8,0.80,"DD")

# Combine datasets into one DataFrame
combined_dataset3 = pd.concat(
    [
        simulate_participant3(subj, sevscore, diagnosis)
        for subj, sevscore, diagnosis in zip(subj_list, sevscore_list, diagnosis_list)
    ],
    ignore_index=True,
)
combined_dataset3

Out[4]:

	rt	response	zSTN	zGPe	participant_id	sevScore	diagnosis
0	1.498441	1.0	0.794418	4.136629	1	0.1	PD
1	0.773924	1.0	5.543321	1.464002	1	0.1	PD
2	1.192634	1.0	1.550020	2.637802	1	0.1	PD
3	3.366054	-1.0	-1.587700	2.577616	1	0.1	PD
4	1.696056	1.0	-0.155725	2.749833	1	0.1	PD
...	...	...	...	...	...	...	...
4795	2.073094	1.0	1.280471	2.647678	8	0.8	DD
4796	2.650885	-1.0	-2.132285	-1.501543	8	0.8	DD
4797	0.852617	-1.0	-2.271584	1.977184	8	0.8	DD
4798	4.889500	1.0	-1.226093	3.485166	8	0.8	DD
4799	1.337355	1.0	-0.775202	3.994313	8	0.8	DD

4800 rows × 7 columns

Model 1: diagnosis as a fixed interaction¶

v ~ 1 + (zSTN + zGPe)*sevScore*C(diagnosis) + ((1 + zSTN + zGPe)|participant_id)

Note idata_kwargs=dict(log_likelihood=True) in the sample() call — without it, az.compare has nothing to work with.

In [5]:

Copied!





model_reg_v_ex4_A1 = hssm.HSSM(
    data=combined_dataset3,
    include=[
        {
            "name": "v",
            "formula": (
                "v ~ 1 + (zSTN + zGPe)*sevScore*C(diagnosis) + "
                "((1 + zSTN + zGPe)|participant_id)"
            ),
            "prior": {
                # All ways to specify priors in the non-regression case
                # work the same way here.
                "Intercept": {"name": "Normal", "mu": 1.5, "sigma": 1.5},
                "zSTN": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zGPe": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zSTN:sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zGPe:sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "1|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
                "zSTN|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
                "zGPe|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
            },
            "link": "identity",
        }
    ],
    noncentered=True,
    p_outlier=0.05,
)
model_reg_v_ex4_A1
model_reg_v_ex4_A1 = hssm.HSSM(
    data=combined_dataset3,
    include=[
        {
            "name": "v",
            "formula": (
                "v ~ 1 + (zSTN + zGPe)*sevScore*C(diagnosis) + "
                "((1 + zSTN + zGPe)|participant_id)"
            ),
            "prior": {
                # All ways to specify priors in the non-regression case
                # work the same way here.
                "Intercept": {"name": "Normal", "mu": 1.5, "sigma": 1.5},
                "zSTN": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zGPe": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zSTN:sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zGPe:sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "1|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
                "zSTN|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
                "zGPe|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
            },
            "link": "identity",
        }
    ],
    noncentered=True,
    p_outlier=0.05,
)
model_reg_v_ex4_A1

Out[5]:

Hierarchical Sequential Sampling Model
Model: ddm

Response variable: rt,response
Likelihood: analytical
Observations: 4800

Parameters:

v:
    Formula: v ~ 1 + (zSTN + zGPe)*sevScore*C(diagnosis) + ((1 + zSTN + zGPe)|participant_id)
    Priors:
        v_Intercept ~ Normal(mu: 1.5, sigma: 1.5)
        v_zSTN ~ Normal(mu: 0.0, sigma: 1.0)
        v_zGPe ~ Normal(mu: 0.0, sigma: 1.0)
        v_sevScore ~ Normal(mu: 0.0, sigma: 1.0)
        v_zSTN:sevScore ~ Normal(mu: 0.0, sigma: 1.0)
        v_zGPe:sevScore ~ Normal(mu: 0.0, sigma: 1.0)
        v_C(diagnosis) ~ Normal(mu: 0.0, sigma: 0.25)
        v_zSTN:C(diagnosis) ~ Normal(mu: 0.0, sigma: 0.25)
        v_zGPe:C(diagnosis) ~ Normal(mu: 0.0, sigma: 0.25)
        v_sevScore:C(diagnosis) ~ Normal(mu: 0.0, sigma: 0.25)
        v_zSTN:sevScore:C(diagnosis) ~ Normal(mu: 0.0, sigma: 0.25)
        v_zGPe:sevScore:C(diagnosis) ~ Normal(mu: 0.0, sigma: 0.25)
        v_1|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 1.0))
        v_zSTN|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 1.0))
        v_zGPe|participant_id ~ Normal(mu: 0.0, sigma: HalfNormal(sigma: 1.0))
    Link: identity
    Explicit bounds: (-inf, inf)

a:
    Prior: HalfNormal(sigma: 2.0)
    Explicit bounds: (0.0, inf)

z:
    Prior: Uniform(lower: 0.0, upper: 1.0)
    Explicit bounds: (0.0, 1.0)

t:
    Prior: HalfNormal(sigma: 2.0)
    Explicit bounds: (0.0, inf)


Lapse probability: 0.05
Lapse distribution: Uniform(lower: 0.0, upper: 20.0)

In [6]:

Copied!





samples_model_reg_v_ex4_A1 = model_reg_v_ex4_A1.sample(
    sampler="numpyro",
    cores=3,
    chains=3,
    draws=500,
    tune=500,
    idata_kwargs=dict(log_likelihood=True),
)
samples_model_reg_v_ex4_A1 = model_reg_v_ex4_A1.sample(
    sampler="numpyro",
    cores=3,
    chains=3,
    draws=500,
    tune=500,
    idata_kwargs=dict(log_likelihood=True),
)

In [7]:

Copied!

az.summary(
    samples_model_reg_v_ex4_A1, var_names=["~_offset", "~_id"], filter_vars="like"
)  # var_names= ["~_offset"])
az.summary(
    samples_model_reg_v_ex4_A1, var_names=["~_offset", "~_id"], filter_vars="like"
)  # var_names= ["~_offset"])

Out[7]:

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
t	0.5101	0.0062	0.5	0.52	2016	1269	1.00	0.00014	9.6e-05
v_zSTN:sevScore	-0.7	0.14	-0.88	-0.52	403	537	1.01	0.0098	0.023
v_zSTN	0.89	0.07	0.8	0.98	423	558	1.01	0.0053	0.012
v_sevScore:C(diagnosis)[PD]	0.045	0.12	-0.15	0.23	1127	937	1.00	0.0036	0.0025
v_zSTN:C(diagnosis)[PD]	-0.047	0.038	-0.11	0.012	1178	1003	1.00	0.0011	0.00077
v_sevScore	0.135	0.15	-0.091	0.35	1101	883	1.00	0.0047	0.0046
v_zGPe:sevScore:C(diagnosis)[PD]	0.062	0.058	-0.033	0.15	854	829	1.01	0.002	0.0014
v_zGPe:sevScore	-0.583	0.065	-0.69	-0.48	687	738	1.01	0.0025	0.0023
v_Intercept	0.578	0.08	0.46	0.7	1091	766	1.01	0.0027	0.0031
v_C(diagnosis)[PD]	-0.342	0.063	-0.44	-0.24	1103	862	1.00	0.0019	0.0014
v_zGPe:C(diagnosis)[PD]	-0.023	0.031	-0.072	0.027	913	924	1.01	0.001	0.00075
v_zGPe	0.395	0.034	0.34	0.45	735	665	1.00	0.0013	0.0012
z	0.5053	0.0063	0.5	0.52	2122	1123	1.01	0.00014	9.6e-05
a	1.4953	0.0158	1.5	1.5	2362	1234	1.00	0.00033	0.00022
v_zSTN:sevScore:C(diagnosis)[PD]	0.06	0.07	-0.052	0.17	1281	934	1.00	0.002	0.0013

Model 2: diagnosis as a grouping level¶

v ~ 1 + (zSTN + zGPe)*sevScore + ((1 + zSTN + zGPe)|participant_id/C(diagnosis))

Same data, same likelihood — the only difference is where diagnosis enters the model structure.

In [8]:

Copied!





model_reg_v_ex4_A2 = hssm.HSSM(
    data=combined_dataset3,
    include=[
        {
            "name": "v",
            "formula": (
                "v ~ 1 + (zSTN + zGPe)*sevScore + "
                "((1 + zSTN + zGPe)|participant_id/C(diagnosis))"
            ),
            "prior": {
                # All ways to specify priors in the non-regression case
                # work the same way here.
                "Intercept": {"name": "Normal", "mu": 1.5, "sigma": 1.0},
                "zSTN": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zGPe": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zSTN:sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zGPe:sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "1|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
                "zSTN|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
                "zGPe|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
            },
            "link": "identity",
        }
    ],
    noncentered=True,
    p_outlier=0.05,
)
model_reg_v_ex4_A2 = hssm.HSSM(
    data=combined_dataset3,
    include=[
        {
            "name": "v",
            "formula": (
                "v ~ 1 + (zSTN + zGPe)*sevScore + "
                "((1 + zSTN + zGPe)|participant_id/C(diagnosis))"
            ),
            "prior": {
                # All ways to specify priors in the non-regression case
                # work the same way here.
                "Intercept": {"name": "Normal", "mu": 1.5, "sigma": 1.0},
                "zSTN": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zGPe": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zSTN:sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "zGPe:sevScore": {"name": "Normal", "mu": 0, "sigma": 1.0},
                "1|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
                "zSTN|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
                "zGPe|participant_id": {
                    "name": "Normal",
                    "mu": 0,
                    "sigma": {"name": "HalfNormal", "sigma": 1},
                },
            },
            "link": "identity",
        }
    ],
    noncentered=True,
    p_outlier=0.05,
)

In [9]:

Copied!





samples_model_reg_v_ex4_A2 = model_reg_v_ex4_A2.sample(
    sampler="numpyro",
    cores=3,
    chains=3,
    draws=500,
    tune=500,
    idata_kwargs=dict(log_likelihood=True),
)
samples_model_reg_v_ex4_A2 = model_reg_v_ex4_A2.sample(
    sampler="numpyro",
    cores=3,
    chains=3,
    draws=500,
    tune=500,
    idata_kwargs=dict(log_likelihood=True),
)

In [10]:

Copied!

az.summary(
    samples_model_reg_v_ex4_A2, var_names=["~_offset", "~_id"], filter_vars="like"
)
az.summary(
    samples_model_reg_v_ex4_A2, var_names=["~_offset", "~_id"], filter_vars="like"
)

Out[10]:

	mean	sd	eti89_lb	eti89_ub	ess_bulk	ess_tail	r_hat	mcse_mean	mcse_sd
t	0.5099	0.0062	0.5	0.52	2252	1146	1.00	0.00013	8.9e-05
v_zSTN:sevScore	-0.682	0.1	-0.84	-0.54	803	768	1.00	0.0037	0.0039
v_zSTN	0.875	0.052	0.8	0.96	877	818	1.00	0.0019	0.0019
v_sevScore	0.15	0.27	-0.29	0.58	836	835	1.00	0.0095	0.0076
v_zGPe:sevScore	-0.55	0.057	-0.64	-0.46	843	571	1.00	0.0021	0.0019
v_Intercept	0.414	0.136	0.2	0.64	941	978	1.00	0.0045	0.0037
v_zGPe	0.381	0.029	0.34	0.43	790	569	1.00	0.0011	0.001
z	0.5051	0.0065	0.49	0.52	2637	1137	1.00	0.00013	9.3e-05
a	1.4966	0.0158	1.5	1.5	2609	1083	1.00	0.00031	0.00022

Compare the models¶

az.compare ranks models by expected log pointwise predictive density (ELPD), estimated via Pareto-smoothed importance-sampling leave-one-out cross-validation (LOO). Higher ELPD = better expected out-of-sample prediction.

In [11]:

Copied!

az.compare(
    {"Model 1": samples_model_reg_v_ex4_A1, "Model 2": samples_model_reg_v_ex4_A2}
)
az.compare(
    {"Model 1": samples_model_reg_v_ex4_A1, "Model 2": samples_model_reg_v_ex4_A2}
)

Out[11]:

	rank	elpd_diff	dse	p_worse	diag_diff	diag_elpd	p	elpd	se	weight
Model 1	0	0.0	0.0	NaN			27.5	-5600.0	91.0	0.72
Model 2	1	-2.0	3.1	0.75	\|elpd_diff\| < 4		35.0	-5600.0	91.0	0.28

Reading the table¶

rank / elpd — Model 1 ranks first, but rank alone is not a verdict.
elpd_diff vs. dse — the difference to the best model (here ≈ 2.0) comes with a standard error (here ≈ 3.1). The difference is smaller than its own standard error, and ArviZ flags |elpd_diff| < 4 as within noise: the data do not distinguish these two models. In that situation, prefer the simpler or more interpretable candidate rather than the nominal winner.
weight — pseudo-Bayesian-model-averaging weights (here ≈ 0.72 / 0.28); useful for averaging predictions, not a posterior probability that a model is "true".
p — the effective number of parameters; a sanity check that a model is not vastly more flexible than its competitor for the same predictive payoff.

Two practical cautions: LOO estimates can be unreliable when Pareto-$k$ diagnostics are high (ArviZ warns when this happens), and ELPD compares predictive performance — a model can predict well while its parameters remain hard to interpret (or vice versa).