[avehtari]

Summary

In cv_varsel(cv_method = "LOO", validate_search = TRUE), the projections of fold i use observation i, although the reference predictions of fold i do not depend on y_i. The submodels then look better than they are, and the larger submodels more so. This PR gives observation i the weight 0 in the projections of its fold. The elpd estimates of validated LOO change, mostly for larger submodels.

This PR was assisted by Claude.

Problem

The search and the evaluation of fold i use the reference draws reweighted with the PSIS LOO weights of observation i, so the reference predictions do not depend on y_i. But loo_varsel() passes refmodel, with refmodel$wobs including observation i, to .select() and perf_eval(). So the submodels of fold i are fitted also at x_i, with the reference model's LOO prediction at x_i as the target, and are then evaluated at the same x_i.

The effect grows with the size of the submodel. From size n − 1 on, a submodel can interpolate the reference predictions at all n inputs. Its fold lppd at x_i is then the reference model's, and the elpd difference to the reference model is the same for every size from n − 1 on. K-fold CV with refits (cv_method = "kfold") does not have this problem, because neither the fold reference model nor the fold projections use the left-out observations.

Example: Gaussian data with n = 30 and P = 40, two relevant predictors, brms reference model with an R2D2 prior, forward search to 35 terms, 5 search clusters, 100 evaluation draws. Elpd difference to the reference model (SE):

size	LOO, master	K-fold, K = n	LOO, this change
0	0.19 (0.70)	0.56 (0.71)	0.15 (0.72)
5	0.29 (0.64)	−3.01 (1.15)	−1.27 (1.24)
10	1.06 (0.51)	−1.12 (1.11)	−0.53 (1.39)
20	0.72 (0.53)	−2.15 (0.84)	−0.21 (1.32)
28	0.57 (0.59)	−4.63 (1.89)	−6.67 (2.47)
29	0.55 (0.60)	−2.70 (1.47)	−2.22 (1.82)
35	0.55 (0.60)	−0.35 (0.32)	0.47 (0.69)

On master, the difference is exactly 0.55 for every size from 29 = n − 1 to 35, and the drop near the interpolation threshold that K-fold CV shows (−4.63 at size 28) is missing. Over sizes 5 to 25, the mean difference is 0.76 with LOO on master, −1.87 with K-fold CV and −0.63 with this change.

In 16 simulated data sets (Gaussian, N 50 and 100, P 50 and 200, R2D2 reference model, forward search, full LOO, 10 MCMC runs each), I compared with a brute-force validated LOO: refit without y_i, then search and evaluation without observation i. The mean bias of the submodel elpd_loo, summed over the observations, was 0.7 nat at P 50 and 6.7 nat at P 200 on master, and 0.6 and 0.3 nat with this change.

Independent test data give a more precise check than brute-force LOO. In 2 400 simulated data sets (Gaussian response, N 50 and 100, P 50, 100 and 200, R2D2 brms reference model, forward search, full validated LOO), I compared validated LOO with and without the change against 1 000 new observations from the same process, scored with the submodels of the full-data path. Bias in nat per observation, LOO minus test, with the Monte Carlo SE over data sets; Δ is the difference to the reference model:

size	elpd, master	elpd, this change	Δ, master	Δ, this change
0	−0.010 (0.002)	−0.012 (0.002)	−0.002 (0.001)	−0.004 (0.001)
≤ 0.1n	+0.002 (0.002)	−0.008 (0.002)	+0.010 (0.001)	−0.001 (0.001)
0.1n to 0.3n	+0.006 (0.002)	−0.004 (0.002)	+0.014 (0.001)	+0.003 (0.001)
0.3n to 0.6n	+0.006 (0.002)	−0.005 (0.002)	+0.014 (0.001)	+0.003 (0.001)
0.6n to 0.9n	+0.013 (0.002)	−0.008 (0.002)	+0.021 (0.001)	−0.000 (0.001)
> 0.9n	+0.089 (0.003)	−0.009 (0.002)	+0.096 (0.003)	−0.001 (0.002)

On master, the optimism grows with the submodel size, to 0.09 nat per observation above 0.9n (0.17 for P 200). With this change, the bias is −0.004 to −0.009 per observation in every size group, like the bias of the intercept-only and the reference model: LOO trains on n − 1 observations and is therefore slightly pessimistic. Δ is within ±0.006 per observation of the test value.

Changes

The PR has five commits: the four code commits below and one that adds the PR number to NEWS.md.

1. Observation i out of the projections of fold i. In one_obs() of loo_varsel(), the search (.select()) and the evaluation (perf_eval()) of fold i get a copy of the reference model object with wobs[i] = 0. perf_eval() gets the unchanged object as refmodel_fulldata, so the test weight and the data of observation i stay the same. The change applies only to validate_search = TRUE. With validate_search = FALSE, the full-data path was selected with y_i, and leaving x_i out of the fold projections made the estimates of the middle sizes more optimistic in the simulation; there, validate_search = TRUE is the remedy.

2. Observation weights in the augmented-data projection. The augmented-data submodel fitters use the probabilities of the response categories as weights and ignore weights, so commit 1 alone did not remove observation i from the augmented-data projection. proj_to_submodl() now multiplies these probabilities by the observation weights if some weight is not 1. Outside the fold copy of commit 1, this does not happen, because init_refmodel() allows only weights 1 for the augmented-data projection. The latent projection uses the Gaussian submodel fitters with weights, so commit 1 covers it.

3. Zero-weight observations out of the data of the multilevel and additive submodel fitters. Not all submodel fitters ignore an observation with weight 0. With weight 0 for one observation, a change of its response still changed the fit of lme4::lmer() (theta by 0.2 to 0.6 in two examples), so commit 1 alone did not remove observation i from Gaussian multilevel submodels, also in the latent projection. mgcv::gam() ignored the response, but its smooth bases still used the predictor values of the observation. For submodels with group-level or additive terms, divmin() now removes the observations with weight 0 from the data before it calls the submodel fitter. The GLM fitters handle weight 0 exactly, so GLM submodels keep all rows (this also keeps the full model matrix that other code can share between the draws). The submodel fits still give predictions for all observations; if the left-out observation is the only one of its group level, the prediction treats that level as new. Commit 4 covers a level of a factor predictor that only the left-out observation has. Observations with weight 0 occur only in the fold projections of commit 1, so other results do not change, and the run time does not change measurably.

4. The same removal for the augmented-data submodels, and unused factor levels. Commit 3 removed the rows with weight 0 only in divmin(). divmin_augdat() passed these rows to its fitters. In my check with 80 observations, the ordinal::clmm() coefficients did not change when these rows were removed. For mclogit::mblogit(), setting the predictor value of the observation with weight 0 to 25 (the other values were standard normal) changed the coefficients by up to 4.8e-3 and the random-effect variance by up to 0.10. MASS::polr() gave no coefficient for a factor level that only the observation with weight 0 has, and the prediction then failed. divmin_augdat() now removes the rows of the observations with weight 0 (in all response categories) before it calls any fitter, with and without group-level terms.

The row removal of commits 3 and 4 can make a level of a factor predictor unused, if only the left-out observation has that level. The fits of lmer(), gam(), clmm(), mblogit() and polr() then have no coefficient for that level. The prediction for the left-out observation failed (for example non-conformable arguments for lmer() submodels). divmin() and divmin_augdat() now store the unused levels and the reference level of each factor in the attribute projpred_lost_lvls of the fits. They also drop the unused levels from the remaining data, so that no fitter builds a column of zeros for them. Before the prediction, mu_fun() replaces the stored levels in newdata with the reference level, the first level that the fit keeps. With treatment contrasts, the unused level then gets the prediction of a level with coefficient 0. The GLM fitter keeps the rows with weight 0 and predicts nearly the same, because its ridge penalty shrinks the coefficient of such a level to almost 0. With other contrasts and for ordered factors, the unused level gets the prediction of the first remaining level. The replacement also covers character and logical predictors, and numeric variables that the formula uses only inside factor() or ordered(). Grouping factors are not replaced; the fit treats an unseen group level as new (see commit 3). proj_linpred() and proj_predict() apply the same replacement to newdata when the projection has observations with weight 0.

Note on the PSIS estimate of the reference model

PSIS-LOO can overestimate the elpd of the reference model when many Pareto k values are high (P 200 in the simulation: by 1.4 nat on average). On master, the bias of the submodels cancels part of that error in the elpd difference. With this change, the default suggest_size() rule selected 28 terms on average at P 200 instead of about 5 (brute force: about 4). A more reliable reference elpd_loo for that case, for example balanced mixture importance sampling, is a separate change.

Results

The change costs no measurable run time: in the simulation, the median time ratio of cv_varsel() with and without it was 1.01. For subsampled LOO (nloo < n), the difference estimator combines the fast estimate without validated search with the validated folds, and it agrees with the full validated LOO of the same version (n = 40, P = 10, 10 subsamples of 10 observations: the mean differs from the full validated LOO by at most 1.3 standard errors per size, as on master).

Tests

The new test_loo_fold_projection.R checks that:

• with n = 8 observations and 9 predictors, the fold lppd of the submodels of sizes n − 1 to 9 is not the same for all these sizes (on master the spread is about 2e-5, because these submodels interpolate the reference predictions at the left-out observation; with the change it is 1.7);
• for one test setup of each combination of model type (GLM, multilevel, additive, additive multilevel), family, projection type (traditional, latent, augmented-data) and observation weights, 15 setups in total, a changed reference target of an observation with weight 0 does not change the submodel predictions, and with its original weight it does; for the multilevel and additive submodels, changed predictor values of the observation with weight 0 do not change the predictions either;
• if an observation with weight 0 is the only one with a factor level, the fits of lmer(), gam(), clmm(), mblogit(), polr() and multinom() store that level, and the prediction uses the reference level for it. The test covers an unused first level, a character predictor, factor(z) with numeric z, and a grouping factor that is also a fixed-effect predictor (the replacement skips it). It also runs the path through proj_to_submodl() and mu_fun() with a reference model.

On master, the first test fails. The second test failed on master for the augmented-data setup in an earlier version of the test; for the Gaussian multilevel case, divmin() on master gives a different fit when the response of the weight-0 observation changes (theta by 0.6). With commits 1 to 3 only, the second test fails for one multilevel augmented-data setup: other predictor values of observation 1 give other values for 144 of 351 predicted category probabilities. The third test also fails. The full suite passes with NOT_CRAN=true: 256 tests, 55 866 expectations, 0 failures, the same 4 warnings as master. With the first two commits merged into #579, the full suite also passed (262 tests), including the comparison of the cross-product and the draw-wise evaluation of validated LOO. #579 changes none of divmin(), divmin_augdat() and mu_fun(). Its cross-product path covers only GLM submodels, which keep all rows, so commits 3 and 4 do not interact with #579. Tested with this change applied to master 73cf0ee (projpred 2.11.0.9000), R 4.6.1, brms 2.23.0, rstanarm 2.32.2, lme4 2.0-6, mgcv 1.9-4, ordinal 2025.12.29, mclogit 0.9.15, cmdstanr 0.9.0.9001, CmdStan 2.39.0, loo 2.10.1.9000, on CSC Roihu (Linux, AMD EPYC).

[avehtari]

Additional information

The problem

For observation $i$, LOO CV with validated search estimates the predictive performance of a submodel $M$ built without observation $i$:

$$ \mathrm{elpd}_i(M) = \log \int p(y_i \mid x_i, \phi) , p_M(\phi \mid D_{-i}) , d\phi , $$

where $D_{-i}$ are the data without $(x_i, y_i)$. A brute-force computation refits the reference model without observation $i$, runs the search on the other observations, and projects each draw $\theta_s$ of $p(\theta \mid D_{-i})$ onto $M$ with the observations $j \ne i$.

projpred approximates $p(\theta \mid D_{-i})$ with PSIS-LOO weights $w_s^{(i)}$ on the draws of the full-data posterior. The reference predictions $\tilde{\mu}_s(x)$ of fold $i$ therefore do not use $y_i$. But master projects each draw with all $n$ observations:

$$ \phi_s^{(i)} = \arg\max_{\phi} \sum_{j=1}^{n} v_j , \mathrm{E}_{\tilde{y}_j \mid \theta_s} \left[ \log p(\tilde{y}_j \mid x_j, \phi) \right] , $$

where $v_j$ are the observation weights. The term $j = i$ fits the submodel at $x_i$ to the reference model's LOO prediction there, and the evaluation then predicts at the same $x_i$.

For a Gaussian submodel, the projection of a draw is a weighted ridge regression of the reference fits $\tilde{\mu}_{js}$ on the model matrix. This is a linear smoother with a hat matrix $H$. The fit at $x_i$ with observation $i$ and the fit without it satisfy

$$ \hat{\mu}_{is} = H_{ii} , \tilde{\mu}_{is} + (1 - H_{ii}) , \hat{\mu}_{is}^{(-i)} , $$

where $\hat{\mu}_{is}^{(-i)}$ is the fit at $x_i$ from the observations $j \ne i$, which the brute force uses. The identity is exact for a fixed penalty and fixed column scales. A numerical check gives differences below 1e-15. glm_ridge() computes the column scales from the weighted data, so in projpred the identity holds up to that small difference.

So master moves the submodel prediction at $x_i$ toward the reference model's LOO prediction, with the leverage $H_{ii}$ as the weight. The leverages sum to the number of model-matrix columns with the intercept, $p + 1$, up to $n$. The shift therefore grows with the submodel size. From $p + 1 \ge n$ on, $H_{ii} \approx 1$ because the ridge penalty is small. The submodel then interpolates the reference fits, its prediction at $x_i$ equals the reference model's LOO prediction, and the elpd difference to the reference model is the same for all these sizes. The brute force and K-fold CV show a drop of the predictive performance near the interpolation threshold, and master does not.

1. Weight 0 for observation $i$ in the projections of its fold

In one_obs() of loo_varsel(), the search and the evaluation of fold $i$ get a copy of the reference model object with $v_i = 0$. The sum in the projection then runs over $j \ne i$, as in the brute force, and the prediction at $x_i$ is $\hat{\mu}_{is}^{(-i)}$. The evaluation still predicts at $x_i$ and uses $y_i$ with its original weight. It takes both from the unchanged object, which it gets as refmodel_fulldata.

The change applies only to validate_search = TRUE. With validate_search = FALSE, the search uses $y_i$, and $x_i$ in the projection partly offsets that selection optimism. In the simulation, leaving $x_i$ out in that case made the estimates for the middle sizes more optimistic.

Subsampled LOO, with $m$ = nloo $< n$, uses the difference estimator

$$ \widehat{\mathrm{elpd}}(M) = \sum_{j=1}^{n} e_j^{\mathrm{fast}} + \frac{n}{m} \sum_{i \in S} \left( e_i^{\mathrm{val}} - e_i^{\mathrm{fast}} \right) , $$

where $e_j^{\mathrm{fast}}$ are the terms of the search on the full data and $e_i^{\mathrm{val}}$ the validated terms for a random subset $S$. The fix changes only the validated terms, so the estimator targets the full validated LOO with the fix.

2. The augmented-data projection

The augmented-data projection expands each observation $j$ into one row per response category $c = 1, \dots, C$. The row for draw $s$ has the reference probability $\pi_{jcs}$ of the category as its weight:

$$ \phi_s = \arg\max_{\phi} \sum_{j=1}^{n} v_j \sum_{c=1}^{C} \pi_{jcs} \log p(c \mid x_j, \phi) . $$

The submodel fitters get only $\pi_{jcs}$ as weights and ignore $v_j$, so commit 1 alone leaves observation $i$ in the projection. proj_to_submodl() now passes $v_j \pi_{jcs}$. Row $(c - 1) n + j$ of the augmented matrix belongs to observation $j$ and category $c$, so the weight vector is $v$ repeated $C$ times. init_refmodel() allows only $v_j = 1$ for the augmented-data projection, so the multiplication changes only the fold copies of commit 1.

3. Zero-weight rows out of the multilevel and additive fits

Not every submodel fitter ignores a row with weight 0. I gave one observation weight 0 and changed its response. The fit of lme4::lmer() changed: the group-level standard deviation moved by 0.2 to 0.6. mgcv::gam() ignored the response but used the predictor value of the row for its smooth bases. glm_ridge(), glmer(), polr(), multinom(), clmm() and mblogit() ignored the row.

For submodels with group-level or additive terms, divmin() now removes the rows with weight 0 before it calls the fitter, so these fits use exactly the observations $j \ne i$. GLM submodels keep all rows: their fitter handles weight 0 exactly, and the draws can share one model matrix.

Check against independent test data

Validated LOO estimates the expected predictive performance of the whole procedure (search, projection, submodel of size $k$) for new data from the same process. Independent test data estimate the same quantity with much less variance than brute-force LOO. For each data set $r$ and size $k$, I compared the mean log predictive density per observation:

$$ b_{rk} = \frac{1}{n} \sum_{i=1}^{n} \log p_{\mathrm{LOO}}(y_i \mid M_k) - \frac{1}{N} \sum_{t=1}^{N} \log p(\tilde{y}_t \mid M_k) , $$

where the second term scores the submodels of the full-data path on $N = 1000$ new observations. The bias is the mean of $b_{rk}$ over data sets, with its Monte Carlo standard error. For the elpd difference to the reference model, I used the same comparison for $\Delta_k = \mathrm{elpd}(M_k) - \mathrm{elpd}(\mathrm{ref})$. LOO trains on $n - 1$ observations and the test comparison on $n$, so LOO should be slightly pessimistic.

Data: 2 400 simulated data sets with a Gaussian response (N 50 and 100, P 50, 100 and 200), an R2D2 brms reference model, forward search and full validated LOO, each with and without the fix on the same data and reference fit. The sizes are grouped relative to $n$; the interpolation threshold is at size $n - 1$.

size	elpd, master	elpd, fix	$\Delta$, master	$\Delta$, fix
0	−0.010 (0.002)	−0.012 (0.002)	−0.002 (0.001)	−0.004 (0.001)
≤ 0.1n	+0.002 (0.002)	−0.008 (0.002)	+0.010 (0.001)	−0.001 (0.001)
0.1n to 0.3n	+0.006 (0.002)	−0.004 (0.002)	+0.014 (0.001)	+0.003 (0.001)
0.3n to 0.6n	+0.006 (0.002)	−0.005 (0.002)	+0.014 (0.001)	+0.003 (0.001)
0.6n to 0.9n	+0.013 (0.002)	−0.008 (0.002)	+0.021 (0.001)	−0.000 (0.001)
> 0.9n	+0.089 (0.003)	−0.009 (0.002)	+0.096 (0.003)	−0.001 (0.002)

The values are LOO minus test in nat per observation, with the Monte Carlo standard error over data sets.

• On master, the submodels look better than the test data show, and more so for larger submodels. Above $0.9n$, the optimism is 0.09 nat per observation, about 5 to 10 nat in total; for P 200 it is 0.17 per observation.
• With the fix, the bias is about −0.004 to −0.009 per observation in every size group. That is the bias of the intercept-only model and of the reference model (−0.004 to −0.013), the expected pessimism of LOO with $n - 1$ training observations. So the fix removes the size-dependent optimism.
• With the fix, $\Delta_k$ is within ±0.006 per observation of the test value for all P, about ±0.3 to ±0.6 nat in total for these $n$.

Interaction with the PSIS estimate of the reference model

suggest_size() uses the difference $\Delta_k = \mathrm{elpd}(M_k) - \mathrm{elpd}(\mathrm{ref})$. With many high Pareto $k$ values, PSIS-LOO overestimates $\mathrm{elpd}(\mathrm{ref})$; in the simulation with $P = 200$, by 1.4 nat on average. On master, the optimism of the submodels from $H_{ii}$ partly cancels this error in $\Delta_k$. The fix removes the optimism of the submodel terms, but the reference term keeps its error. $\Delta_k$ can then become too low, and suggest_size() selects larger submodels. A more accurate $\mathrm{elpd}(\mathrm{ref})$ for this case, for example from balanced mixture importance sampling, is a separate change.