fix: correct elpd_diff sign and return a data.frame in loo_compare() for subsampled loo objects

Author: Florence Bockting

NEWS.md

@@ -1,5 +1,7 @@
 # loo (development version)
 
+* Fix `loo_compare()` when used with subsampling: compute model comparison by comparison-model minus reference-model and change output structure from matrix
+to data.frame by @florence-bockting in #TODO
 * Update user messages in `print()` by @ishaan-arora-1, @florence-bockting in 
 #328.
 

R/loo_compare.psis_loo_ss_list.R

@@ -24,8 +24,12 @@ loo_compare.psis_loo_ss_list <- function(x, ...) {
   for(i in 2:length(ord)){
     elpd_diff_mat[i,] <- loo_compare_ss(ref_loo = x[ord[1]], compare_loo = x[ord[i]])
   }
-  comp <- cbind(elpd_diff_mat, comp)
-  rownames(comp) <- rnms
+  comp <- cbind(
+    data.frame(model = rnms, stringsAsFactors = FALSE),
+    as.data.frame(elpd_diff_mat),
+    as.data.frame(comp)
+  )
+  rownames(comp) <- NULL
 
   class(comp) <- c("compare.loo_ss", "compare.loo", class(comp))
   return(comp)
@@ -86,13 +90,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 +116,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
@@ -174,15 +178,16 @@ loo_compare_checks.psis_loo_ss_list <- function(loos) {
 #' @rdname loo_compare
 #' @export
 print.compare.loo_ss <- function(x, ..., digits = 1) {
-  xcopy <- x
-  if (NCOL(xcopy) >= 2) {
-    xcopy <- xcopy[, c("elpd_diff", "se_diff", "subsampling_se_diff")]
-  }
-  print(.fr(xcopy, digits), quote = FALSE)
+  cols <- c("model", "elpd_diff", "se_diff", "subsampling_se_diff")
+  cols <- intersect(cols, colnames(x))
+
+  x_sub <- x[, cols, drop = FALSE]
+  x_sub[setdiff(cols, "model")] <- .fr(x_sub[setdiff(cols, "model")], digits)
+  print(as.data.frame(x_sub), quote = FALSE, row.names = FALSE)
+
   invisible(x)
 }
 
-
 #' Compute comparison matrix for `psis_loo_ss` objects
 #' @noRd
 #' @keywords internal

tests/testthat/_snaps/loo_subsampling_cases.md

@@ -104,34 +104,34 @@
     Code
       print(comp)
     Output
-             elpd_diff se_diff subsampling_se_diff
-      model2  0.0       0.0     0.0               
-      model1 16.5      22.5     0.4               
+        model elpd_diff se_diff subsampling_se_diff
+       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               
+        model elpd_diff se_diff subsampling_se_diff
+       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               
+        model elpd_diff se_diff subsampling_se_diff
+       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               
+        model elpd_diff se_diff subsampling_se_diff
+       model2       0.0     0.0                 0.0
+       model1     -16.5     4.4                 0.3
 

tests/testthat/test_loo_subsampling.R

@@ -1047,30 +1047,39 @@ test_that("loo_compare_subsample", {
     )
   )
 
-  expect_equal(lcss[, 1], lcsso[, 1], tolerance = 1)
-  expect_equal(lcss2[, 1], lcsso[, 1], tolerance = 1)
-  expect_equal(lcssohh[, 1], lcsso[, 1], tolerance = 1)
-  expect_equal(lcssf1[, 1], lcsso[, 1], tolerance = 1)
-  expect_equal(lcssf2[, 1], lcsso[, 1], tolerance = 1)
+  expect_equal(lcss$elpd_diff, lcsso$elpd_diff, tolerance = 1)
+  expect_equal(lcss2$elpd_diff, lcsso$elpd_diff, tolerance = 1)
+  expect_equal(lcssohh$elpd_diff, lcsso$elpd_diff, tolerance = 1)
+  expect_equal(lcssf1$elpd_diff, lcsso$elpd_diff, tolerance = 1)
+  expect_equal(lcssf2$elpd_diff, lcsso$elpd_diff, tolerance = 1)
 
-  expect_gt(lcss[, 2][2], lcsso[, 2][2])
-  expect_gt(lcss[, 2][3], lcsso[, 2][3])
-  expect_gt(lcss2[, 2][2], lcsso[, 2][2])
-  expect_equal(lcss2[, 2][3], lcsso[, 2][3])
-  expect_gt(lcssohh[, 2][2], lcsso[, 2][2])
-  expect_equal(lcssohh[, 2][3], lcsso[, 2][3])
+  expect_gt(lcss$se_diff[2], lcsso$se_diff[2])
+  expect_gt(lcss$se_diff[3], lcsso$se_diff[3])
+  expect_gt(lcss2$se_diff[2], lcsso$se_diff[2])
+  expect_equal(lcss2$se_diff[3], lcsso$se_diff[3])
+  expect_gt(lcssohh$se_diff[2], lcsso$se_diff[2])
+  expect_equal(lcssohh$se_diff[3], lcsso$se_diff[3])
 
   expect_silent(
     lcss2m <- loo:::loo_compare.psis_loo_ss_list(x = list(lss2o1, lss3o1))
   )
-  expect_equal(unname(lcss2m[,]), unname(lcsso[1:2, ]))
+  expect_equal(
+    lcss2m[, setdiff(colnames(lcss2m), "model")],
+    lcsso[1:2, setdiff(colnames(lcsso), "model")]
+  )
 
   expect_snapshot(lcssapi <- loo_compare(lss1, lss2, lss3))
   expect_equal(lcssapi, lcss)
   expect_warning(lcssohhapi <- loo_compare(lss1, lss2hh, lss3o1))
   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$elpd_diff[2], 0)
+    expect_lt(m$elpd_diff[3], 0)
+    expect_true("data.frame" %in% class(m))
+  }
 })
 
 test_that("Test 'tis' and 'sis'", {

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               
+model   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               
+model   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   
+model   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