diff --git a/tests/testthat/test_0_helpers.R b/tests/testthat/test_0_helpers.R
index d9523b21..e17f3ec1 100644
--- a/tests/testthat/test_0_helpers.R
+++ b/tests/testthat/test_0_helpers.R
@@ -1,5 +1,3 @@
-library(loo)
-
 LLarr <- example_loglik_array()
 LLmat <- example_loglik_matrix()
 
@@ -24,18 +22,26 @@ test_that("reshaping functions result in correct dimensions", {
 })
 
 test_that("reshaping functions throw correct errors", {
-  expect_error(llmatrix_to_array(LLmat, chain_id = rep(1:2, times = c(400, 600))),
-               regexp = "Not all chains have same number of iterations",
-               fixed = TRUE)
-  expect_error(llmatrix_to_array(LLmat, chain_id = rep(1:2, each = 400)),
-               regexp = "Number of rows in matrix not equal to length(chain_id)",
-               fixed = TRUE)
-  expect_error(llmatrix_to_array(LLmat, chain_id = rep(2:3, each = 500)),
-               regexp = "max(chain_id) not equal to the number of chains",
-               fixed = TRUE)
-  expect_error(llmatrix_to_array(LLmat, chain_id = rnorm(1000)),
-               regexp = "all(chain_id == as.integer(chain_id)) is not TRUE",
-               fixed = TRUE)
+  expect_error(
+    llmatrix_to_array(LLmat, chain_id = rep(1:2, times = c(400, 600))),
+    regexp = "Not all chains have same number of iterations",
+    fixed = TRUE
+  )
+  expect_error(
+    llmatrix_to_array(LLmat, chain_id = rep(1:2, each = 400)),
+    regexp = "Number of rows in matrix not equal to length(chain_id)",
+    fixed = TRUE
+  )
+  expect_error(
+    llmatrix_to_array(LLmat, chain_id = rep(2:3, each = 500)),
+    regexp = "max(chain_id) not equal to the number of chains",
+    fixed = TRUE
+  )
+  expect_error(
+    llmatrix_to_array(LLmat, chain_id = rnorm(1000)),
+    regexp = "all(chain_id == as.integer(chain_id)) is not TRUE",
+    fixed = TRUE
+  )
 })
 
 test_that("colLogMeanExps(x) = log(colMeans(exp(x))) ", {
@@ -54,9 +60,14 @@ test_that("validating log-lik objects and functions works", {
 })
 
 test_that("nlist works", {
-  a <- 1; b <- 2; c <- 3;
+  a <- 1
+  b <- 2
+  c <- 3
   nlist_val <- list(nlist(a, b, c), nlist(a, b, c = "tornado"))
-  nlist_ans <- list(list(a = 1, b = 2, c = 3), list(a = 1, b = 2, c = "tornado"))
+  nlist_ans <- list(
+    list(a = 1, b = 2, c = 3),
+    list(a = 1, b = 2, c = "tornado")
+  )
   expect_equal(nlist_val, nlist_ans)
   expect_equal(nlist(a = 1, b = 2, c = 3), list(a = 1, b = 2, c = 3))
 })
@@ -69,6 +80,5 @@ test_that("loo_cores works", {
 
   options(loo.cores = 2)
   expect_warning(expect_equal(loo_cores(10), 2), "deprecated")
-  options(loo.cores=NULL)
+  options(loo.cores = NULL)
 })
-
diff --git a/tests/testthat/test_E_loo.R b/tests/testthat/test_E_loo.R
index cfbacf1e..f7acaf04 100644
--- a/tests/testthat/test_E_loo.R
+++ b/tests/testthat/test_E_loo.R
@@ -1,5 +1,3 @@
-library(loo)
-
 LLarr <- example_loglik_array()
 LLmat <- example_loglik_matrix()
 LLvec <- LLmat[, 1]
@@ -17,15 +15,54 @@ log_rats <- -LLmat
 E_test_mean <- E_loo(x, psis_mat, type = "mean", log_ratios = log_rats)
 E_test_var <- E_loo(x, psis_mat, type = "var", log_ratios = log_rats)
 E_test_sd <- E_loo(x, psis_mat, type = "sd", log_ratios = log_rats)
-E_test_quant <- E_loo(x, psis_mat, type = "quantile", probs = 0.5, log_ratios = log_rats)
-E_test_quant2 <- E_loo(x, psis_mat, type = "quantile", probs = c(0.1, 0.9), log_ratios = log_rats)
+E_test_quant <- E_loo(
+  x,
+  psis_mat,
+  type = "quantile",
+  probs = 0.5,
+  log_ratios = log_rats
+)
+E_test_quant2 <- E_loo(
+  x,
+  psis_mat,
+  type = "quantile",
+  probs = c(0.1, 0.9),
+  log_ratios = log_rats
+)
 
 # vector method
-E_test_mean_vec <- E_loo(x[, 1], psis_vec, type = "mean", log_ratios = log_rats[,1])
-E_test_var_vec <- E_loo(x[, 1], psis_vec, type = "var", log_ratios = log_rats[,1])
-E_test_sd_vec <- E_loo(x[, 1], psis_vec, type = "sd", log_ratios = log_rats[,1])
-E_test_quant_vec <- E_loo(x[, 1], psis_vec, type = "quant", probs = 0.5, log_ratios = log_rats[,1])
-E_test_quant_vec2 <- E_loo(x[, 1], psis_vec, type = "quant", probs = c(0.1, 0.5, 0.9), log_ratios = log_rats[,1])
+E_test_mean_vec <- E_loo(
+  x[, 1],
+  psis_vec,
+  type = "mean",
+  log_ratios = log_rats[, 1]
+)
+E_test_var_vec <- E_loo(
+  x[, 1],
+  psis_vec,
+  type = "var",
+  log_ratios = log_rats[, 1]
+)
+E_test_sd_vec <- E_loo(
+  x[, 1],
+  psis_vec,
+  type = "sd",
+  log_ratios = log_rats[, 1]
+)
+E_test_quant_vec <- E_loo(
+  x[, 1],
+  psis_vec,
+  type = "quant",
+  probs = 0.5,
+  log_ratios = log_rats[, 1]
+)
+E_test_quant_vec2 <- E_loo(
+  x[, 1],
+  psis_vec,
+  type = "quant",
+  probs = c(0.1, 0.5, 0.9),
+  log_ratios = log_rats[, 1]
+)
 
 # E_loo_khat
 khat <- loo:::E_loo_khat.matrix(x, psis_mat, log_rats)
@@ -114,11 +151,11 @@ test_that("E_loo throws correct errors and warnings", {
   # warnings
   expect_no_warning(E_loo.matrix(x, psis_mat))
   # no warnings if x is constant, binary, NA, NaN, Inf
-  expect_no_warning(E_loo.matrix(x*0, psis_mat))
-  expect_no_warning(E_loo.matrix(0+(x>0), psis_mat))
-  expect_no_warning(E_loo.matrix(x+NA, psis_mat))   
-  expect_no_warning(E_loo.matrix(x*NaN, psis_mat))
-  expect_no_warning(E_loo.matrix(x*Inf, psis_mat))
+  expect_no_warning(E_loo.matrix(x * 0, psis_mat))
+  expect_no_warning(E_loo.matrix(0 + (x > 0), psis_mat))
+  expect_no_warning(E_loo.matrix(x + NA, psis_mat))
+  expect_no_warning(E_loo.matrix(x * NaN, psis_mat))
+  expect_no_warning(E_loo.matrix(x * Inf, psis_mat))
   expect_no_warning(E_test <- E_loo.default(x[, 1], psis_vec))
   expect_length(E_test$pareto_k, 1)
 
@@ -161,7 +198,6 @@ test_that("weighted quantiles work", {
     quantile(xx, probs, names = FALSE)
   }
 
-
   set.seed(123)
   pr <- seq(0.025, 0.975, 0.025)
 
diff --git a/tests/testthat/test_compare.R b/tests/testthat/test_compare.R
index e70a25f9..256c263c 100644
--- a/tests/testthat/test_compare.R
+++ b/tests/testthat/test_compare.R
@@ -1,4 +1,3 @@
-library(loo)
 set.seed(123)
 
 LLarr <- example_loglik_array()
@@ -12,9 +11,11 @@ test_that("loo_compare throws appropriate errors", {
   w4 <- suppressWarnings(waic(LLarr[,, -(1:2)]))
 
   expect_error(loo_compare(2, 3), "must be a list if not a 'loo' object")
-  expect_error(loo_compare(w1, w2, x = list(w1, w2)),
-               "If 'x' is a list then '...' should not be specified")
-  expect_error(loo_compare(w1, list(1,2,3)), "class 'loo'")
+  expect_error(
+    loo_compare(w1, w2, x = list(w1, w2)),
+    "If 'x' is a list then '...' should not be specified"
+  )
+  expect_error(loo_compare(w1, list(1, 2, 3)), "class 'loo'")
   expect_error(loo_compare(w1), "requires at least two models")
   expect_error(loo_compare(x = list(w1)), "requires at least two models")
   expect_error(loo_compare(w1, w3), "same number of data points")
@@ -22,35 +23,51 @@ test_that("loo_compare throws appropriate errors", {
 })
 
 test_that("loo_compare throws appropriate warnings", {
-  w3 <- w1; w4 <- w2
+  w3 <- w1
+  w4 <- w2
   class(w3) <- class(w4) <- c("kfold", "loo")
   attr(w3, "K") <- 2
   attr(w4, "K") <- 3
-  expect_warning(loo_compare(w3, w4), "Not all kfold objects have the same K value")
+  expect_warning(
+    loo_compare(w3, w4),
+    "Not all kfold objects have the same K value"
+  )
 
   class(w4) <- c("psis_loo", "loo")
   attr(w4, "K") <- NULL
   expect_warning(loo_compare(w3, w4), "Comparing LOO-CV to K-fold-CV")
 
-  w3 <- w1; w4 <- w2
+  w3 <- w1
+  w4 <- w2
   attr(w3, "yhash") <- "a"
   attr(w4, "yhash") <- "b"
   expect_warning(loo_compare(w3, w4), "Not all models have the same y variable")
 
   set.seed(123)
-  w_list <- lapply(1:25, function(x) suppressWarnings(waic(LLarr + rnorm(1, 0, 0.1))))
-  expect_warning(loo_compare(w_list),
-                 "Difference in performance potentially due to chance")
-
-  w_list_short <- lapply(1:4, function(x) suppressWarnings(waic(LLarr + rnorm(1, 0, 0.1))))
+  w_list <- lapply(1:25, function(x) {
+    suppressWarnings(waic(LLarr + rnorm(1, 0, 0.1)))
+  })
+  expect_warning(
+    loo_compare(w_list),
+    "Difference in performance potentially due to chance"
+  )
+
+  w_list_short <- lapply(1:4, function(x) {
+    suppressWarnings(waic(LLarr + rnorm(1, 0, 0.1)))
+  })
   expect_no_warning(loo_compare(w_list_short))
 })
 
 
-
 comp_colnames <- c(
-  "elpd_diff", "se_diff", "elpd_waic", "se_elpd_waic",
-  "p_waic", "se_p_waic", "waic", "se_waic"
+  "elpd_diff",
+  "se_diff",
+  "elpd_waic",
+  "se_elpd_waic",
+  "p_waic",
+  "se_p_waic",
+  "waic",
+  "se_waic"
 )
 
 test_that("loo_compare returns expected results (2 models)", {
@@ -59,15 +76,15 @@ test_that("loo_compare returns expected results (2 models)", {
   expect_equal(colnames(comp1), comp_colnames)
   expect_equal(rownames(comp1), c("model1", "model2"))
   expect_output(print(comp1), "elpd_diff")
-  expect_equal(comp1[1:2,1], c(0, 0), ignore_attr = TRUE)
-  expect_equal(comp1[1:2,2], c(0, 0), ignore_attr = TRUE)
+  expect_equal(comp1[1:2, 1], c(0, 0), ignore_attr = TRUE)
+  expect_equal(comp1[1:2, 2], c(0, 0), ignore_attr = TRUE)
 
   comp2 <- loo_compare(w1, w2)
   expect_s3_class(comp2, "compare.loo")
   expect_equal(colnames(comp2), comp_colnames)
-  
+
   expect_snapshot_value(comp2, style = "serialize")
-  
+
   # specifying objects via ... and via arg x gives equal results
   expect_equal(comp2, loo_compare(x = list(w1, w2)))
 })
@@ -79,7 +96,7 @@ test_that("loo_compare returns expected result (3 models)", {
 
   expect_equal(colnames(comp1), comp_colnames)
   expect_equal(rownames(comp1), c("model1", "model2", "model3"))
-  expect_equal(comp1[1,1], 0)
+  expect_equal(comp1[1, 1], 0)
   expect_s3_class(comp1, "compare.loo")
   expect_s3_class(comp1, "matrix")
 
@@ -119,34 +136,53 @@ test_that("compare returns expected result (3 models)", {
   expect_equal(
     colnames(comp1),
     c(
-      "elpd_diff", "se_diff", "elpd_waic", "se_elpd_waic",
-      "p_waic", "se_p_waic", "waic", "se_waic"
-    ))
+      "elpd_diff",
+      "se_diff",
+      "elpd_waic",
+      "se_elpd_waic",
+      "p_waic",
+      "se_p_waic",
+      "waic",
+      "se_waic"
+    )
+  )
   expect_equal(rownames(comp1), c("w1", "w2", "w3"))
-  expect_equal(comp1[1,1], 0)
+  expect_equal(comp1[1, 1], 0)
   expect_s3_class(comp1, "compare.loo")
   expect_s3_class(comp1, "matrix")
   expect_snapshot_value(comp1, style = "serialize")
 
   # specifying objects via '...' gives equivalent results (equal
   # except rownames) to using 'x' argument
-  expect_warning(comp_via_list <- loo::compare(x = list(w1, w2, w3)), "Deprecated")
+  expect_warning(
+    comp_via_list <- loo::compare(x = list(w1, w2, w3)),
+    "Deprecated"
+  )
   expect_equal(comp1, comp_via_list, ignore_attr = TRUE)
 })
 
 test_that("compare throws appropriate errors", {
-  expect_error(suppressWarnings(loo::compare(w1, w2, x = list(w1, w2))),
-               "should not be specified")
-  expect_error(suppressWarnings(loo::compare(x = 2)),
-               "must be a list")
-  expect_error(suppressWarnings(loo::compare(x = list(2))),
-               "should have class 'loo'")
-  expect_error(suppressWarnings(loo::compare(x = list(w1))),
-               "requires at least two models")
-
-  w3 <- suppressWarnings(waic(LLarr2[,,-1]))
-  expect_error(suppressWarnings(loo::compare(x = list(w1, w3))),
-               "same number of data points")
-  expect_error(suppressWarnings(loo::compare(x = list(w1, w2, w3))),
-               "same number of data points")
+  expect_error(
+    suppressWarnings(loo::compare(w1, w2, x = list(w1, w2))),
+    "should not be specified"
+  )
+  expect_error(suppressWarnings(loo::compare(x = 2)), "must be a list")
+  expect_error(
+    suppressWarnings(loo::compare(x = list(2))),
+    "should have class 'loo'"
+  )
+  expect_error(
+    suppressWarnings(loo::compare(x = list(w1))),
+    "requires at least two models"
+  )
+
+  w3 <- suppressWarnings(waic(LLarr2[,, -1]))
+  expect_error(
+    suppressWarnings(loo::compare(x = list(w1, w3))),
+    "same number of data points"
+  )
+  expect_error(
+    suppressWarnings(loo::compare(x = list(w1, w2, w3))),
+    "same number of data points"
+  )
 })
diff --git a/tests/testthat/test_deprecated_extractors.R b/tests/testthat/test_deprecated_extractors.R
index 6f27767e..4cfda2f9 100644
--- a/tests/testthat/test_deprecated_extractors.R
+++ b/tests/testthat/test_deprecated_extractors.R
@@ -1,4 +1,3 @@
-library(loo)
 options(mc.cores = 1)
 set.seed(123)
 
diff --git a/tests/testthat/test_extract_log_lik.R b/tests/testthat/test_extract_log_lik.R
index bd55033e..4efddb39 100644
--- a/tests/testthat/test_extract_log_lik.R
+++ b/tests/testthat/test_extract_log_lik.R
@@ -1,5 +1,3 @@
-library(loo)
-
 test_that("extract_log_lik throws appropriate errors", {
   x1 <- rnorm(100)
   expect_error(extract_log_lik(x1), regexp = "Not a stanfit object")
diff --git a/tests/testthat/test_gpdfit.R b/tests/testthat/test_gpdfit.R
index a1b1c902..df40a914 100644
--- a/tests/testthat/test_gpdfit.R
+++ b/tests/testthat/test_gpdfit.R
@@ -1,22 +1,20 @@
-library(loo)
-
 test_that("gpdfit returns correct result", {
   set.seed(123)
   x <- rexp(100)
-  gpdfit_val_old <- unlist(gpdfit(x, wip=FALSE, min_grid_pts = 80))
+  gpdfit_val_old <- unlist(gpdfit(x, wip = FALSE, min_grid_pts = 80))
   expect_snapshot_value(gpdfit_val_old, style = "serialize")
 
-  gpdfit_val_wip <- unlist(gpdfit(x, wip=TRUE, min_grid_pts = 80))
+  gpdfit_val_wip <- unlist(gpdfit(x, wip = TRUE, min_grid_pts = 80))
   expect_snapshot_value(gpdfit_val_wip, style = "serialize")
 
-  gpdfit_val_wip_default_grid <- unlist(gpdfit(x, wip=TRUE))
+  gpdfit_val_wip_default_grid <- unlist(gpdfit(x, wip = TRUE))
   expect_snapshot_value(gpdfit_val_wip_default_grid, style = "serialize")
 })
 
 test_that("qgpd returns the correct result ", {
   probs <- seq(from = 0, to = 1, by = 0.25)
   q1 <- qgpd(probs, k = 1, sigma = 1)
-  expect_equal(q1, c(0, 1/3, 1, 3, Inf))
+  expect_equal(q1, c(0, 1 / 3, 1, 3, Inf))
 
   q2 <- qgpd(probs, k = 1, sigma = 0)
   expect_true(all(is.nan(q2)))
diff --git a/tests/testthat/test_kfold_helpers.R b/tests/testthat/test_kfold_helpers.R
index c35b192e..fd95de3e 100644
--- a/tests/testthat/test_kfold_helpers.R
+++ b/tests/testthat/test_kfold_helpers.R
@@ -1,4 +1,3 @@
-library(loo)
 set.seed(14014)
 
 test_that("kfold_split_random works", {
@@ -27,7 +26,9 @@ test_that("kfold_split_stratified works", {
   expect_silent(fold_strat <- kfold_split_stratified(5, y)) # used to be a warning before fixing issue #277
   tab <- table(fold_strat, y)
   expect_equal(tab[1, ], c("1" = 4, "2" = 8, "3" = 1))
-  for (i in 2:nrow(tab)) expect_equal(tab[i, ], c("1" = 4, "2" = 8, "3" = 0))
+  for (i in 2:nrow(tab)) {
+    expect_equal(tab[i, ], c("1" = 4, "2" = 8, "3" = 0))
+  }
 })
 
 test_that("kfold_split_grouped works", {
@@ -45,33 +46,93 @@ test_that("kfold_split_grouped works", {
   fold_group <- kfold_split_grouped(K = 10, x = grp)
   expect_equal(sum(table(fold_group)), length(grp) - 4)
 
-  grp <- rep(c("A","B"), each = 20)
+  grp <- rep(c("A", "B"), each = 20)
   fold_group <- kfold_split_grouped(K = 2, x = grp)
   expect_equal(fold_group, as.integer(as.factor(grp)))
 })
 
 test_that("kfold helpers throw correct errors", {
   expect_error(kfold_split_random(10), "!is.null(N) is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_random(10.5, 100), "K == as.integer(K) is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_random(10, 100.5), "N == as.integer(N) is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_random(K = c(1,1), N = 100), "length(K) == 1 is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_random(N = c(100, 100)), "length(N) == 1 is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_random(K = 5, N = 4), "K <= N is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_random(K = 1, N = 4), "K > 1 is not TRUE", fixed = TRUE)
+  expect_error(
+    kfold_split_random(10.5, 100),
+    "K == as.integer(K) is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_random(10, 100.5),
+    "N == as.integer(N) is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_random(K = c(1, 1), N = 100),
+    "length(K) == 1 is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_random(N = c(100, 100)),
+    "length(N) == 1 is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_random(K = 5, N = 4),
+    "K <= N is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_random(K = 1, N = 4),
+    "K > 1 is not TRUE",
+    fixed = TRUE
+  )
 
   y <- sample(c(0, 1), size = 200, replace = TRUE, prob = c(0.05, 0.95))
-  expect_error(kfold_split_stratified(10), "!is.null(x) is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_stratified(10.5, y), "K == as.integer(K) is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_stratified(K = c(1,1), y), "length(K) == 1 is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_stratified(K = 201, y), "K <= length(x) is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_stratified(K = 1, y), "K > 1 is not TRUE", fixed = TRUE)
+  expect_error(
+    kfold_split_stratified(10),
+    "!is.null(x) is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_stratified(10.5, y),
+    "K == as.integer(K) is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_stratified(K = c(1, 1), y),
+    "length(K) == 1 is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_stratified(K = 201, y),
+    "K <= length(x) is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_stratified(K = 1, y),
+    "K > 1 is not TRUE",
+    fixed = TRUE
+  )
 
   grp <- gl(n = 50, k = 15)
   expect_error(kfold_split_grouped(10), "!is.null(x) is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_grouped(3, c(1,1,1)), "'K' must not be bigger than the number of levels/groups in 'x'", fixed = TRUE)
-  expect_error(kfold_split_grouped(10.5, grp), "K == as.integer(K) is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_grouped(K = c(1,1), grp), "length(K) == 1 is not TRUE", fixed = TRUE)
-  expect_error(kfold_split_grouped(K = 1, grp), "K > 1 is not TRUE", fixed = TRUE)
+  expect_error(
+    kfold_split_grouped(3, c(1, 1, 1)),
+    "'K' must not be bigger than the number of levels/groups in 'x'",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_grouped(10.5, grp),
+    "K == as.integer(K) is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_grouped(K = c(1, 1), grp),
+    "length(K) == 1 is not TRUE",
+    fixed = TRUE
+  )
+  expect_error(
+    kfold_split_grouped(K = 1, grp),
+    "K > 1 is not TRUE",
+    fixed = TRUE
+  )
 })
 
 
@@ -82,4 +143,3 @@ test_that("print_dims.kfold works", {
   attr(xx, "K") <- NULL
   expect_silent(print_dims(xx))
 })
-
diff --git a/tests/testthat/test_loo_and_waic.R b/tests/testthat/test_loo_and_waic.R
index 663d9858..5cdc823e 100644
--- a/tests/testthat/test_loo_and_waic.R
+++ b/tests/testthat/test_loo_and_waic.R
@@ -1,4 +1,3 @@
-library(loo)
 options(mc.cores = 1)
 set.seed(123)
 
@@ -82,7 +81,10 @@ test_that("loo returns object with correct structure", {
   expect_named(loo1$diagnostics, c("pareto_k", "n_eff", "r_eff"))
   expect_equal(dimnames(loo1$estimates)[[1]], c("elpd_loo", "p_loo", "looic"))
   expect_equal(dimnames(loo1$estimates)[[2]], c("Estimate", "SE"))
-  expect_equal(colnames(loo1$pointwise), c("elpd_loo", "mcse_elpd_loo", "p_loo", "looic", "influence_pareto_k"))
+  expect_equal(
+    colnames(loo1$pointwise),
+    c("elpd_loo", "mcse_elpd_loo", "p_loo", "looic", "influence_pareto_k")
+  )
   expect_equal(dim(loo1), dim(LLmat))
 })
 
@@ -105,7 +107,10 @@ test_that("elpd returns object with correct structure", {
 
 
 test_that("two pareto k values are equal", {
-  expect_identical(loo1$pointwise[,"influence_pareto_k"], loo1$diagnostics$pareto_k)
+  expect_identical(
+    loo1$pointwise[, "influence_pareto_k"],
+    loo1$diagnostics$pareto_k
+  )
 })
 
 test_that("loo.array and loo.matrix give same result", {
@@ -141,32 +146,55 @@ test_that("loo, waic, and elpd error with vector input", {
 })
 
 
-
 # testing function methods ------------------------------------------------
 source(test_path("data-for-tests/function_method_stuff.R"))
 
 waic_with_fn <- waic(llfun, data = data, draws = draws)
 waic_with_mat <- waic(llmat_from_fn)
 
-loo_with_fn <- loo(llfun, data = data, draws = draws,
-                   r_eff = rep(1, nrow(data)))
-loo_with_mat <- loo(llmat_from_fn, r_eff = rep(1, ncol(llmat_from_fn)),
-                    save_psis = TRUE)
+loo_with_fn <- loo(
+  llfun,
+  data = data,
+  draws = draws,
+  r_eff = rep(1, nrow(data))
+)
+loo_with_mat <- loo(
+  llmat_from_fn,
+  r_eff = rep(1, ncol(llmat_from_fn)),
+  save_psis = TRUE
+)
 
 test_that("loo.cores deprecation warning works with function method", {
   options(loo.cores = 1)
-  expect_warning(loo(llfun, cores = 2, data = data, draws = draws, r_eff = rep(1, nrow(data))),
-                 "loo.cores")
-  options(loo.cores=NULL)
+  expect_warning(
+    loo(
+      llfun,
+      cores = 2,
+      data = data,
+      draws = draws,
+      r_eff = rep(1, nrow(data))
+    ),
+    "loo.cores"
+  )
+  options(loo.cores = NULL)
 })
 
 test_that("loo_i results match loo results for ith data point", {
   expect_no_warning(
     loo_i_val <- loo_i(i = 2, llfun = llfun, data = data, draws = draws),
   )
-  expect_equal(loo_i_val$pointwise[, "elpd_loo"], loo_with_fn$pointwise[2, "elpd_loo"])
-  expect_equal(loo_i_val$pointwise[, "p_loo"], loo_with_fn$pointwise[2, "p_loo"])
-  expect_equal(loo_i_val$diagnostics$pareto_k, loo_with_fn$diagnostics$pareto_k[2])
+  expect_equal(
+    loo_i_val$pointwise[, "elpd_loo"],
+    loo_with_fn$pointwise[2, "elpd_loo"]
+  )
+  expect_equal(
+    loo_i_val$pointwise[, "p_loo"],
+    loo_with_fn$pointwise[2, "p_loo"]
+  )
+  expect_equal(
+    loo_i_val$diagnostics$pareto_k,
+    loo_with_fn$diagnostics$pareto_k[2]
+  )
   expect_equal(loo_i_val$diagnostics$n_eff, loo_with_fn$diagnostics$n_eff[2])
 })
 
@@ -178,16 +206,31 @@ test_that("function and matrix methods return same result", {
 })
 
 test_that("loo.function runs with multiple cores", {
-  loo_with_fn1 <- loo(llfun, data = data, draws = draws,
-                     r_eff = rep(1, nrow(data)), cores = 1)
-  loo_with_fn2 <- loo(llfun, data = data, draws = draws,
-                      r_eff = rep(1, nrow(data)), cores = 2)
+  loo_with_fn1 <- loo(
+    llfun,
+    data = data,
+    draws = draws,
+    r_eff = rep(1, nrow(data)),
+    cores = 1
+  )
+  loo_with_fn2 <- loo(
+    llfun,
+    data = data,
+    draws = draws,
+    r_eff = rep(1, nrow(data)),
+    cores = 2
+  )
   expect_identical(loo_with_fn2$estimates, loo_with_fn1$estimates)
 })
 
 test_that("save_psis option to loo.function makes correct psis object", {
-  loo_with_fn2 <- loo.function(llfun, data = data, draws = draws,
-                      r_eff = rep(1, nrow(data)), save_psis = TRUE)
+  loo_with_fn2 <- loo.function(
+    llfun,
+    data = data,
+    draws = draws,
+    r_eff = rep(1, nrow(data)),
+    save_psis = TRUE
+  )
   expect_identical(loo_with_fn2$psis_object, loo_with_mat$psis_object)
 })
 
@@ -196,4 +239,3 @@ test_that("loo doesn't throw r_eff warnings", {
   expect_no_warning(loo(-LLmat))
   expect_no_warning(loo(llfun, data = data, draws = draws))
 })
-
diff --git a/tests/testthat/test_loo_approximate_posterior.R b/tests/testthat/test_loo_approximate_posterior.R
index f241ba2f..06fd503c 100644
--- a/tests/testthat/test_loo_approximate_posterior.R
+++ b/tests/testthat/test_loo_approximate_posterior.R
@@ -1,5 +1,3 @@
-suppressPackageStartupMessages(library(loo))
-
 # Create test data
 # Checked by Mans M and Paul B 24th of June 2019
 set.seed(123)
@@ -10,58 +8,95 @@ a0 <- 1
 b0 <- 1
 p <- 0.5
 y <- rbinom(N, size = K, prob = p)
-fake_data <- data.frame(y,K)
+fake_data <- data.frame(y, K)
 
 # The log posterior
 log_post <- function(p, y, a0, b0, K) {
-  log_lik <- sum(dbinom(x = y, size = K, prob = p, log = TRUE))  # the log likelihood
-  log_post <- log_lik + dbeta(x = p, shape1 = a0, shape2 = b0, log = TRUE)  # the log prior
+  log_lik <- sum(dbinom(x = y, size = K, prob = p, log = TRUE)) # the log likelihood
+  log_post <- log_lik + dbeta(x = p, shape1 = a0, shape2 = b0, log = TRUE) # the log prior
   log_post
 }
-it <- optim(par = 0.5, fn = log_post, control = list(fnscale = -1),
-            hessian = TRUE,  y = y, a0 = a0, b0 = b0, K = K,
-            lower = 0.01, upper = 0.99, method = "Brent")
+it <- optim(
+  par = 0.5,
+  fn = log_post,
+  control = list(fnscale = -1),
+  hessian = TRUE,
+  y = y,
+  a0 = a0,
+  b0 = b0,
+  K = K,
+  lower = 0.01,
+  upper = 0.99,
+  method = "Brent"
+)
 lap_params <- c(mu = it$par, sd = sqrt(solve(-it$hessian)))
 
 a <- a0 + sum(y)
 b <- b0 + N * K - sum(y)
 fake_true_posterior <- as.matrix(rbeta(S, a, b))
-fake_laplace_posterior <- as.matrix(rnorm(n = S, mean = lap_params["mu"], sd = lap_params["sd"]))
+fake_laplace_posterior <- as.matrix(rnorm(
+  n = S,
+  mean = lap_params["mu"],
+  sd = lap_params["sd"]
+))
 # mean(fake_laplace_posterior); sd(fake_laplace_posterior)
 
 p_draws <- as.vector(fake_laplace_posterior)
 log_p <- numeric(S)
-for(s in 1:S){
+for (s in 1:S) {
   log_p[s] <- log_post(p_draws[s], y = y, a0 = a0, b0 = b0, K = K)
 }
-log_g <- as.vector(dnorm(as.vector(fake_laplace_posterior), mean = lap_params["mu"], sd = lap_params["sd"], log = TRUE))
+log_g <- as.vector(dnorm(
+  as.vector(fake_laplace_posterior),
+  mean = lap_params["mu"],
+  sd = lap_params["sd"],
+  log = TRUE
+))
 
 llfun <- function(data_i, draws) {
   dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE)
 }
 
 ll <- matrix(0, nrow = S, ncol = N)
-for(j in 1:N){
-  ll[, j] <- llfun(data_i = fake_data[j, , drop=FALSE], draws = fake_laplace_posterior)
+for (j in 1:N) {
+  ll[, j] <- llfun(
+    data_i = fake_data[j, , drop = FALSE],
+    draws = fake_laplace_posterior
+  )
 }
 
 
 test_that("loo_approximate_posterior.array works as loo_approximate_posterior.matrix", {
-
   # Create array with two "chains"
-  log_p_mat <- matrix(log_p, nrow = (S/2), ncol = 2)
-  log_g_mat <- matrix(log_g, nrow = (S/2), ncol = 2)
-  ll_array <- array(0, dim = c((S/2), 2 , ncol(ll)))
-  ll_array[,1,] <- ll[1:(S/2),]
-  ll_array[,2,] <- ll[(S/2 + 1):S,]
+  log_p_mat <- matrix(log_p, nrow = (S / 2), ncol = 2)
+  log_g_mat <- matrix(log_g, nrow = (S / 2), ncol = 2)
+  ll_array <- array(0, dim = c((S / 2), 2, ncol(ll)))
+  ll_array[, 1, ] <- ll[1:(S / 2), ]
+  ll_array[, 2, ] <- ll[(S / 2 + 1):S, ]
 
   # Assert that they are ok
-  expect_equal(ll_array[1:2,1,1:2], ll[1:2,1:2], ignore_attr = TRUE)
-  expect_equal(ll_array[1:2,2,1:2], ll[(S/2+1):((S/2)+2),1:2], ignore_attr = TRUE)
+  expect_equal(ll_array[1:2, 1, 1:2], ll[1:2, 1:2], ignore_attr = TRUE)
+  expect_equal(
+    ll_array[1:2, 2, 1:2],
+    ll[(S / 2 + 1):((S / 2) + 2), 1:2],
+    ignore_attr = TRUE
+  )
 
   # Compute aploo
-  expect_silent(aploo1 <- loo_approximate_posterior.matrix(x = ll, log_p = log_p, log_g = log_g))
-  expect_silent(aploo2 <- loo_approximate_posterior.array(x = ll_array, log_p = log_p_mat, log_g = log_g_mat))
+  expect_silent(
+    aploo1 <- loo_approximate_posterior.matrix(
+      x = ll,
+      log_p = log_p,
+      log_g = log_g
+    )
+  )
+  expect_silent(
+    aploo2 <- loo_approximate_posterior.array(
+      x = ll_array,
+      log_p = log_p_mat,
+      log_g = log_g_mat
+    )
+  )
   expect_silent(aploo1b <- loo.matrix(x = ll, r_eff = rep(1, N)))
 
   # Check equivalence
@@ -71,24 +106,53 @@ test_that("loo_approximate_posterior.array works as loo_approximate_posterior.ma
   expect_failure(expect_equal(class(aploo1), class(aploo1b)))
 
   # Should fail with matrix
-  expect_error(aploo2 <- loo_approximate_posterior.matrix(x = ll, log_p = as.matrix(log_p), log_g = log_g))
-  expect_error(aploo2 <- loo_approximate_posterior.matrix(x = ll, log_p = as.matrix(log_p), log_g = as.matrix(log_g)))
+  expect_error(
+    aploo2 <- loo_approximate_posterior.matrix(
+      x = ll,
+      log_p = as.matrix(log_p),
+      log_g = log_g
+    )
+  )
+  expect_error(
+    aploo2 <- loo_approximate_posterior.matrix(
+      x = ll,
+      log_p = as.matrix(log_p),
+      log_g = as.matrix(log_g)
+    )
+  )
 
   # Expect log_p and log_g be stored in the approximate_posterior in the same way
   expect_length(aploo1$approximate_posterior$log_p, nrow(ll))
   expect_length(aploo1$approximate_posterior$log_g, nrow(ll))
-  expect_equal(aploo1$approximate_posterior$log_p, aploo2$approximate_posterior$log_p)
-  expect_equal(aploo1$approximate_posterior$log_g, aploo2$approximate_posterior$log_g)
-
+  expect_equal(
+    aploo1$approximate_posterior$log_p,
+    aploo2$approximate_posterior$log_p
+  )
+  expect_equal(
+    aploo1$approximate_posterior$log_g,
+    aploo2$approximate_posterior$log_g
+  )
 })
 
 test_that("loo_approximate_posterior.function works as loo_approximate_posterior.matrix", {
-
-
   # Compute aploo
-  expect_silent(aploo1 <- loo_approximate_posterior.matrix(x = ll, log_p = log_p, log_g = log_g))
+  expect_silent(
+    aploo1 <- loo_approximate_posterior.matrix(
+      x = ll,
+      log_p = log_p,
+      log_g = log_g
+    )
+  )
   expect_silent(aploo1b <- loo.matrix(x = ll, r_eff = rep(1, N)))
-  expect_silent(aploo2 <- loo_approximate_posterior.function(x = llfun, log_p = log_p, log_g = log_g, data = fake_data, draws = fake_laplace_posterior))
+  expect_silent(
+    aploo2 <- loo_approximate_posterior.function(
+      x = llfun,
+      log_p = log_p,
+      log_g = log_g,
+      data = fake_data,
+      draws = fake_laplace_posterior
+    )
+  )
 
   # Check equivalence
   expect_equal(aploo1$estimates, aploo2$estimates)
@@ -97,11 +161,20 @@ test_that("loo_approximate_posterior.function works as loo_approximate_posterior
 
   # Check equivalence
   # Expect log_p and log_g be stored in the approximate_posterior in the same way
-  expect_length(aploo2$approximate_posterior$log_p, nrow(fake_laplace_posterior))
-  expect_length(aploo2$approximate_posterior$log_g, nrow(fake_laplace_posterior))
-  expect_equal(aploo1$approximate_posterior$log_p, aploo2$approximate_posterior$log_p)
-  expect_equal(aploo1$approximate_posterior$log_g, aploo2$approximate_posterior$log_g)
-
+  expect_length(
+    aploo2$approximate_posterior$log_p,
+    nrow(fake_laplace_posterior)
+  )
+  expect_length(
+    aploo2$approximate_posterior$log_g,
+    nrow(fake_laplace_posterior)
+  )
+  expect_equal(
+    aploo1$approximate_posterior$log_p,
+    aploo2$approximate_posterior$log_p
+  )
+  expect_equal(
+    aploo1$approximate_posterior$log_g,
+    aploo2$approximate_posterior$log_g
+  )
 })
-
-
diff --git a/tests/testthat/test_loo_moment_matching.R b/tests/testthat/test_loo_moment_matching.R
index 039c3b35..f96ce70b 100644
--- a/tests/testthat/test_loo_moment_matching.R
+++ b/tests/testthat/test_loo_moment_matching.R
@@ -1,35 +1,51 @@
-library(loo)
 options(mc.cores = 1)
 
 set.seed(123)
 S <- 4000
 
 # helper functions for sampling from the posterior distribution
-rinvchisq <- function(n, df, scale = 1/df, ...)
-{
-  if ((length(scale) != 1) & (length(scale) != n))
+rinvchisq <- function(n, df, scale = 1 / df, ...) {
+  if ((length(scale) != 1) & (length(scale) != n)) {
     stop("scale should be a scalar or a vector of the same length as x")
-  if (df <= 0)
+  }
+  if (df <= 0) {
     stop("df must be greater than zero")
-  if (any(scale <= 0))
+  }
+  if (any(scale <= 0)) {
     stop("scale must be greater than zero")
-  return((df*scale)/rchisq(n, df = df))
+  }
+  return((df * scale) / rchisq(n, df = df))
 }
 
-dinvchisq <- function(x, df, scale=1/df, log = FALSE, ...)
-{
-  if (df <= 0)
+dinvchisq <- function(x, df, scale = 1 / df, log = FALSE, ...) {
+  if (df <= 0) {
     stop("df must be greater than zero")
-  if (scale <= 0)
+  }
+  if (scale <= 0) {
     stop("scale must be greater than zero")
-  nu <- df/2
-  if (log)
-    return(ifelse(x > 0, nu*log(nu) - log(gamma(nu)) + nu*log(scale) -
-                    (nu + 1)*log(x) - (nu*scale/x), NA))
-  else
-    return(ifelse(x > 0,
-                  (((nu)^(nu))/gamma(nu)) * (scale^nu) *
-                    (x^(-(nu + 1))) * exp(-nu*scale/x), NA))
+  }
+  nu <- df / 2
+  if (log) {
+    return(ifelse(
+      x > 0,
+      nu *
+        log(nu) -
+        log(gamma(nu)) +
+        nu * log(scale) -
+        (nu + 1) * log(x) -
+        (nu * scale / x),
+      NA
+    ))
+  } else {
+    return(ifelse(
+      x > 0,
+      (((nu)^(nu)) / gamma(nu)) *
+        (scale^nu) *
+        (x^(-(nu + 1))) *
+        exp(-nu * scale / x),
+      NA
+    ))
+  }
 }
 
 
@@ -43,11 +59,15 @@ y_tilde <- 11
 y[1] <- y_tilde
 
 ymean <- mean(y)
-s2 <- sum((y - ymean)^2)/(n - 1)
+s2 <- sum((y - ymean)^2) / (n - 1)
 
 # draws from the posterior distribution when including all observations
 draws_full_posterior_sigma2 <- rinvchisq(S, n - 1, s2)
-draws_full_posterior_mu <- rnorm(S, ymean, sqrt(draws_full_posterior_sigma2/n))
+draws_full_posterior_mu <- rnorm(
+  S,
+  ymean,
+  sqrt(draws_full_posterior_sigma2 / n)
+)
 
 
 # create a dummy model object
@@ -64,15 +84,8 @@ x$draws <- data.frame(
 )
 
 
-
-
-
-
-
-
 # implement functions for moment matching loo
 
-
 # extract original posterior draws
 post_draws_test <- function(x, ...) {
   as.matrix(x$draws)
@@ -80,42 +93,50 @@ post_draws_test <- function(x, ...) {
 
 # extract original log lik draws
 log_lik_i_test <- function(x, i, ...) {
-  -0.5*log(2*pi) - log(x$draws$sigma) - 1.0/(2*x$draws$sigma^2)*(x$data$y[i] - x$draws$mu)^2
+  -0.5 *
+    log(2 * pi) -
+    log(x$draws$sigma) -
+    1.0 / (2 * x$draws$sigma^2) * (x$data$y[i] - x$draws$mu)^2
 }
 
 
-loglik <- matrix(0,S,n)
+loglik <- matrix(0, S, n)
 for (j in seq(n)) {
-  loglik[,j] <- log_lik_i_test(x, j)
+  loglik[, j] <- log_lik_i_test(x, j)
 }
 
 
-
 # mu, log(sigma)
 unconstrain_pars_test <- function(x, pars, ...) {
   upars <- as.matrix(pars)
-  upars[,2] <- log(upars[,2])
+  upars[, 2] <- log(upars[, 2])
   upars
 }
 
 log_prob_upars_test <- function(x, upars, ...) {
-  dinvchisq(exp(upars[,2])^2,x$data$n - 1,x$data$s2, log = TRUE) +
-    dnorm(upars[,1],x$data$ymean,exp(upars[,2])/sqrt(x$data$n), log = TRUE)
+  dinvchisq(exp(upars[, 2])^2, x$data$n - 1, x$data$s2, log = TRUE) +
+    dnorm(
+      upars[, 1],
+      x$data$ymean,
+      exp(upars[, 2]) / sqrt(x$data$n),
+      log = TRUE
+    )
 }
 
 # compute log_lik_i values based on the unconstrained parameters
 log_lik_i_upars_test <- function(x, upars, i, ...) {
-  -0.5*log(2*pi) - upars[,2] - 1.0/(2*exp(upars[,2])^2)*(x$data$y[i] - upars[,1])^2
+  -0.5 *
+    log(2 * pi) -
+    upars[, 2] -
+    1.0 / (2 * exp(upars[, 2])^2) * (x$data$y[i] - upars[, 1])^2
 }
 
 
-
 upars <- unconstrain_pars_test(x, x$draws)
-lwi_1 <- -loglik[,1]
+lwi_1 <- -loglik[, 1]
 lwi_1 <- lwi_1 - matrixStats::logSumExp(lwi_1)
 
 
-
 test_that("log_prob_upars_test works", {
   upars <- unconstrain_pars_test(x, x$draws)
   xloo <- list()
@@ -123,13 +144,13 @@ test_that("log_prob_upars_test works", {
   xloo$data$y <- y[-1]
   xloo$data$n <- n - 1
   xloo$data$ymean <- mean(y[-1])
-  xloo$data$s2 <- sum((y[-1] - mean(y[-1]))^2)/(n - 2)
+  xloo$data$s2 <- sum((y[-1] - mean(y[-1]))^2) / (n - 2)
 
-  post1 <- log_prob_upars_test(x,upars)
+  post1 <- log_prob_upars_test(x, upars)
   post1 <- post1 - matrixStats::logSumExp(post1)
-  post2 <- log_prob_upars_test(xloo,upars) + loglik[,1]
+  post2 <- log_prob_upars_test(xloo, upars) + loglik[, 1]
   post2 <- post2 - matrixStats::logSumExp(post2)
-  expect_equal(post1,post2)
+  expect_equal(post1, post2)
 })
 
 
@@ -138,98 +159,207 @@ test_that("loo_moment_match.default warnings work", {
   loo_manual <- suppressWarnings(loo(loglik))
   loo_manual_tis <- suppressWarnings(loo(loglik, is_method = "tis"))
 
+  expect_warning(
+    loo_moment_match(
+      x,
+      loo_manual,
+      post_draws_test,
+      log_lik_i_test,
+      unconstrain_pars_test,
+      log_prob_upars_test,
+      log_lik_i_upars_test,
+      max_iters = 30L,
+      k_thres = 0.5,
+      split = FALSE,
+      cov = TRUE,
+      cores = 1
+    ),
+    "The accuracy of self-normalized importance sampling"
+  )
 
-  expect_warning(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                              unconstrain_pars_test, log_prob_upars_test,
-                              log_lik_i_upars_test, max_iters = 30L,
-                              k_thres = 0.5, split = FALSE,
-                              cov = TRUE, cores = 1), "The accuracy of self-normalized importance sampling")
-
-  expect_warning(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                              unconstrain_pars_test, log_prob_upars_test,
-                              log_lik_i_upars_test, max_iters = 30L,
-                              split = FALSE,
-                              cov = TRUE, cores = 1), "The accuracy of self-normalized importance sampling")
-
-  expect_no_warning(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                              unconstrain_pars_test, log_prob_upars_test,
-                              log_lik_i_upars_test, max_iters = 30L,
-                              k_thres = 100, split = TRUE,
-                              cov = TRUE, cores = 1))
-
-  expect_snapshot(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                              unconstrain_pars_test, log_prob_upars_test,
-                              log_lik_i_upars_test, max_iters = 1,
-                              k_thres = 0.5, split = TRUE,
-                              cov = TRUE, cores = 1))
-
-  expect_error(loo_moment_match(x, loo_manual_tis, post_draws_test, log_lik_i_test,
-                       unconstrain_pars_test, log_prob_upars_test,
-                       log_lik_i_upars_test, max_iters = 30L,
-                       k_thres = 0.5, split = TRUE,
-                       cov = TRUE, cores = 1), "loo_moment_match currently supports only")
-})
-
+  expect_warning(
+    loo_moment_match(
+      x,
+      loo_manual,
+      post_draws_test,
+      log_lik_i_test,
+      unconstrain_pars_test,
+      log_prob_upars_test,
+      log_lik_i_upars_test,
+      max_iters = 30L,
+      split = FALSE,
+      cov = TRUE,
+      cores = 1
+    ),
+    "The accuracy of self-normalized importance sampling"
+  )
 
+  expect_no_warning(loo_moment_match(
+    x,
+    loo_manual,
+    post_draws_test,
+    log_lik_i_test,
+    unconstrain_pars_test,
+    log_prob_upars_test,
+    log_lik_i_upars_test,
+    max_iters = 30L,
+    k_thres = 100,
+    split = TRUE,
+    cov = TRUE,
+    cores = 1
+  ))
+
+  expect_snapshot(loo_moment_match(
+    x,
+    loo_manual,
+    post_draws_test,
+    log_lik_i_test,
+    unconstrain_pars_test,
+    log_prob_upars_test,
+    log_lik_i_upars_test,
+    max_iters = 1,
+    k_thres = 0.5,
+    split = TRUE,
+    cov = TRUE,
+    cores = 1
+  ))
+
+  expect_error(
+    loo_moment_match(
+      x,
+      loo_manual_tis,
+      post_draws_test,
+      log_lik_i_test,
+      unconstrain_pars_test,
+      log_prob_upars_test,
+      log_lik_i_upars_test,
+      max_iters = 30L,
+      k_thres = 0.5,
+      split = TRUE,
+      cov = TRUE,
+      cores = 1
+    ),
+    "loo_moment_match currently supports only"
+  )
+})
 
 
 test_that("loo_moment_match.default works", {
-
   # allow -Inf
   lwi_x <- lwi_1
   lwi_x[which.min(lwi_1)] <- -Inf
-  expect_no_error(suppressWarnings(importance_sampling.default(lwi_1, method = "psis", r_eff = 1, cores = 1)))
+  expect_no_error(suppressWarnings(importance_sampling.default(
+    lwi_1,
+    method = "psis",
+    r_eff = 1,
+    cores = 1
+  )))
 
   # loo object
   loo_manual <- suppressWarnings(loo(loglik))
 
-  loo_moment_match_object <- suppressWarnings(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                                                unconstrain_pars_test, log_prob_upars_test,
-                                                log_lik_i_upars_test, max_iters = 30L,
-                                                k_thres = 0.8, split = FALSE,
-                                                cov = TRUE, cores = 1))
+  loo_moment_match_object <- suppressWarnings(loo_moment_match(
+    x,
+    loo_manual,
+    post_draws_test,
+    log_lik_i_test,
+    unconstrain_pars_test,
+    log_prob_upars_test,
+    log_lik_i_upars_test,
+    max_iters = 30L,
+    k_thres = 0.8,
+    split = FALSE,
+    cov = TRUE,
+    cores = 1
+  ))
 
   # diagnostic Pareto k decreases but influence pareto k stays the same
-  expect_lt(loo_moment_match_object$diagnostics$pareto_k[1], loo_moment_match_object$pointwise[1,"influence_pareto_k"])
-  expect_equal(loo_moment_match_object$pointwise[,"influence_pareto_k"],loo_manual$pointwise[,"influence_pareto_k"])
-  expect_equal(loo_moment_match_object$pointwise[,"influence_pareto_k"],loo_manual$diagnostics$pareto_k)
+  expect_lt(
+    loo_moment_match_object$diagnostics$pareto_k[1],
+    loo_moment_match_object$pointwise[1, "influence_pareto_k"]
+  )
+  expect_equal(
+    loo_moment_match_object$pointwise[, "influence_pareto_k"],
+    loo_manual$pointwise[, "influence_pareto_k"]
+  )
+  expect_equal(
+    loo_moment_match_object$pointwise[, "influence_pareto_k"],
+    loo_manual$diagnostics$pareto_k
+  )
 
   expect_snapshot_value(loo_moment_match_object, style = "serialize")
-  
-  loo_moment_match_object2 <- suppressWarnings(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                                                unconstrain_pars_test, log_prob_upars_test,
-                                                log_lik_i_upars_test, max_iters = 30L,
-                                                k_thres = 0.5, split = FALSE,
-                                                cov = TRUE, cores = 1))
+
+  loo_moment_match_object2 <- suppressWarnings(loo_moment_match(
+    x,
+    loo_manual,
+    post_draws_test,
+    log_lik_i_test,
+    unconstrain_pars_test,
+    log_prob_upars_test,
+    log_lik_i_upars_test,
+    max_iters = 30L,
+    k_thres = 0.5,
+    split = FALSE,
+    cov = TRUE,
+    cores = 1
+  ))
 
   expect_snapshot_value(loo_moment_match_object2, style = "serialize")
 
-  loo_moment_match_object3 <- suppressWarnings(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                                                unconstrain_pars_test, log_prob_upars_test,
-                                                log_lik_i_upars_test, max_iters = 30L,
-                                                k_thres = 0.5, split = TRUE,
-                                                cov = TRUE, cores = 1))
+  loo_moment_match_object3 <- suppressWarnings(loo_moment_match(
+    x,
+    loo_manual,
+    post_draws_test,
+    log_lik_i_test,
+    unconstrain_pars_test,
+    log_prob_upars_test,
+    log_lik_i_upars_test,
+    max_iters = 30L,
+    k_thres = 0.5,
+    split = TRUE,
+    cov = TRUE,
+    cores = 1
+  ))
 
   expect_snapshot_value(loo_moment_match_object3, style = "serialize")
 
-  loo_moment_match_object4 <- suppressWarnings(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                                                unconstrain_pars_test, log_prob_upars_test,
-                                                log_lik_i_upars_test, max_iters = 30L,
-                                                k_thres = 100, split = FALSE,
-                                                cov = TRUE, cores = 1))
-
-  expect_equal(loo_manual,loo_moment_match_object4)
+  loo_moment_match_object4 <- suppressWarnings(loo_moment_match(
+    x,
+    loo_manual,
+    post_draws_test,
+    log_lik_i_test,
+    unconstrain_pars_test,
+    log_prob_upars_test,
+    log_lik_i_upars_test,
+    max_iters = 30L,
+    k_thres = 100,
+    split = FALSE,
+    cov = TRUE,
+    cores = 1
+  ))
+
+  expect_equal(loo_manual, loo_moment_match_object4)
 
   loo_manual_with_psis <- suppressWarnings(loo(loglik, save_psis = TRUE))
-  loo_moment_match_object5 <- suppressWarnings(loo_moment_match(x, loo_manual_with_psis, post_draws_test, log_lik_i_test,
-                                          unconstrain_pars_test, log_prob_upars_test,
-                                          log_lik_i_upars_test, max_iters = 30L,
-                                          k_thres = 0.8, split = FALSE,
-                                          cov = TRUE, cores = 1))
-
-  expect_equal(loo_moment_match_object5$diagnostics,loo_moment_match_object5$psis_object$diagnostics)
-
-
+  loo_moment_match_object5 <- suppressWarnings(loo_moment_match(
+    x,
+    loo_manual_with_psis,
+    post_draws_test,
+    log_lik_i_test,
+    unconstrain_pars_test,
+    log_prob_upars_test,
+    log_lik_i_upars_test,
+    max_iters = 30L,
+    k_thres = 0.8,
+    split = FALSE,
+    cov = TRUE,
+    cores = 1
+  ))
+
+  expect_equal(
+    loo_moment_match_object5$diagnostics,
+    loo_moment_match_object5$psis_object$diagnostics
+  )
 })
 
 test_that("variance and covariance transformations work", {
@@ -237,57 +367,99 @@ test_that("variance and covariance transformations work", {
 
   set.seed(8493874)
   draws_full_posterior_sigma2 <- rinvchisq(S, n - 1, s2)
-  draws_full_posterior_mu <- rnorm(S, ymean, sqrt(draws_full_posterior_sigma2/n))
+  draws_full_posterior_mu <- rnorm(
+    S,
+    ymean,
+    sqrt(draws_full_posterior_sigma2 / n)
+  )
 
   x$draws <- data.frame(
     mu = draws_full_posterior_mu,
     sigma = sqrt(draws_full_posterior_sigma2)
   )
-  loglik <- matrix(0,S,n)
+  loglik <- matrix(0, S, n)
   for (j in seq(n)) {
-    loglik[,j] <- log_lik_i_test(x, j)
+    loglik[, j] <- log_lik_i_test(x, j)
   }
 
   upars <- unconstrain_pars_test(x, x$draws)
-  lwi_1 <- -loglik[,1]
+  lwi_1 <- -loglik[, 1]
   lwi_1 <- lwi_1 - matrixStats::logSumExp(lwi_1)
 
   loo_manual <- suppressWarnings(loo(loglik))
 
-  loo_moment_match_object <- suppressWarnings(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                                         unconstrain_pars_test, log_prob_upars_test,
-                                         log_lik_i_upars_test, max_iters = 30L,
-                                         k_thres = 0.0, split = FALSE,
-                                         cov = TRUE, cores = 1))
+  loo_moment_match_object <- suppressWarnings(loo_moment_match(
+    x,
+    loo_manual,
+    post_draws_test,
+    log_lik_i_test,
+    unconstrain_pars_test,
+    log_prob_upars_test,
+    log_lik_i_upars_test,
+    max_iters = 30L,
+    k_thres = 0.0,
+    split = FALSE,
+    cov = TRUE,
+    cores = 1
+  ))
 
   expect_snapshot_value(loo_moment_match_object, style = "serialize")
-
 })
 
 test_that("loo_moment_match.default works with multiple cores", {
-
   # loo object
   loo_manual <- suppressWarnings(loo(loglik))
 
-  loo_moment_match_manual3 <- suppressWarnings(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                                                 unconstrain_pars_test, log_prob_upars_test,
-                                                 log_lik_i_upars_test, max_iters = 30L,
-                                                 k_thres = 0.5, split = FALSE,
-                                                 cov = TRUE, cores = 1))
-
-  loo_moment_match_manual4 <- suppressWarnings(loo_moment_match(x, loo_manual, post_draws_test, log_lik_i_test,
-                                                 unconstrain_pars_test, log_prob_upars_test,
-                                                 log_lik_i_upars_test, max_iters = 30L,
-                                                 k_thres = 0.5, split = FALSE,
-                                                 cov = TRUE, cores = 2))
-
-  expect_equal(loo_moment_match_manual3$diagnostics$pareto_k, loo_moment_match_manual4$diagnostics$pareto_k)
-  expect_equal(loo_moment_match_manual3$diagnostics$n_eff, loo_moment_match_manual4$diagnostics$n_eff, tolerance = 5e-4)
-
-  expect_equal(loo_moment_match_manual3$estimates, loo_moment_match_manual4$estimates)
+  loo_moment_match_manual3 <- suppressWarnings(loo_moment_match(
+    x,
+    loo_manual,
+    post_draws_test,
+    log_lik_i_test,
+    unconstrain_pars_test,
+    log_prob_upars_test,
+    log_lik_i_upars_test,
+    max_iters = 30L,
+    k_thres = 0.5,
+    split = FALSE,
+    cov = TRUE,
+    cores = 1
+  ))
+
+  loo_moment_match_manual4 <- suppressWarnings(loo_moment_match(
+    x,
+    loo_manual,
+    post_draws_test,
+    log_lik_i_test,
+    unconstrain_pars_test,
+    log_prob_upars_test,
+    log_lik_i_upars_test,
+    max_iters = 30L,
+    k_thres = 0.5,
+    split = FALSE,
+    cov = TRUE,
+    cores = 2
+  ))
+
+  expect_equal(
+    loo_moment_match_manual3$diagnostics$pareto_k,
+    loo_moment_match_manual4$diagnostics$pareto_k
+  )
+  expect_equal(
+    loo_moment_match_manual3$diagnostics$n_eff,
+    loo_moment_match_manual4$diagnostics$n_eff,
+    tolerance = 5e-4
+  )
 
-  expect_equal(loo_moment_match_manual3$pointwise, loo_moment_match_manual4$pointwise, tolerance=5e-4)
+  expect_equal(
+    loo_moment_match_manual3$estimates,
+    loo_moment_match_manual4$estimates
+  )
 
+  expect_equal(
+    loo_moment_match_manual3$pointwise,
+    loo_moment_match_manual4$pointwise,
+    tolerance = 5e-4
+  )
 })
 
 
@@ -295,50 +467,97 @@ test_that("loo_moment_match_split works", {
   # skip on M1 Mac until we figure out why this test fails only on M1 Mac
   skip_if(Sys.info()[["sysname"]] == "Darwin" && R.version$arch == "aarch64")
 
-  is_obj_1 <- suppressWarnings(importance_sampling.default(lwi_1, method = "psis", r_eff = 1, cores = 1))
+  is_obj_1 <- suppressWarnings(importance_sampling.default(
+    lwi_1,
+    method = "psis",
+    r_eff = 1,
+    cores = 1
+  ))
   lwi_1_ps <- as.vector(weights(is_obj_1))
 
   split <- loo_moment_match_split(
-    x, upars, cov = FALSE, total_shift = c(0,0), total_scaling = c(1,1), total_mapping = diag(c(1,1)), i = 1,
-    log_prob_upars = log_prob_upars_test, log_lik_i_upars = log_lik_i_upars_test,
-    cores = 1, r_eff_i = 1, is_method = "psis")
+    x,
+    upars,
+    cov = FALSE,
+    total_shift = c(0, 0),
+    total_scaling = c(1, 1),
+    total_mapping = diag(c(1, 1)),
+    i = 1,
+    log_prob_upars = log_prob_upars_test,
+    log_lik_i_upars = log_lik_i_upars_test,
+    cores = 1,
+    r_eff_i = 1,
+    is_method = "psis"
+  )
 
-  expect_named(split,c("lwi", "lwfi", "log_liki", "r_eff_i"))
+  expect_named(split, c("lwi", "lwfi", "log_liki", "r_eff_i"))
 
-  expect_equal(lwi_1_ps,split$lwi)
+  expect_equal(lwi_1_ps, split$lwi)
 
   split2 <- loo_moment_match_split(
-    x, upars, cov = FALSE, total_shift = c(-0.1,-0.2), total_scaling = c(0.7,0.7),
-    total_mapping = matrix(c(1,0.1,0.1,1),2,2), i = 1,
-    log_prob_upars = log_prob_upars_test, log_lik_i_upars = log_lik_i_upars_test,
-    cores = 1, r_eff_i = 1, is_method = "psis")
+    x,
+    upars,
+    cov = FALSE,
+    total_shift = c(-0.1, -0.2),
+    total_scaling = c(0.7, 0.7),
+    total_mapping = matrix(c(1, 0.1, 0.1, 1), 2, 2),
+    i = 1,
+    log_prob_upars = log_prob_upars_test,
+    log_lik_i_upars = log_lik_i_upars_test,
+    cores = 1,
+    r_eff_i = 1,
+    is_method = "psis"
+  )
 
   expect_snapshot_value(split2, style = "serialize")
 })
 
 test_that("passing arguments works", {
-  log_lik_i_upars_test_additional_argument <- function(x, upars, i, passed_arg = FALSE, ...) {
+  log_lik_i_upars_test_additional_argument <- function(
+    x,
+    upars,
+    i,
+    passed_arg = FALSE,
+    ...
+  ) {
     if (!passed_arg) {
       warning("passed_arg was not passed here")
     }
-    -0.5*log(2*pi) - upars[,2] - 1.0/(2*exp(upars[,2])^2)*(x$data$y[i] - upars[,1])^2
-
+    -0.5 *
+      log(2 * pi) -
+      upars[, 2] -
+      1.0 / (2 * exp(upars[, 2])^2) * (x$data$y[i] - upars[, 1])^2
   }
-  unconstrain_pars_test_additional_argument <- function(x, pars, passed_arg = FALSE, ...) {
+  unconstrain_pars_test_additional_argument <- function(
+    x,
+    pars,
+    passed_arg = FALSE,
+    ...
+  ) {
     if (!passed_arg) {
       warning("passed_arg was not passed here")
     }
     upars <- as.matrix(pars)
-    upars[,2] <- log(upars[,2])
+    upars[, 2] <- log(upars[, 2])
     upars
   }
 
-  log_prob_upars_test_additional_argument <- function(x, upars, passed_arg = FALSE, ...) {
+  log_prob_upars_test_additional_argument <- function(
+    x,
+    upars,
+    passed_arg = FALSE,
+    ...
+  ) {
     if (!passed_arg) {
       warning("passed_arg was not passed here")
     }
-    dinvchisq(exp(upars[,2])^2,x$data$n - 1,x$data$s2, log = TRUE) +
-      dnorm(upars[,1],x$data$ymean,exp(upars[,2])/sqrt(x$data$n), log = TRUE)
+    dinvchisq(exp(upars[, 2])^2, x$data$n - 1, x$data$s2, log = TRUE) +
+      dnorm(
+        upars[, 1],
+        x$data$ymean,
+        exp(upars[, 2]) / sqrt(x$data$n),
+        log = TRUE
+      )
   }
   post_draws_test_additional_argument <- function(x, passed_arg = FALSE, ...) {
     if (!passed_arg) {
@@ -346,18 +565,36 @@ test_that("passing arguments works", {
     }
     as.matrix(x$draws)
   }
-  log_lik_i_test_additional_argument <- function(x, i, passed_arg = FALSE, ...) {
+  log_lik_i_test_additional_argument <- function(
+    x,
+    i,
+    passed_arg = FALSE,
+    ...
+  ) {
     if (!passed_arg) {
       warning("passed_arg was not passed here")
     }
-    -0.5*log(2*pi) - log(x$draws$sigma) - 1.0/(2*x$draws$sigma^2)*(x$data$y[i] - x$draws$mu)^2
+    -0.5 *
+      log(2 * pi) -
+      log(x$draws$sigma) -
+      1.0 / (2 * x$draws$sigma^2) * (x$data$y[i] - x$draws$mu)^2
   }
 
   # loo object
   loo_manual <- suppressWarnings(loo(loglik))
-  expect_silent(loo_moment_match(x, loo_manual, post_draws_test_additional_argument, log_lik_i_test_additional_argument,
-                                 unconstrain_pars_test_additional_argument, log_prob_upars_test_additional_argument,
-                                 log_lik_i_upars_test_additional_argument, max_iters = 30L,
-                                 k_thres = 0.5, split = TRUE,
-                                 cov = TRUE, cores = 1, passed_arg = TRUE))
+  expect_silent(loo_moment_match(
+    x,
+    loo_manual,
+    post_draws_test_additional_argument,
+    log_lik_i_test_additional_argument,
+    unconstrain_pars_test_additional_argument,
+    log_prob_upars_test_additional_argument,
+    log_lik_i_upars_test_additional_argument,
+    max_iters = 30L,
+    k_thres = 0.5,
+    split = TRUE,
+    cov = TRUE,
+    cores = 1,
+    passed_arg = TRUE
+  ))
 })
diff --git a/tests/testthat/test_loo_subsampling.R b/tests/testthat/test_loo_subsampling.R
index 3bd38dee..84e8d785 100644
--- a/tests/testthat/test_loo_subsampling.R
+++ b/tests/testthat/test_loo_subsampling.R
@@ -1,14 +1,18 @@
-library(loo)
 options(mc.cores = 1)
 
 test_that("overall loo_subampling works as expected (compared with loo) for diff_est", {
   set.seed(123)
-  N <- 1000; K <- 10; S <- 1000; a0 <- 3; b0 <- 2
+  N <- 1000
+  K <- 10
+  S <- 1000
+  a0 <- 3
+  b0 <- 2
   p <- 0.7
   y <- rbinom(N, size = K, prob = p)
-  a <- a0 + sum(y); b <- b0 + N * K - sum(y)
+  a <- a0 + sum(y)
+  b <- b0 + N * K - sum(y)
   fake_posterior <- as.matrix(rbeta(S, a, b))
-  fake_data <- data.frame(y,K)
+  fake_data <- data.frame(y, K)
   rm(N, K, S, a0, b0, p, y, a, b)
   llfun_test <- function(data_i, draws) {
     # each time called internally within loo the arguments will be equal to:
@@ -17,63 +21,227 @@ test_that("overall loo_subampling works as expected (compared with loo) for diff
     dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE)
   }
 
-  expect_silent(true_loo <- loo(llfun_test, draws = fake_posterior, data = fake_data, r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    true_loo <- loo(
+      llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(true_loo, "psis_loo")
-  expect_silent(loo_ss <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 500, loo_approximation = "plpd", r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 500,
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(loo_ss, "psis_loo_ss")
 
   # Check consistency
-  expect_equal(loo_ss$pointwise[, "elpd_loo_approx"],
-                    loo_ss$loo_subsampling$elpd_loo_approx[loo_ss$pointwise[, "idx"]], ignore_attr = TRUE)
+  expect_equal(
+    loo_ss$pointwise[, "elpd_loo_approx"],
+    loo_ss$loo_subsampling$elpd_loo_approx[loo_ss$pointwise[, "idx"]],
+    ignore_attr = TRUE
+  )
 
   # Expect values
   z <- 2
-  expect_lte(loo_ss$estimates["elpd_loo", "Estimate"] - z * loo_ss$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_gte(loo_ss$estimates["elpd_loo", "Estimate"] + z * loo_ss$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_lte(loo_ss$estimates["p_loo", "Estimate"] - z * loo_ss$estimates["p_loo", "subsampling SE"], true_loo$estimates["p_loo", "Estimate"])
-  expect_gte(loo_ss$estimates["p_loo", "Estimate"] + z * loo_ss$estimates["p_loo", "subsampling SE"], true_loo$estimates["p_loo", "Estimate"])
-  expect_lte(loo_ss$estimates["looic", "Estimate"] - z * loo_ss$estimates["looic", "subsampling SE"], true_loo$estimates["looic", "Estimate"])
-  expect_gte(loo_ss$estimates["looic", "Estimate"] + z * loo_ss$estimates["looic", "subsampling SE"], true_loo$estimates["looic", "Estimate"])
-
-  expect_failure(expect_equal(true_loo$estimates["elpd_loo", "Estimate"], loo_ss$estimates["elpd_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(true_loo$estimates["p_loo", "Estimate"], loo_ss$estimates["p_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(true_loo$estimates["looic", "Estimate"], loo_ss$estimates["looic", "Estimate"], tolerance = 0.00000001))
+  expect_lte(
+    loo_ss$estimates["elpd_loo", "Estimate"] -
+      z * loo_ss$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["elpd_loo", "Estimate"] +
+      z * loo_ss$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss$estimates["p_loo", "Estimate"] -
+      z * loo_ss$estimates["p_loo", "subsampling SE"],
+    true_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["p_loo", "Estimate"] +
+      z * loo_ss$estimates["p_loo", "subsampling SE"],
+    true_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss$estimates["looic", "Estimate"] -
+      z * loo_ss$estimates["looic", "subsampling SE"],
+    true_loo$estimates["looic", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["looic", "Estimate"] +
+      z * loo_ss$estimates["looic", "subsampling SE"],
+    true_loo$estimates["looic", "Estimate"]
+  )
+
+  expect_failure(expect_equal(
+    true_loo$estimates["elpd_loo", "Estimate"],
+    loo_ss$estimates["elpd_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    true_loo$estimates["p_loo", "Estimate"],
+    loo_ss$estimates["p_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    true_loo$estimates["looic", "Estimate"],
+    loo_ss$estimates["looic", "Estimate"],
+    tolerance = 0.00000001
+  ))
 
   # Test that observations works as expected
-  expect_message(loo_ss2 <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = obs_idx(loo_ss), loo_approximation = "plpd", r_eff = rep(1, nrow(fake_data))))
+  expect_message(
+    loo_ss2 <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = obs_idx(loo_ss),
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_equal(loo_ss2$estimates, loo_ss$estimates, tolerance = 0.00000001)
-  expect_silent(loo_ss2 <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = loo_ss, loo_approximation = "plpd", r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss2 <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = loo_ss,
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_equal(loo_ss2$estimates, loo_ss$estimates, tolerance = 0.00000001)
 
   # Test lpd
-  expect_silent(loo_ss_lpd <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 500, loo_approximation = "lpd", r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss_lpd <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 500,
+      loo_approximation = "lpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(loo_ss_lpd, "psis_loo_ss")
   z <- 2
-  expect_lte(loo_ss_lpd$estimates["elpd_loo", "Estimate"] - z * loo_ss_lpd$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_gte(loo_ss_lpd$estimates["elpd_loo", "Estimate"] + z * loo_ss_lpd$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_lte(loo_ss_lpd$estimates["p_loo", "Estimate"] - z * loo_ss_lpd$estimates["p_loo", "subsampling SE"], true_loo$estimates["p_loo", "Estimate"])
-  expect_gte(loo_ss_lpd$estimates["p_loo", "Estimate"] + z * loo_ss_lpd$estimates["p_loo", "subsampling SE"], true_loo$estimates["p_loo", "Estimate"])
-  expect_lte(loo_ss_lpd$estimates["looic", "Estimate"] - z * loo_ss_lpd$estimates["looic", "subsampling SE"], true_loo$estimates["looic", "Estimate"])
-  expect_gte(loo_ss_lpd$estimates["looic", "Estimate"] + z * loo_ss_lpd$estimates["looic", "subsampling SE"], true_loo$estimates["looic", "Estimate"])
-
-  expect_failure(expect_equal(true_loo$estimates["elpd_loo", "Estimate"], loo_ss_lpd$estimates["elpd_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(true_loo$estimates["p_loo", "Estimate"], loo_ss_lpd$estimates["p_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(true_loo$estimates["looic", "Estimate"], loo_ss_lpd$estimates["looic", "Estimate"], tolerance = 0.00000001))
-
-  expect_silent(loo_ss_lpd10 <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 500, loo_approximation = "lpd", loo_approximation_draws = 10, r_eff = rep(1, nrow(fake_data))))
+  expect_lte(
+    loo_ss_lpd$estimates["elpd_loo", "Estimate"] -
+      z * loo_ss_lpd$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss_lpd$estimates["elpd_loo", "Estimate"] +
+      z * loo_ss_lpd$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss_lpd$estimates["p_loo", "Estimate"] -
+      z * loo_ss_lpd$estimates["p_loo", "subsampling SE"],
+    true_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss_lpd$estimates["p_loo", "Estimate"] +
+      z * loo_ss_lpd$estimates["p_loo", "subsampling SE"],
+    true_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss_lpd$estimates["looic", "Estimate"] -
+      z * loo_ss_lpd$estimates["looic", "subsampling SE"],
+    true_loo$estimates["looic", "Estimate"]
+  )
+  expect_gte(
+    loo_ss_lpd$estimates["looic", "Estimate"] +
+      z * loo_ss_lpd$estimates["looic", "subsampling SE"],
+    true_loo$estimates["looic", "Estimate"]
+  )
+
+  expect_failure(expect_equal(
+    true_loo$estimates["elpd_loo", "Estimate"],
+    loo_ss_lpd$estimates["elpd_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    true_loo$estimates["p_loo", "Estimate"],
+    loo_ss_lpd$estimates["p_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    true_loo$estimates["looic", "Estimate"],
+    loo_ss_lpd$estimates["looic", "Estimate"],
+    tolerance = 0.00000001
+  ))
+
+  expect_silent(
+    loo_ss_lpd10 <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 500,
+      loo_approximation = "lpd",
+      loo_approximation_draws = 10,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(loo_ss_lpd10, "psis_loo_ss")
 
   z <- 2
-  expect_lte(loo_ss_lpd10$estimates["elpd_loo", "Estimate"] - z * loo_ss_lpd10$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_gte(loo_ss_lpd10$estimates["elpd_loo", "Estimate"] + z * loo_ss_lpd10$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_lte(loo_ss_lpd10$estimates["p_loo", "Estimate"] - z * loo_ss_lpd10$estimates["p_loo", "subsampling SE"], true_loo$estimates["p_loo", "Estimate"])
-  expect_gte(loo_ss_lpd10$estimates["p_loo", "Estimate"] + z * loo_ss_lpd10$estimates["p_loo", "subsampling SE"], true_loo$estimates["p_loo", "Estimate"])
-  expect_lte(loo_ss_lpd10$estimates["looic", "Estimate"] - z * loo_ss_lpd10$estimates["looic", "subsampling SE"], true_loo$estimates["looic", "Estimate"])
-  expect_gte(loo_ss_lpd10$estimates["looic", "Estimate"] + z * loo_ss_lpd10$estimates["looic", "subsampling SE"], true_loo$estimates["looic", "Estimate"])
-
-  expect_failure(expect_equal(true_loo$estimates["elpd_loo", "Estimate"], loo_ss_lpd10$estimates["elpd_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(true_loo$estimates["p_loo", "Estimate"], loo_ss_lpd10$estimates["p_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(true_loo$estimates["looic", "Estimate"], loo_ss_lpd10$estimates["looic", "Estimate"], tolerance = 0.00000001))
+  expect_lte(
+    loo_ss_lpd10$estimates["elpd_loo", "Estimate"] -
+      z * loo_ss_lpd10$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss_lpd10$estimates["elpd_loo", "Estimate"] +
+      z * loo_ss_lpd10$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss_lpd10$estimates["p_loo", "Estimate"] -
+      z * loo_ss_lpd10$estimates["p_loo", "subsampling SE"],
+    true_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss_lpd10$estimates["p_loo", "Estimate"] +
+      z * loo_ss_lpd10$estimates["p_loo", "subsampling SE"],
+    true_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss_lpd10$estimates["looic", "Estimate"] -
+      z * loo_ss_lpd10$estimates["looic", "subsampling SE"],
+    true_loo$estimates["looic", "Estimate"]
+  )
+  expect_gte(
+    loo_ss_lpd10$estimates["looic", "Estimate"] +
+      z * loo_ss_lpd10$estimates["looic", "subsampling SE"],
+    true_loo$estimates["looic", "Estimate"]
+  )
+
+  expect_failure(expect_equal(
+    true_loo$estimates["elpd_loo", "Estimate"],
+    loo_ss_lpd10$estimates["elpd_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    true_loo$estimates["p_loo", "Estimate"],
+    loo_ss_lpd10$estimates["p_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    true_loo$estimates["looic", "Estimate"],
+    loo_ss_lpd10$estimates["looic", "Estimate"],
+    tolerance = 0.00000001
+  ))
 
   # Test conversion of objects
   expect_silent(true_loo_2 <- loo:::as.psis_loo.psis_loo(true_loo))
@@ -83,244 +251,645 @@ test_that("overall loo_subampling works as expected (compared with loo) for diff
   expect_failure(expect_s3_class(true_loo_conv, "psis_loo_ss"))
   expect_equal(true_loo_conv, true_loo)
   expect_error(loo:::as.psis_loo.psis_loo_ss(loo_ss))
-
 })
 
 test_that("loo with subsampling of all observations works as ordinary loo.", {
   set.seed(123)
-  N <- 1000; K <- 10; S <- 1000; a0 <- 3; b0 <- 2
+  N <- 1000
+  K <- 10
+  S <- 1000
+  a0 <- 3
+  b0 <- 2
   p <- 0.7
   y <- rbinom(N, size = K, prob = p)
-  a <- a0 + sum(y); b <- b0 + N * K - sum(y)
+  a <- a0 + sum(y)
+  b <- b0 + N * K - sum(y)
   fake_posterior <- as.matrix(rbeta(S, a, b))
-  fake_data <- data.frame(y,K)
+  fake_data <- data.frame(y, K)
   rm(N, K, S, a0, b0, p, y, a, b)
   llfun_test <- function(data_i, draws) {
     dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE)
   }
 
-  expect_silent(true_loo <- loo(llfun_test, draws = fake_posterior, data = fake_data, r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    true_loo <- loo(
+      llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(true_loo, "psis_loo")
-  expect_silent(loo_ss <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 1000, loo_approximation = "plpd", r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 1000,
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(loo_ss, "psis_loo_ss")
-  expect_error(loo_ss <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 1001, loo_approximation = "plpd", r_eff = rep(1, nrow(fake_data))))
-
-  expect_equal(true_loo$estimates["elpd_loo", "Estimate"], loo_ss$estimates["elpd_loo", "Estimate"], tolerance = 0.00000001)
-  expect_equal(true_loo$estimates["p_loo", "Estimate"], loo_ss$estimates["p_loo", "Estimate"], tolerance = 0.00000001)
-  expect_equal(true_loo$estimates["looic", "Estimate"], loo_ss$estimates["looic", "Estimate"], tolerance = 0.00000001)
-
-  expect_equal(dim(true_loo),dim(loo_ss))
+  expect_error(
+    loo_ss <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 1001,
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+
+  expect_equal(
+    true_loo$estimates["elpd_loo", "Estimate"],
+    loo_ss$estimates["elpd_loo", "Estimate"],
+    tolerance = 0.00000001
+  )
+  expect_equal(
+    true_loo$estimates["p_loo", "Estimate"],
+    loo_ss$estimates["p_loo", "Estimate"],
+    tolerance = 0.00000001
+  )
+  expect_equal(
+    true_loo$estimates["looic", "Estimate"],
+    loo_ss$estimates["looic", "Estimate"],
+    tolerance = 0.00000001
+  )
+
+  expect_equal(dim(true_loo), dim(loo_ss))
   expect_equal(true_loo$diagnostics, loo_ss$diagnostics)
   expect_equal(max(loo_ss$pointwise[, "m_i"]), 1)
-
 })
 
 test_that("overall loo_subsample works with diff_srs as expected (compared with loo)", {
   set.seed(123)
-  N <- 1000; K <- 10; S <- 1000; a0 <- 3; b0 <- 2
+  N <- 1000
+  K <- 10
+  S <- 1000
+  a0 <- 3
+  b0 <- 2
   p <- 0.7
   y <- rbinom(N, size = K, prob = p)
-  a <- a0 + sum(y); b <- b0 + N * K - sum(y)
+  a <- a0 + sum(y)
+  b <- b0 + N * K - sum(y)
   fake_posterior <- as.matrix(rbeta(S, a, b))
-  fake_data <- data.frame(y,K)
+  fake_data <- data.frame(y, K)
   rm(N, K, S, a0, b0, p, y, a, b)
   llfun_test <- function(data_i, draws) {
     dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE)
   }
 
-  expect_silent(true_loo <- loo(x = llfun_test, draws = fake_posterior, data = fake_data, r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 200, loo_approximation = "plpd", estimator = "diff_srs", r_eff  = rep(1, nrow(fake_data))))
-  expect_equal(true_loo$estimates[1,1], loo_ss$estimates[1,1], tolerance = 0.1)
-
+  expect_silent(
+    true_loo <- loo(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_silent(
+    loo_ss <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 200,
+      loo_approximation = "plpd",
+      estimator = "diff_srs",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_equal(
+    true_loo$estimates[1, 1],
+    loo_ss$estimates[1, 1],
+    tolerance = 0.1
+  )
 })
 
 test_that("Test the srs estimator with 'none' approximation", {
   set.seed(123)
-  N <- 1000; K <- 10; S <- 1000; a0 <- 3; b0 <- 2
+  N <- 1000
+  K <- 10
+  S <- 1000
+  a0 <- 3
+  b0 <- 2
   p <- 0.7
   y <- rbinom(N, size = K, prob = p)
-  a <- a0 + sum(y); b <- b0 + N * K - sum(y)
+  a <- a0 + sum(y)
+  b <- b0 + N * K - sum(y)
   fake_posterior <- as.matrix(rbeta(S, a, b))
-  fake_data <- data.frame(y,K)
+  fake_data <- data.frame(y, K)
   rm(N, K, S, a0, b0, p, y, a, b)
   llfun_test <- function(data_i, draws) {
     dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE)
   }
 
-  expect_silent(true_loo <- loo(llfun_test, draws = fake_posterior, data = fake_data, r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    true_loo <- loo(
+      llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(true_loo, "psis_loo")
-  expect_silent(loo_ss <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 200, loo_approximation = "none", estimator = "srs", r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 200,
+      loo_approximation = "none",
+      estimator = "srs",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(loo_ss, "psis_loo_ss")
-  expect_error(loo_ss <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 1100, loo_approximation = "none", estimator = "srs", r_eff = rep(1, nrow(fake_data))))
+  expect_error(
+    loo_ss <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 1100,
+      loo_approximation = "none",
+      estimator = "srs",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
 
   expect_equal(length(obs_idx(loo_ss)), nobs(loo_ss))
 
-
   # Check consistency
-  expect_equal(loo_ss$pointwise[, "elpd_loo_approx"],
-                    loo_ss$loo_subsampling$elpd_loo_approx[loo_ss$pointwise[, "idx"]], ignore_attr = TRUE)
+  expect_equal(
+    loo_ss$pointwise[, "elpd_loo_approx"],
+    loo_ss$loo_subsampling$elpd_loo_approx[loo_ss$pointwise[, "idx"]],
+    ignore_attr = TRUE
+  )
 
   # Expect values
   z <- 2
-  expect_lte(loo_ss$estimates["elpd_loo", "Estimate"] - z * loo_ss$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_gte(loo_ss$estimates["elpd_loo", "Estimate"] + z * loo_ss$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_lte(loo_ss$estimates["p_loo", "Estimate"] - z * loo_ss$estimates["p_loo", "subsampling SE"], true_loo$estimates["p_loo", "Estimate"])
-  expect_gte(loo_ss$estimates["p_loo", "Estimate"] + z * loo_ss$estimates["p_loo", "subsampling SE"], true_loo$estimates["p_loo", "Estimate"])
-  expect_lte(loo_ss$estimates["looic", "Estimate"] - z * loo_ss$estimates["looic", "subsampling SE"], true_loo$estimates["looic", "Estimate"])
-  expect_gte(loo_ss$estimates["looic", "Estimate"] + z * loo_ss$estimates["looic", "subsampling SE"], true_loo$estimates["looic", "Estimate"])
-
-  expect_failure(expect_equal(true_loo$estimates["elpd_loo", "Estimate"], loo_ss$estimates["elpd_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(true_loo$estimates["p_loo", "Estimate"], loo_ss$estimates["p_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(true_loo$estimates["looic", "Estimate"], loo_ss$estimates["looic", "Estimate"], tolerance = 0.00000001))
-
+  expect_lte(
+    loo_ss$estimates["elpd_loo", "Estimate"] -
+      z * loo_ss$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["elpd_loo", "Estimate"] +
+      z * loo_ss$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss$estimates["p_loo", "Estimate"] -
+      z * loo_ss$estimates["p_loo", "subsampling SE"],
+    true_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["p_loo", "Estimate"] +
+      z * loo_ss$estimates["p_loo", "subsampling SE"],
+    true_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss$estimates["looic", "Estimate"] -
+      z * loo_ss$estimates["looic", "subsampling SE"],
+    true_loo$estimates["looic", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["looic", "Estimate"] +
+      z * loo_ss$estimates["looic", "subsampling SE"],
+    true_loo$estimates["looic", "Estimate"]
+  )
+
+  expect_failure(expect_equal(
+    true_loo$estimates["elpd_loo", "Estimate"],
+    loo_ss$estimates["elpd_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    true_loo$estimates["p_loo", "Estimate"],
+    loo_ss$estimates["p_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    true_loo$estimates["looic", "Estimate"],
+    loo_ss$estimates["looic", "Estimate"],
+    tolerance = 0.00000001
+  ))
 })
 
 test_that("Test the Hansen-Hurwitz estimator", {
   set.seed(123)
-  N <- 1000; K <- 10; S <- 1000; a0 <- 3; b0 <- 2
+  N <- 1000
+  K <- 10
+  S <- 1000
+  a0 <- 3
+  b0 <- 2
   p <- 0.7
   y <- rbinom(N, size = K, prob = p)
-  a <- a0 + sum(y); b <- b0 + N * K - sum(y)
+  a <- a0 + sum(y)
+  b <- b0 + N * K - sum(y)
   fake_posterior <- as.matrix(rbeta(S, a, b))
-  fake_data <- data.frame(y,K)
+  fake_data <- data.frame(y, K)
   rm(N, K, S, a0, b0, p, y, a, b)
   llfun_test <- function(data_i, draws) {
     dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE)
   }
 
-  expect_silent(true_loo <- loo(llfun_test, draws = fake_posterior, data = fake_data, r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    true_loo <- loo(
+      llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(true_loo, "psis_loo")
-  expect_silent(loo_ss <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 300, loo_approximation = "plpd", estimator = "hh_pps", r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 300,
+      loo_approximation = "plpd",
+      estimator = "hh_pps",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(loo_ss, "psis_loo_ss")
-  expect_silent(loo_ss_max <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 1100, loo_approximation = "plpd", estimator = "hh_pps", r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss_max <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 1100,
+      loo_approximation = "plpd",
+      estimator = "hh_pps",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(loo_ss_max, "psis_loo_ss")
-  expect_silent(loo_ss_max2 <- update(loo_ss, draws = fake_posterior, data = fake_data, observations = 1100, r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss_max2 <- update(
+      loo_ss,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 1100,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_equal(nobs(loo_ss_max2), 1100)
   expect_gt(max(loo_ss_max2$pointwise[, "m_i"]), 1)
-  expect_error(loo_ss_max2 <- update(loo_ss_max2, draws = fake_posterior, data = fake_data, observations = 300, r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss_max3 <- update(loo_ss, draws = fake_posterior, data = fake_data, observations = 1500, r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss2 <- update(loo_ss, draws = fake_posterior, data = fake_data, observations = loo_ss, r_eff = rep(1, nrow(fake_data))))
-  expect_error(loo_ss2 <- update(loo_ss, draws = fake_posterior, data = fake_data, observations = loo_ss, loo_approximation = "lpd", r_eff = rep(1, nrow(fake_data))))
+  expect_error(
+    loo_ss_max2 <- update(
+      loo_ss_max2,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 300,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_silent(
+    loo_ss_max3 <- update(
+      loo_ss,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 1500,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_silent(
+    loo_ss2 <- update(
+      loo_ss,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = loo_ss,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_error(
+    loo_ss2 <- update(
+      loo_ss,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = loo_ss,
+      loo_approximation = "lpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_equal(loo_ss$estimates, loo_ss2$estimates)
   expect_equal(length(obs_idx(loo_ss_max)), length(obs_idx(loo_ss_max2)))
   expect_equal(length(obs_idx(loo_ss_max)), nobs(loo_ss_max))
 
-
   # Check consistency
-  expect_equal(loo_ss$pointwise[, "elpd_loo_approx"],
-                    loo_ss$loo_subsampling$elpd_loo_approx[loo_ss$pointwise[, "idx"]], ignore_attr = TRUE)
+  expect_equal(
+    loo_ss$pointwise[, "elpd_loo_approx"],
+    loo_ss$loo_subsampling$elpd_loo_approx[loo_ss$pointwise[, "idx"]],
+    ignore_attr = TRUE
+  )
   # Check consistency
-  expect_equal(loo_ss_max$pointwise[, "elpd_loo_approx"],
-                    loo_ss_max$loo_subsampling$elpd_loo_approx[loo_ss_max$pointwise[, "idx"]], ignore_attr = TRUE)
+  expect_equal(
+    loo_ss_max$pointwise[, "elpd_loo_approx"],
+    loo_ss_max$loo_subsampling$elpd_loo_approx[loo_ss_max$pointwise[, "idx"]],
+    ignore_attr = TRUE
+  )
 
   # Expect values
   z <- 2
-  expect_lte(loo_ss$estimates["elpd_loo", "Estimate"] - z * loo_ss$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_gte(loo_ss$estimates["elpd_loo", "Estimate"] + z * loo_ss$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_lte(loo_ss$estimates["p_loo", "Estimate"] - z * loo_ss$estimates["p_loo", "subsampling SE"], true_loo$estimates["p_loo", "Estimate"])
-  expect_gte(loo_ss$estimates["p_loo", "Estimate"] + z * loo_ss$estimates["p_loo", "subsampling SE"], true_loo$estimates["p_loo", "Estimate"])
-  expect_lte(loo_ss$estimates["looic", "Estimate"] - z * loo_ss$estimates["looic", "subsampling SE"], true_loo$estimates["looic", "Estimate"])
-  expect_gte(loo_ss$estimates["looic", "Estimate"] + z * loo_ss$estimates["looic", "subsampling SE"], true_loo$estimates["looic", "Estimate"])
-
-  expect_failure(expect_equal(true_loo$estimates["elpd_loo", "Estimate"], loo_ss$estimates["elpd_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(true_loo$estimates["p_loo", "Estimate"], loo_ss$estimates["p_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(true_loo$estimates["looic", "Estimate"], loo_ss$estimates["looic", "Estimate"], tolerance = 0.00000001))
-
-  expect_lte(loo_ss_max$estimates["elpd_loo", "Estimate"] - z * loo_ss_max$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-  expect_gte(loo_ss_max$estimates["elpd_loo", "Estimate"] + z * loo_ss_max$estimates["elpd_loo", "subsampling SE"], true_loo$estimates["elpd_loo", "Estimate"])
-
+  expect_lte(
+    loo_ss$estimates["elpd_loo", "Estimate"] -
+      z * loo_ss$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["elpd_loo", "Estimate"] +
+      z * loo_ss$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss$estimates["p_loo", "Estimate"] -
+      z * loo_ss$estimates["p_loo", "subsampling SE"],
+    true_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["p_loo", "Estimate"] +
+      z * loo_ss$estimates["p_loo", "subsampling SE"],
+    true_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss$estimates["looic", "Estimate"] -
+      z * loo_ss$estimates["looic", "subsampling SE"],
+    true_loo$estimates["looic", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["looic", "Estimate"] +
+      z * loo_ss$estimates["looic", "subsampling SE"],
+    true_loo$estimates["looic", "Estimate"]
+  )
+
+  expect_failure(expect_equal(
+    true_loo$estimates["elpd_loo", "Estimate"],
+    loo_ss$estimates["elpd_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    true_loo$estimates["p_loo", "Estimate"],
+    loo_ss$estimates["p_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    true_loo$estimates["looic", "Estimate"],
+    loo_ss$estimates["looic", "Estimate"],
+    tolerance = 0.00000001
+  ))
+
+  expect_lte(
+    loo_ss_max$estimates["elpd_loo", "Estimate"] -
+      z * loo_ss_max$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss_max$estimates["elpd_loo", "Estimate"] +
+      z * loo_ss_max$estimates["elpd_loo", "subsampling SE"],
+    true_loo$estimates["elpd_loo", "Estimate"]
+  )
 })
 
 
 test_that("update.psis_loo_ss works as expected (compared with loo)", {
-
-
   set.seed(123)
-  N <- 1000; K <- 10; S <- 1000; a0 <- 3; b0 <- 2
+  N <- 1000
+  K <- 10
+  S <- 1000
+  a0 <- 3
+  b0 <- 2
   p <- 0.7
   y <- rbinom(N, size = K, prob = p)
-  a <- a0 + sum(y); b <- b0 + N * K - sum(y)
+  a <- a0 + sum(y)
+  b <- b0 + N * K - sum(y)
   fake_posterior <- as.matrix(rbeta(S, a, b))
-  fake_data <- data.frame(y,K)
+  fake_data <- data.frame(y, K)
   rm(N, K, S, a0, b0, p, y, a, b)
   llfun_test <- function(data_i, draws) {
     dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE)
   }
 
-  expect_silent(true_loo <- loo(llfun_test, draws = fake_posterior, data = fake_data, r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    true_loo <- loo(
+      llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(true_loo, "psis_loo")
-  expect_silent(loo_ss <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 500, loo_approximation = "plpd", r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 500,
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_s3_class(loo_ss, "psis_loo_ss")
 
   # Check error when draws and data dimensions differ
-  expect_error(loo_ss2 <- update(object = loo_ss, draws = cbind(fake_posterior, 1), data = fake_data, observations = 600, r_eff = rep(1, nrow(fake_data))))
-  expect_error(loo_ss2 <- update(object = loo_ss, draws = fake_posterior, data = fake_data[-1,], observations = 600, r_eff = rep(1, nrow(fake_data))))
+  expect_error(
+    loo_ss2 <- update(
+      object = loo_ss,
+      draws = cbind(fake_posterior, 1),
+      data = fake_data,
+      observations = 600,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_error(
+    loo_ss2 <- update(
+      object = loo_ss,
+      draws = fake_posterior,
+      data = fake_data[-1, ],
+      observations = 600,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
 
   # Add tests for adding observations
-  expect_silent(loo_ss2 <- update(object = loo_ss, draws = fake_posterior, data = fake_data, observations = 600, r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss2 <- update(
+      object = loo_ss,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 600,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_equal(dim(loo_ss2)[2] - dim(loo_ss)[2], expected = 100)
   expect_equal(dim(loo_ss2)[2], expected = dim(loo_ss2$pointwise)[1])
   expect_length(loo_ss2$diagnostics$pareto_k, 600)
   expect_length(loo_ss2$diagnostics$n_eff, 600)
-  for(i in 1:nrow(loo_ss2$estimates)) {
-    expect_lt(loo_ss2$estimates[i, "subsampling SE"],
-              loo_ss$estimates[i, "subsampling SE"])
+  for (i in 1:nrow(loo_ss2$estimates)) {
+    expect_lt(
+      loo_ss2$estimates[i, "subsampling SE"],
+      loo_ss$estimates[i, "subsampling SE"]
+    )
   }
 
-  expect_silent(loo_ss2b <- update(object = loo_ss, draws = fake_posterior, data = fake_data))
+  expect_silent(
+    loo_ss2b <- update(
+      object = loo_ss,
+      draws = fake_posterior,
+      data = fake_data
+    )
+  )
   expect_equal(loo_ss2b$estimates, loo_ss$estimates)
   expect_equal(loo_ss2b$pointwise, loo_ss$pointwise)
   expect_equal(loo_ss2b$diagnostics$pareto_k, loo_ss$diagnostics$pareto_k)
   expect_equal(loo_ss2b$diagnostics$n_eff, loo_ss$diagnostics$n_eff)
 
-  expect_silent(loo_ss3 <- update(object = loo_ss2, draws = fake_posterior, data = fake_data, observations = loo_ss))
+  expect_silent(
+    loo_ss3 <- update(
+      object = loo_ss2,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = loo_ss
+    )
+  )
   expect_equal(loo_ss3$estimates, loo_ss$estimates)
   expect_equal(loo_ss3$pointwise, loo_ss$pointwise)
   expect_equal(loo_ss3$diagnostics$pareto_k, loo_ss$diagnostics$pareto_k)
   expect_equal(loo_ss3$diagnostics$n_eff, loo_ss$diagnostics$n_eff)
 
-  expect_silent(loo_ss4 <- update(object = loo_ss, draws = fake_posterior, data = fake_data, observations = 1000, r_eff = rep(1, nrow(fake_data))))
-  expect_equal(loo_ss4$estimates[,1], true_loo$estimates[,1])
-  expect_equal(loo_ss4$estimates[,2], true_loo$estimates[,2], tolerance = 0.001)
-
-  expect_silent(loo_ss5 <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 1000, loo_approximation = "plpd", r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss4 <- update(
+      object = loo_ss,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 1000,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_equal(loo_ss4$estimates[, 1], true_loo$estimates[, 1])
+  expect_equal(
+    loo_ss4$estimates[, 2],
+    true_loo$estimates[, 2],
+    tolerance = 0.001
+  )
+
+  expect_silent(
+    loo_ss5 <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 1000,
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
 
   ss4_order <- order(loo_ss4$pointwise[, "idx"])
-  expect_equal(loo_ss4$pointwise[ss4_order,c(1,3,4)], loo_ss5$pointwise[,c(1,3,4)])
-  expect_equal(loo_ss4$diagnostics$pareto_k[ss4_order], loo_ss5$diagnostics$pareto_k)
+  expect_equal(
+    loo_ss4$pointwise[ss4_order, c(1, 3, 4)],
+    loo_ss5$pointwise[, c(1, 3, 4)]
+  )
+  expect_equal(
+    loo_ss4$diagnostics$pareto_k[ss4_order],
+    loo_ss5$diagnostics$pareto_k
+  )
   expect_equal(loo_ss4$diagnostics$n_eff[ss4_order], loo_ss5$diagnostics$n_eff)
-  expect_equal(loo_ss4$pointwise[ss4_order,c(1,3,4)], true_loo$pointwise[,c(1,3,4)])
-  expect_equal(loo_ss4$diagnostics$pareto_k[ss4_order], true_loo$diagnostics$pareto_k)
+  expect_equal(
+    loo_ss4$pointwise[ss4_order, c(1, 3, 4)],
+    true_loo$pointwise[, c(1, 3, 4)]
+  )
+  expect_equal(
+    loo_ss4$diagnostics$pareto_k[ss4_order],
+    true_loo$diagnostics$pareto_k
+  )
   expect_equal(loo_ss4$diagnostics$n_eff[ss4_order], true_loo$diagnostics$n_eff)
 
-  expect_error(loo_ss_min <- update(object = loo_ss, draws = fake_posterior, data = fake_data, observations = 50, r_eff = rep(1, nrow(fake_data))))
+  expect_error(
+    loo_ss_min <- update(
+      object = loo_ss,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 50,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
 
   expect_silent(true_loo_ss <- loo:::as.psis_loo_ss.psis_loo(true_loo))
-  expect_silent(loo_ss_subset0 <- update(true_loo_ss, observations = loo_ss, r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss_subset0 <- update(
+      true_loo_ss,
+      observations = loo_ss,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_true(identical(obs_idx(loo_ss_subset0), obs_idx(loo_ss)))
-  expect_silent(loo_ss_subset1 <- update(object = loo_ss, observations = loo_ss, r_eff = rep(1, nrow(fake_data))))
-  expect_message(loo_ss_subset2 <- update(object = loo_ss, observations = obs_idx(loo_ss)[1:10], r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss_subset1 <- update(
+      object = loo_ss,
+      observations = loo_ss,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_message(
+    loo_ss_subset2 <- update(
+      object = loo_ss,
+      observations = obs_idx(loo_ss)[1:10],
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_equal(nobs(loo_ss_subset2), 10)
 
-
   expect_silent(true_loo_ss <- loo:::as.psis_loo_ss.psis_loo(true_loo))
   set.seed(4711)
-  expect_silent(loo_ss2 <- update(object = loo_ss, draws = fake_posterior, data = fake_data, observations = 600, r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss2_subset0 <- update(object = true_loo_ss, observations = loo_ss2, r_eff = rep(1, nrow(fake_data))))
+  expect_silent(
+    loo_ss2 <- update(
+      object = loo_ss,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 600,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_silent(
+    loo_ss2_subset0 <- update(
+      object = true_loo_ss,
+      observations = loo_ss2,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
   expect_true(setequal(obs_idx(loo_ss2), obs_idx(loo_ss2_subset0)))
   expect_true(identical(obs_idx(loo_ss2), obs_idx(loo_ss2_subset0)))
   expect_true(identical(loo_ss2$diagnostic, loo_ss2_subset0$diagnostic))
 
   # Add tests for changing approx variable
-  expect_silent(loo_ss_lpd <- update(object = loo_ss, draws = fake_posterior, data = fake_data, loo_approximation = "lpd", r_eff = rep(1, nrow(fake_data))))
-  expect_failure(expect_equal(loo_ss_lpd$loo_subsampling$elpd_loo_approx, loo_ss$loo_subsampling$elpd_loo_approx))
+  expect_silent(
+    loo_ss_lpd <- update(
+      object = loo_ss,
+      draws = fake_posterior,
+      data = fake_data,
+      loo_approximation = "lpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_failure(expect_equal(
+    loo_ss_lpd$loo_subsampling$elpd_loo_approx,
+    loo_ss$loo_subsampling$elpd_loo_approx
+  ))
   expect_equal(dim(loo_ss_lpd)[2], dim(loo_ss)[2])
   expect_equal(dim(loo_ss_lpd)[2], dim(loo_ss_lpd$pointwise)[1])
   expect_length(loo_ss_lpd$diagnostics$pareto_k, 500)
   expect_length(loo_ss_lpd$diagnostics$n_eff, 500)
-  expect_failure(expect_equal(loo_ss_lpd$estimates[1, "subsampling SE"], loo_ss$estimates[1, "subsampling SE"]))
-  expect_failure(expect_equal(loo_ss_lpd$estimates[3, "subsampling SE"], loo_ss$estimates[3, "subsampling SE"]))
-
+  expect_failure(expect_equal(
+    loo_ss_lpd$estimates[1, "subsampling SE"],
+    loo_ss$estimates[1, "subsampling SE"]
+  ))
+  expect_failure(expect_equal(
+    loo_ss_lpd$estimates[3, "subsampling SE"],
+    loo_ss$estimates[3, "subsampling SE"]
+  ))
 })
 
 test_that("loo_compare_subsample", {
@@ -332,69 +901,169 @@ test_that("loo_compare_subsample", {
   x2 <- rnorm(N)
   x3 <- rnorm(N)
   sigma <- 2
-  y <- rnorm(N, 1 + 2*x1 - 2*x2 - 1*x3, sd = sigma)
+  y <- rnorm(N, 1 + 2 * x1 - 2 * x2 - 1 * x3, sd = sigma)
   X <- cbind("x0" = rep(1, N), x1, x2, x3)
 
   # Generate samples from posterior
-  samples_blin <- function(X, y, sigma, draws = 1000){
-    XtX <- t(X)%*%X
+  samples_blin <- function(X, y, sigma, draws = 1000) {
+    XtX <- t(X) %*% X
     b_hat <- solve(XtX) %*% (t(X) %*% y)
     Lambda_n = XtX + diag(ncol(X))
-    mu_n <- solve(Lambda_n) %*% (XtX %*% b_hat + diag(ncol(X)) %*% rep(0,ncol(X)))
+    mu_n <- solve(Lambda_n) %*%
+      (XtX %*% b_hat + diag(ncol(X)) %*% rep(0, ncol(X)))
     L <- t(chol(sigma^2 * solve(Lambda_n)))
-    draws_mat <- matrix(0, ncol=ncol(X), nrow = draws)
-    for(i in 1:draws){
+    draws_mat <- matrix(0, ncol = ncol(X), nrow = draws)
+    for (i in 1:draws) {
       z <- rnorm(length(mu_n))
-      draws_mat[i,] <- L %*% z + mu_n
+      draws_mat[i, ] <- L %*% z + mu_n
     }
     draws_mat
   }
 
-  fake_posterior1 <- samples_blin(X[,1:2], y, sigma, draws = 1000)
-  fake_posterior2 <- samples_blin(X[,1:3], y, sigma, draws = 1000)
+  fake_posterior1 <- samples_blin(X[, 1:2], y, sigma, draws = 1000)
+  fake_posterior2 <- samples_blin(X[, 1:3], y, sigma, draws = 1000)
   fake_posterior3 <- samples_blin(X, y, sigma, draws = 1000)
 
-  fake_data1 <- data.frame(y, X[,1:2])
-  fake_data2 <- data.frame(y, X[,1:3])
+  fake_data1 <- data.frame(y, X[, 1:2])
+  fake_data2 <- data.frame(y, X[, 1:3])
   fake_data3 <- data.frame(y, X)
 
   llfun_test <- function(data_i, draws) {
-    dnorm(x = data_i$y, mean = draws %*% t(data_i[,-1, drop=FALSE]), sd = sigma, log = TRUE)
+    dnorm(
+      x = data_i$y,
+      mean = draws %*% t(data_i[, -1, drop = FALSE]),
+      sd = sigma,
+      log = TRUE
+    )
   }
 
-  expect_silent(l1 <- loo(llfun_test, data = fake_data1, draws = fake_posterior1, r_eff = rep(1, N)))
-  expect_silent(l2 <- loo(llfun_test, data = fake_data2, draws = fake_posterior2, r_eff = rep(1, N)))
-  expect_silent(l3 <- loo(llfun_test, data = fake_data3, draws = fake_posterior3, r_eff = rep(1, N)))
-
-  expect_silent(lss1 <- loo_subsample(llfun_test, data = fake_data1, draws = fake_posterior1, observations = 100, r_eff = rep(1, N)))
-  expect_silent(lss2 <- loo_subsample(llfun_test, data = fake_data2, draws = fake_posterior2, observations = 100, r_eff = rep(1, N)))
-  expect_silent(lss3 <- loo_subsample(llfun_test, data = fake_data3, draws = fake_posterior3, observations = 100, r_eff = rep(1, N)))
-  expect_silent(lss2o1 <- loo_subsample(llfun_test, data = fake_data2, draws = fake_posterior2, observations = lss1, r_eff = rep(1, N)))
-  expect_silent(lss3o1 <- loo_subsample(llfun_test, data = fake_data3, draws = fake_posterior3, observations = lss1, r_eff = rep(1, N)))
-  expect_silent(lss2hh <- loo_subsample(llfun_test, data = fake_data2, draws = fake_posterior2, observations = 100, estimator = "hh_pps", r_eff = rep(1, N)))
-
-  expect_snapshot(lcss <- loo:::loo_compare.psis_loo_ss_list(x = list(lss1, lss2, lss3)))
-  expect_warning(lcss2 <- loo:::loo_compare.psis_loo_ss_list(x = list(lss1, lss2, lss3o1)))
-  expect_silent(lcsso <- loo:::loo_compare.psis_loo_ss_list(x = list(lss1, lss2o1, lss3o1)))
-  expect_warning(lcssohh <- loo:::loo_compare.psis_loo_ss_list(x = list(lss1, lss2hh, lss3o1)))
-  expect_message(lcssf1 <- loo:::loo_compare.psis_loo_ss_list(x = list(loo:::as.psis_loo_ss.psis_loo(l1), lss2o1, lss3o1)))
-  expect_message(lcssf2 <- loo:::loo_compare.psis_loo_ss_list(x = list(loo:::as.psis_loo_ss.psis_loo(l1), lss2o1, loo:::as.psis_loo_ss.psis_loo(l3))))
-
-  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_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_silent(lcss2m <- loo:::loo_compare.psis_loo_ss_list(x = list(lss2o1, lss3o1)))
-  expect_equal(unname(lcss2m[,]), unname(lcsso[1:2,]))
+  expect_silent(
+    l1 <- loo(
+      llfun_test,
+      data = fake_data1,
+      draws = fake_posterior1,
+      r_eff = rep(1, N)
+    )
+  )
+  expect_silent(
+    l2 <- loo(
+      llfun_test,
+      data = fake_data2,
+      draws = fake_posterior2,
+      r_eff = rep(1, N)
+    )
+  )
+  expect_silent(
+    l3 <- loo(
+      llfun_test,
+      data = fake_data3,
+      draws = fake_posterior3,
+      r_eff = rep(1, N)
+    )
+  )
+
+  expect_silent(
+    lss1 <- loo_subsample(
+      llfun_test,
+      data = fake_data1,
+      draws = fake_posterior1,
+      observations = 100,
+      r_eff = rep(1, N)
+    )
+  )
+  expect_silent(
+    lss2 <- loo_subsample(
+      llfun_test,
+      data = fake_data2,
+      draws = fake_posterior2,
+      observations = 100,
+      r_eff = rep(1, N)
+    )
+  )
+  expect_silent(
+    lss3 <- loo_subsample(
+      llfun_test,
+      data = fake_data3,
+      draws = fake_posterior3,
+      observations = 100,
+      r_eff = rep(1, N)
+    )
+  )
+  expect_silent(
+    lss2o1 <- loo_subsample(
+      llfun_test,
+      data = fake_data2,
+      draws = fake_posterior2,
+      observations = lss1,
+      r_eff = rep(1, N)
+    )
+  )
+  expect_silent(
+    lss3o1 <- loo_subsample(
+      llfun_test,
+      data = fake_data3,
+      draws = fake_posterior3,
+      observations = lss1,
+      r_eff = rep(1, N)
+    )
+  )
+  expect_silent(
+    lss2hh <- loo_subsample(
+      llfun_test,
+      data = fake_data2,
+      draws = fake_posterior2,
+      observations = 100,
+      estimator = "hh_pps",
+      r_eff = rep(1, N)
+    )
+  )
+
+  expect_snapshot(
+    lcss <- loo:::loo_compare.psis_loo_ss_list(x = list(lss1, lss2, lss3))
+  )
+  expect_warning(
+    lcss2 <- loo:::loo_compare.psis_loo_ss_list(x = list(lss1, lss2, lss3o1))
+  )
+  expect_silent(
+    lcsso <- loo:::loo_compare.psis_loo_ss_list(x = list(lss1, lss2o1, lss3o1))
+  )
+  expect_warning(
+    lcssohh <- loo:::loo_compare.psis_loo_ss_list(
+      x = list(lss1, lss2hh, lss3o1)
+    )
+  )
+  expect_message(
+    lcssf1 <- loo:::loo_compare.psis_loo_ss_list(
+      x = list(loo:::as.psis_loo_ss.psis_loo(l1), lss2o1, lss3o1)
+    )
+  )
+  expect_message(
+    lcssf2 <- loo:::loo_compare.psis_loo_ss_list(
+      x = list(
+        loo:::as.psis_loo_ss.psis_loo(l1),
+        lss2o1,
+        loo:::as.psis_loo_ss.psis_loo(l3)
+      )
+    )
+  )
+
+  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_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_silent(
+    lcss2m <- loo:::loo_compare.psis_loo_ss_list(x = list(lss2o1, lss3o1))
+  )
+  expect_equal(unname(lcss2m[,]), unname(lcsso[1:2, ]))
 
   expect_snapshot(lcssapi <- loo_compare(lss1, lss2, lss3))
   expect_equal(lcssapi, lcss)
@@ -402,98 +1071,181 @@ test_that("loo_compare_subsample", {
   expect_equal(lcssohhapi, lcssohh)
   expect_silent(lcss2mapi <- loo_compare(lss2o1, lss3o1))
   expect_equal(lcss2mapi, lcss2m)
-
 })
 
 test_that("Test 'tis' and 'sis'", {
   skip_on_cran()
 
   set.seed(123)
-  N <- 1000; K <- 10; S <- 1000; a0 <- 3; b0 <- 2
+  N <- 1000
+  K <- 10
+  S <- 1000
+  a0 <- 3
+  b0 <- 2
   p <- 0.7
   y <- rbinom(N, size = K, prob = p)
-  a <- a0 + sum(y); b <- b0 + N * K - sum(y)
+  a <- a0 + sum(y)
+  b <- b0 + N * K - sum(y)
   fake_posterior <- draws <- as.matrix(rbeta(S, a, b))
-  fake_data <- data.frame(y,K)
+  fake_data <- data.frame(y, K)
   rm(N, K, S, a0, b0, p, y, a, b)
   llfun_test <- function(data_i, draws) {
     dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE)
   }
 
-  expect_silent(loo_ss_full <-
-                  loo_subsample(x = llfun_test,
-                                draws = fake_posterior,
-                                data = fake_data,
-                                observations = 1000,
-                                loo_approximation = "plpd",
-                                r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss_plpd <-
-                  loo_subsample(x = llfun_test,
-                                draws = fake_posterior,
-                                data = fake_data,
-                                observations = 100,
-                                loo_approximation = "plpd",
-                                r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss_tis_S1000 <-
-                  loo_subsample(x = llfun_test,
-                                draws = fake_posterior,
-                                data = fake_data,
-                                observations = 100,
-                                loo_approximation = "tis",
-                                r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss_tis_S100 <-
-                  loo_subsample(x = llfun_test,
-                                draws = fake_posterior,
-                                data = fake_data,
-                                observations = 100,
-                                loo_approximation = "tis",
-                                loo_approximation_draws = 100,
-                                r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss_tis_S10 <-
-                  loo_subsample(x = llfun_test,
-                                draws = fake_posterior,
-                                data = fake_data,
-                                observations = 100,
-                                loo_approximation = "tis",
-                                loo_approximation_draws = 10,
-                                r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss_sis_S1000 <-
-                  loo_subsample(x = llfun_test,
-                                draws = fake_posterior,
-                                data = fake_data,
-                                observations = 100,
-                                loo_approximation = "sis",
-                                r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss_sis_S100 <-
-                  loo_subsample(x = llfun_test,
-                                draws = fake_posterior,
-                                data = fake_data,
-                                observations = 100,
-                                loo_approximation = "sis",
-                                loo_approximation_draws = 100,
-                                r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss_sis_S10 <-
-                  loo_subsample(x = llfun_test,
-                                draws = fake_posterior,
-                                data = fake_data,
-                                observations = 100,
-                                loo_approximation = "sis",
-                                loo_approximation_draws = 10,
-                                r_eff = rep(1, nrow(fake_data))))
-
+  expect_silent(
+    loo_ss_full <-
+      loo_subsample(
+        x = llfun_test,
+        draws = fake_posterior,
+        data = fake_data,
+        observations = 1000,
+        loo_approximation = "plpd",
+        r_eff = rep(1, nrow(fake_data))
+      )
+  )
+  expect_silent(
+    loo_ss_plpd <-
+      loo_subsample(
+        x = llfun_test,
+        draws = fake_posterior,
+        data = fake_data,
+        observations = 100,
+        loo_approximation = "plpd",
+        r_eff = rep(1, nrow(fake_data))
+      )
+  )
+  expect_silent(
+    loo_ss_tis_S1000 <-
+      loo_subsample(
+        x = llfun_test,
+        draws = fake_posterior,
+        data = fake_data,
+        observations = 100,
+        loo_approximation = "tis",
+        r_eff = rep(1, nrow(fake_data))
+      )
+  )
+  expect_silent(
+    loo_ss_tis_S100 <-
+      loo_subsample(
+        x = llfun_test,
+        draws = fake_posterior,
+        data = fake_data,
+        observations = 100,
+        loo_approximation = "tis",
+        loo_approximation_draws = 100,
+        r_eff = rep(1, nrow(fake_data))
+      )
+  )
+  expect_silent(
+    loo_ss_tis_S10 <-
+      loo_subsample(
+        x = llfun_test,
+        draws = fake_posterior,
+        data = fake_data,
+        observations = 100,
+        loo_approximation = "tis",
+        loo_approximation_draws = 10,
+        r_eff = rep(1, nrow(fake_data))
+      )
+  )
+  expect_silent(
+    loo_ss_sis_S1000 <-
+      loo_subsample(
+        x = llfun_test,
+        draws = fake_posterior,
+        data = fake_data,
+        observations = 100,
+        loo_approximation = "sis",
+        r_eff = rep(1, nrow(fake_data))
+      )
+  )
+  expect_silent(
+    loo_ss_sis_S100 <-
+      loo_subsample(
+        x = llfun_test,
+        draws = fake_posterior,
+        data = fake_data,
+        observations = 100,
+        loo_approximation = "sis",
+        loo_approximation_draws = 100,
+        r_eff = rep(1, nrow(fake_data))
+      )
+  )
+  expect_silent(
+    loo_ss_sis_S10 <-
+      loo_subsample(
+        x = llfun_test,
+        draws = fake_posterior,
+        data = fake_data,
+        observations = 100,
+        loo_approximation = "sis",
+        loo_approximation_draws = 10,
+        r_eff = rep(1, nrow(fake_data))
+      )
+  )
 
   SEs <- 4
-  expect_gt(loo_ss_tis_S1000$estimates["elpd_loo", "Estimate"] + SEs*loo_ss_tis_S1000$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_lt(loo_ss_tis_S1000$estimates["elpd_loo", "Estimate"] - SEs*loo_ss_tis_S1000$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_gt(loo_ss_tis_S100$estimates["elpd_loo", "Estimate"] + SEs*loo_ss_tis_S100$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_lt(loo_ss_tis_S100$estimates["elpd_loo", "Estimate"] - SEs*loo_ss_tis_S100$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_gt(loo_ss_tis_S10$estimates["elpd_loo", "Estimate"] + SEs*loo_ss_tis_S10$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_lt(loo_ss_tis_S10$estimates["elpd_loo", "Estimate"] - SEs*loo_ss_tis_S10$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-
-  expect_gt(loo_ss_sis_S1000$estimates["elpd_loo", "Estimate"] + SEs*loo_ss_sis_S1000$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_lt(loo_ss_sis_S1000$estimates["elpd_loo", "Estimate"] - SEs*loo_ss_sis_S1000$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_gt(loo_ss_sis_S100$estimates["elpd_loo", "Estimate"] + SEs*loo_ss_sis_S100$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_lt(loo_ss_sis_S100$estimates["elpd_loo", "Estimate"] - SEs*loo_ss_sis_S100$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_gt(loo_ss_sis_S10$estimates["elpd_loo", "Estimate"] + SEs*loo_ss_sis_S10$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_lt(loo_ss_sis_S10$estimates["elpd_loo", "Estimate"] - SEs*loo_ss_sis_S10$estimates["elpd_loo", "subsampling SE"], loo_ss_full$estimates["elpd_loo", "Estimate"])
+  expect_gt(
+    loo_ss_tis_S1000$estimates["elpd_loo", "Estimate"] +
+      SEs * loo_ss_tis_S1000$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lt(
+    loo_ss_tis_S1000$estimates["elpd_loo", "Estimate"] -
+      SEs * loo_ss_tis_S1000$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gt(
+    loo_ss_tis_S100$estimates["elpd_loo", "Estimate"] +
+      SEs * loo_ss_tis_S100$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lt(
+    loo_ss_tis_S100$estimates["elpd_loo", "Estimate"] -
+      SEs * loo_ss_tis_S100$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gt(
+    loo_ss_tis_S10$estimates["elpd_loo", "Estimate"] +
+      SEs * loo_ss_tis_S10$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lt(
+    loo_ss_tis_S10$estimates["elpd_loo", "Estimate"] -
+      SEs * loo_ss_tis_S10$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+
+  expect_gt(
+    loo_ss_sis_S1000$estimates["elpd_loo", "Estimate"] +
+      SEs * loo_ss_sis_S1000$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lt(
+    loo_ss_sis_S1000$estimates["elpd_loo", "Estimate"] -
+      SEs * loo_ss_sis_S1000$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gt(
+    loo_ss_sis_S100$estimates["elpd_loo", "Estimate"] +
+      SEs * loo_ss_sis_S100$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lt(
+    loo_ss_sis_S100$estimates["elpd_loo", "Estimate"] -
+      SEs * loo_ss_sis_S100$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gt(
+    loo_ss_sis_S10$estimates["elpd_loo", "Estimate"] +
+      SEs * loo_ss_sis_S10$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lt(
+    loo_ss_sis_S10$estimates["elpd_loo", "Estimate"] -
+      SEs * loo_ss_sis_S10$estimates["elpd_loo", "subsampling SE"],
+    loo_ss_full$estimates["elpd_loo", "Estimate"]
+  )
 })
diff --git a/tests/testthat/test_loo_subsampling_approximations.R b/tests/testthat/test_loo_subsampling_approximations.R
index 3f25ab2c..e5a02d58 100644
--- a/tests/testthat/test_loo_subsampling_approximations.R
+++ b/tests/testthat/test_loo_subsampling_approximations.R
@@ -1,13 +1,17 @@
-library(loo)
 options(mc.cores = 1)
 
 generate_test_elpd_dataset <- function() {
-  N <- 10; K <- 10; S <- 1000; a0 <- 3; b0 <- 2
+  N <- 10
+  K <- 10
+  S <- 1000
+  a0 <- 3
+  b0 <- 2
   p <- 0.7
   y <- rbinom(N, size = K, prob = p)
-  a <- a0 + sum(y); b <- b0 + N * K - sum(y)
+  a <- a0 + sum(y)
+  b <- b0 + N * K - sum(y)
   fake_posterior <- draws <- as.matrix(rbeta(S, a, b))
-  fake_data <- data.frame(y,K)
+  fake_data <- data.frame(y, K)
   rm(N, K, S, a0, b0, p, y, a, b)
 
   list(fake_posterior = fake_posterior, fake_data = fake_data)
@@ -24,37 +28,78 @@ test_elpd_loo_approximation <- function(cores) {
   }
 
   # Compute plpd approximation
-  expect_silent(pi_vals <- loo:::elpd_loo_approximation(.llfun = llfun_test, data = fake_data, draws = fake_posterior, loo_approximation = "plpd", cores = cores))
+  expect_silent(
+    pi_vals <- loo:::elpd_loo_approximation(
+      .llfun = llfun_test,
+      data = fake_data,
+      draws = fake_posterior,
+      loo_approximation = "plpd",
+      cores = cores
+    )
+  )
   # Compute it manually
   point <- mean(fake_posterior)
   llik <- dbinom(fake_data$y, size = fake_data$K, prob = point, log = TRUE)
   abs_lliks <- abs(llik)
-  man_elpd_loo_approximation <- abs_lliks/sum(abs_lliks)
-  expect_equal(abs(pi_vals)/sum(abs(pi_vals)), man_elpd_loo_approximation, tolerance = 0.00001)
+  man_elpd_loo_approximation <- abs_lliks / sum(abs_lliks)
+  expect_equal(
+    abs(pi_vals) / sum(abs(pi_vals)),
+    man_elpd_loo_approximation,
+    tolerance = 0.00001
+  )
 
   # Compute lpd approximation
-  expect_silent(pi_vals <- loo:::elpd_loo_approximation(.llfun = llfun_test, data = fake_data, draws = fake_posterior, loo_approximation = "lpd", cores = cores))
+  expect_silent(
+    pi_vals <- loo:::elpd_loo_approximation(
+      .llfun = llfun_test,
+      data = fake_data,
+      draws = fake_posterior,
+      loo_approximation = "lpd",
+      cores = cores
+    )
+  )
   # Compute it manually
   llik <- numeric(10)
-  for(i in seq_along(fake_data$y)){
-    llik[i] <- loo:::logMeanExp(dbinom(fake_data$y[i], size = fake_data$K, prob = fake_posterior, log = TRUE))
+  for (i in seq_along(fake_data$y)) {
+    llik[i] <- loo:::logMeanExp(dbinom(
+      fake_data$y[i],
+      size = fake_data$K,
+      prob = fake_posterior,
+      log = TRUE
+    ))
   }
   abs_lliks <- abs(llik)
-  man_approx_loo_variable <- abs_lliks/sum(abs_lliks)
-  expect_equal(abs(pi_vals)/sum(abs(pi_vals)), man_approx_loo_variable, tolerance = 0.00001)
+  man_approx_loo_variable <- abs_lliks / sum(abs_lliks)
+  expect_equal(
+    abs(pi_vals) / sum(abs(pi_vals)),
+    man_approx_loo_variable,
+    tolerance = 0.00001
+  )
 
   # Compute waic approximation
-  expect_silent(pi_vals_waic <- loo:::elpd_loo_approximation(.llfun = llfun_test, data = fake_data, draws = fake_posterior, loo_approximation = "waic", cores = cores))
+  expect_silent(
+    pi_vals_waic <- loo:::elpd_loo_approximation(
+      .llfun = llfun_test,
+      data = fake_data,
+      draws = fake_posterior,
+      loo_approximation = "waic",
+      cores = cores
+    )
+  )
   expect_true(all(pi_vals > pi_vals_waic))
   expect_true(sum(pi_vals) - sum(pi_vals_waic) < 1)
 
   # Compute tis approximation
-  expect_silent(pi_vals_tis <- loo:::elpd_loo_approximation(.llfun = llfun_test,
-                                                            data = fake_data,
-                                                            draws = fake_posterior,
-                                                            loo_approximation = "tis",
-                                                            loo_approximation_draws = 100,
-                                                            cores = cores))
+  expect_silent(
+    pi_vals_tis <- loo:::elpd_loo_approximation(
+      .llfun = llfun_test,
+      data = fake_data,
+      draws = fake_posterior,
+      loo_approximation = "tis",
+      loo_approximation_draws = 100,
+      cores = cores
+    )
+  )
   expect_true(all(pi_vals > pi_vals_tis))
   expect_true(sum(pi_vals) - sum(pi_vals_tis) < 1)
 }
@@ -69,55 +114,155 @@ test_that("elpd_loo_approximation with multiple cores", {
 
 test_that("Test loo_approximation_draws", {
   set.seed(123)
-  N <- 1000; K <- 10; S <- 1000; a0 <- 3; b0 <- 2
+  N <- 1000
+  K <- 10
+  S <- 1000
+  a0 <- 3
+  b0 <- 2
   p <- 0.7
   y <- rbinom(N, size = K, prob = p)
-  a <- a0 + sum(y); b <- b0 + N * K - sum(y)
+  a <- a0 + sum(y)
+  b <- b0 + N * K - sum(y)
   fake_posterior <- draws <- as.matrix(rbeta(S, a, b))
-  fake_data <- data.frame(y,K)
+  fake_data <- data.frame(y, K)
   rm(N, K, S, a0, b0, p, y, a, b)
   llfun_test <- function(data_i, draws) {
     dbinom(data_i$y, size = data_i$K, prob = draws, log = TRUE)
   }
 
-  expect_silent(res1 <- loo:::elpd_loo_approximation(.llfun = llfun_test, data = fake_data, draws = fake_posterior, loo_approximation = "waic", loo_approximation_draws = NULL, cores = 1))
-  expect_silent(res2 <- loo:::elpd_loo_approximation(.llfun = llfun_test, data = fake_data, draws = fake_posterior, loo_approximation = "waic", loo_approximation_draws = 10, cores = 1))
-  expect_silent(res3 <- loo:::elpd_loo_approximation(.llfun = llfun_test, data = fake_data, draws = fake_posterior[1:10*100,], loo_approximation = "waic", loo_approximation_draws = NULL, cores = 1))
-  expect_silent(res4 <- loo:::elpd_loo_approximation(.llfun = llfun_test, data = fake_data, draws = fake_posterior[1:10*100,, drop = FALSE], loo_approximation = "waic", loo_approximation_draws = NULL, cores = 1))
+  expect_silent(
+    res1 <- loo:::elpd_loo_approximation(
+      .llfun = llfun_test,
+      data = fake_data,
+      draws = fake_posterior,
+      loo_approximation = "waic",
+      loo_approximation_draws = NULL,
+      cores = 1
+    )
+  )
+  expect_silent(
+    res2 <- loo:::elpd_loo_approximation(
+      .llfun = llfun_test,
+      data = fake_data,
+      draws = fake_posterior,
+      loo_approximation = "waic",
+      loo_approximation_draws = 10,
+      cores = 1
+    )
+  )
+  expect_silent(
+    res3 <- loo:::elpd_loo_approximation(
+      .llfun = llfun_test,
+      data = fake_data,
+      draws = fake_posterior[1:10 * 100, ],
+      loo_approximation = "waic",
+      loo_approximation_draws = NULL,
+      cores = 1
+    )
+  )
+  expect_silent(
+    res4 <- loo:::elpd_loo_approximation(
+      .llfun = llfun_test,
+      data = fake_data,
+      draws = fake_posterior[1:10 * 100, , drop = FALSE],
+      loo_approximation = "waic",
+      loo_approximation_draws = NULL,
+      cores = 1
+    )
+  )
   expect_failure(expect_equal(res1, res3))
   expect_equal(res2, res3)
 
-  expect_silent(loo_ss1 <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 100, loo_approximation = "plpd", r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss2 <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 100, loo_approximation = "plpd", loo_approximation_draws = 10, r_eff = rep(1, nrow(fake_data))))
-  expect_silent(loo_ss3 <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 100, loo_approximation = "plpd", loo_approximation_draws = 31, r_eff = rep(1, nrow(fake_data))))
-  expect_error(loo_ss4 <- loo_subsample(x = llfun_test, draws = fake_posterior, data = fake_data, observations = 100, loo_approximation = "plpd", loo_approximation_draws = 3100, r_eff = rep(1, nrow(fake_data))))
-
-  expect_equal(names(loo_ss1$loo_subsampling), c("elpd_loo_approx", "loo_approximation", "loo_approximation_draws", "estimator", ".llfun", ".llgrad", ".llhess", "data_dim", "ndraws"))
+  expect_silent(
+    loo_ss1 <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 100,
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_silent(
+    loo_ss2 <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 100,
+      loo_approximation = "plpd",
+      loo_approximation_draws = 10,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_silent(
+    loo_ss3 <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 100,
+      loo_approximation = "plpd",
+      loo_approximation_draws = 31,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+  expect_error(
+    loo_ss4 <- loo_subsample(
+      x = llfun_test,
+      draws = fake_posterior,
+      data = fake_data,
+      observations = 100,
+      loo_approximation = "plpd",
+      loo_approximation_draws = 3100,
+      r_eff = rep(1, nrow(fake_data))
+    )
+  )
+
+  expect_equal(
+    names(loo_ss1$loo_subsampling),
+    c(
+      "elpd_loo_approx",
+      "loo_approximation",
+      "loo_approximation_draws",
+      "estimator",
+      ".llfun",
+      ".llgrad",
+      ".llhess",
+      "data_dim",
+      "ndraws"
+    )
+  )
   expect_null(loo_ss1$loo_subsampling$loo_approximation_draws)
   expect_equal(loo_ss2$loo_subsampling$loo_approximation_draws, 10L)
   expect_equal(loo_ss3$loo_subsampling$loo_approximation_draws, 31L)
 })
 
 
-
 test_that("waic using delta method and gradient", {
-
-
-  if (FALSE){
+  if (FALSE) {
     # Code to generate testdata - saved and loaded to avoid dependency of mvtnorm
     set.seed(123)
-    N <- 400; beta <- c(1,2); X_full <- matrix(rep(1,N), ncol = 1); X_full <- cbind(X_full, runif(N)); S <- 1000
-    y_full <- rnorm(n = N, mean = X_full%*%beta, sd = 1)
-    X <- X_full; y <- y_full
-    Lambda_0 <- diag(length(beta)); mu_0 <- c(0,0)
-    b_hat <- solve(t(X)%*%X)%*%t(X)%*%y
-    mu_n <- solve(t(X)%*%X)%*%(t(X)%*%X%*%b_hat + Lambda_0%*%mu_0)
-    Lambda_n <- t(X)%*%X + Lambda_0
+    N <- 400
+    beta <- c(1, 2)
+    X_full <- matrix(rep(1, N), ncol = 1)
+    X_full <- cbind(X_full, runif(N))
+    S <- 1000
+    y_full <- rnorm(n = N, mean = X_full %*% beta, sd = 1)
+    X <- X_full
+    y <- y_full
+    Lambda_0 <- diag(length(beta))
+    mu_0 <- c(0, 0)
+    b_hat <- solve(t(X) %*% X) %*% t(X) %*% y
+    mu_n <- solve(t(X) %*% X) %*% (t(X) %*% X %*% b_hat + Lambda_0 %*% mu_0)
+    Lambda_n <- t(X) %*% X + Lambda_0
     # Uncomment row below when running. Commented out to remove CHECK warnings
     # fake_posterior <- mvtnorm::rmvnorm(n = S, mean = mu_n, sigma = solve(Lambda_n))
     colnames(fake_posterior) <- c("a", "b")
     fake_data <- data.frame(y, X)
-    save(fake_posterior, fake_data, file = test_path("data-for-tests/normal_reg_waic_test_example.rda"))
+    save(
+      fake_posterior,
+      fake_data,
+      file = test_path("data-for-tests/normal_reg_waic_test_example.rda")
+    )
   } else {
     load(file = test_path("data-for-tests/normal_reg_waic_test_example.rda"))
   }
@@ -125,53 +270,148 @@ test_that("waic using delta method and gradient", {
   .llfun <- function(data_i, draws) {
     # data_i: ith row of fdata (fake_data[i,, drop=FALSE])
     # draws: entire fake_posterior matrix
-    dnorm(data_i$y, mean = draws[, c("a", "b")] %*% t(as.matrix(data_i[, c("X1", "X2")])), sd = 1, log = TRUE)
+    dnorm(
+      data_i$y,
+      mean = draws[, c("a", "b")] %*% t(as.matrix(data_i[, c("X1", "X2")])),
+      sd = 1,
+      log = TRUE
+    )
   }
 
   .llgrad <- function(data_i, draws) {
     x_i <- data_i[, "X2"]
-    gr <- cbind(data_i$y - draws[,"a"] - draws[,"b"]*x_i,
-                (data_i$y - draws[,"a"] - draws[,"b"]*x_i) * x_i)
+    gr <- cbind(
+      data_i$y - draws[, "a"] - draws[, "b"] * x_i,
+      (data_i$y - draws[, "a"] - draws[, "b"] * x_i) * x_i
+    )
     colnames(gr) <- c("a", "b")
     gr
   }
 
   fake_posterior <- cbind(fake_posterior, runif(nrow(fake_posterior)))
 
-  expect_silent(approx_loo_waic <- loo:::elpd_loo_approximation(.llfun, data = fake_data, draws = fake_posterior, cores = 1, loo_approximation = "waic"))
-  expect_silent(approx_loo_waic_delta <- loo:::elpd_loo_approximation(.llfun, data = fake_data, draws = fake_posterior, cores = 1, loo_approximation = "waic_grad", .llgrad = .llgrad))
-  expect_silent(approx_loo_waic_delta_diag <- loo:::elpd_loo_approximation(.llfun, data = fake_data, draws = fake_posterior, cores = 1, loo_approximation = "waic_grad_marginal", .llgrad = .llgrad))
+  expect_silent(
+    approx_loo_waic <- loo:::elpd_loo_approximation(
+      .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      loo_approximation = "waic"
+    )
+  )
+  expect_silent(
+    approx_loo_waic_delta <- loo:::elpd_loo_approximation(
+      .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      loo_approximation = "waic_grad",
+      .llgrad = .llgrad
+    )
+  )
+  expect_silent(
+    approx_loo_waic_delta_diag <- loo:::elpd_loo_approximation(
+      .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      loo_approximation = "waic_grad_marginal",
+      .llgrad = .llgrad
+    )
+  )
 
   # Test that the approaches should not deviate too much
   diff_waic_delta <- mean(approx_loo_waic - approx_loo_waic_delta)
   diff_waic_delta_diag <- mean(approx_loo_waic - approx_loo_waic_delta_diag)
-  expect_equal(approx_loo_waic,approx_loo_waic_delta_diag, tolerance = 0.1)
-  expect_equal(approx_loo_waic,approx_loo_waic_delta, tolerance = 0.01)
+  expect_equal(approx_loo_waic, approx_loo_waic_delta_diag, tolerance = 0.1)
+  expect_equal(approx_loo_waic, approx_loo_waic_delta, tolerance = 0.01)
 
   # Test usage in subsampling_loo
-  expect_silent(loo_ss_waic <- loo_subsample(x = .llfun, data = fake_data, draws = fake_posterior, cores = 1, r_eff = rep(1, nrow(fake_data)), loo_approximation = "waic", observations = 50, llgrad = .llgrad))
-  expect_silent(loo_ss_waic_delta <- loo_subsample(x = .llfun, data = fake_data, draws = fake_posterior, cores = 1, r_eff = rep(1, nrow(fake_data)), loo_approximation = "waic_grad", observations = 50, llgrad = .llgrad))
-  expect_silent(loo_ss_waic_delta_marginal <- loo_subsample(x = .llfun, data = fake_data, draws = fake_posterior, cores = 1, r_eff = rep(1, nrow(fake_data)), loo_approximation = "waic_grad_marginal", observations = 50, llgrad = .llgrad))
-  expect_silent(loo_ss_plpd <- loo_subsample(x = .llfun, data = fake_data, draws = fake_posterior, cores = 1, r_eff = rep(1, nrow(fake_data)), loo_approximation = "plpd", observations = 50, llgrad = .llgrad))
-  expect_error(loo_ss_waic_delta <- loo_subsample(x = .llfun, data = fake_data, draws = fake_posterior, cores = 1, r_eff = rep(1, nrow(fake_data)), loo_approximation = "waic_grad", observations = 50))
+  expect_silent(
+    loo_ss_waic <- loo_subsample(
+      x = .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      r_eff = rep(1, nrow(fake_data)),
+      loo_approximation = "waic",
+      observations = 50,
+      llgrad = .llgrad
+    )
+  )
+  expect_silent(
+    loo_ss_waic_delta <- loo_subsample(
+      x = .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      r_eff = rep(1, nrow(fake_data)),
+      loo_approximation = "waic_grad",
+      observations = 50,
+      llgrad = .llgrad
+    )
+  )
+  expect_silent(
+    loo_ss_waic_delta_marginal <- loo_subsample(
+      x = .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      r_eff = rep(1, nrow(fake_data)),
+      loo_approximation = "waic_grad_marginal",
+      observations = 50,
+      llgrad = .llgrad
+    )
+  )
+  expect_silent(
+    loo_ss_plpd <- loo_subsample(
+      x = .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      r_eff = rep(1, nrow(fake_data)),
+      loo_approximation = "plpd",
+      observations = 50,
+      llgrad = .llgrad
+    )
+  )
+  expect_error(
+    loo_ss_waic_delta <- loo_subsample(
+      x = .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      r_eff = rep(1, nrow(fake_data)),
+      loo_approximation = "waic_grad",
+      observations = 50
+    )
+  )
 })
 
 test_that("waic using delta 2nd order method", {
-
-
-  if (FALSE){
+  if (FALSE) {
     # Code to generate testdata - saved and loaded to avoid dependency of MCMCPack
     set.seed(123)
-    N <- 100; beta <- c(1,2); X_full <- matrix(rep(1,N), ncol = 1); X_full <- cbind(X_full, runif(N)); S <- 1000
-    y_full <- rnorm(n = N, mean = X_full%*%beta, sd = 0.5)
-    X <- X_full; y <- y_full
+    N <- 100
+    beta <- c(1, 2)
+    X_full <- matrix(rep(1, N), ncol = 1)
+    X_full <- cbind(X_full, runif(N))
+    S <- 1000
+    y_full <- rnorm(n = N, mean = X_full %*% beta, sd = 0.5)
+    X <- X_full
+    y <- y_full
     # Uncomment row below when running. Commented out to remove CHECK warnings
     # fake_posterior <- MCMCpack::MCMCregress(y~x, data = data.frame(y = y,x=X[,2]), thin = 10, mcmc = 10000) # Because Im lazy
     fake_posterior <- as.matrix(fake_posterior)
-    fake_posterior[,"sigma2"] <- sqrt(fake_posterior[,"sigma2"])
+    fake_posterior[, "sigma2"] <- sqrt(fake_posterior[, "sigma2"])
     colnames(fake_posterior) <- c("a", "b", "sigma")
     fake_data <- data.frame(y, X)
-    save(fake_posterior, fake_data, file = test_path("data-for-tests/normal_reg_waic_test_example2.rda"), compression_level = 9)
+    save(
+      fake_posterior,
+      fake_data,
+      file = test_path("data-for-tests/normal_reg_waic_test_example2.rda"),
+      compression_level = 9
+    )
   } else {
     load(file = test_path("data-for-tests/normal_reg_waic_test_example2.rda"))
   }
@@ -179,40 +419,51 @@ test_that("waic using delta 2nd order method", {
   .llfun <- function(data_i, draws) {
     # data_i: ith row of fdata (data_i <- fake_data[i,, drop=FALSE])
     # draws: entire fake_posterior matrix
-    dnorm(data_i$y, mean = draws[, c("a", "b")] %*% t(as.matrix(data_i[, c("X1", "X2")])), sd = draws[, c("sigma")], log = TRUE)
+    dnorm(
+      data_i$y,
+      mean = draws[, c("a", "b")] %*% t(as.matrix(data_i[, c("X1", "X2")])),
+      sd = draws[, c("sigma")],
+      log = TRUE
+    )
   }
 
   .llgrad <- function(data_i, draws) {
-    sigma <- draws[,"sigma"]
+    sigma <- draws[, "sigma"]
     sigma2 <- sigma^2
-    b <- draws[,"b"]
-    a <- draws[,"a"]
+    b <- draws[, "b"]
+    a <- draws[, "a"]
     x_i <- unlist(data_i[, c("X1", "X2")])
-    e <- (data_i$y - draws[,"a"] * x_i[1] - draws[,"b"] * x_i[2])
+    e <- (data_i$y - draws[, "a"] * x_i[1] - draws[, "b"] * x_i[2])
 
-    gr <- cbind(e * x_i[1] / sigma2,
-                e * x_i[2] / sigma2,
-                - 1 / sigma + e^2 / (sigma2 * sigma))
+    gr <- cbind(
+      e * x_i[1] / sigma2,
+      e * x_i[2] / sigma2,
+      -1 / sigma + e^2 / (sigma2 * sigma)
+    )
     colnames(gr) <- c("a", "b", "sigma")
     gr
   }
 
   .llhess <- function(data_i, draws) {
-    hess_array <- array(0, dim = c(ncol(draws), ncol(draws), nrow(draws)), dimnames = list(colnames(draws),colnames(draws),NULL))
-    sigma <- draws[,"sigma"]
+    hess_array <- array(
+      0,
+      dim = c(ncol(draws), ncol(draws), nrow(draws)),
+      dimnames = list(colnames(draws), colnames(draws), NULL)
+    )
+    sigma <- draws[, "sigma"]
     sigma2 <- sigma^2
-    sigma3 <- sigma2*sigma
-    b <- draws[,"b"]
-    a <- draws[,"a"]
+    sigma3 <- sigma2 * sigma
+    b <- draws[, "b"]
+    a <- draws[, "a"]
     x_i <- unlist(data_i[, c("X1", "X2")])
-    e <- (data_i$y - draws[,"a"] * x_i[1] - draws[,"b"] * x_i[2])
-
-    hess_array[1,1,] <- - x_i[1]^2 / sigma2
-    hess_array[1,2,] <- hess_array[2,1,] <- - x_i[1] * x_i[2] / sigma2
-    hess_array[2,2,] <- - x_i[2]^2 / sigma2
-    hess_array[3,1,] <- hess_array[1,3,] <- -2 * x_i[1] * e / sigma3
-    hess_array[3,2,] <- hess_array[2,3,] <- -2 * x_i[2] * e / sigma3
-    hess_array[3,3,] <- 1 / sigma2 - 3 * e^2 / (sigma2^2)
+    e <- (data_i$y - draws[, "a"] * x_i[1] - draws[, "b"] * x_i[2])
+
+    hess_array[1, 1, ] <- -x_i[1]^2 / sigma2
+    hess_array[1, 2, ] <- hess_array[2, 1, ] <- -x_i[1] * x_i[2] / sigma2
+    hess_array[2, 2, ] <- -x_i[2]^2 / sigma2
+    hess_array[3, 1, ] <- hess_array[1, 3, ] <- -2 * x_i[1] * e / sigma3
+    hess_array[3, 2, ] <- hess_array[2, 3, ] <- -2 * x_i[2] * e / sigma3
+    hess_array[3, 3, ] <- 1 / sigma2 - 3 * e^2 / (sigma2^2)
     hess_array
   }
 
@@ -220,52 +471,133 @@ test_that("waic using delta 2nd order method", {
   fake_posterior <- cbind(fake_posterior, runif(nrow(fake_posterior)))
   #draws <- fake_posterior <- cbind(fake_posterior, runif(nrow(fake_posterior)))
 
-  expect_silent(approx_loo_waic <- loo:::elpd_loo_approximation(.llfun, data = fake_data, draws = fake_posterior, cores = 1, loo_approximation = "waic"))
-  expect_silent(approx_loo_waic_delta <- loo:::elpd_loo_approximation(.llfun, data = fake_data, draws = fake_posterior, cores = 1, loo_approximation = "waic_grad", .llgrad = .llgrad))
-  expect_silent(approx_loo_waic_delta2 <- loo:::elpd_loo_approximation(.llfun, data = fake_data, draws = fake_posterior, cores = 1, loo_approximation = "waic_hess", .llgrad = .llgrad, .llhess = .llhess))
+  expect_silent(
+    approx_loo_waic <- loo:::elpd_loo_approximation(
+      .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      loo_approximation = "waic"
+    )
+  )
+  expect_silent(
+    approx_loo_waic_delta <- loo:::elpd_loo_approximation(
+      .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      loo_approximation = "waic_grad",
+      .llgrad = .llgrad
+    )
+  )
+  expect_silent(
+    approx_loo_waic_delta2 <- loo:::elpd_loo_approximation(
+      .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      loo_approximation = "waic_hess",
+      .llgrad = .llgrad,
+      .llhess = .llhess
+    )
+  )
 
   # Test that the approaches should not deviate too much
-  expect_equal(approx_loo_waic,approx_loo_waic_delta2, tolerance = 0.01)
-  expect_equal(approx_loo_waic,approx_loo_waic_delta, tolerance = 0.01)
-
-  expect_silent(test_loo_ss_waic <- loo_subsample(x = .llfun, data = fake_data, draws = fake_posterior, cores = 1, r_eff = rep(1, nrow(fake_data)), loo_approximation = "waic", observations = 50, llgrad = .llgrad))
-  expect_error(test_loo_ss_delta2 <- loo_subsample(x = .llfun, data = fake_data, draws = fake_posterior, cores = 1, r_eff = rep(1, nrow(fake_data)), loo_approximation = "waic_hess", observations = 50, llgrad = .llgrad))
-  expect_silent(test_loo_ss_delta2 <- loo_subsample(x = .llfun, data = fake_data, draws = fake_posterior, cores = 1, r_eff = rep(1, nrow(fake_data)), loo_approximation = "waic_hess", observations = 50, llgrad = .llgrad, llhess = .llhess))
-  expect_silent(test_loo_ss_delta <- loo_subsample(x = .llfun, data = fake_data, draws = fake_posterior, cores = 1, r_eff = rep(1, nrow(fake_data)), loo_approximation = "waic_grad", observations = 50, llgrad = .llgrad))
-  expect_silent(test_loo_ss_point <- loo_subsample(x = .llfun, data = fake_data, draws = fake_posterior, cores = 1, r_eff = rep(1, nrow(fake_data)), loo_approximation = "plpd", observations = 50, llgrad = .llgrad))
+  expect_equal(approx_loo_waic, approx_loo_waic_delta2, tolerance = 0.01)
+  expect_equal(approx_loo_waic, approx_loo_waic_delta, tolerance = 0.01)
+
+  expect_silent(
+    test_loo_ss_waic <- loo_subsample(
+      x = .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      r_eff = rep(1, nrow(fake_data)),
+      loo_approximation = "waic",
+      observations = 50,
+      llgrad = .llgrad
+    )
+  )
+  expect_error(
+    test_loo_ss_delta2 <- loo_subsample(
+      x = .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      r_eff = rep(1, nrow(fake_data)),
+      loo_approximation = "waic_hess",
+      observations = 50,
+      llgrad = .llgrad
+    )
+  )
+  expect_silent(
+    test_loo_ss_delta2 <- loo_subsample(
+      x = .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      r_eff = rep(1, nrow(fake_data)),
+      loo_approximation = "waic_hess",
+      observations = 50,
+      llgrad = .llgrad,
+      llhess = .llhess
+    )
+  )
+  expect_silent(
+    test_loo_ss_delta <- loo_subsample(
+      x = .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      r_eff = rep(1, nrow(fake_data)),
+      loo_approximation = "waic_grad",
+      observations = 50,
+      llgrad = .llgrad
+    )
+  )
+  expect_silent(
+    test_loo_ss_point <- loo_subsample(
+      x = .llfun,
+      data = fake_data,
+      draws = fake_posterior,
+      cores = 1,
+      r_eff = rep(1, nrow(fake_data)),
+      loo_approximation = "plpd",
+      observations = 50,
+      llgrad = .llgrad
+    )
+  )
 })
 
 
 test_that("whhest works as expected", {
-
-
   N <- 100
   m <- 10
-  z <- rep(1/N, m)
+  z <- rep(1 / N, m)
   y <- 1:10
-  m_i <- rep(1,m)
+  m_i <- rep(1, m)
   expect_silent(whe <- loo:::whhest(z = z, m_i = m_i, y = y, N = N))
   expect_equal(whe$y_hat_ppz, 550)
-  man_var <- (sum((whe$y_hat_ppz - y/z)^2)/(m-1))/m
+  man_var <- (sum((whe$y_hat_ppz - y / z)^2) / (m - 1)) / m
   expect_equal(whe$v_hat_y_ppz, man_var)
-  z <- 1:10/(sum(1:10)*10)
+  z <- 1:10 / (sum(1:10) * 10)
   expect_silent(whe <- loo:::whhest(z = z, m_i = m_i, y = y, N = N))
   expect_equal(whe$y_hat_ppz, 550)
   expect_equal(whe$v_hat_y_ppz, 0)
 
   # School book example
   # https://newonlinecourses.science.psu.edu/stat506/node/15/
-  z <- c(650/15650, 2840/15650, 3200/15650)
+  z <- c(650 / 15650, 2840 / 15650, 3200 / 15650)
   y <- c(420, 1785, 2198)
-  m_i <- c(1,1,1)
+  m_i <- c(1, 1, 1)
   N <- 10
   expect_silent(whe <- loo:::whhest(z = z, m_i = m_i, y = y, N = N))
   expect_equal(round(whe$y_hat_ppz, 2), 10232.75, tolerance = 0)
   expect_equal(whe$v_hat_y_ppz, 73125.74, tolerance = 0.01)
   # Double check that it is rounding error
-  man_var_round <- (sum((round(y/z,2) - 10232.75)^2)) * (1/2) * (1/3)
+  man_var_round <- (sum((round(y / z, 2) - 10232.75)^2)) * (1 / 2) * (1 / 3)
   expect_equal(man_var_round, 73125.74, tolerance = 0.001)
-  man_var_exact <- (sum((y/z - 10232.75)^2)) * (1/2) * (1/3)
+  man_var_exact <- (sum((y / z - 10232.75)^2)) * (1 / 2) * (1 / 3)
   expect_equal(whe$v_hat_y_ppz, man_var_exact, tolerance = 0.001)
 
   # Add test for variance estimation
@@ -273,7 +605,7 @@ test_that("whhest works as expected", {
   m <- 10
   y <- rep(1:10, 1)
   true_var <- var(rep(y, 10)) * (99)
-  z <- rep(1/N, m)
+  z <- rep(1 / N, m)
   m_i <- rep(100000, m)
   expect_silent(whe <- loo:::whhest(z = z, m_i = m_i, y = y, N = N))
   expect_equal(true_var, whe$hat_v_y_ppz, tolerance = 0.01)
@@ -282,24 +614,21 @@ test_that("whhest works as expected", {
   N <- 100
   y <- rep(1:10, 2)
   m <- length(y)
-  z <- rep(1/N, m)
-  m_i <- rep(1,m)
+  z <- rep(1 / N, m)
+  m_i <- rep(1, m)
   expect_silent(whe1 <- loo:::whhest(z = z, m_i = m_i, y = y, N = N))
   y <- rep(1:10)
   m <- length(y)
-  z <- rep(1/N, m)
-  m_i <- rep(2,m)
+  z <- rep(1 / N, m)
+  m_i <- rep(2, m)
   expect_silent(whe2 <- loo:::whhest(z = z, m_i = m_i, y = y, N = N))
   expect_equal(whe1$y_hat_ppz, whe2$y_hat_ppz)
   expect_equal(whe1$v_hat_y_ppz, whe2$v_hat_y_ppz)
   expect_equal(whe1$hat_v_y_ppz, whe1$hat_v_y_ppz)
-
 })
 
 
 test_that("srs_diff_est works as expected", {
-
-
   set.seed(1234)
   N <- 1000
   y_true <- 1:N
@@ -307,7 +636,7 @@ test_that("srs_diff_est works as expected", {
   y_approx <- rnorm(N, y_true, 0.1)
   m <- 100
   sigma_hat <- y_hat <- se_y_hat <- numeric(10000)
-  for(i in 1:10000){
+  for (i in 1:10000) {
     y_idx <- sample(1:N, size = m)
     y <- y_true[y_idx]
     res <- loo:::srs_diff_est(y_approx, y, y_idx)
@@ -317,7 +646,8 @@ test_that("srs_diff_est works as expected", {
   }
   expect_equal(mean(y_hat), sum(y_true), tolerance = 0.1)
 
-  in_ki <- y_hat + 2 * se_y_hat > sum(y_true) & y_hat - 2*se_y_hat < sum(y_true)
+  in_ki <- y_hat + 2 * se_y_hat > sum(y_true) &
+    y_hat - 2 * se_y_hat < sum(y_true)
   expect_equal(mean(in_ki), 0.95, tolerance = 0.01)
 
   # Should be  unbiased
@@ -330,22 +660,45 @@ test_that("srs_diff_est works as expected", {
   expect_equal(res$y_hat, 500500, tolerance = 0.0001)
   expect_equal(res$v_y_hat, 0, tolerance = 0.0001)
   expect_equal(sqrt(res$hat_v_y), sigma_hat_true, tolerance = 0.1)
-
 })
 
 test_that("srs_est works as expected", {
-
-
   set.seed(1234)
   # Cochran 1976 example Table 2.2
 
-  y <- c(rep(42,23),rep(41,4), 36, 32, 29, 27, 27, 23, 19, 16, 16, 15, 15, 14, 11, 10, 9, 7, 6, 6, 6, 5, 5, 4, 3)
+  y <- c(
+    rep(42, 23),
+    rep(41, 4),
+    36,
+    32,
+    29,
+    27,
+    27,
+    23,
+    19,
+    16,
+    16,
+    15,
+    15,
+    14,
+    11,
+    10,
+    9,
+    7,
+    6,
+    6,
+    6,
+    5,
+    5,
+    4,
+    3
+  )
   expect_equal(sum(y), 1471)
   approx_loo <- rep(0L, 676)
   expect_equal(sum(y^2), 54497)
   res <- loo:::srs_est(y = y, approx_loo)
   expect_equal(res$y_hat, 19888, tolerance = 0.0001)
-  expect_equal(res$v_y_hat, 676^2*229*(1-0.074)/50, tolerance = 0.0001)
+  expect_equal(res$v_y_hat, 676^2 * 229 * (1 - 0.074) / 50, tolerance = 0.0001)
   expect_equal(res$hat_v_y, 676 * var(y), tolerance = 0.0001)
 
   # Simulation example
@@ -356,7 +709,7 @@ test_that("srs_est works as expected", {
 
   m <- 100
   y_hat <- se_y_hat <- sigma_hat <- numeric(10000)
-  for(i in 1:10000){
+  for (i in 1:10000) {
     y_idx <- sample(1:N, size = m)
     y <- y_true[y_idx]
     res <- loo:::srs_est(y = y, y_approx = y_true)
@@ -366,7 +719,8 @@ test_that("srs_est works as expected", {
   }
   expect_equal(mean(y_hat), sum(y_true), tolerance = 0.1)
 
-  in_ki <- y_hat + 2*se_y_hat > sum(y_true) & y_hat - 2*se_y_hat < sum(y_true)
+  in_ki <- y_hat + 2 * se_y_hat > sum(y_true) &
+    y_hat - 2 * se_y_hat < sum(y_true)
   expect_equal(mean(in_ki), 0.95, tolerance = 0.01)
 
   # Should be  unbiased
@@ -378,5 +732,4 @@ test_that("srs_est works as expected", {
   res <- loo:::srs_est(y, y_true)
   expect_equal(res$y_hat, 500500, tolerance = 0.0001)
   expect_equal(res$v_y_hat, 0, tolerance = 0.0001)
-
 })
diff --git a/tests/testthat/test_loo_subsampling_cases.R b/tests/testthat/test_loo_subsampling_cases.R
index 26f517d2..1f53192b 100644
--- a/tests/testthat/test_loo_subsampling_cases.R
+++ b/tests/testthat/test_loo_subsampling_cases.R
@@ -1,4 +1,3 @@
-library(loo)
 options(mc.cores = 1)
 
 test_that("Test loo_subsampling and loo_approx with radon data", {
@@ -11,62 +10,187 @@ test_that("Test loo_subsampling and loo_approx with radon data", {
   log_g_test <- log_q
   draws_test <- draws
   data_test <- data
-  rm(llfun, log_p,log_q, draws, data)
+  rm(llfun, log_p, log_q, draws, data)
 
   set.seed(134)
-  expect_silent(full_loo <- loo(llfun_test, draws = draws_test, data = data_test, r_eff = rep(1, nrow(data_test))))
+  expect_silent(
+    full_loo <- loo(
+      llfun_test,
+      draws = draws_test,
+      data = data_test,
+      r_eff = rep(1, nrow(data_test))
+    )
+  )
   expect_s3_class(full_loo, "psis_loo")
 
   set.seed(134)
-  expect_silent(loo_ss <- loo_subsample(x = llfun_test, draws = draws_test, data = data_test, observations = 200, loo_approximation = "plpd", r_eff = rep(1, nrow(data_test))))
+  expect_silent(
+    loo_ss <- loo_subsample(
+      x = llfun_test,
+      draws = draws_test,
+      data = data_test,
+      observations = 200,
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(data_test))
+    )
+  )
   expect_s3_class(loo_ss, "psis_loo_ss")
 
   set.seed(134)
-  expect_silent(loo_ap_ss <- loo_subsample(x = llfun_test, draws = draws_test, data = data_test, log_p = log_p_test, log_g = log_g_test, observations = 200, loo_approximation = "plpd", r_eff = rep(1, nrow(data_test))))
+  expect_silent(
+    loo_ap_ss <- loo_subsample(
+      x = llfun_test,
+      draws = draws_test,
+      data = data_test,
+      log_p = log_p_test,
+      log_g = log_g_test,
+      observations = 200,
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(data_test))
+    )
+  )
   expect_s3_class(loo_ap_ss, "psis_loo_ss")
   expect_s3_class(loo_ap_ss, "psis_loo_ap")
 
-  expect_silent(loo_ap_ss_full <- loo_subsample(x = llfun_test, log_p = log_p_test, log_g = log_g_test, draws = draws_test, data = data_test, observations = NULL, loo_approximation = "plpd", r_eff = rep(1, nrow(data_test))))
+  expect_silent(
+    loo_ap_ss_full <- loo_subsample(
+      x = llfun_test,
+      log_p = log_p_test,
+      log_g = log_g_test,
+      draws = draws_test,
+      data = data_test,
+      observations = NULL,
+      loo_approximation = "plpd",
+      r_eff = rep(1, nrow(data_test))
+    )
+  )
   expect_failure(expect_s3_class(loo_ap_ss_full, "psis_loo_ss"))
   expect_s3_class(loo_ap_ss_full, "psis_loo_ap")
 
   # Expect similar results
   z <- 2
-  expect_lte(loo_ss$estimates["elpd_loo", "Estimate"] - z * loo_ss$estimates["elpd_loo", "subsampling SE"], full_loo$estimates["elpd_loo", "Estimate"])
-  expect_gte(loo_ss$estimates["elpd_loo", "Estimate"] + z * loo_ss$estimates["elpd_loo", "subsampling SE"], full_loo$estimates["elpd_loo", "Estimate"])
-  expect_lte(loo_ss$estimates["p_loo", "Estimate"] - z * loo_ss$estimates["p_loo", "subsampling SE"], full_loo$estimates["p_loo", "Estimate"])
-  expect_gte(loo_ss$estimates["p_loo", "Estimate"] + z * loo_ss$estimates["p_loo", "subsampling SE"], full_loo$estimates["p_loo", "Estimate"])
-  expect_lte(loo_ss$estimates["looic", "Estimate"] - z * loo_ss$estimates["looic", "subsampling SE"], full_loo$estimates["looic", "Estimate"])
-  expect_gte(loo_ss$estimates["looic", "Estimate"] + z * loo_ss$estimates["looic", "subsampling SE"], full_loo$estimates["looic", "Estimate"])
-
-  expect_failure(expect_equal(full_loo$estimates["elpd_loo", "Estimate"], loo_ss$estimates["elpd_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(full_loo$estimates["p_loo", "Estimate"], loo_ss$estimates["p_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(full_loo$estimates["looic", "Estimate"], loo_ss$estimates["looic", "Estimate"], tolerance = 0.00000001))
+  expect_lte(
+    loo_ss$estimates["elpd_loo", "Estimate"] -
+      z * loo_ss$estimates["elpd_loo", "subsampling SE"],
+    full_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["elpd_loo", "Estimate"] +
+      z * loo_ss$estimates["elpd_loo", "subsampling SE"],
+    full_loo$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss$estimates["p_loo", "Estimate"] -
+      z * loo_ss$estimates["p_loo", "subsampling SE"],
+    full_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["p_loo", "Estimate"] +
+      z * loo_ss$estimates["p_loo", "subsampling SE"],
+    full_loo$estimates["p_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ss$estimates["looic", "Estimate"] -
+      z * loo_ss$estimates["looic", "subsampling SE"],
+    full_loo$estimates["looic", "Estimate"]
+  )
+  expect_gte(
+    loo_ss$estimates["looic", "Estimate"] +
+      z * loo_ss$estimates["looic", "subsampling SE"],
+    full_loo$estimates["looic", "Estimate"]
+  )
+
+  expect_failure(expect_equal(
+    full_loo$estimates["elpd_loo", "Estimate"],
+    loo_ss$estimates["elpd_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    full_loo$estimates["p_loo", "Estimate"],
+    loo_ss$estimates["p_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    full_loo$estimates["looic", "Estimate"],
+    loo_ss$estimates["looic", "Estimate"],
+    tolerance = 0.00000001
+  ))
 
   z <- 2
-  expect_lte(loo_ap_ss$estimates["elpd_loo", "Estimate"] - z * loo_ap_ss$estimates["elpd_loo", "subsampling SE"], loo_ap_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_gte(loo_ap_ss$estimates["elpd_loo", "Estimate"] + z * loo_ap_ss$estimates["elpd_loo", "subsampling SE"], loo_ap_ss_full$estimates["elpd_loo", "Estimate"])
-  expect_lte(loo_ap_ss$estimates["p_loo", "Estimate"] - z * loo_ap_ss$estimates["p_loo", "subsampling SE"], loo_ap_ss_full$estimates["p_loo", "Estimate"])
-  expect_gte(loo_ap_ss$estimates["p_loo", "Estimate"] + z * loo_ap_ss$estimates["p_loo", "subsampling SE"], loo_ap_ss_full$estimates["p_loo", "Estimate"])
-  expect_lte(loo_ap_ss$estimates["looic", "Estimate"] - z * loo_ap_ss$estimates["looic", "subsampling SE"], loo_ap_ss_full$estimates["looic", "Estimate"])
-  expect_gte(loo_ap_ss$estimates["looic", "Estimate"] + z * loo_ap_ss$estimates["looic", "subsampling SE"], loo_ap_ss_full$estimates["looic", "Estimate"])
-
-  expect_failure(expect_equal(loo_ap_ss_full$estimates["elpd_loo", "Estimate"], loo_ap_ss$estimates["elpd_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(loo_ap_ss_full$estimates["p_loo", "Estimate"], loo_ap_ss$estimates["p_loo", "Estimate"], tolerance = 0.00000001))
-  expect_failure(expect_equal(loo_ap_ss_full$estimates["looic", "Estimate"], loo_ap_ss$estimates["looic", "Estimate"], tolerance = 0.00000001))
+  expect_lte(
+    loo_ap_ss$estimates["elpd_loo", "Estimate"] -
+      z * loo_ap_ss$estimates["elpd_loo", "subsampling SE"],
+    loo_ap_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ap_ss$estimates["elpd_loo", "Estimate"] +
+      z * loo_ap_ss$estimates["elpd_loo", "subsampling SE"],
+    loo_ap_ss_full$estimates["elpd_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ap_ss$estimates["p_loo", "Estimate"] -
+      z * loo_ap_ss$estimates["p_loo", "subsampling SE"],
+    loo_ap_ss_full$estimates["p_loo", "Estimate"]
+  )
+  expect_gte(
+    loo_ap_ss$estimates["p_loo", "Estimate"] +
+      z * loo_ap_ss$estimates["p_loo", "subsampling SE"],
+    loo_ap_ss_full$estimates["p_loo", "Estimate"]
+  )
+  expect_lte(
+    loo_ap_ss$estimates["looic", "Estimate"] -
+      z * loo_ap_ss$estimates["looic", "subsampling SE"],
+    loo_ap_ss_full$estimates["looic", "Estimate"]
+  )
+  expect_gte(
+    loo_ap_ss$estimates["looic", "Estimate"] +
+      z * loo_ap_ss$estimates["looic", "subsampling SE"],
+    loo_ap_ss_full$estimates["looic", "Estimate"]
+  )
+
+  expect_failure(expect_equal(
+    loo_ap_ss_full$estimates["elpd_loo", "Estimate"],
+    loo_ap_ss$estimates["elpd_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    loo_ap_ss_full$estimates["p_loo", "Estimate"],
+    loo_ap_ss$estimates["p_loo", "Estimate"],
+    tolerance = 0.00000001
+  ))
+  expect_failure(expect_equal(
+    loo_ap_ss_full$estimates["looic", "Estimate"],
+    loo_ap_ss$estimates["looic", "Estimate"],
+    tolerance = 0.00000001
+  ))
 
   # Correct printout
-  expect_failure(expect_output(print(full_loo), "Posterior approximation correction used\\."))
-  expect_failure(expect_output(print(full_loo), "subsampled log-likelihood\nvalues"))
-
-  expect_failure(expect_output(print(loo_ss), "Posterior approximation correction used\\."))
+  expect_failure(expect_output(
+    print(full_loo),
+    "Posterior approximation correction used\\."
+  ))
+  expect_failure(expect_output(
+    print(full_loo),
+    "subsampled log-likelihood\nvalues"
+  ))
+
+  expect_failure(expect_output(
+    print(loo_ss),
+    "Posterior approximation correction used\\."
+  ))
   expect_output(print(loo_ss), "subsampled log-likelihood\nvalues")
 
   expect_output(print(loo_ap_ss), "Posterior approximation correction used\\.")
   expect_output(print(loo_ap_ss), "subsampled log-likelihood\nvalues")
 
-  expect_output(print(loo_ap_ss_full), "Posterior approximation correction used\\.")
-  expect_failure(expect_output(print(loo_ap_ss_full), "subsampled log-likelihood\nvalues"))
+  expect_output(
+    print(loo_ap_ss_full),
+    "Posterior approximation correction used\\."
+  )
+  expect_failure(expect_output(
+    print(loo_ap_ss_full),
+    "subsampled log-likelihood\nvalues"
+  ))
 
   # Test conversion of objects
   expect_silent(loo_ap_full <- loo:::as.psis_loo.psis_loo(loo_ap_ss_full))
@@ -77,43 +201,60 @@ test_that("Test loo_subsampling and loo_approx with radon data", {
   expect_silent(loo_ap_full2 <- loo:::as.psis_loo.psis_loo_ss(loo_ap_full_ss))
   expect_s3_class(loo_ap_full2, "psis_loo_ap")
   expect_failure(expect_s3_class(loo_ap_full2, "psis_loo_ss"))
-  expect_equal(loo_ap_full2,loo_ap_full)
+  expect_equal(loo_ap_full2, loo_ap_full)
 
   # Test update
   set.seed(4712)
-  expect_silent(loo_ss2 <- update(loo_ss, draws = draws_test, data = data_test, observations = 1000, r_eff = rep(1, nrow(data_test))))
+  expect_silent(
+    loo_ss2 <- update(
+      loo_ss,
+      draws = draws_test,
+      data = data_test,
+      observations = 1000,
+      r_eff = rep(1, nrow(data_test))
+    )
+  )
   expect_gt(dim(loo_ss2)[2], dim(loo_ss)[2])
   expect_gt(dim(loo_ss2$pointwise)[1], dim(loo_ss$pointwise)[1])
   expect_equal(nobs(loo_ss), 200)
   expect_equal(nobs(loo_ss2), 1000)
-  for(i in 1:nrow(loo_ss2$estimates)) {
-    expect_lt(loo_ss2$estimates[i, "subsampling SE"],
-              loo_ss$estimates[i, "subsampling SE"])
+  for (i in 1:nrow(loo_ss2$estimates)) {
+    expect_lt(
+      loo_ss2$estimates[i, "subsampling SE"],
+      loo_ss$estimates[i, "subsampling SE"]
+    )
   }
 
   set.seed(4712)
-  expect_silent(loo_ap_ss2 <- update(object = loo_ap_ss, draws = draws_test, data = data_test, observations = 2000))
+  expect_silent(
+    loo_ap_ss2 <- update(
+      object = loo_ap_ss,
+      draws = draws_test,
+      data = data_test,
+      observations = 2000
+    )
+  )
   expect_gt(dim(loo_ap_ss2)[2], dim(loo_ap_ss)[2])
   expect_gt(dim(loo_ap_ss2$pointwise)[1], dim(loo_ap_ss$pointwise)[1])
   expect_equal(nobs(loo_ap_ss), 200)
   expect_equal(nobs(loo_ap_ss2), 2000)
-  for(i in 1:nrow(loo_ap_ss2$estimates)) {
-    expect_lt(loo_ap_ss2$estimates[i, "subsampling SE"],
-              loo_ap_ss$estimates[i, "subsampling SE"])
+  for (i in 1:nrow(loo_ap_ss2$estimates)) {
+    expect_lt(
+      loo_ap_ss2$estimates[i, "subsampling SE"],
+      loo_ap_ss$estimates[i, "subsampling SE"]
+    )
   }
 
   expect_equal(round(full_loo$estimates), round(loo_ap_ss_full$estimates))
   expect_failure(expect_equal(full_loo$estimates, loo_ap_ss_full$estimates))
   expect_equal(dim(full_loo), dim(loo_ap_ss_full))
   expect_s3_class(loo_ap_ss_full, "psis_loo_ap")
-
 })
 
 
 test_that("Test the vignette", {
   skip_on_cran()
 
-
   # NOTE
   # If any of these test fails, the vignette probably needs to be updated
 
@@ -138,9 +279,12 @@ test_that("Test the vignette", {
     # logistic <- function(x) {1  / (1 + exp(-x))}
     # logit <- function(x) {log(x) - log(1-x)}
     llfun_logistic <- function(data_i, draws) {
-      x_i <- as.matrix(data_i[, which(grepl(colnames(data_i), pattern = "X")), drop=FALSE])
+      x_i <- as.matrix(data_i[,
+        which(grepl(colnames(data_i), pattern = "X")),
+        drop = FALSE
+      ])
       y_pred <- draws %*% t(x_i)
-      dbinom(x = data_i$y, size = 1, prob = 1  / (1 + exp(-y_pred)), log = TRUE)
+      dbinom(x = data_i$y, size = 1, prob = 1 / (1 + exp(-y_pred)), log = TRUE)
     }
 
     # Prepare data
@@ -173,96 +317,231 @@ test_that("Test the vignette", {
     fit_2 <- stan(fit = fit_1, data = standata, seed = 4711)
     parameter_draws_2 <- extract(fit_2)$beta
 
-    save(llfun_logistic,
-         stan_df, stan_df2,
-         parameter_draws, parameter_draws_laplace, parameter_draws_2,
-         log_p, log_g,
-         file = test_path("data-for-tests/loo_subsample_vignette.rda"), compression_level = 9)
-
+    save(
+      llfun_logistic,
+      stan_df,
+      stan_df2,
+      parameter_draws,
+      parameter_draws_laplace,
+      parameter_draws_2,
+      log_p,
+      log_g,
+      file = test_path("data-for-tests/loo_subsample_vignette.rda"),
+      compression_level = 9
+    )
   } else {
     load(test_path("data-for-tests/loo_subsample_vignette.rda"))
   }
 
   set.seed(4711)
-  expect_no_warning(looss_1 <- loo_subsample(llfun_logistic, draws = parameter_draws, data = stan_df, observations = 100))
-  expect_output(print(looss_1), "Computed from 4000 by 100 subsampled log-likelihood")
+  expect_no_warning(
+    looss_1 <- loo_subsample(
+      llfun_logistic,
+      draws = parameter_draws,
+      data = stan_df,
+      observations = 100
+    )
+  )
+  expect_output(
+    print(looss_1),
+    "Computed from 4000 by 100 subsampled log-likelihood"
+  )
   expect_output(print(looss_1), "values from 3020 total observations.")
-  expect_output(print(looss_1), "MCSE and ESS estimates assume independent draws")
+  expect_output(
+    print(looss_1),
+    "MCSE and ESS estimates assume independent draws"
+  )
   expect_output(print(looss_1), "elpd_loo  -1968.5 15.6            0.3")
   expect_output(print(looss_1), "p_loo         3.1  0.1            0.4")
   expect_s3_class(looss_1, c("psis_loo_ss", "psis_loo", "loo"))
 
   set.seed(4711)
-  expect_no_warning(looss_1b <- update(looss_1, draws = parameter_draws, data = stan_df, observations = 200))
-  expect_output(print(looss_1b), "Computed from 4000 by 200 subsampled log-likelihood")
+  expect_no_warning(
+    looss_1b <- update(
+      looss_1,
+      draws = parameter_draws,
+      data = stan_df,
+      observations = 200
+    )
+  )
+  expect_output(
+    print(looss_1b),
+    "Computed from 4000 by 200 subsampled log-likelihood"
+  )
   expect_output(print(looss_1b), "values from 3020 total observations.")
-  expect_output(print(looss_1b), "MCSE and ESS estimates assume independent draws")
+  expect_output(
+    print(looss_1b),
+    "MCSE and ESS estimates assume independent draws"
+  )
   expect_output(print(looss_1b), "elpd_loo  -1968.3 15.6            0.2")
   expect_output(print(looss_1b), "p_loo         3.2  0.1            0.4")
   expect_s3_class(looss_1b, c("psis_loo_ss", "psis_loo", "loo"))
 
   set.seed(4711)
-  expect_no_warning(looss_2 <- loo_subsample(x = llfun_logistic, draws = parameter_draws, data = stan_df, observations = 100, estimator = "hh_pps", loo_approximation = "lpd", loo_approximation_draws = 100))
-  expect_output(print(looss_2), "Computed from 4000 by 100 subsampled log-likelihood")
+  expect_no_warning(
+    looss_2 <- loo_subsample(
+      x = llfun_logistic,
+      draws = parameter_draws,
+      data = stan_df,
+      observations = 100,
+      estimator = "hh_pps",
+      loo_approximation = "lpd",
+      loo_approximation_draws = 100
+    )
+  )
+  expect_output(
+    print(looss_2),
+    "Computed from 4000 by 100 subsampled log-likelihood"
+  )
   expect_output(print(looss_2), "values from 3020 total observations.")
-  expect_output(print(looss_2), "MCSE and ESS estimates assume independent draws")
+  expect_output(
+    print(looss_2),
+    "MCSE and ESS estimates assume independent draws"
+  )
   # Currently failing
   # expect_output(print(looss_2), "elpd_loo  -1968.9 15.4            0.5")
   # expect_output(print(looss_2), "p_loo         3.5  0.2            0.5")
   expect_s3_class(looss_2, c("psis_loo_ss", "psis_loo", "loo"))
 
   set.seed(4711)
-  expect_no_warning(aploo_1 <- loo_approximate_posterior(llfun_logistic, draws = parameter_draws_laplace, data = stan_df, log_p = log_p, log_g = log_g))
-  expect_output(print(aploo_1), "Computed from 2000 by 3020 log-likelihood matrix")
-  expect_output(print(aploo_1), "MCSE and ESS estimates assume independent draws")
+  expect_no_warning(
+    aploo_1 <- loo_approximate_posterior(
+      llfun_logistic,
+      draws = parameter_draws_laplace,
+      data = stan_df,
+      log_p = log_p,
+      log_g = log_g
+    )
+  )
+  expect_output(
+    print(aploo_1),
+    "Computed from 2000 by 3020 log-likelihood matrix"
+  )
+  expect_output(
+    print(aploo_1),
+    "MCSE and ESS estimates assume independent draws"
+  )
   expect_output(print(aploo_1), "elpd_loo  -1968.4 15.6")
   expect_output(print(aploo_1), "p_loo         3.2  0.2")
   expect_output(print(aploo_1), "Posterior approximation correction used.")
   expect_output(print(aploo_1), "All Pareto k estimates are good")
-  expect_equal(length(pareto_k_ids(aploo_1,threshold=0.5)), 31)
+  expect_equal(length(pareto_k_ids(aploo_1, threshold = 0.5)), 31)
   expect_s3_class(aploo_1, c("psis_loo_ap", "psis_loo", "loo"))
 
   set.seed(4711)
-  expect_no_warning(looapss_1 <- loo_subsample(llfun_logistic, draws = parameter_draws_laplace, data = stan_df, log_p = log_p, log_g = log_g, observations = 100))
-  expect_output(print(looapss_1), "Computed from 2000 by 100 subsampled log-likelihood")
-  expect_output(print(looapss_1), "MCSE and ESS estimates assume independent draws")
+  expect_no_warning(
+    looapss_1 <- loo_subsample(
+      llfun_logistic,
+      draws = parameter_draws_laplace,
+      data = stan_df,
+      log_p = log_p,
+      log_g = log_g,
+      observations = 100
+    )
+  )
+  expect_output(
+    print(looapss_1),
+    "Computed from 2000 by 100 subsampled log-likelihood"
+  )
+  expect_output(
+    print(looapss_1),
+    "MCSE and ESS estimates assume independent draws"
+  )
   expect_output(print(looapss_1), "values from 3020 total observations.")
   expect_output(print(looapss_1), "elpd_loo  -1968.2 15.6            0.4")
   expect_output(print(looapss_1), "p_loo         2.9  0.1            0.5")
   expect_output(print(looapss_1), "All Pareto k estimates are good")
-  expect_equal(length(pareto_k_ids(looapss_1,threshold=0.5)), 3)
+  expect_equal(length(pareto_k_ids(looapss_1, threshold = 0.5)), 3)
 
   # Loo compare
   set.seed(4711)
-  expect_no_warning(looss_1 <- loo_subsample(llfun_logistic, draws = parameter_draws, data = stan_df, observations = 100))
+  expect_no_warning(
+    looss_1 <- loo_subsample(
+      llfun_logistic,
+      draws = parameter_draws,
+      data = stan_df,
+      observations = 100
+    )
+  )
   set.seed(4712)
-  expect_no_warning(looss_2 <- loo_subsample(x = llfun_logistic, draws = parameter_draws_2, data = stan_df2, observations = 100))
-  expect_output(print(looss_2), "Computed from 4000 by 100 subsampled log-likelihood")
-  expect_output(print(looss_2), "MCSE and ESS estimates assume independent draws")
+  expect_no_warning(
+    looss_2 <- loo_subsample(
+      x = llfun_logistic,
+      draws = parameter_draws_2,
+      data = stan_df2,
+      observations = 100
+    )
+  )
+  expect_output(
+    print(looss_2),
+    "Computed from 4000 by 100 subsampled log-likelihood"
+  )
+  expect_output(
+    print(looss_2),
+    "MCSE and ESS estimates assume independent draws"
+  )
   expect_output(print(looss_2), "values from 3020 total observations.")
   expect_output(print(looss_2), "elpd_loo  -1952.0 16.2            0.2")
   expect_output(print(looss_2), "p_loo         2.6  0.1            0.3")
 
-  expect_warning(comp <- loo_compare(looss_1, looss_2), "Different subsamples in 'model2' and 'model1'. Naive diff SE is used.")
+  expect_warning(
+    comp <- loo_compare(looss_1, looss_2),
+    "Different subsamples in 'model2' and 'model1'. Naive diff SE is used."
+  )
   expect_output(print(comp), "model1 16.5      22.5     0.4")
 
   set.seed(4712)
-  expect_no_warning(looss_2_m <- loo_subsample(x = llfun_logistic, draws = parameter_draws_2, data = stan_df2, observations = looss_1))
-  expect_message(looss_2_m <- suppressWarnings(loo_subsample(x = llfun_logistic, draws = parameter_draws_2, data = stan_df2, observations = obs_idx(looss_1))),
-                 "Simple random sampling with replacement assumed.")
+  expect_no_warning(
+    looss_2_m <- loo_subsample(
+      x = llfun_logistic,
+      draws = parameter_draws_2,
+      data = stan_df2,
+      observations = looss_1
+    )
+  )
+  expect_message(
+    looss_2_m <- suppressWarnings(loo_subsample(
+      x = llfun_logistic,
+      draws = parameter_draws_2,
+      data = stan_df2,
+      observations = obs_idx(looss_1)
+    )),
+    "Simple random sampling with replacement assumed."
+  )
 
   expect_silent(comp <- loo_compare(looss_1, looss_2_m))
   expect_output(print(comp), "model1 16.1       4.4     0.1")
 
   set.seed(4712)
-  expect_no_warning(looss_1 <- update(looss_1, draws = parameter_draws, data = stan_df, observations = 200))
-  expect_no_warning(looss_2_m <- update(looss_2_m, draws = parameter_draws_2, data = stan_df2, observations = looss_1))
+  expect_no_warning(
+    looss_1 <- update(
+      looss_1,
+      draws = parameter_draws,
+      data = stan_df,
+      observations = 200
+    )
+  )
+  expect_no_warning(
+    looss_2_m <- update(
+      looss_2_m,
+      draws = parameter_draws_2,
+      data = stan_df2,
+      observations = looss_1
+    )
+  )
   expect_silent(comp2 <- loo_compare(looss_1, looss_2_m))
   expect_output(print(comp2), "model1 16.3       4.4     0.1")
 
-  expect_no_warning(looss_2_full <- loo(x = llfun_logistic, draws = parameter_draws_2, data = stan_df2))
-  expect_message(comp3 <- loo_compare(x = list(looss_1, looss_2_full)),
-                 "Estimated elpd_diff using observations included in loo calculations for all models.")
+  expect_no_warning(
+    looss_2_full <- loo(
+      x = llfun_logistic,
+      draws = parameter_draws_2,
+      data = stan_df2
+    )
+  )
+  expect_message(
+    comp3 <- loo_compare(x = list(looss_1, looss_2_full)),
+    "Estimated elpd_diff using observations included in loo calculations for all models."
+  )
   expect_output(print(comp3), "model1 16.5       4.4     0.3")
-
 })
diff --git a/tests/testthat/test_model_weighting.R b/tests/testthat/test_model_weighting.R
index 6397d7b3..b15730ba 100644
--- a/tests/testthat/test_model_weighting.R
+++ b/tests/testthat/test_model_weighting.R
@@ -1,20 +1,18 @@
-library(loo)
-
 # generate fake data
 set.seed(123)
-y<-rnorm(50,0,1)
-sd_sim1<- abs(rnorm(500,1.5, 0.1))
-sd_sim2<- abs(rnorm(500,1.2, 0.1))
-sd_sim3<- abs(rnorm(500,1, 0.05))
+y <- rnorm(50, 0, 1)
+sd_sim1 <- abs(rnorm(500, 1.5, 0.1))
+sd_sim2 <- abs(rnorm(500, 1.2, 0.1))
+sd_sim3 <- abs(rnorm(500, 1, 0.05))
 log_lik1 <- log_lik2 <- log_lik3 <- matrix(NA, 500, 50)
-for(s in 1:500) {
-  log_lik1[s,] <- dnorm(y,-1,sd_sim1[s], log=T)
-  log_lik2[s,] <-  dnorm(y,0.7,sd_sim2[s], log=T)
-  log_lik3[s,] <-  dnorm(y,1,sd_sim3[s], log=T)
+for (s in 1:500) {
+  log_lik1[s, ] <- dnorm(y, -1, sd_sim1[s], log = T)
+  log_lik2[s, ] <- dnorm(y, 0.7, sd_sim2[s], log = T)
+  log_lik3[s, ] <- dnorm(y, 1, sd_sim3[s], log = T)
 }
 
-ll_list <- list(log_lik1, log_lik2,log_lik3)
-r_eff_list <- list(rep(0.9,50), rep(0.9,50), rep(0.9,50))
+ll_list <- list(log_lik1, log_lik2, log_lik3)
+r_eff_list <- list(rep(0.9, 50), rep(0.9, 50), rep(0.9, 50))
 
 loo_list <- lapply(1:length(ll_list), function(j) {
   loo(ll_list[[j]], r_eff = r_eff_list[[j]])
@@ -23,28 +21,46 @@ loo_list <- lapply(1:length(ll_list), function(j) {
 tol <- 0.01 # absolute tolerance of weights
 
 test_that("loo_model_weights throws correct errors and warnings", {
-  expect_error(loo_model_weights(log_lik1), "list of matrices or a list of 'psis_loo' objects")
+  expect_error(
+    loo_model_weights(log_lik1),
+    "list of matrices or a list of 'psis_loo' objects"
+  )
   expect_error(loo_model_weights(list(log_lik1)), "At least two models")
   expect_error(loo_model_weights(list(loo_list[[1]])), "At least two models")
-  expect_error(loo_model_weights(list(log_lik1), method = "pseudobma"), "At least two models")
+  expect_error(
+    loo_model_weights(list(log_lik1), method = "pseudobma"),
+    "At least two models"
+  )
 
-  expect_error(loo_model_weights(list(log_lik1, log_lik2[-1, ])), "same dimensions")
-  expect_error(loo_model_weights(list(log_lik1, log_lik2, log_lik3[, -1])), "same dimensions")
+  expect_error(
+    loo_model_weights(list(log_lik1, log_lik2[-1, ])),
+    "same dimensions"
+  )
+  expect_error(
+    loo_model_weights(list(log_lik1, log_lik2, log_lik3[, -1])),
+    "same dimensions"
+  )
 
   loo_list2 <- loo_list
   attr(loo_list2[[3]], "dims") <- c(10, 10)
   expect_error(loo_model_weights(loo_list2), "same dimensions")
 
-  expect_error(loo_model_weights(ll_list, r_eff_list = r_eff_list[-1]),
-               "one component for each model")
+  expect_error(
+    loo_model_weights(ll_list, r_eff_list = r_eff_list[-1]),
+    "one component for each model"
+  )
 
   r_eff_list[[3]] <- rep(0.9, 51)
-  expect_error(loo_model_weights(ll_list, r_eff_list = r_eff_list),
-               "same length as the number of columns")
+  expect_error(
+    loo_model_weights(ll_list, r_eff_list = r_eff_list),
+    "same length as the number of columns"
+  )
 
-  expect_error(loo_model_weights(list(loo_list[[1]], 2)),
-               "List elements must all be 'psis_loo' objects or log-likelihood matrices",
-               fixed = TRUE)
+  expect_error(
+    loo_model_weights(list(loo_list[[1]], 2)),
+    "List elements must all be 'psis_loo' objects or log-likelihood matrices",
+    fixed = TRUE
+  )
 
   expect_no_warning(loo_model_weights(ll_list))
 })
@@ -52,32 +68,43 @@ test_that("loo_model_weights throws correct errors and warnings", {
 
 test_that("loo_model_weights (stacking and pseudo-BMA) gives expected result", {
   w1 <- loo_model_weights(ll_list, method = "stacking", r_eff_list = r_eff_list)
-  expect_type(w1,"double")
+  expect_type(w1, "double")
   expect_s3_class(w1, "stacking_weights")
   expect_length(w1, 3)
-  expect_named(w1, paste0("model"  ,c(1:3)))
+  expect_named(w1, paste0("model", c(1:3)))
   expect_snapshot_value(as.numeric(w1), style = "serialize")
   expect_output(print(w1), "Method: stacking")
-  
+
   w1_b <- loo_model_weights(loo_list)
   expect_identical(w1, w1_b)
-  
-  w2 <- loo_model_weights(ll_list, r_eff_list=r_eff_list,
-    method = "pseudobma", BB = TRUE)
-    expect_type(w2, "double")
-    expect_s3_class(w2, "pseudobma_bb_weights")
-    expect_length(w2, 3)
-    expect_named(w2, paste0("model", c(1:3)))
-    expect_snapshot_value(as.numeric(w2), style = "serialize")
+
+  w2 <- loo_model_weights(
+    ll_list,
+    r_eff_list = r_eff_list,
+    method = "pseudobma",
+    BB = TRUE
+  )
+  expect_type(w2, "double")
+  expect_s3_class(w2, "pseudobma_bb_weights")
+  expect_length(w2, 3)
+  expect_named(w2, paste0("model", c(1:3)))
+  expect_snapshot_value(as.numeric(w2), style = "serialize")
   expect_output(print(w2), "Method: pseudo-BMA+")
 
-  w3 <- loo_model_weights(ll_list, r_eff_list=r_eff_list,
-                      method = "pseudobma", BB = FALSE)
-  expect_type(w3,"double")
+  w3 <- loo_model_weights(
+    ll_list,
+    r_eff_list = r_eff_list,
+    method = "pseudobma",
+    BB = FALSE
+  )
+  expect_type(w3, "double")
   expect_length(w3, 3)
-  expect_named(w3, paste0("model"  ,c(1:3)))
-  expect_equal(as.numeric(w3), c(5.365279e-05, 9.999436e-01, 2.707028e-06),
-               tolerance  = tol)
+  expect_named(w3, paste0("model", c(1:3)))
+  expect_equal(
+    as.numeric(w3),
+    c(5.365279e-05, 9.999436e-01, 2.707028e-06),
+    tolerance = tol
+  )
   expect_output(print(w3), "Method: pseudo-BMA")
 
   w3_b <- loo_model_weights(loo_list, method = "pseudobma", BB = FALSE)
@@ -108,4 +135,3 @@ test_that("loo_model_weights uses correct names for list of loo objects", {
     c("a", "b", "c")
   )
 })
-
diff --git a/tests/testthat/test_pointwise.R b/tests/testthat/test_pointwise.R
index 5707692c..7157b986 100644
--- a/tests/testthat/test_pointwise.R
+++ b/tests/testthat/test_pointwise.R
@@ -1,5 +1,3 @@
-library(loo)
-
 loo1 <- suppressWarnings(loo(example_loglik_matrix()))
 
 test_that("pointwise throws the right errors", {
@@ -27,8 +25,14 @@ test_that("pointwise throws the right errors", {
 
 test_that("pointwise returns correct estimate", {
   expect_equal(pointwise(loo1, "elpd_loo"), loo1$pointwise[, "elpd_loo"])
-  expect_equal(pointwise(loo1, "mcse_elpd_loo"), loo1$pointwise[, "mcse_elpd_loo"])
+  expect_equal(
+    pointwise(loo1, "mcse_elpd_loo"),
+    loo1$pointwise[, "mcse_elpd_loo"]
+  )
   expect_equal(pointwise(loo1, "p_loo"), loo1$pointwise[, "p_loo"])
   expect_equal(pointwise(loo1, "looic"), loo1$pointwise[, "looic"])
-  expect_equal(pointwise(loo1, "influence_pareto_k"), loo1$pointwise[, "influence_pareto_k"])
+  expect_equal(
+    pointwise(loo1, "influence_pareto_k"),
+    loo1$pointwise[, "influence_pareto_k"]
+  )
 })
diff --git a/tests/testthat/test_print_plot.R b/tests/testthat/test_print_plot.R
index d59b52fe..f69ed6f0 100644
--- a/tests/testthat/test_print_plot.R
+++ b/tests/testthat/test_print_plot.R
@@ -1,4 +1,3 @@
-library(loo)
 set.seed(1414)
 
 LLarr <- example_loglik_array()
@@ -27,26 +26,42 @@ test_that("plot methods throw appropriate errors/warnings", {
 
   loo1$diagnostics$pareto_k[1:5] <- Inf
   psis1$diagnostics$pareto_k[1:5] <- Inf
-  expect_warning(plot(loo1), regexp = "estimates are Inf/NA/NaN and not plotted.")
-  expect_warning(plot(psis1), regexp = "estimates are Inf/NA/NaN and not plotted.")
+  expect_warning(
+    plot(loo1),
+    regexp = "estimates are Inf/NA/NaN and not plotted."
+  )
+  expect_warning(
+    plot(psis1),
+    regexp = "estimates are Inf/NA/NaN and not plotted."
+  )
 })
 
 
-
 # printing ----------------------------------------------------------------
-lldim_msg <- paste0("Computed from ", prod(dim(LLarr)[1:2]) , " by ",
-                    dim(LLarr)[3], " log-likelihood matrix")
-lwdim_msg <- paste0("Computed from ", prod(dim(LLarr)[1:2]) , " by ",
-                    dim(LLarr)[3], " log-weights matrix")
-
-test_that("print.waic output is ok",{
+lldim_msg <- paste0(
+  "Computed from ",
+  prod(dim(LLarr)[1:2]),
+  " by ",
+  dim(LLarr)[3],
+  " log-likelihood matrix"
+)
+lwdim_msg <- paste0(
+  "Computed from ",
+  prod(dim(LLarr)[1:2]),
+  " by ",
+  dim(LLarr)[3],
+  " log-weights matrix"
+)
+
+test_that("print.waic output is ok", {
   expect_output(print(waic1), lldim_msg)
-  expect_output(print(waic1),
+  expect_output(
+    print(waic1),
     "p_waic estimates greater than 0.4. We recommend trying loo instead."
   )
 })
 
-test_that("print.psis_loo and print.psis output ok",{
+test_that("print.psis_loo and print.psis output ok", {
   expect_output(print(psis1), lwdim_msg)
   expect_output(print(psis1), "Pareto k estimates are good")
   expect_output(print(loo1), lldim_msg)
@@ -80,13 +95,23 @@ test_that("pareto_k_influence_values works for psis_loo objects, errors for psis
   kloo2 <- pareto_k_values(loo1)
   expect_identical(kloo, kloo2)
 
-  expect_error(pareto_k_influence_values(psis1), "No Pareto k influence estimates found")
-  expect_error(pareto_k_influence_values(waic1), "No Pareto k influence estimates found")
+  expect_error(
+    pareto_k_influence_values(psis1),
+    "No Pareto k influence estimates found"
+  )
+  expect_error(
+    pareto_k_influence_values(waic1),
+    "No Pareto k influence estimates found"
+  )
 })
 
 test_that("pareto_k_ids identifies correct observations", {
   for (j in 1:5) {
-    loo1$diagnostics$pareto_k <- psis1$diagnostics$pareto_k <- runif(32, .25, 1.25)
+    loo1$diagnostics$pareto_k <- psis1$diagnostics$pareto_k <- runif(
+      32,
+      .25,
+      1.25
+    )
     expect_identical(
       pareto_k_ids(loo1, threshold = 0.5),
       pareto_k_ids(psis1, threshold = 0.5)
@@ -105,7 +130,7 @@ test_that("pareto_k_ids identifies correct observations", {
 test_that("pareto_k_table gives correct output", {
   threshold <- ps_khat_threshold(dim(psis1)[1])
   psis1$diagnostics$pareto_k[1:10] <- runif(10, 0, threshold)
-  psis1$diagnostics$pareto_k[11:20] <- runif(10, threshold+0.01, 0.99)
+  psis1$diagnostics$pareto_k[11:20] <- runif(10, threshold + 0.01, 0.99)
   psis1$diagnostics$pareto_k[21:32] <- runif(12, 1, 10)
   k <- pareto_k_values(psis1)
   tab <- pareto_k_table(psis1)
@@ -115,9 +140,9 @@ test_that("pareto_k_table gives correct output", {
   expect_equal(sum(tab[, "Count"]), length(k))
   expect_equal(sum(tab[, "Proportion"]), 1)
 
-  expect_equal(sum(k <= threshold), tab[1,1])
-  expect_equal(sum(k > threshold & k <= 1), tab[2,1])
-  expect_equal(sum(k > 1), tab[3,1])
+  expect_equal(sum(k <= threshold), tab[1, 1])
+  expect_equal(sum(k > threshold & k <= 1), tab[2, 1])
+  expect_equal(sum(k > 1), tab[3, 1])
 
   # if n_eff is NULL
   psis1$diagnostics$n_eff <- NULL
@@ -126,13 +151,14 @@ test_that("pareto_k_table gives correct output", {
   expect_equal(unname(tab2[, "Min. n_eff"]), rep(NA_real_, 3))
 
   psis1$diagnostics$pareto_k[1:32] <- 0.4
-  expect_output(print(pareto_k_table(psis1)),
-                paste0("All Pareto k estimates are good (k < ", round(threshold,2), ")"),
-                fixed = TRUE)
+  expect_output(
+    print(pareto_k_table(psis1)),
+    paste0("All Pareto k estimates are good (k < ", round(threshold, 2), ")"),
+    fixed = TRUE
+  )
 })
 
 
-
 # psis_neff and mcse_loo --------------------------------------------------
 test_that("psis_n_eff_values extractor works", {
   n_eff_psis <- psis1$diagnostics$n_eff
@@ -159,4 +185,3 @@ test_that("mcse_loo returns NA when it should", {
 test_that("mcse_loo errors if not psis_loo object", {
   expect_error(mcse_loo(psis1), "psis_loo")
 })
-
diff --git a/tests/testthat/test_psis.R b/tests/testthat/test_psis.R
index 246eae04..de0a5512 100644
--- a/tests/testthat/test_psis.R
+++ b/tests/testthat/test_psis.R
@@ -1,6 +1,5 @@
-library(loo)
-options(mc.cores=1)
-options(loo.cores=NULL)
+options(mc.cores = 1)
+options(loo.cores = NULL)
 set.seed(123)
 
 LLarr <- example_loglik_array()
@@ -12,7 +11,7 @@ r_eff_vec <- relative_eff(exp(LLvec), chain_id = chain_id)
 psis1 <- psis(log_ratios = -LLarr, r_eff = r_eff_arr)
 
 test_that("psis results haven't changed", {
-    expect_snapshot_value(psis1, style = "serialize")
+  expect_snapshot_value(psis1, style = "serialize")
 })
 
 test_that("psis returns object with correct structure", {
@@ -46,12 +45,12 @@ test_that("psis throws correct errors and warnings", {
   # r_eff=NULL no warnings
   expect_silent(psis(-LLarr, r_eff = NULL))
   expect_silent(psis(-LLmat, r_eff = NULL))
-  expect_silent(psis(-LLmat[,1], r_eff = NULL))
+  expect_silent(psis(-LLmat[, 1], r_eff = NULL))
 
   # r_eff=NA disables warnings
   expect_silent(psis(-LLarr, r_eff = NA))
   expect_silent(psis(-LLmat, r_eff = NA))
-  expect_silent(psis(-LLmat[,1], r_eff = NA))
+  expect_silent(psis(-LLmat[, 1], r_eff = NA))
 
   # r_eff default and r_eff=NA give same answer
   expect_equal(
@@ -61,7 +60,7 @@ test_that("psis throws correct errors and warnings", {
 
   # r_eff=NULL and r_eff=NA give same answer
   expect_equal(
-    suppressWarnings(psis(-LLarr, r_eff=NULL)),
+    suppressWarnings(psis(-LLarr, r_eff = NULL)),
     psis(-LLarr, r_eff = NA)
   )
 
@@ -76,13 +75,13 @@ test_that("psis throws correct errors and warnings", {
   expect_error(psis(-LLarr, r_eff = r_eff_arr), "mix NA and not NA values")
 
   # tail length warnings
-  expect_snapshot(psis(-LLarr[1:5,, ]))
+  expect_snapshot(psis(-LLarr[1:5, , ]))
 
   # no NAs or non-finite values allowed
-  LLmat[1,1] <- NA
+  LLmat[1, 1] <- NA
   expect_error(psis(-LLmat), "NAs not allowed in input")
 
-  LLmat[1,1] <- 1
+  LLmat[1, 1] <- 1
   LLmat[10, 2] <- -Inf
   expect_error(psis(-LLmat), "All input values must be finite or -Inf")
   # log ratio of -Inf is allowed
@@ -90,7 +89,10 @@ test_that("psis throws correct errors and warnings", {
   expect_no_error(psis(-LLmat))
 
   # no lists allowed
-  expect_error(expect_warning(psis(as.list(-LLvec))), "List not allowed as input")
+  expect_error(
+    expect_warning(psis(as.list(-LLvec))),
+    "List not allowed as input"
+  )
 
   # if array, must be 3-D array
   dim(LLarr) <- c(2, 250, 2, 32)
@@ -105,10 +107,11 @@ test_that("throw_tail_length_warnings gives correct output", {
   expect_silent(throw_tail_length_warnings(10))
   expect_equal(throw_tail_length_warnings(10), 10)
   expect_warning(throw_tail_length_warnings(1), "Not enough tail samples")
-  expect_warning(throw_tail_length_warnings(c(1, 10, 2)),
-                 "Skipping the following columns: 1, 3")
-  expect_warning(throw_tail_length_warnings(rep(1, 21)),
-                 "11 more not printed")
+  expect_warning(
+    throw_tail_length_warnings(c(1, 10, 2)),
+    "Skipping the following columns: 1, 3"
+  )
+  expect_warning(throw_tail_length_warnings(rep(1, 21)), "11 more not printed")
 })
 
 
@@ -141,8 +144,11 @@ test_that("psis_n_eff methods works properly", {
 
 
 test_that("do_psis_i throws warning if all tail values the same", {
-  xx <- c(1,2,3,4,4,4,4,4,4,4,4)
-  expect_warning(val <- do_psis_i(xx, tail_len_i = 6), "all tail values are the same")
+  xx <- c(1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4)
+  expect_warning(
+    val <- do_psis_i(xx, tail_len_i = 6),
+    "all tail values are the same"
+  )
   expect_equal(val$pareto_k, Inf)
 })
 
@@ -150,9 +156,8 @@ test_that("psis_smooth_tail returns original tail values if k is infinite", {
   # skip on M1 Mac until we figure out why this test fails only on M1 Mac
   skip_if(Sys.info()[["sysname"]] == "Darwin" && R.version$arch == "aarch64")
 
-  xx <- c(1,2,3,4,4,4,4,4,4,4,4)
+  xx <- c(1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4)
   val <- suppressWarnings(psis_smooth_tail(xx, 3))
   expect_equal(val$tail, xx)
   expect_equal(val$k, Inf)
 })
-
diff --git a/tests/testthat/test_psis_approximate_posterior.R b/tests/testthat/test_psis_approximate_posterior.R
index c4d3f5e5..8d0ffe70 100644
--- a/tests/testthat/test_psis_approximate_posterior.R
+++ b/tests/testthat/test_psis_approximate_posterior.R
@@ -1,5 +1,3 @@
-library(loo)
-
 load(test_path("data-for-tests/test_data_psis_approximate_posterior.rda"))
 
 test_that("Laplace approximation, independent posterior", {
@@ -8,14 +6,25 @@ test_that("Laplace approximation, independent posterior", {
   ll <- test_data_psis_approximate_posterior$laplace_independent$log_liks
   expect_silent(
     psis_lap <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_lap, "psis")
   expect_lt(pareto_k_values(psis_lap), 0.7)
 
   expect_silent(
     psis_lap_ll <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, log_liks = ll , cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        log_liks = ll,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_lap_ll, "loo")
   expect_true(all(pareto_k_values(psis_lap_ll) < 0.7))
@@ -28,14 +37,25 @@ test_that("Laplace approximation, correlated posterior", {
   ll <- test_data_psis_approximate_posterior$laplace_correlated$log_liks
   expect_silent(
     psis_lap <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_lap, "psis")
   expect_lt(pareto_k_values(psis_lap), 0.7)
 
   expect_silent(
     psis_lap_ll <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, log_liks = ll , cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        log_liks = ll,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_lap_ll, "loo")
   expect_true(all(pareto_k_values(psis_lap_ll) < 0.7))
@@ -47,35 +67,56 @@ test_that("Laplace approximation, normal model", {
   ll <- test_data_psis_approximate_posterior$laplace_normal$log_liks
   expect_no_warning(
     psis_lap <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_lap, "psis")
   expect_gt(pareto_k_values(psis_lap), 0.5)
 
   expect_warning(
     psis_lap_ll <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, log_liks = ll , cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        log_liks = ll,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_lap_ll, "loo")
   expect_true(all(pareto_k_values(psis_lap_ll) > 0.5))
 })
 
 
-
 test_that("ADVI fullrank approximation, independent posterior", {
   log_p <- test_data_psis_approximate_posterior$fullrank_independent$log_p
   log_g <- test_data_psis_approximate_posterior$fullrank_independent$log_q
   ll <- test_data_psis_approximate_posterior$fullrank_independent$log_liks
   expect_silent(
     psis_advi <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi, "psis")
   expect_lt(pareto_k_values(psis_advi), 0.7)
 
   expect_silent(
     psis_advi_ll <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, log_liks = ll , cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        log_liks = ll,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi_ll, "loo")
   expect_true(all(pareto_k_values(psis_advi_ll) < 0.7))
@@ -88,14 +129,25 @@ test_that("ADVI fullrank approximation, correlated posterior", {
   ll <- test_data_psis_approximate_posterior$fullrank_correlated$log_liks
   expect_silent(
     psis_advi <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi, "psis")
   expect_lt(pareto_k_values(psis_advi), 0.7)
 
   expect_silent(
     psis_advi_ll <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, log_liks = ll , cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        log_liks = ll,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi_ll, "loo")
   expect_true(all(pareto_k_values(psis_advi_ll) < 0.7))
@@ -107,14 +159,25 @@ test_that("ADVI fullrank approximation, correlated posterior", {
   ll <- test_data_psis_approximate_posterior$fullrank_normal$log_liks
   expect_warning(
     psis_advi <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi, "psis")
   expect_gt(pareto_k_values(psis_advi), 0.7)
 
   expect_warning(
     psis_advi_ll <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, log_liks = ll , cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        log_liks = ll,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi_ll, "loo")
   expect_true(all(pareto_k_values(psis_advi_ll) > 0.7))
@@ -127,14 +190,25 @@ test_that("ADVI meanfield approximation, independent posterior", {
   ll <- test_data_psis_approximate_posterior$meanfield_independent$log_liks
   expect_silent(
     psis_advi <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi, "psis")
   expect_lt(pareto_k_values(psis_advi), 0.7)
 
   expect_silent(
     psis_advi_ll <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, log_liks = ll , cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        log_liks = ll,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi_ll, "loo")
   expect_true(all(pareto_k_values(psis_advi_ll) < 0.7))
@@ -147,14 +221,25 @@ test_that("ADVI meanfield approximation, correlated posterior", {
   ll <- test_data_psis_approximate_posterior$meanfield_correlated$log_liks
   expect_warning(
     psis_advi <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi, "psis")
   expect_gt(pareto_k_values(psis_advi), 0.7)
 
   expect_warning(
     psis_advi_ll <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, log_liks = ll , cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        log_liks = ll,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi_ll, "loo")
   expect_true(all(pareto_k_values(psis_advi_ll) > 0.5))
@@ -167,14 +252,25 @@ test_that("ADVI meanfield approximation, normal model", {
   ll <- test_data_psis_approximate_posterior$meanfield_normal$log_liks
   expect_warning(
     psis_advi <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi, "psis")
   expect_gt(pareto_k_values(psis_advi), 0.7)
 
   expect_warning(
     psis_advi_ll <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, log_liks = ll , cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        log_liks = ll,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi_ll, "loo")
   expect_true(all(pareto_k_values(psis_advi_ll) > 0.7))
@@ -182,20 +278,30 @@ test_that("ADVI meanfield approximation, normal model", {
 
 
 test_that("ADVI meanfield approximation, normal model", {
-
   log_p <- test_data_psis_approximate_posterior$meanfield_normal$log_p
   log_g <- test_data_psis_approximate_posterior$meanfield_normal$log_q
   ll <- test_data_psis_approximate_posterior$meanfield_normal$log_liks
   expect_warning(
     psis_advi <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi, "psis")
   expect_gt(pareto_k_values(psis_advi), 0.7)
 
   expect_warning(
     psis_advi_ll <-
-      psis_approximate_posterior(log_p = log_p, log_g = log_g, log_liks = ll , cores = 1, save_psis = FALSE)
+      psis_approximate_posterior(
+        log_p = log_p,
+        log_g = log_g,
+        log_liks = ll,
+        cores = 1,
+        save_psis = FALSE
+      )
   )
   expect_s3_class(psis_advi_ll, "loo")
   expect_true(all(pareto_k_values(psis_advi_ll) > 0.7))
@@ -203,14 +309,18 @@ test_that("ADVI meanfield approximation, normal model", {
 
 
 test_that("Deprecation of log_q argument", {
-
   log_p <- test_data_psis_approximate_posterior$laplace_independent$log_p
   log_g <- test_data_psis_approximate_posterior$laplace_independent$log_q
   ll <- test_data_psis_approximate_posterior$laplace_independent$log_liks
   expect_warning(
     psis_lap <-
-      loo:::psis_approximate_posterior(log_p = log_p, log_q = log_g, cores = 1, save_psis = FALSE)
-    , regexp = "argument log_q has been changed to log_g"
+      loo:::psis_approximate_posterior(
+        log_p = log_p,
+        log_q = log_g,
+        cores = 1,
+        save_psis = FALSE
+      ),
+    regexp = "argument log_q has been changed to log_g"
   )
   expect_s3_class(psis_lap, "psis")
   expect_lt(pareto_k_values(psis_lap), 0.7)
diff --git a/tests/testthat/test_psislw.R b/tests/testthat/test_psislw.R
index c65b1f9b..d272779b 100644
--- a/tests/testthat/test_psislw.R
+++ b/tests/testthat/test_psislw.R
@@ -1,4 +1,3 @@
-library(loo)
 SW <- suppressWarnings
 
 set.seed(123)
@@ -38,19 +37,24 @@ test_that("psislw function and matrix methods return same result", {
   set.seed(024)
 
   # fake data and posterior draws
-  N <- 50; K <- 10; S <- 100; a0 <- 3; b0 <- 2
+  N <- 50
+  K <- 10
+  S <- 100
+  a0 <- 3
+  b0 <- 2
   p <- rbeta(1, a0, b0)
   y <- rbinom(N, size = K, prob = p)
-  a <- a0 + sum(y); b <- b0 + N * K - sum(y)
+  a <- a0 + sum(y)
+  b <- b0 + N * K - sum(y)
   draws <- rbeta(S, a, b)
-  data <- data.frame(y,K)
+  data <- data.frame(y, K)
   llfun <- function(i, data, draws) {
     dbinom(data$y, size = data$K, prob = draws, log = TRUE)
   }
   psislw_with_fn <- SW(psislw(llfun = llfun, llargs = nlist(data, draws, N, S)))
 
   # Check that we get same answer if using log-likelihood matrix
-  ll <- sapply(1:N, function(i) llfun(i, data[i,, drop=FALSE], draws))
+  ll <- sapply(1:N, function(i) llfun(i, data[i, , drop = FALSE], draws))
   psislw_with_mat <- SW(psislw(-ll))
   expect_equal(psislw_with_fn, psislw_with_mat)
 })
@@ -58,8 +62,12 @@ test_that("psislw function and matrix methods return same result", {
 test_that("psislw_warnings helper works properly", {
   k <- c(0, 0.1, 0.55, 0.75)
   expect_silent(psislw_warnings(k[1:2]))
-  expect_warning(psislw_warnings(k[1:3]),
-                 "Some Pareto k diagnostic values are slightly high")
-  expect_warning(psislw_warnings(k),
-                 "Some Pareto k diagnostic values are too high")
+  expect_warning(
+    psislw_warnings(k[1:3]),
+    "Some Pareto k diagnostic values are slightly high"
+  )
+  expect_warning(
+    psislw_warnings(k),
+    "Some Pareto k diagnostic values are too high"
+  )
 })
diff --git a/tests/testthat/test_relative_eff.R b/tests/testthat/test_relative_eff.R
index 6def1df5..195a2c47 100644
--- a/tests/testthat/test_relative_eff.R
+++ b/tests/testthat/test_relative_eff.R
@@ -1,4 +1,3 @@
-library(loo)
 options(mc.cores = 1)
 set.seed(123)
 
@@ -13,8 +12,8 @@ test_that("relative_eff is equal to ESS / S", {
   dims <- dim(LLarr)
   ess <- r_eff <- rep(NA, dims[3])
   for (j in 1:dims[3]) {
-    r_eff[j] <- relative_eff(LLarr[,,1, drop=FALSE])
-    ess[j] <- posterior::ess_mean(LLarr[,,1])
+    r_eff[j] <- relative_eff(LLarr[,, 1, drop = FALSE])
+    ess[j] <- posterior::ess_mean(LLarr[,, 1])
   }
   S <- prod(dim(LLarr)[1:2])
   expect_equal(r_eff, ess / S)
@@ -30,7 +29,13 @@ test_that("relative_eff matrix and function methods return identical output", {
   source(test_path("data-for-tests/function_method_stuff.R"))
   chain <- rep(1, nrow(draws))
   r_eff_mat <- relative_eff(llmat_from_fn, chain_id = chain)
-  r_eff_fn <- relative_eff(llfun, chain_id = chain, data = data, draws = draws, cores = 1)
+  r_eff_fn <- relative_eff(
+    llfun,
+    chain_id = chain,
+    data = data,
+    draws = draws,
+    cores = 1
+  )
   expect_identical(r_eff_mat, r_eff_fn)
 })
 
diff --git a/tests/testthat/test_tisis.R b/tests/testthat/test_tisis.R
index 0d41c953..8616ef28 100644
--- a/tests/testthat/test_tisis.R
+++ b/tests/testthat/test_tisis.R
@@ -1,6 +1,5 @@
-library(loo)
-options(mc.cores=1)
-options(loo.cores=NULL)
+options(mc.cores = 1)
+options(loo.cores = NULL)
 set.seed(123)
 
 LLarr <- example_loglik_array()
@@ -23,7 +22,6 @@ test_that("tis and is runs", {
 })
 
 test_that("tis() and sis() returns object with correct structure for tis/sis", {
-
   expect_false(is.psis(tis1))
   expect_false(is.psis(is1))
   expect_true(is.tis(tis1))
@@ -75,7 +73,6 @@ test_that("psis methods give same results", {
 })
 
 
-
 test_that("tis throws correct errors and warnings", {
   # r_eff default no warnings
   expect_silent(tis(-LLarr))
@@ -85,23 +82,23 @@ test_that("tis throws correct errors and warnings", {
   # r_eff=NULL no warnings
   expect_silent(tis(-LLarr, r_eff = NULL))
   expect_silent(tis(-LLmat, r_eff = NULL))
-  expect_silent(tis(-LLmat[,1], r_eff = NULL))
+  expect_silent(tis(-LLmat[, 1], r_eff = NULL))
 
   # r_eff=NA no warnings
   expect_silent(tis(-LLarr, r_eff = NA))
   expect_silent(tis(-LLmat, r_eff = NA))
-  expect_silent(tis(-LLmat[,1], r_eff = NA))
+  expect_silent(tis(-LLmat[, 1], r_eff = NA))
 
   # r_eff default and r_eff=NA give same answer
   expect_equal(
     suppressWarnings(tis(-LLarr)),
-                     tis(-LLarr, r_eff = NA)
+    tis(-LLarr, r_eff = NA)
   )
 
   # r_eff=NULL and r_eff=NA give same answer
   expect_equal(
     suppressWarnings(tis(-LLarr, r_eff = NULL)),
-                     tis(-LLarr, r_eff = NA)
+    tis(-LLarr, r_eff = NA)
   )
 
   # r_eff scalar is fine
@@ -115,17 +112,20 @@ test_that("tis throws correct errors and warnings", {
   expect_error(tis(-LLarr, r_eff = r_eff_arr), "mix NA and not NA values")
 
   # no NAs or non-finite values allowed
-  LLmat[1,1] <- NA
+  LLmat[1, 1] <- NA
   expect_error(tis(-LLmat), "NAs not allowed in input")
 
-  LLmat[1,1] <- 1
+  LLmat[1, 1] <- 1
   LLmat[10, 2] <- -Inf
   expect_error(tis(-LLmat), "All input values must be finite or -Inf")
   LLmat[10, 2] <- Inf
   expect_no_error(tis(-LLmat))
 
   # no lists allowed
-  expect_error(expect_warning(tis(as.list(-LLvec)), "List not allowed as input"))
+  expect_error(expect_warning(
+    tis(as.list(-LLvec)),
+    "List not allowed as input"
+  ))
 
   # if array, must be 3-D array
   dim(LLarr) <- c(2, 250, 2, 32)
@@ -140,21 +140,73 @@ test_that("tis throws correct errors and warnings", {
 test_that("explict test of values for 'sis' and 'tis'", {
   lw <- 1:16
   expect_silent(tis_true <- tis(log_ratios = lw, r_eff = NA))
-  expect_equal(as.vector(weights(tis_true, log = TRUE, normalize = FALSE)),
-               c(-14.0723, -13.0723, -12.0723, -11.0723, -10.0723, -9.0723, -8.0723, -7.0723, -6.0723, -5.0723, -4.0723, -3.0723, -2.0723, -1.0723, -0.0723, 0.) + 15.07238, tolerance = 0.001)
+  expect_equal(
+    as.vector(weights(tis_true, log = TRUE, normalize = FALSE)),
+    c(
+      -14.0723,
+      -13.0723,
+      -12.0723,
+      -11.0723,
+      -10.0723,
+      -9.0723,
+      -8.0723,
+      -7.0723,
+      -6.0723,
+      -5.0723,
+      -4.0723,
+      -3.0723,
+      -2.0723,
+      -1.0723,
+      -0.0723,
+      0.
+    ) +
+      15.07238,
+    tolerance = 0.001
+  )
   expect_silent(is_true <- sis(log_ratios = lw, r_eff = NA))
-  expect_equal(as.vector(weights(is_true, log = TRUE, normalize = FALSE)),
-               lw, tolerance = 0.00001)
+  expect_equal(
+    as.vector(weights(is_true, log = TRUE, normalize = FALSE)),
+    lw,
+    tolerance = 0.00001
+  )
 
-  lw <- c(0.7609420, 1.3894140, 0.4158346, 2.5307927, 4.3379119, 2.4159240, 2.2462172, 0.8057697, 0.9333107, 1.5599302)
+  lw <- c(
+    0.7609420,
+    1.3894140,
+    0.4158346,
+    2.5307927,
+    4.3379119,
+    2.4159240,
+    2.2462172,
+    0.8057697,
+    0.9333107,
+    1.5599302
+  )
 
   expect_silent(tis_true <- tis(log_ratios = lw, r_eff = NA))
-  expect_equal(as.vector(weights(tis_true, log = TRUE, normalize = FALSE)),
-               c(-2.931, -2.303, -3.276, -1.161,  0, -1.276, -1.446, -2.886, -2.759, -2.132) + 3.692668,
-               tolerance = 0.001)
+  expect_equal(
+    as.vector(weights(tis_true, log = TRUE, normalize = FALSE)),
+    c(
+      -2.931,
+      -2.303,
+      -3.276,
+      -1.161,
+      0,
+      -1.276,
+      -1.446,
+      -2.886,
+      -2.759,
+      -2.132
+    ) +
+      3.692668,
+    tolerance = 0.001
+  )
   expect_silent(is_true <- sis(log_ratios = lw, r_eff = NA))
-  expect_equal(as.vector(weights(is_true, log = TRUE, normalize = FALSE)),
-               lw, tolerance = 0.00001)
+  expect_equal(
+    as.vector(weights(is_true, log = TRUE, normalize = FALSE)),
+    lw,
+    tolerance = 0.00001
+  )
 })