8  Model Comparison for Large Data

When full PSIS-LOO-CV is costly, subsampling lets us estimate predictive accuracy from a smaller set of observations (Vehtari, Gelman, and Gabry 2017; Magnusson et al. 2020). We will use the wells data to see how this works in ArviZ, first for one model and then for a model comparison. We also show how the same workflow can be combined with a Laplace approximation. The theory behind the method is discussed in Section 7.9. This example is adapted from a vignette in the loo package documentation [https://mc-stan.org/loo/articles/loo2-large-data.html].

8.1 Workflow for model comparison

Subsampled model comparison follows six steps.

1. Compute posteriors \(p_A(\theta \mid y)\) and \(p_B(\theta \mid y)\) for models \(A\) and \(B\)
2. Calculate \(\tilde\pi_i\) for all \(n\) observations using an appropriate approximation
3. Use simple random sampling to select \(m\) observations (start with \(m \approx 100\))
4. Calculate \(\pi_j\) for the subsample using PSIS-LOO-CV, checking Pareto-\(k\) diagnostics for unreliable importance sampling approximations
5. Estimate \(\operatorname{elpd}_A\), \(\operatorname{elpd}_B\), and their difference
6. If the subsampling uncertainty in the ELPD difference is too large, use update_subsample() to add the same observations to both models

Using the same observations for each model is especially important. It lets ArviZ account for the correlation between the models’ predictions when estimating \(\operatorname{elpd}_D = \operatorname{elpd}_A - \operatorname{elpd}_B\). Separate subsamples lose this benefit and can make the difference look much more uncertain (Magnusson et al. 2020).

8.2 Wells data and logistic model

To demonstrate the PSIS-LOO-CV subsampling method, we use data from a household survey in a region of Bangladesh where drinking water was contaminated with arsenic. Households with unsafe arsenic levels in their wells were asked whether they would switch to using a neighbors’ safe well. The goal is to predict this binary switching decision using household characteristics (Gelman and Hill 2007).

We focus on a simple logistic regression with two predictors, the arsenic concentration in the household’s well and the distance to the nearest safe well. The dataset contains \(n = 3020\) observations. It is small enough to compute full PSIS-LOO for comparison while still providing a useful example of the subsampling workflow.

We model the indicator of switching water wells using a logistic regression with intercept, arsenic concentration, and distance to the nearest safe well in hundreds of meters. The design matrix is \(X = [1, \text{dist}/100, \text{arsenic}]\) and the likelihood is given by

\[ \Pr(y_i = 1 \mid \beta) = \operatorname{logit}^{-1}(X_i^{\mathsf{T}} \beta), \qquad \beta \sim \mathcal{N}(0, I_3). \]

import pandas as pd
import numpy as np

wells = pd.read_csv("../data/wells.csv")
wells["dist100"] = wells["dist"] / 100

X = np.column_stack([
    np.ones(len(wells)),        
    wells["dist100"].values,     
    wells["arsenic"].values 
])
y = wells["switch"].values

We can fit the model with PyMC, Bambi, or CmdStanPy and store the posterior draws in a DataTree object.

PyMC

import pymc as pm
import xarray as xr

with pm.Model() as model:
    beta = pm.Normal("beta", mu=0, sigma=1, shape=3)
    logit_p = pm.math.dot(X, beta)
    pm.Bernoulli("y", logit_p=logit_p, observed=y)

    data = pm.sample(
        draws=1000,
        tune=1000,
        chains=4,
        random_seed=4711,
        progressbar=True
    )
    pm.compute_log_likelihood(data)

data["constant_data"] = xr.Dataset({"X": (["obs_id", "coef"], X)})

Bambi

import bambi as bmb

wells["y"] = wells["switch"]

model = bmb.Model(
    'y ~ dist100 + arsenic',
    data=wells,
    family='bernoulli',
    center_predictors=False,
    priors={
        'Intercept': bmb.Prior('Normal', mu=0, sigma=1),
        'dist100': bmb.Prior('Normal', mu=0, sigma=1),
        'arsenic': bmb.Prior('Normal', mu=0, sigma=1)
    }
)

data = model.fit(
    draws=1000,
    tune=1000,
    chains=4,
    random_seed=4711
)
model.compute_log_likelihood(data)

CmdStanPy

from cmdstanpy import CmdStanModel

stan_code = """
data {
  int<lower=0> N;
  matrix[N, 3] X;
  array[N] int<lower=0, upper=1> y;
}
parameters {
  vector[3] beta;
}
model {
  beta ~ std_normal();
  y ~ bernoulli_logit_glm(X, 0, beta);
}
generated quantities {
  vector[N] log_lik;
  {
    vector[N] eta = X * beta;
    for (n in 1:N)
      log_lik[n] = bernoulli_logit_lpmf(y[n] | eta[n]);
  }
}
"""

with open("wells_logistic.stan", "w") as f:
    f.write(stan_code)

stan_model = CmdStanModel(stan_file="wells_logistic.stan")

fit = stan_model.sample(
    data={"N": len(y), "X": X, "y": y},
    chains=4,
    iter_warmup=1000,
    iter_sampling=1000,
    seed=4711,
    show_progress=False,
)

data = az.from_cmdstanpy(
    posterior=fit,
    log_likelihood={"y": "log_lik"},
    observed_data={"y": y},
    coords={"obs_id": np.arange(len(y))},
    dims={"y": ["obs_id"]},
)
data["constant_data"] = xr.Dataset({"X": (["obs_id", "coef"], X)})

Passing a fitted model directly is supported for PyMC and Bambi. For other backends, we can provide a custom log-likelihood function as shown below.

8.3 Subsampled PSIS-LOO-CV with the PLPD surrogate

The point log predictive density (PLPD) is a fast surrogate for each leave-one-out term \(\pi_i\). Magnusson et al. (2020) define it as

\[ \tilde\pi_i^{\text{plpd}} = \log p(y_i \mid X_i, \hat{\beta}), \quad \hat{\beta} = \mathbb{E}[\beta \mid y]. \]

The approximation evaluates each observation once at the posterior mean. ArviZ then computes PSIS-LOO using all posterior draws for the sampled observations and uses those values to correct the surrogate estimate.

Here we select the PLPD surrogate with method="plpd". By passing the fitted model, ArviZ builds the function used to evaluate the surrogate. The log_likelihood group in data provides the pointwise log likelihood values used by PSIS-LOO for the selected observations.

loo_plpd = az.loo_subsample(
    data=data,
    observations=100,
    var_name="y",
    method="plpd",
    model=model,
    pointwise=True,
    seed=SEED,
)

loo_plpd

Computed from 4000 by 100 subsampled log-likelihood
values from 3020 total observations.

         Estimate   SE subsampling SE
elpd_loo   -1968.4 15.6            0.3
p_loo          3.2

------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)      100  100.0%
   (0.70, 1]   (bad)         0    0.0%
    (1, Inf)   (very bad)    0    0.0%

The result includes subsampling_se, which measures the uncertainty added by using only \(100\) of the \(3020\) observations. In this run we obtained \(\widehat{\operatorname{elpd}}_{\text{loo}}\approx-1968.4\) with \(\text{SE}\approx15.6\) and \(\text{subsampling SE} \approx0.3\). All Pareto-\(k\) values in the subsample were below \(0.7\).

8.3.1 Custom log-likelihood functions

If automatic model support is not available, we can write the PLPD function ourselves. loo_subsample() calls log_lik_fn(observed, data) after replacing the posterior variables listed in param_names with their posterior means. The function must return one log likelihood value for each observation.

For logistic regression, the PLPD value is

\[ \log p(y_i \mid \hat{\beta}, X_i) = y_i \eta_i - \log(1 + \exp(\eta_i)), \qquad \eta_i = X_i^{\mathsf{T}}\hat{\beta}. \]

We stored the design matrix in the constant_data group, so the function can use it directly.

def log_lik_fn(obs, data):
    beta = data.posterior["beta"].values
    X = data.constant_data["X"].values
    logit_pred = X @ beta
    return obs * logit_pred - np.logaddexp(0, logit_pred)

The function evaluates the Bernoulli log probability from the logits and returns one value for each row of X. The stored log_likelihood group is still used for the full PSIS-LOO calculations on the subsample. When passing log_lik_fn, we use param_names to tell ArviZ which posterior variables it needs for the PLPD surrogate.

loo_plpd_custom = az.loo_subsample(
    data=data,
    observations=100,
    var_name="y",
    method="plpd",
    log_lik_fn=log_lik_fn,
    param_names=["beta"],
    pointwise=True,
    seed=SEED,
)

loo_plpd_custom

Computed from 4000 by 100 subsampled log-likelihood
values from 3020 total observations.

         Estimate   SE subsampling SE
elpd_loo   -1968.4 15.6            0.3
p_loo          3.2

------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)      100  100.0%
   (0.70, 1]   (bad)         0    0.0%
    (1, Inf)   (very bad)    0    0.0%

Using the same seed selects the same observations, and the custom and model-based routes give the same estimate.

8.4 Subsampled PSIS-LOO-CV with the LPD surrogate

The log predictive density (LPD) offers another surrogate. Instead of evaluating the likelihood only at the posterior mean, it averages over the posterior draws.

\[ \tilde\pi_i^{\text{lpd}} = \log \left( \frac{1}{S} \sum_{s=1}^S p(y_i \mid X_i, \beta_s) \right), \]

Here \(\beta_s\) are posterior draws. Using all draws costs more than PLPD but can provide a better surrogate for some models. We use the same subsample size and random seed so both calculations select the same observations.

Note

ArviZ can use the fitted model to build the PLPD surrogate. The LPD surrogate can be calculated directly from the posterior log likelihood values stored in data, so the fitted model is not needed in the call below.

loo_lpd = az.loo_subsample(
    data=data,
    observations=100,
    var_name="y",
    method="lpd",
    pointwise=True,
    seed=SEED,
)

loo_lpd

Computed from 4000 by 100 subsampled log-likelihood
values from 3020 total observations.

         Estimate   SE subsampling SE
elpd_loo   -1968.5 15.6            0.4
p_loo          3.2

------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)      100  100.0%
   (0.70, 1]   (bad)         0    0.0%
    (1, Inf)   (very bad)    0    0.0%

The LPD estimate is about \(-1968.5\), close to the PLPD estimate of about \(-1968.4\). Both have \(\text{SE}\approx15.6\). In this example, PLPD has a slightly smaller subsampling SE than LPD. The choice of surrogate depends on the model, so it is useful to compare the reported subsampling uncertainty when both are available.

8.5 Incremental sampling updates

Magnusson et al. (2020) recommend increasing \(m\) until the subsampling uncertainty is small enough for the task at hand. update_subsample() extends an existing subsample without discarding its selected observations. Here we add \(200\) observations, taking the subsample from \(100\) to \(300\).

loo_update = az.update_subsample(
    loo_plpd,
    data,
    observations=200,
    var_name="y",
    method="plpd",
    model=model,
    seed=SEED,
)

loo_update

Computed from 4000 by 300 subsampled log-likelihood
values from 3020 total observations.

         Estimate   SE subsampling SE
elpd_loo   -1968.4 15.6            0.2
p_loo          3.1

------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)      300  100.0%
   (0.70, 1]   (bad)         0    0.0%
    (1, Inf)   (very bad)    0    0.0%

The subsampling SE decreases from about \(0.3\) to \(0.2\), while the estimated ELPD changes very little.

8.6 PSIS-LOO-CV with an approximate posterior

PSIS-LOO-CV can also be used with an approximate posterior. This can reduce the time spent fitting the model before LOO-CV, though the result still depends on the quality of the posterior approximation.

We use a Laplace approximation, which represents the posterior with a multivariate normal distribution centered at the posterior mode \(\hat{\theta}\). Its covariance is the inverse Hessian of the negative log posterior evaluated at the mode (Tierney and Kadane 1986). For a model with likelihood \(p(y \mid \theta)\) and prior \(p(\theta)\), the approximation is

\[ q(\theta \mid y) = \mathcal{N}(\theta \mid \hat{\theta}, \Sigma), \]

where

\[ \hat{\theta} = \arg\max_\theta \log p(\theta \mid y) \quad \text{ and } \quad \Sigma = \left[-\nabla_\theta^2 \log p(\theta \mid y)\big|_{\theta=\hat{\theta}}\right]^{-1}. \]

Draws from \(q(\theta \mid y)\) can be obtained efficiently by sampling from the fitted normal distribution.

The following uses pymc_extras.fit_laplace to fit the wells model. It returns draws from the approximate posterior together with estimates of its mode and covariance.

import pymc as pm
from pymc_extras import fit_laplace

with pm.Model() as model_laplace:
    beta = pm.Normal("beta", mu=0, sigma=1, shape=3)
    logit_p = pm.math.dot(X, beta)
    pm.Bernoulli("y", logit_p=logit_p, observed=y)

    data_laplace = fit_laplace(
        draws=2000,
        random_seed=SEED,
        progressbar=True
    )

The fit_laplace function stores the posterior mode in the fit group as mean_vector and the covariance in covariance_matrix. Posterior draws are stored in the posterior group as usual. We attach the observed data and design matrix for use with loo_approximate_posterior().

data_laplace["observed_data"] = xr.Dataset({
    "y": ("obs_id", y)
}, coords={"obs_id": np.arange(len(y))})
data_laplace["constant_data"] = xr.Dataset({"X": (["obs_id", "coef"], X)})

8.6.1 Computing log density corrections

To use loo_approximate_posterior(), we supply log_p and log_q. These are the target and approximate posterior log densities evaluated at draws from the approximate posterior. The target density only needs to be known up to an additive constant. The importance ratios are

\[ r_{is} \propto \frac{p(\theta_s \mid y)}{q(\theta_s \mid y)\,p(y_i \mid \theta_s)}. \]

For the Laplace approximation, log_q is the log density of the fitted normal distribution. We compute log_p as the log prior plus the log likelihood, which gives the target log density up to an additive constant.

from scipy.stats import multivariate_normal

beta_samples = data_laplace.posterior["beta"]
X_xr = data_laplace.constant_data["X"]
y_xr = data_laplace.observed_data["y"]

param_dim = [d for d in beta_samples.dims if d not in ["chain", "draw"]][0]
beta_renamed = beta_samples.rename({param_dim: "coef"})

logit_pred = xr.dot(X_xr, beta_renamed, dims=["coef"])
log_lik_laplace = y_xr * logit_pred - np.logaddexp(0, logit_pred)

data_laplace["log_likelihood"] = xr.Dataset({"y": log_lik_laplace})

mean_vals = data_laplace.fit["mean_vector"].values
cov_matrix = data_laplace.fit["covariance_matrix"].values

log_q_vals = multivariate_normal.logpdf(beta_renamed.values, mean=mean_vals, cov=cov_matrix)
log_q_da = xr.DataArray(
    log_q_vals,
    dims=["chain", "draw"],
    coords={"chain": beta_renamed.chain, "draw": beta_renamed.draw}
)

log_prior = -0.5 * (beta_renamed ** 2).sum("coef") - 1.5 * np.log(2 * np.pi)
log_p_da = log_prior + log_lik_laplace.sum("obs_id")

The log prior uses the standard normal density \(\mathcal{N}(0, I_3)\) specified in the model, and the log likelihood is summed over all observations. We then pass log_p_da and log_q_da to loo_approximate_posterior().

8.6.2 Running LOO-CV with an approximate posterior

We can now compute LOO-CV from the approximate posterior. loo_approximate_posterior() adjusts the importance ratios using the difference between the target posterior and the Laplace approximation.

loo_ap = az.loo_approximate_posterior(
    data=data_laplace,
    log_p=log_p_da,
    log_q=log_q_da,
    var_name="y",
    pointwise=True
)

loo_ap

Computed from 8000 posterior samples and 3020 observations log-likelihood matrix.
Posterior approximation correction used.

         Estimate       SE
elpd_loo -1968.42    15.59
p_loo        3.17        -
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)     3020  100.0%
   (0.70, 1]   (bad)         0    0.0%
    (1, Inf)   (very bad)    0    0.0%

The Pareto-\(k\) diagnostics remain available and help us assess the importance sampling correction. The output also records that the posterior approximation correction was used.

8.7 Combining approximate posterior with subsampling

We can combine the approximate posterior correction with subsampling by passing log_p and log_q to loo_subsample(). The importance weights account for draws from the approximate posterior, while the difference estimator uses the corrected LOO values for the sampled observations to adjust the PLPD surrogate. The reported subsampling_se measures the remaining uncertainty from sampling observations.

loo_ap_ss = az.loo_subsample(
    data=data_laplace,
    observations=100,
    var_name="y",
    method="plpd",
    model=model_laplace,
    log_p=log_p_da,
    log_q=log_q_da,
    pointwise=True,
    seed=SEED
)

loo_ap_ss

Computed from 8000 by 100 subsampled log-likelihood
values from 3020 total observations. Posterior approximation correction used.

         Estimate   SE subsampling SE
elpd_loo   -1968.6 15.6            0.4
p_loo          3.4

------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)      100  100.0%
   (0.70, 1]   (bad)         0    0.0%
    (1, Inf)   (very bad)    0    0.0%

The Laplace approximation reduces the fitting cost, while subsampling reduces the number of observations that require PSIS-LOO calculations. The Pareto-\(k\) and subsampling diagnostics remain important checks.

In this example, the subsampled result is close to the full approximate posterior result and has a small subsampling SE.

8.8 Comparing models with subsampled LOO

We can compare subsampled LOO estimates with compare() just as we compare regular LOO estimates. We will fit a second model using log(arsenic) instead of arsenic as a predictor and compare the two models.

PyMC

X_log = X.copy()
X_log[:, 2] = np.log(X[:, 2])

with pm.Model() as model2:
    beta2 = pm.Normal("beta", mu=0, sigma=1, shape=3)
    logit_p2 = pm.math.dot(X_log, beta2)
    pm.Bernoulli("y", logit_p=logit_p2, observed=y)

    data2 = pm.sample(
        draws=1000,
        tune=1000,
        chains=4,
        random_seed=SEED,
        progressbar=True
    )
    pm.compute_log_likelihood(data2)

data2["constant_data"] = xr.Dataset({"X": (["obs_id", "coef"], X_log)})

Bambi

wells["log_arsenic"] = np.log(wells["arsenic"])

model2 = bmb.Model(
    'y ~ dist100 + log_arsenic',
    data=wells,
    family='bernoulli',
    center_predictors=False,
    priors={
        'Intercept': bmb.Prior('Normal', mu=0, sigma=1),
        'dist100': bmb.Prior('Normal', mu=0, sigma=1),
        'log_arsenic': bmb.Prior('Normal', mu=0, sigma=1)
    }
)

data2 = model2.fit(
    draws=1000,
    tune=1000,
    chains=4,
    random_seed=SEED
)
model2.compute_log_likelihood(data2)

CmdStanPy

X_log = X.copy()
X_log[:, 2] = np.log(X[:, 2])

fit2 = stan_model.sample(
    data={"N": len(y), "X": X_log, "y": y},
    chains=4,
    iter_warmup=1000,
    iter_sampling=1000,
    seed=SEED,
    show_progress=False,
)

data2 = az.from_cmdstanpy(
    posterior=fit2,
    log_likelihood={"y": "log_lik"},
    observed_data={"y": y},
    coords={"obs_id": np.arange(len(y))},
    dims={"y": ["obs_id"]},
)
data2["constant_data"] = xr.Dataset({"X": (["obs_id", "coef"], X_log)})

8.8.1 Comparison with different subsamples

We first compute subsampled LOO for the second model using a different random seed than the first model. This means the two models will be evaluated on different subsets of observations.

loo_ss_2_diff = az.loo_subsample(
    data=data2,
    observations=100,
    var_name="y",
    method="plpd",
    model=model2,
    pointwise=True,
    seed=315
)

loo_ss_2_diff

Computed from 4000 by 100 subsampled log-likelihood
values from 3020 total observations.

         Estimate   SE subsampling SE
elpd_loo   -1952.3 16.2            0.3
p_loo          3.2

------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)      100  100.0%
   (0.70, 1]   (bad)         0    0.0%
    (1, Inf)   (very bad)    0    0.0%

Stacking needs pointwise LOO values for every observation, so we use method="pseudo-bma" for the subsampled results.

comparison_diff = az.compare(
    {
        "model1_arsenic": loo_plpd,
        "model2_log_arsenic": loo_ss_2_diff,
    },
    method="pseudo-bma",
    round_to=1,
)

comparison_diff

/opt/hostedtoolcache/Python/3.13.16/x64/lib/python3.13/site-packages/arviz_stats/loo/compare.py:378: UserWarning: Different subsamples used in 'model2_log_arsenic' and 'model1_arsenic'. Naive diff SE is used.
  warnings.warn(

	rank	elpd_diff	dse	p_worse	diag_diff	diag_elpd	p	elpd	se	weight	subsampling_dse
model2_log_arsenic	0	0.0	0.0	NaN			3.2	-1952.3	16.2	1.0	0.0
model1_arsenic	1	-16.1	22.5	0.8			3.2	-1968.4	15.6	0.0	0.4

NoteNote on dse and subsampling_dse

dse measures the usual uncertainty in the ELPD difference across possible datasets. subsampling_dse measures the additional uncertainty from estimating that difference with a subsample. We inspect subsampling_dse separately to decide whether more observations are needed.

Because the subsamples differ, ArviZ warns that it cannot use the correlation between the models’ predictions. The resulting dse treats the two estimates as independent and is much larger.

8.8.2 Comparison with the same subsample

We now reuse the same observations by passing the indices stored in loo_plpd.loo_subsample_observations to the second loo_subsample() call.

loo_ss_2 = az.loo_subsample(
    data=data2,
    observations=loo_plpd.loo_subsample_observations,
    var_name="y",
    method="plpd",
    model=model2,
    pointwise=True,
)

loo_ss_2

Computed from 4000 by 100 subsampled log-likelihood
values from 3020 total observations.

         Estimate   SE subsampling SE
elpd_loo   -1952.2 16.2            0.2
p_loo          3.2

------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)      100  100.0%
   (0.70, 1]   (bad)         0    0.0%
    (1, Inf)   (very bad)    0    0.0%

comparison_same = az.compare(
    {
        "model1_arsenic": loo_plpd,
        "model2_log_arsenic": loo_ss_2,
    },
    method="pseudo-bma",
    round_to=1,
)

comparison_same

	rank	elpd_diff	dse	p_worse	diag_diff	diag_elpd	p	elpd	se	weight	subsampling_dse
model2_log_arsenic	0	0.0	0.0	NaN			3.2	-1952.2	16.2	1.0	0.0
model1_arsenic	1	-16.2	4.4	1.0			3.2	-1968.4	15.6	0.0	0.1

With shared observations, dse falls from about \(22\) to \(4.4\) because ArviZ can account for the correlation between the models. The second model, which uses log-transformed arsenic, has the higher ELPD. The small subsampling_dse suggests that \(100\) observations are enough for this comparison. If it were larger, we could use update_subsample() to add the same observations to both models.

8.9 Comparing full LOO to subsampled LOO

We can compare a subsampled LOO result with a full result computed over all \(N = 3020\) observations. Here we compute full PSIS-LOO-CV for the first model using loo().

loo_full = az.loo(data, var_name="y", pointwise=True)
loo_full

Computed from 4000 posterior samples and 3020 observations log-likelihood matrix.

         Estimate       SE
elpd_loo -1968.47    15.60
p_loo        3.24        -
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.70]   (good)     3020  100.0%
   (0.70, 1]   (bad)         0    0.0%
    (1, Inf)   (very bad)    0    0.0%

ArviZ can also estimate the model difference when only one of the LOO results uses subsampling.

comparison_mixed = az.compare(
    {
        "model1_full_loo": loo_full,
        "model2_subsampled_loo": loo_ss_2,
    },
    method="pseudo-bma",
    round_to=1,
)

comparison_mixed

/opt/hostedtoolcache/Python/3.13.16/x64/lib/python3.13/site-packages/arviz_stats/loo/compare.py:445: UserWarning: Estimated elpd_diff using observations included in loo calculations for all models.
  warnings.warn(

	rank	elpd_diff	dse	p_worse	diag_diff	diag_elpd	p	elpd	se	weight	subsampling_dse
model2_subsampled_loo	0	0.0	0.0	NaN			3.2	-1952.2	16.2	1.0	0.0
model1_full_loo	1	-16.3	4.4	1.0			3.2	-1968.5	15.6	0.0	0.2

ArviZ uses the \(100\) shared observations for the paired correction, while the surrogate values still cover all \(3020\) observations. The reported elpd_diff therefore estimates the difference across the full dataset. The warning tells us that the paired correction is limited to the shared observations.

8.10 Practical considerations

Start with a modest number of observations and inspect subsampling_se together with the Pareto-\(k\) diagnostics. If the subsampling uncertainty is too large for the question, add observations with update_subsample(). PLPD is quick to compute, while LPD can be useful when the posterior mean is not a good summary of the predictive distribution.

When comparing models, reuse the same subsample. This lets ArviZ account for the relationship between the models’ predictions and usually gives a more useful estimate of the uncertainty in their ELPD difference.

Gelman, Andrew, and Jennifer Hill. 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge: Cambridge University Press. https://sites.stat.columbia.edu/gelman/arm/.

Magnusson, Måns, Michael Riis Andersen, Johan Jonasson, and Aki Vehtari. 2020. “Leave-One-Out Cross-Validation for Model Comparison in Large Data.” In Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics. Vol. 108. Proceedings of Machine Learning Research. PMLR. https://arxiv.org/abs/2001.00980.

Tierney, Luke, and Joseph B. Kadane. 1986. “Accurate Approximations for Posterior Moments and Marginal Densities.” Journal of the American Statistical Association 81 (393): 82–86. https://doi.org/10.2307/2287970.

Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32. https://doi.org/10.1007/s11222-016-9696-4.