diff --git a/NEWS.md b/NEWS.md
index 291db72b..854c05ce 100644
--- a/NEWS.md
+++ b/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.
 
diff --git a/R/loo_compare.psis_loo_ss_list.R b/R/loo_compare.psis_loo_ss_list.R
index 5ef5f202..194eb69f 100644
--- a/R/loo_compare.psis_loo_ss_list.R
+++ b/R/loo_compare.psis_loo_ss_list.R
@@ -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
diff --git a/tests/testthat/_snaps/loo_subsampling_cases.md b/tests/testthat/_snaps/loo_subsampling_cases.md
index 677e4fa1..4cff4eb9 100644
--- a/tests/testthat/_snaps/loo_subsampling_cases.md
+++ b/tests/testthat/_snaps/loo_subsampling_cases.md
@@ -105,8 +105,8 @@
       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              
 
 ---
 
@@ -114,8 +114,8 @@
       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              
 
 ---
 
@@ -123,8 +123,8 @@
       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              
 
 ---
 
@@ -132,6 +132,6 @@
       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              
 
diff --git a/tests/testthat/test_loo_subsampling.R b/tests/testthat/test_loo_subsampling.R
index 84e8d785..bf566ff7 100644
--- a/tests/testthat/test_loo_subsampling.R
+++ b/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'", {
diff --git a/vignettes/loo2-large-data.Rmd b/vignettes/loo2-large-data.Rmd
index fbaa5827..7f7083f0 100644
--- a/vignettes/loo2-large-data.Rmd
+++ b/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