Subsampled LOO

[fweber144]

In cv_varsel.refmodel(), the observation-specific weights are hard-coded to be all equal (line

projpred/R/cv_varsel.R

Line 221 in 0cabd73

w <- rep(1, NCOL(solution_terms_cv_ch))

), but in case of subsampling LOO-CV, they might be different (as the comment says there, too). Thus, shouldn't that line be modified?

[AlejandroCatalina]

Yes. This is a bit less straightforward to fix as we can compute both approximate LOO or full LOO, so I think the right weights should be returned by loo_varsel and kfold_varsel themselves. I'll work on this after resolving the easier issues.

[fweber144]

Related comments/discussions may be found in issue #173 (enumeration point 4) and PR #188. Now that I understand the purpose of w better, I think constant weights for LOO-CV are ok (no matter whether subsampled or not) but they might be inappropriate for K-fold CV with folds of very different sizes.

[fweber144]

I have now read Magnusson et al. (2019) and Magnusson et al. (2020), but I'm still confused how projpred wants to perform the subsampled LOO CV. More importantly, I have the strong feeling that the current implementation is not correct. What I learned from Magnusson et al. (2019) and Magnusson et al. (2020) makes me think the problems with the implementation in projpred get deeper than I initially thought.

What is important for projpred is that:

1. Magnusson et al. (2019) propose to perform probability-proportional-to-size sampling (PPS) combined with the Hansen-Hurwitz (HH) estimator.
2. Magnusson et al. (2020) propose to perform simple random sampling (SRS) combined with the difference estimator. The reason for this new proposal is that the approach by Magnusson et al. (2019) is not optimal for model comparisons (as explained by Magnusson et al., 2020).

However, in my opinion, neither the approach by Magnusson et al. (2019) nor that by Magnusson et al. (2020) is currently used in projpred. The key equation from Magnusson et al. (2019) is equation (4). The key equation from Magnusson et al. (2020) is equation (7). However, I can't find any of these two equations in projpred. The appropriate place would probably be in lines

projpred/R/summary_funs.R

Lines 171 to 179 in dd4bc36

	if (!is.null(lppd.bs)) {
	value <- sum((lppd - lppd.bs) * weights, na.rm = TRUE)
	value.se <- weighted.sd(lppd - lppd.bs, weights,
	na.rm = TRUE) / sqrt(n_notna) * n_notna
	} else {
	value <- sum(lppd * weights, na.rm = TRUE)
	value.se <- weighted.sd(lppd, weights,
	na.rm = TRUE) / sqrt(n_notna) * n_notna
	}

(these lines only refer to the ELPD as performance statistic; projpred has more performance statistics, of course), but this does not correspond to equation (4) from Magnusson et al. (2019) and neither to equation (7) from Magnusson et al. (2020). The SE estimators from these lines also don't correspond to the SE estimators from Magnusson et al. (2019, equation (6)) and Magnusson et al. (2020, equation (9)). The subsampling itself (performed in lines

projpred/R/cv_varsel.R

Lines 839 to 840 in dd4bc36

	w <- exp(lppd - max(lppd))
	inds <- sample(seq_along(lppd), size = nloo, prob = w)

) also confuses me. For example, for the approach by Magnusson et al. (2019), I would have expected the logarithms of the predictive densities (not the predictive densities themselves) in the definition of w, and then sampling with replacement.

So if the current implementation of subsampled LOO CV in projpred is not correct, then I think projpred should not support subsampled LOO CV until this is fixed.

Of course, I could be totally wrong and neither Magnusson et al. (2019) nor Magnusson et al. (2020) are relevant for subsampled LOO CV in projpred. But then, I would be even more confused because there is this comment:

projpred/R/cv_varsel.R

Lines 821 to 825 in dd4bc36

	## decide which points to go through in the validation based on
	## proportional-to-size subsampling as implemented in Magnusson, M., Riis
	## Andersen, M., Jonasson, J. and Vehtari, A. (2019). Leave-One-Out
	## Cross-Validation for Large Data. In International Conference on Machine
	## Learning.

which indicates that projpred wants to use the approach by Magnusson et al. (2019).

If my assumption is correct and projpred needs a fix, then a way to implement either the approach by Magnusson et al. (2019) or that by Magnusson et al. (2020) might be to rely on loo::loo_subsample() (see this vignette). This is just a rough idea, though, and I don't know if it is really feasible.

Of course, given these new findings, my statement

[...] I think constant weights for LOO-CV are ok (no matter whether subsampled or not) [...]

from the comment above (see also enumeration point 4 in #173) is not correct. At that time, I hadn't realized yet that projpred (probably) wants to use the approach by Magnusson et al. (2019) which requires sampling with replacement, which in turn means that each observation can occur more than once in the subsample (in contrast to my statement

[...] after all, they all represent a single observation [...]

from #188) so that when aggregating across the subsampled LOO CV folds, non-constant weights make sense (but note that these would have to be inverse-probability weights, in contrast to the currently used non-inversed probability weights).

References

Magnusson, M., Andersen, M., Jonasson, J., and Vehtari, A. (2019). Bayesian leave-one-out cross-validation for large data. In Proceedings of the 36th International Conference on Machine Learning, 4244-4253. URL: https://proceedings.mlr.press/v97/magnusson19a.html.

Magnusson, M., Vehtari, A., Jonasson, J., and Andersen, M. (2020). Leave-one-out cross-validation for Bayesian model comparison in large data. In Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, 341-351. URL: https://proceedings.mlr.press/v108/magnusson20a.html.

[fweber144]

An alternative to not supporting subsampled LOO CV at all anymore might be to throw a warning concerning a potentially incorrect implementation.

[fweber144]

See also

I also noticed that <vsel_object>$summaries$ref$w (so w for the reference model) is always NULL which is strange because it looks like lines

projpred/R/summary_funs.R

Lines 76 to 78 in 33047d6

	summ <- summ_ref
	res <- get_stat(summ$mu, summ$lppd, varsel$d_test, stat, mu.bs = mu.bs,
	lppd.bs = lppd.bs, weights = summ$w, alpha = alpha, ...)

expect a non-NULL element w also for the reference model (which would also make sense if w is supposed to contain the weights of the CV folds the observations are in because the reference model's performance is also cross-validated if the submodels' performances are).

in this comment. When fixing subsampled LOO CV, I guess that also needs to be fixed.

[fweber144]

More issues/questions concerning subsampled LOO CV:

1. NAs in get_stat()'s baseline objects mu.bs and lppd.bs are not handled consistently across the different stats (I think the handling from RMSE and AUC also needs to be used for all other stats). Furthermore, I think n_notna needs to be updated according to the NAs in the baseline objects mu.bs and lppd.bs if baseline = "best". I have now (PR Fix #349 #350) also added n_notna.bs which could perhaps be re-used.
2. For the AUC, the modified mu doesn't even seem to be used (see lines

   projpred/R/summary_funs.R

   Lines 295 to 307 in e7080db

   	auc.data <- cbind(y, mu, wcv)
   	if (!is.null(mu.bs)) {
   	mu.bs[is.na(mu)] <- NA # compute the AUCs using only those points
   	mu[is.na(mu.bs)] <- NA # for which both mu and mu.bs are non-NA
   	auc.data.bs <- cbind(y, mu.bs, wcv)
   	value <- auc(auc.data) - auc(auc.data.bs)
   	value.bootstrap1 <- bootstrap(auc.data, auc, ...)
   	value.bootstrap2 <- bootstrap(auc.data.bs, auc, ...)
   	value.se <- sd(value.bootstrap1 - value.bootstrap2, na.rm = TRUE)
   	lq_uq <- quantile(value.bootstrap1 - value.bootstrap2,
   	probs = c(alpha_half, one_minus_alpha_half),
   	names = FALSE, na.rm = TRUE)
   	} else {

   ).
3. It seems like loo_varsel()$d_test always contains all observations as test points, but line

   projpred/R/summary_funs.R

   Line 103 in e7080db

   res_diff <- get_stat(summ$mu, summ$lppd, varsel$d_test, stat,

   doesn't reduce this down to the subsampled observations. Is this on purpose?