Fix loo_compare() for subsampling results

NEWS.md

@@ -1,5 +1,6 @@
 # loo (development version)
 
+* Fix `loo_compare()` when used with subsampling: compute model comparison by comparison-model minus reference-model by @florence-bockting in #391
 * Update user messages in `print()` by @ishaan-arora-1, @florence-bockting in 
 #328.
 

R/loo_compare.psis_loo_ss_list.R

@@ -16,7 +16,7 @@ loo_compare.psis_loo_ss_list <- function(x, ...) {
 
   comp <- loo_compare_matrix.psis_loo_ss_list(x)
   ord <- loo_compare_order(x)
-  names(x) <- rownames(comp)[ord]
+  names(x) <- rownames(comp)[order(ord)]
 
   rnms <- rownames(comp)
   elpd_diff_mat <- matrix(0, nrow = nrow(comp), ncol = 3,
@@ -86,13 +86,13 @@ loo_compare_ss_naive <- function(ref_loo, compare_loo){
   checkmate::assert_class(ref_loo[[1]], "psis_loo_ss")
   checkmate::assert_class(compare_loo[[1]], "psis_loo_ss")
 
-  elpd_loo_diff <- ref_loo[[1]]$estimates["elpd_loo","Estimate"] - compare_loo[[1]]$estimates["elpd_loo","Estimate"]
+  elpd_loo_diff <- compare_loo[[1]]$estimates["elpd_loo", "Estimate"] - ref_loo[[1]]$estimates["elpd_loo", "Estimate"]
   elpd_loo_diff_se <- sqrt(
-    (ref_loo[[1]]$estimates["elpd_loo","SE"])^2 +
-    (compare_loo[[1]]$estimates["elpd_loo","SE"])^2)
+    (ref_loo[[1]]$estimates["elpd_loo", "SE"])^2 +
+    (compare_loo[[1]]$estimates["elpd_loo", "SE"])^2)
   elpd_loo_diff_subsampling_se <- sqrt(
-      (ref_loo[[1]]$estimates["elpd_loo","subsampling SE"])^2 +
-      (compare_loo[[1]]$estimates["elpd_loo","subsampling SE"])^2)
+      (ref_loo[[1]]$estimates["elpd_loo", "subsampling SE"])^2 +
+      (compare_loo[[1]]$estimates["elpd_loo", "subsampling SE"])^2)
 
   c(elpd_loo_diff, elpd_loo_diff_se, elpd_loo_diff_subsampling_se)
 }
@@ -112,8 +112,8 @@ loo_compare_ss_diff <- function(ref_loo, compare_loo){
   checkmate::assert_true(ref_loo[[1]]$loo_subsampling$loo_approximation != "none")
   checkmate::assert_true(compare_loo[[1]]$loo_subsampling$loo_approximation != "none")
 
-  diff_approx <- ref_loo[[1]]$loo_subsampling$elpd_loo_approx - compare_loo[[1]]$loo_subsampling$elpd_loo_approx
-  diff_sample <- ref_loo[[1]]$pointwise[,"elpd_loo"] - compare_loo[[1]]$pointwise[,"elpd_loo"]
+  diff_approx <- compare_loo[[1]]$loo_subsampling$elpd_loo_approx - ref_loo[[1]]$loo_subsampling$elpd_loo_approx
+  diff_sample <- compare_loo[[1]]$pointwise[,"elpd_loo"] - ref_loo[[1]]$pointwise[,"elpd_loo"]
   est <- srs_diff_est(diff_approx, y = diff_sample, y_idx = ref_loo[[1]]$pointwise[,"idx"])
 
   elpd_loo_diff <- est$y_hat
@@ -182,7 +182,6 @@ print.compare.loo_ss <- function(x, ..., digits = 1) {
   invisible(x)
 }
 
-
 #' Compute comparison matrix for `psis_loo_ss` objects
 #' @noRd
 #' @keywords internal

tests/testthat/_snaps/loo_subsampling_cases.md

@@ -105,33 +105,33 @@
       print(comp)
     Output
              elpd_diff se_diff subsampling_se_diff
-      model2  0.0       0.0     0.0               
-      model1 16.5      22.5     0.4               
+      model2   0.0       0.0     0.0              
+      model1 -16.5      22.5     0.4              
 
 ---
 
     Code
       print(comp)
     Output
              elpd_diff se_diff subsampling_se_diff
-      model2  0.0       0.0     0.0               
-      model1 16.1       4.4     0.1               
+      model2   0.0       0.0     0.0              
+      model1 -16.1       4.4     0.1              
 
 ---
 
     Code
       print(comp2)
     Output
              elpd_diff se_diff subsampling_se_diff
-      model2  0.0       0.0     0.0               
-      model1 16.3       4.4     0.1               
+      model2   0.0       0.0     0.0              
+      model1 -16.3       4.4     0.1              
 
 ---
 
     Code
       print(comp3)
     Output
              elpd_diff se_diff subsampling_se_diff
-      model2  0.0       0.0     0.0               
-      model1 16.5       4.4     0.3               
+      model2   0.0       0.0     0.0              
+      model1 -16.5       4.4     0.3              
 

tests/testthat/test_loo_subsampling.R

@@ -1071,6 +1071,11 @@ test_that("loo_compare_subsample", {
   expect_equal(lcssohhapi, lcssohh)
   expect_silent(lcss2mapi <- loo_compare(lss2o1, lss3o1))
   expect_equal(lcss2mapi, lcss2m)
+  # check that comparison is comp - ref model (i.e., elpd_diff is neg.)
+  for (m in list(lcss, lcss2, lcssohh)) {
+    expect_lt(m[2, "elpd_diff"], 0)
+    expect_lt(m[3, "elpd_diff"], 0)
+  }
 })
 
 test_that("Test 'tis' and 'sis'", {
@@ -1249,3 +1254,17 @@ test_that("Test 'tis' and 'sis'", {
     loo_ss_full$estimates["elpd_loo", "Estimate"]
   )
 })
+
+# subsampling should agree with non-subsampling when nothing
+# is actually subsampled (pins the direction *and* the magnitude)
+LL <- example_loglik_array()
+l1 <- loo(LL)
+l2 <- loo(LL + 1)
+l3 <- loo(LL + 2)
+l1ss <- loo:::as.psis_loo_ss.psis_loo(l1)
+l2ss <- loo:::as.psis_loo_ss.psis_loo(l2)
+l3ss <- loo:::as.psis_loo_ss.psis_loo(l3)
+expect_equal(
+  unname(loo_compare(l1ss, l2ss, l3ss)[, "elpd_diff"]),
+  loo_compare(l1, l2, l3)$elpd_diff
+)

vignettes/loo2-large-data.Rmd

@@ -481,9 +481,9 @@ print(comp)
 ```
 Warning: Different subsamples in 'model2' and 'model1'. Naive diff SE is used.
 
-       elpd_diff se_diff subsampling_se_diff
-model2  0.0       0.0     0.0               
-model1 16.5      22.5     0.4               
+        elpd_diff se_diff subsampling_se_diff
+model2    0.0      0.0     0.0               
+model1  -16.5     22.5     0.4               
 ```
 
 This new object `comp` contains the estimated difference of expected
@@ -538,23 +538,23 @@ subsampled observations.
 
 ```{r, eval=FALSE}
 comp <- loo_compare(loo_ss_1, loo_ss_2)
-print(comp) 
+print(comp)
 ```
 
 ```
-       elpd_diff se_diff subsampling_se_diff
-model2  0.0       0.0     0.0               
-model1 16.1       4.4     0.1               
+        elpd_diff se_diff subsampling_se_diff
+model2    0.0      0.0     0.0               
+model1  -16.1      4.4     0.1               
 ```
 
 First, notice that now the `se_diff` is now around 4 (as opposed to 20 when using
-different subsamples). The first column shows the difference in ELPD relative to
-the model with the largest ELPD. In this case, the difference in `elpd` and its
-scale relative to the approximate standard error of the difference) indicates a
-preference for the second model (`model2`). Since the subsampling uncertainty is
-so small in this case it can effectively be ignored. If we need larger
-subsamples we can simply add samples using the `update()` method demonstrated
-earlier.
+different subsamples). The `elpd_diff` column shows the difference in ELPD
+relative to the model with the largest ELPD. In this case, the difference in
+`elpd` and its scale relative to the approximate standard error of the
+difference) indicates a preference for the second model (`model2`).
+Since the subsampling uncertainty is so small in this case it can effectively
+be ignored. If we need larger subsamples we can simply add samples using the
+`update()` method demonstrated earlier.
 
 It is also possible to compare a subsampled loo computation with a full loo object.
 
@@ -580,9 +580,9 @@ the loo calculations for both `model1` and `model2` are included in the
 computations for the comparison.
 
 ```
-       elpd_diff se_diff subsampling_se_diff
-model2  0.0       0.0     0.0               
-model1 16.3       4.4     0.3   
+        elpd_diff se_diff subsampling_se_diff
+model2    0.0     0.0     0.0               
+model1  -16.3     4.4     0.3   
 ```
 
 Here we actually see an increase in `subsampling_se_diff`, but this is due to a