diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index c9ec15c..562fe0f 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -1,46 +1,51 @@
-# Workflow derived from https://github.com/r-lib/actions/tree/master/examples
-# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
-on:
-  push:
-    branches: [main, master]
-  pull_request:
-    branches: [main, master]
-
-name: R-CMD-check
-
-jobs:
-  R-CMD-check:
-    runs-on: ${{ matrix.config.os }}
-
-    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
-
-    strategy:
-      fail-fast: false
-      matrix:
-        config:
-          - {os: macOS-latest,   r: 'release'}
-          - {os: windows-latest, r: 'release'}
-          - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
-          - {os: ubuntu-latest,   r: 'release'}
-          - {os: ubuntu-latest,   r: 'oldrel-1'}
-
-    env:
-      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
-      R_KEEP_PKG_SOURCE: yes
-
-    steps:
-      - uses: actions/checkout@v4
-
-      - uses: r-lib/actions/setup-pandoc@v2
-
-      - uses: r-lib/actions/setup-r@v2
-        with:
-          r-version: ${{ matrix.config.r }}
-          http-user-agent: ${{ matrix.config.http-user-agent }}
-          use-public-rspm: true
-
-      - uses: r-lib/actions/setup-r-dependencies@v2
-        with:
-          extra-packages: rcmdcheck
-
-      - uses: r-lib/actions/check-r-package@v2
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+
+name: R-CMD-check.yaml
+
+permissions: read-all
+
+jobs:
+  R-CMD-check:
+    runs-on: ${{ matrix.config.os }}
+
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: macos-latest,   r: 'release'}
+          - {os: windows-latest, r: 'release'}
+          - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
+          - {os: ubuntu-latest,   r: 'release'}
+          - {os: ubuntu-latest,   r: 'oldrel-1'}
+
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          r-version: ${{ matrix.config.r }}
+          http-user-agent: ${{ matrix.config.http-user-agent }}
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::rcmdcheck
+          needs: check
+
+      - uses: r-lib/actions/check-r-package@v2
+        with:
+          upload-snapshots: true
+          build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")'
diff --git a/.github/workflows/linter.yaml b/.github/workflows/linter.yaml
index ae78a83..27f2b44 100644
--- a/.github/workflows/linter.yaml
+++ b/.github/workflows/linter.yaml
@@ -1,26 +1,30 @@
-name: Super-Linter
-
-on:
-  pull_request:
-    branches: [ master, main ]
-  push:
-    branches: [ master, main ]
-
-jobs:
-  lint:
-    name: Lint Code Base
-    runs-on: ubuntu-latest
-
-    steps:
-      - name: Checkout Code
-        uses: actions/checkout@v3
-
-      - name: Lint Code Base
-        uses: github/super-linter@v4
-        env:
-          VALIDATE_ALL_CODEBASE: true
-          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-
-          VALIDATE_R_LINTR: true
-
-          FILTER_REGEX_INCLUDE: R/.*|tests/.*\.R
+name: Lint Code Base
+
+on:
+  push:
+    branches: [main]
+  pull_request:
+    branches: [main]
+
+jobs:
+  run-lint:
+    name: Lint Code Base
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout Code
+        uses: actions/checkout@v4
+        with:
+          # Full git history is needed for Super-Linter to detect changed files
+          fetch-depth: 0
+
+      - name: Super-Linter
+        uses: super-linter/super-linter@v7
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          # Enable only R linting (optional; defaults to all languages if omitted)
+          VALIDATE_R: true
+          # Lint the entire codebase (true) or just changed files (false)
+          VALIDATE_ALL_CODEBASE: false
+          # Set the default branch (required for diff-based linting)
+          DEFAULT_BRANCH: main
+          FILTER_REGEX_INCLUDE: R/.*|tests/.*\.R
diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
index 38ddfe5..0ab748d 100644
--- a/.github/workflows/test-coverage.yaml
+++ b/.github/workflows/test-coverage.yaml
@@ -1,30 +1,62 @@
-# Workflow derived from https://github.com/r-lib/actions/tree/master/examples
-# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
-on:
-  push:
-    branches: [main, master]
-  pull_request:
-    branches: [main, master]
-
-name: test-coverage
-
-jobs:
-  test-coverage:
-    runs-on: ubuntu-latest
-    env:
-      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
-
-    steps:
-      - uses: actions/checkout@v4
-
-      - uses: r-lib/actions/setup-r@v2
-        with:
-          use-public-rspm: true
-
-      - uses: r-lib/actions/setup-r-dependencies@v2
-        with:
-          extra-packages: covr
-
-      - name: Test coverage
-        run: covr::codecov()
-        shell: Rscript {0}
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+
+name: test-coverage.yaml
+
+permissions: read-all
+
+jobs:
+  test-coverage:
+    runs-on: ubuntu-latest
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::covr, any::xml2
+          needs: coverage
+
+      - name: Test coverage
+        run: |
+          cov <- covr::package_coverage(
+            quiet = FALSE,
+            clean = FALSE,
+            install_path = file.path(normalizePath(Sys.getenv("RUNNER_TEMP"), winslash = "/"), "package")
+          )
+          print(cov)
+          covr::to_cobertura(cov)
+        shell: Rscript {0}
+
+      - uses: codecov/codecov-action@v5
+        with:
+          # Fail if error if not on PR, or if on PR and token is given
+          fail_ci_if_error: ${{ github.event_name != 'pull_request' || secrets.CODECOV_TOKEN }}
+          files: ./cobertura.xml
+          plugins: noop
+          disable_search: true
+          token: ${{ secrets.CODECOV_TOKEN }}
+
+      - name: Show testthat output
+        if: always()
+        run: |
+          ## --------------------------------------------------------------------
+          find '${{ runner.temp }}/package' -name 'testthat.Rout*' -exec cat '{}' \; || true
+        shell: bash
+
+      - name: Upload test results
+        if: failure()
+        uses: actions/upload-artifact@v4
+        with:
+          name: coverage-test-failures
+          path: ${{ runner.temp }}/package
diff --git a/DESCRIPTION b/DESCRIPTION
index 902b1d4..f0bd4e4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -16,7 +16,7 @@ Description: Contains a tidy method for adonis results, custom permutation
 License: GPL-3 + file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.3
 Imports:
     generics (>= 0.0.2),
     rlang,
@@ -25,11 +25,10 @@ Imports:
     dplyr (>= 1.0.0),
     tibble (>= 3.0.0),
     ape,
-    magrittr,
     ggplot2,
     stringr
 Depends:
-    R (>= 2.10)
+    R (>= 4.1.0)
 Suggests:
     covr,
     testthat
diff --git a/NAMESPACE b/NAMESPACE
index e08d9b4..3f8817e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,7 +4,6 @@ S3method(plot,nmdsplus)
 S3method(plot,pcoaplus)
 S3method(tidy,adonis)
 S3method(tidy,anova.cca)
-export("%>%")
 export(adonisplus)
 export(adonispost)
 export(nmdsplus)
@@ -13,4 +12,3 @@ export(shuffle_between_groups)
 export(shuffle_within_groups)
 export(tidy)
 importFrom(generics,tidy)
-importFrom(magrittr,"%>%")
diff --git a/R/adonis-custom.R b/R/adonis-custom.R
index 2e727fb..c4f61b3 100644
--- a/R/adonis-custom.R
+++ b/R/adonis-custom.R
@@ -1,3 +1,5 @@
+utils::globalVariables(c("SampleID", "subject_id", "study_group", "term"))
+
 #' Permutational multivariate analysis of variance, plus
 #'
 #' @param data Data to use in the test.
@@ -31,11 +33,18 @@
 #' subjects, preserving the value within each subject. Conversely, the values
 #' of the time point will be shuffled only within each subject.
 #' @export
-adonisplus <- function(data, distmat, formula, sample_id_var = SampleID,
-                       rep_meas_var = subject_id, shuffle = NULL,
-                       permutations = 999, seed = 42) {
-  sample_ids <- data %>%
-    dplyr::pull({{ sample_id_var }}) %>%
+adonisplus <- function(
+  data,
+  distmat,
+  formula,
+  sample_id_var = SampleID,
+  rep_meas_var = subject_id,
+  shuffle = NULL,
+  permutations = 999,
+  seed = 42
+) {
+  sample_ids <- data |>
+    dplyr::pull({{ sample_id_var }}) |>
     as.character()
   distmat <- usedist::dist_subset(distmat, sample_ids)
   formula <- stats::as.formula(formula)
@@ -44,12 +53,15 @@ adonisplus <- function(data, distmat, formula, sample_id_var = SampleID,
 
   set.seed(seed)
   a_observed <- vegan::adonis2(
-    formula = formula, data = data, permutations = permutations, by = "terms"
+    formula = formula,
+    data = data,
+    permutations = permutations,
+    by = "terms"
   )
   result <- tidy.anova.cca(a_observed)
 
   if (!is.null(shuffle)) {
-    rep_meas_vals <- data %>%
+    rep_meas_vals <- data |>
       dplyr::pull({{ rep_meas_var }})
 
     # Here, vars is the variables that will undergo custom permutations
@@ -82,8 +94,11 @@ adonisplus <- function(data, distmat, formula, sample_id_var = SampleID,
         trial_data[[var]] <- new_vals
       }
       trial_a <- vegan::adonis2(
-        formula, trial_data, permutations = 4, by = "terms"
-        )
+        formula,
+        trial_data,
+        permutations = 4,
+        by = "terms"
+      )
       trial_result <- tidy.anova.cca(trial_a)
       trial_result$statistic[term_idx]
     })
@@ -109,29 +124,43 @@ adonisplus <- function(data, distmat, formula, sample_id_var = SampleID,
 #'   comparison.
 #' @return The results in tidy format.
 #' @export
-adonispost <- function(data, ..., which = study_group, alpha = 0.05) {
+adonispost <- function(
+  data,
+  ...,
+  which = study_group,
+  alpha = 0.05
+) {
   var_name <- rlang::as_name(rlang::ensym(which))
 
-  result_main <- adonisplus(data, ...) %>%
-    dplyr::mutate(comparison = paste("All", var_name)) %>%
-    dplyr::select(comparison, term, dplyr::everything()) %>%
+  result_cols <- c(
+    "comparison",
+    "term",
+    "df",
+    "sumsq",
+    "r.squared",
+    "statistic",
+    "p.value"
+  )
+  result_main <- adonisplus(data, ...) |>
+    dplyr::mutate(comparison = paste("All", var_name)) |>
+    dplyr::select(dplyr::all_of(result_cols)) |>
     dplyr::filter(!(term %in% c("Residual", "Total")))
 
-  var_levels <- data %>%
-    dplyr::pull({{ which }}) %>%
-    as.factor() %>%
+  var_levels <- data |>
+    dplyr::pull({{ which }}) |>
+    as.factor() |>
     levels()
   pairs <- utils::combn(var_levels, 2, simplify = FALSE)
 
   make_pairwise_comparison <- function(pair) {
-    pair_data <- data %>%
+    pair_data <- data |>
       dplyr::filter({{ which }} %in% pair)
-    adonisplus(pair_data, ...) %>%
-      dplyr::mutate(comparison = paste(pair, collapse = " - ")) %>%
-      dplyr::select(comparison, term, dplyr::everything()) %>%
+    adonisplus(pair_data, ...) |>
+      dplyr::mutate(comparison = paste(pair, collapse = " - ")) |>
+      dplyr::select(dplyr::all_of(result_cols)) |>
       dplyr::filter(!(term %in% c("Residual", "Total")))
   }
-  result_posthoc <- lapply(pairs, make_pairwise_comparison) %>%
+  result_posthoc <- lapply(pairs, make_pairwise_comparison) |>
     dplyr::bind_rows()
   dplyr::bind_rows(result_main, result_posthoc)
 }
diff --git a/R/ordination.R b/R/ordination.R
index 6272e13..303df3d 100644
--- a/R/ordination.R
+++ b/R/ordination.R
@@ -1,3 +1,5 @@
+utils::globalVariables(c("Axis.1", "Axis.2", "MDS1", "MDS2"))
+
 #' Principal coordinates analysis (PCoA), plus
 #'
 #' @param data A data frame giving information on the objects for the PCoA.
@@ -19,25 +21,30 @@
 #'
 #'   For \code{plot.pcoaplus}, a ggplot object.
 #' @export
-pcoaplus <- function(data, distmat, sample_id_var = SampleID, num_axes = 2) {
+pcoaplus <- function(
+  data,
+  distmat,
+  sample_id_var = SampleID,
+  num_axes = 2
+) {
   num_axes <- as.integer(num_axes)
   stopifnot("num_axes must be 2 or more" = num_axes >= 2)
   axis_numbers <- 1:num_axes
   axis_names <- paste0("Axis.", axis_numbers)
 
-  sample_ids <- data %>%
-    dplyr::pull({{ sample_id_var }}) %>%
+  sample_ids <- data |>
+    dplyr::pull({{ sample_id_var }}) |>
     as.character()
   stopifnot("Duplicated sample IDs" = anyDuplicated(sample_ids) == 0)
 
   distmat <- usedist::dist_subset(distmat, sample_ids)
   sample_id_var_name <- rlang::as_name(rlang::ensym(sample_id_var))
   pcoa_obj <- ape::pcoa(distmat)
-  pcoa_df <- pcoa_obj$vectors %>%
-    `[`(sample_ids, axis_names) %>%
-    as.data.frame() %>%
-    tibble::rownames_to_column(sample_id_var_name) %>%
-    tibble::as_tibble() %>%
+  pcoa_df <- pcoa_obj$vectors |>
+    (\(x) x[sample_ids, axis_names])() |>
+    as.data.frame() |>
+    tibble::rownames_to_column(sample_id_var_name) |>
+    tibble::as_tibble() |>
     dplyr::left_join(data, by = sample_id_var_name)
 
   pctvar <- (pcoa_obj$values$Relative_eig * 100)[axis_numbers]
@@ -56,7 +63,7 @@ pcoaplus <- function(data, distmat, sample_id_var = SampleID, num_axes = 2) {
 #' @describeIn pcoaplus Make a principal coordinates scatter plot
 #' @export
 plot.pcoaplus <- function(x, ...) {
-  x %>%
+  x |>
     ggplot2::ggplot() +
     ggplot2::geom_point(ggplot2::aes(x = Axis.1, y = Axis.2, ...)) +
     ggplot2::coord_equal() +
@@ -83,23 +90,28 @@ plot.pcoaplus <- function(x, ...) {
 #'
 #'   For \code{plot.nmdsplus}, a ggplot object.
 #' @export
-nmdsplus <- function(data, distmat, sample_id_var = SampleID, num_axes = 2) {
+nmdsplus <- function(
+  data,
+  distmat,
+  sample_id_var = SampleID,
+  num_axes = 2
+) {
   num_axes <- as.integer(num_axes)
   stopifnot("num_axes must be 2 or more" = num_axes >= 2)
 
-  sample_ids <- data %>%
-    dplyr::pull({{ sample_id_var }}) %>%
+  sample_ids <- data |>
+    dplyr::pull({{ sample_id_var }}) |>
     as.character()
   stopifnot("Duplicated sample IDs" = anyDuplicated(sample_ids) == 0)
 
   distmat <- usedist::dist_subset(distmat, sample_ids)
   sample_id_var_name <- rlang::as_name(rlang::ensym(sample_id_var))
   nmds_obj <- vegan::monoMDS(distmat, k = num_axes)
-  nmds_df <- nmds_obj$points %>%
-    `[`(sample_ids, ) %>%
-    as.data.frame() %>%
-    tibble::rownames_to_column(sample_id_var_name) %>%
-    tibble::as_tibble() %>%
+  nmds_df <- nmds_obj$points |>
+    (\(x) x[sample_ids, ])() |>
+    as.data.frame() |>
+    tibble::rownames_to_column(sample_id_var_name) |>
+    tibble::as_tibble() |>
     dplyr::left_join(data, by = sample_id_var_name)
 
   attr(nmds_df, "stress") <- nmds_obj$stress
@@ -111,7 +123,7 @@ nmdsplus <- function(data, distmat, sample_id_var = SampleID, num_axes = 2) {
 #' @describeIn nmdsplus Make an NMDS scatter plot
 #' @export
 plot.nmdsplus <- function(x, ...) {
-  x %>%
+  x |>
     ggplot2::ggplot() +
     ggplot2::geom_point(ggplot2::aes(x = MDS1, y = MDS2, ...))
 }
diff --git a/R/shuffle.R b/R/shuffle.R
index e2f23fd..4c0986b 100644
--- a/R/shuffle.R
+++ b/R/shuffle.R
@@ -24,7 +24,7 @@ shuffle_between_groups <- function(x, g) {
   g <- as.factor(g)
   x <- as.factor(x)
   # What is the value of x for every unique value of g?
-  x_vals_per_g <- tapply(x, g, unique, simplify = F)
+  x_vals_per_g <- tapply(x, g, unique, simplify = FALSE)
   # Check that x has one unique value for each value of g. If not, raise an
   # error with an informative message.
   if (!all(sapply(x_vals_per_g, length) == 1)) {
diff --git a/R/tidy.R b/R/tidy.R
index 4189cb7..e6084b7 100644
--- a/R/tidy.R
+++ b/R/tidy.R
@@ -22,7 +22,13 @@ generics::tidy
 tidy.adonis <- function(x, ...) {
   ret <- tibble::as_tibble(x$aov.tab, rownames = "term")
   colnames(ret) <- c(
-    "term", "df", "sumsq", "meansq", "statistic", "r.squared", "p.value"
+    "term",
+    "df",
+    "sumsq",
+    "meansq",
+    "statistic",
+    "r.squared",
+    "p.value"
   )
   attr(ret, "heading") <- NULL
   ret
@@ -45,7 +51,12 @@ tidy.adonis <- function(x, ...) {
 tidy.anova.cca <- function(x, ...) {
   ret <- tibble::as_tibble(x, rownames = "term")
   colnames(ret) <- c(
-    "term", "df", "sumsq", "r.squared", "statistic", "p.value"
+    "term",
+    "df",
+    "sumsq",
+    "r.squared",
+    "statistic",
+    "p.value"
   )
   attr(ret, "Random.seed") <- NULL
   attr(ret, "F.perm") <- NULL
diff --git a/R/util.R b/R/util.R
index 8805b15..66257ed 100644
--- a/R/util.R
+++ b/R/util.R
@@ -3,8 +3,12 @@ check_lhs <- function(formula, expected) {
   expected_lhs <- expected[[2]]
   if (formula_lhs != expected_lhs) {
     stop(
-      "LHS of formula must be exactly '", expected_lhs, "'. ",
-      "Saw '", formula_lhs, "'."
+      "LHS of formula must be exactly '",
+      expected_lhs,
+      "'. ",
+      "Saw '",
+      formula_lhs,
+      "'."
     )
   }
 }
diff --git a/R/utils-pipe.R b/R/utils-pipe.R
deleted file mode 100644
index fd0b1d1..0000000
--- a/R/utils-pipe.R
+++ /dev/null
@@ -1,14 +0,0 @@
-#' Pipe operator
-#'
-#' See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details.
-#'
-#' @name %>%
-#' @rdname pipe
-#' @keywords internal
-#' @export
-#' @importFrom magrittr %>%
-#' @usage lhs \%>\% rhs
-#' @param lhs A value or the magrittr placeholder.
-#' @param rhs A function call using the magrittr semantics.
-#' @return The result of calling `rhs(lhs)`.
-NULL
diff --git a/README.Rmd b/README.Rmd
index 4c4ad4b..15f328f 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -6,7 +6,7 @@ output: github_document
 
 ```{r, echo = FALSE}
 knitr::opts_chunk$set(
-  fig.path = "tools/readme/"
+  fig.path = "man/figures/README-"
 )
 ```
 
@@ -17,8 +17,8 @@ devtools::load_all()
 # adonisplus
 
 <!-- badges: start -->
-  [![R-CMD-check](https://github.com/PennChopMicrobiomeProgram/adonisplus/workflows/R-CMD-check/badge.svg)](https://github.com/PennChopMicrobiomeProgram/adonisplus/actions)
-  [![Codecov test coverage](https://codecov.io/gh/PennChopMicrobiomeProgram/adonisplus/branch/master/graph/badge.svg)](https://app.codecov.io/gh/PennChopMicrobiomeProgram/adonisplus?branch=master)
+[![R-CMD-check](https://github.com/PennChopMicrobiomeProgram/adonisplus/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/PennChopMicrobiomeProgram/adonisplus/actions/workflows/R-CMD-check.yaml)
+[![Codecov test coverage](https://codecov.io/gh/PennChopMicrobiomeProgram/adonisplus/graph/badge.svg)](https://app.codecov.io/gh/PennChopMicrobiomeProgram/adonisplus)
 <!-- badges: end -->
 
 The `adonisplus` package provides some utilities for permutational multivariate
@@ -67,9 +67,9 @@ the last time point obtained for each participant. We'll make a new data frame,
 `farmm_final`, containing just the samples we need.
 
 ```{r}
-farmm_final <- farmm_samples %>%
-  group_by(subject_id) %>%
-  filter(study_day == max(study_day)) %>%
+farmm_final <- farmm_samples |>
+  group_by(subject_id) |>
+  filter(study_day == max(study_day)) |>
   ungroup()
 ```
 
@@ -83,7 +83,7 @@ distance matrix over and over as you filter and re-arrange the samples. Here,
 we specify this column as `sample_id`.
 
 ```{r}
-farmm_final %>%
+farmm_final |>
   pcoaplus(distmat = farmm_bc, sample_id_var = sample_id)
 ```
 
@@ -95,8 +95,8 @@ default mapping to color the samples by diet. We also add a few elements to the
 `ggplot` object, to set the color scale and change the theme.
 
 ```{r farmm_final_pcoa}
-farmm_final %>%
-  pcoaplus(distmat = farmm_bc, sample_id_var = sample_id) %>%
+farmm_final |>
+  pcoaplus(distmat = farmm_bc, sample_id_var = sample_id) |>
   plot(color = diet) +
   scale_color_brewer(palette = "Set1") +
   theme_bw()
@@ -113,7 +113,7 @@ an additional argument to let the function know which column of the data frame
 corresponds to the IDs in the distance matrix.
 
 ```{r}
-farmm_final %>%
+farmm_final |>
   adonisplus(
     distmat = farmm_bc, formula = distmat ~ diet,
     sample_id_var = sample_id)
@@ -132,7 +132,7 @@ are the same as `adonisplus()`, but we add an additional argument, `which`, to
 specify the variable on which we want to carry out pairwise comparisons.
 
 ```{r}
-farmm_final %>%
+farmm_final |>
   adonispost(
     distmat = farmm_bc, formula = distmat ~ diet,
     sample_id_var = sample_id, which = diet)
@@ -148,8 +148,8 @@ time for each diet. To get an overview, we'll re-generate Figure 2A from the
 paper.
 
 ```{r farmm_pcoa}
-farmm_samples %>%
-  pcoaplus(distmat = farmm_bc, sample_id_var = "sample_id") %>%
+farmm_samples |>
+  pcoaplus(distmat = farmm_bc, sample_id_var = "sample_id") |>
   plot(color = study_day) +
   facet_grid(~ diet) +
   scale_color_viridis_c(direction = -1) +
@@ -166,7 +166,7 @@ To investigate both claims, we'll limit the data set to the pre-antibiotics
 period, and call it `farmm_preabx`.
 
 ```{r}
-farmm_preabx <- farmm_samples %>%
+farmm_preabx <- farmm_samples |>
   filter(antibiotics %in% "pre")
 ```
 
@@ -178,7 +178,7 @@ tell it to randomly re-assign diet between subjects, and randomly shuffle the
 study days within each subject.
 
 ```{r}
-farmm_preabx %>%
+farmm_preabx |>
   adonisplus(
     distmat = farmm_bc, formula = distmat ~ diet * study_day,
     sample_id_var = sample_id, rep_meas_var = subject_id,
@@ -202,7 +202,7 @@ diet pairs are different. As before, we use `adonispost()` and tell it to run
 pairwise comparisons of diet by setting `which = diet`.
 
 ```{r}
-farmm_preabx %>%
+farmm_preabx |>
   adonispost(
     distmat = farmm_bc, formula = distmat ~ diet * study_day,
     sample_id_var = sample_id, rep_meas_var = subject_id,
diff --git a/README.md b/README.md
index fe0af2f..8b67e53 100644
--- a/README.md
+++ b/README.md
@@ -5,9 +5,9 @@
 
 <!-- badges: start -->
 
-[![R-CMD-check](https://github.com/PennChopMicrobiomeProgram/adonisplus/workflows/R-CMD-check/badge.svg)](https://github.com/PennChopMicrobiomeProgram/adonisplus/actions)
+[![R-CMD-check](https://github.com/PennChopMicrobiomeProgram/adonisplus/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/PennChopMicrobiomeProgram/adonisplus/actions/workflows/R-CMD-check.yaml)
 [![Codecov test
-coverage](https://codecov.io/gh/PennChopMicrobiomeProgram/adonisplus/branch/master/graph/badge.svg)](https://app.codecov.io/gh/PennChopMicrobiomeProgram/adonisplus?branch=master)
+coverage](https://codecov.io/gh/PennChopMicrobiomeProgram/adonisplus/graph/badge.svg)](https://app.codecov.io/gh/PennChopMicrobiomeProgram/adonisplus)
 <!-- badges: end -->
 
 The `adonisplus` package provides some utilities for permutational
@@ -76,9 +76,9 @@ We’ll make a new data frame, `farmm_final`, containing just the samples
 we need.
 
 ``` r
-farmm_final <- farmm_samples %>%
-  group_by(subject_id) %>%
-  filter(study_day == max(study_day)) %>%
+farmm_final <- farmm_samples |>
+  group_by(subject_id) |>
+  filter(study_day == max(study_day)) |>
   ungroup()
 ```
 
@@ -93,7 +93,7 @@ over and over as you filter and re-arrange the samples. Here, we specify
 this column as `sample_id`.
 
 ``` r
-farmm_final %>%
+farmm_final |>
   pcoaplus(distmat = farmm_bc, sample_id_var = sample_id)
 ```
 
@@ -123,14 +123,14 @@ by diet. We also add a few elements to the `ggplot` object, to set the
 color scale and change the theme.
 
 ``` r
-farmm_final %>%
-  pcoaplus(distmat = farmm_bc, sample_id_var = sample_id) %>%
+farmm_final |>
+  pcoaplus(distmat = farmm_bc, sample_id_var = sample_id) |>
   plot(color = diet) +
   scale_color_brewer(palette = "Set1") +
   theme_bw()
 ```
 
-![](tools/readme/farmm_final_pcoa-1.png)<!-- -->
+![](man/figures/README-farmm_final_pcoa-1.png)<!-- -->
 
 It looks like the microbiome of the EEN diet group is different from
 that in the omnivore and vegan groups. Let’s test for differences
@@ -144,7 +144,7 @@ additional argument to let the function know which column of the data
 frame corresponds to the IDs in the distance matrix.
 
 ``` r
-farmm_final %>%
+farmm_final |>
   adonisplus(
     distmat = farmm_bc, formula = distmat ~ diet,
     sample_id_var = sample_id)
@@ -153,7 +153,7 @@ farmm_final %>%
     ## # A tibble: 3 × 6
     ##   term        df sumsq r.squared statistic p.value
     ##   <chr>    <dbl> <dbl>     <dbl>     <dbl>   <dbl>
-    ## 1 Model        2  1.84     0.168      2.73   0.001
+    ## 1 diet         2  1.84     0.168      2.73   0.001
     ## 2 Residual    27  9.07     0.832     NA     NA    
     ## 3 Total       29 10.9      1         NA     NA
 
@@ -171,7 +171,7 @@ from each other. To run the pairwise comparisons, we’ll use
 the variable on which we want to carry out pairwise comparisons.
 
 ``` r
-farmm_final %>%
+farmm_final |>
   adonispost(
     distmat = farmm_bc, formula = distmat ~ diet,
     sample_id_var = sample_id, which = diet)
@@ -180,10 +180,10 @@ farmm_final %>%
     ## # A tibble: 4 × 7
     ##   comparison       term     df sumsq r.squared statistic p.value
     ##   <chr>            <chr> <dbl> <dbl>     <dbl>     <dbl>   <dbl>
-    ## 1 All diet         Model     2 1.84     0.168       2.73   0.001
-    ## 2 Omnivore - Vegan Model     1 0.445    0.0671      1.30   0.159
-    ## 3 Omnivore - EEN   Model     1 1.15     0.159       3.39   0.001
-    ## 4 Vegan - EEN      Model     1 1.16     0.165       3.56   0.001
+    ## 1 All diet         diet      2 1.84     0.168       2.73   0.001
+    ## 2 Omnivore - Vegan diet      1 0.445    0.0671      1.30   0.159
+    ## 3 Omnivore - EEN   diet      1 1.15     0.159       3.39   0.001
+    ## 4 Vegan - EEN      diet      1 1.16     0.165       3.56   0.001
 
 The results are as expected. So, far we could have done all of this work
 using the functions in `vegan` and `ape`, without much difficulty.
@@ -195,15 +195,15 @@ over time for each diet. To get an overview, we’ll re-generate Figure 2A
 from the paper.
 
 ``` r
-farmm_samples %>%
-  pcoaplus(distmat = farmm_bc, sample_id_var = "sample_id") %>%
+farmm_samples |>
+  pcoaplus(distmat = farmm_bc, sample_id_var = "sample_id") |>
   plot(color = study_day) +
   facet_grid(~ diet) +
   scale_color_viridis_c(direction = -1) +
   theme_bw()
 ```
 
-![](tools/readme/farmm_pcoa-1.png)<!-- -->
+![](man/figures/README-farmm_pcoa-1.png)<!-- -->
 
 The paper reports that “EEN led to a significant change in the
 microbiota composition within 3 days of the dietary phase relative to
@@ -216,7 +216,7 @@ To investigate both claims, we’ll limit the data set to the
 pre-antibiotics period, and call it `farmm_preabx`.
 
 ``` r
-farmm_preabx <- farmm_samples %>%
+farmm_preabx <- farmm_samples |>
   filter(antibiotics %in% "pre")
 ```
 
@@ -229,19 +229,21 @@ between subjects, and randomly shuffle the study days within each
 subject.
 
 ``` r
-farmm_preabx %>%
+farmm_preabx |>
   adonisplus(
     distmat = farmm_bc, formula = distmat ~ diet * study_day,
     sample_id_var = sample_id, rep_meas_var = subject_id,
     shuffle = c(diet = "between", study_day = "within"))
 ```
 
-    ## # A tibble: 3 × 6
-    ##   term        df sumsq r.squared statistic p.value
-    ##   <chr>    <dbl> <dbl>     <dbl>     <dbl>   <dbl>
-    ## 1 Model        5  5.51     0.132      4.38   0.001
-    ## 2 Residual   144 36.2      0.868     NA     NA    
-    ## 3 Total      149 41.7      1         NA     NA
+    ## # A tibble: 5 × 6
+    ##   term              df  sumsq r.squared statistic p.value
+    ##   <chr>          <dbl>  <dbl>     <dbl>     <dbl>   <dbl>
+    ## 1 diet               2  4.34     0.104       8.62   0.001
+    ## 2 study_day          1  0.504    0.0121      2.00   0.001
+    ## 3 diet:study_day     2  0.663    0.0159      1.32   0.001
+    ## 4 Residual         144 36.2      0.868      NA     NA    
+    ## 5 Total            149 41.7      1          NA     NA
 
 If you run this function yourself, you’ll notice that it takes a lot
 longer than it takes to run `adonis()`. As we randomly re-assign diets
@@ -256,7 +258,7 @@ which diet pairs are different. As before, we use `adonispost()` and
 tell it to run pairwise comparisons of diet by setting `which = diet`.
 
 ``` r
-farmm_preabx %>%
+farmm_preabx |>
   adonispost(
     distmat = farmm_bc, formula = distmat ~ diet * study_day,
     sample_id_var = sample_id, rep_meas_var = subject_id,
@@ -264,13 +266,21 @@ farmm_preabx %>%
     which = diet)
 ```
 
-    ## # A tibble: 4 × 7
-    ##   comparison       term     df sumsq r.squared statistic p.value
-    ##   <chr>            <chr> <dbl> <dbl>     <dbl>     <dbl>   <dbl>
-    ## 1 All diet         Model     5  5.51    0.132       4.38   0.001
-    ## 2 Omnivore - Vegan Model     3  1.79    0.0679      2.36   0.15 
-    ## 3 Omnivore - EEN   Model     3  3.04    0.110       4.07   0.001
-    ## 4 Vegan - EEN      Model     3  4.06    0.149       5.35   0.001
+    ## # A tibble: 12 × 7
+    ##    comparison       term              df  sumsq r.squared statistic p.value
+    ##    <chr>            <chr>          <dbl>  <dbl>     <dbl>     <dbl>   <dbl>
+    ##  1 All diet         diet               2 4.34     0.104       8.62    0.001
+    ##  2 All diet         study_day          1 0.504    0.0121      2.00    0.001
+    ##  3 All diet         diet:study_day     2 0.663    0.0159      1.32    0.001
+    ##  4 Omnivore - Vegan diet               1 1.68     0.0635      6.61    0.119
+    ##  5 Omnivore - Vegan study_day          1 0.0804   0.00305     0.317   0.726
+    ##  6 Omnivore - Vegan diet:study_day     1 0.0358   0.00136     0.141   0.988
+    ##  7 Omnivore - EEN   diet               1 1.89     0.0683      7.59    0.025
+    ##  8 Omnivore - EEN   study_day          1 0.628    0.0227      2.53    0.001
+    ##  9 Omnivore - EEN   diet:study_day     1 0.523    0.0189      2.10    0.001
+    ## 10 Vegan - EEN      diet               1 2.99     0.109      11.8     0.001
+    ## 11 Vegan - EEN      study_day          1 0.668    0.0245      2.64    0.001
+    ## 12 Vegan - EEN      diet:study_day     1 0.399    0.0146      1.58    0.001
 
 In the pairwise comparisons, we find that the microbiome of the omnivore
 and vegan groups was not different during the pre-antibiotics period,
diff --git a/data-raw/farmm.R b/data-raw/farmm.R
index b73abda..afd35ea 100644
--- a/data-raw/farmm.R
+++ b/data-raw/farmm.R
@@ -2,35 +2,56 @@ library(tidyverse)
 library(usedist)
 
 farmm_samples <- read_tsv(
-  "data-raw/farmm_samples.tsv", show_col_types = FALSE) %>%
-  filter(SampleType %in% "Feces") %>%
-  filter(Keep) %>%
-  filter(!(study_day %in% "PS")) %>%
-  mutate(study_day = as.integer(study_day)) %>%
-  mutate(study_group = fct_relevel(study_group, "Omnivore", "Vegan", "EEN")) %>%
-  mutate(current_antibiotics = fct_recode(
-    current_antibiotics, pre = "Pre Antibiotics",
-    current = "Antibiotics Treatment", post = "Post Antibiotics")) %>%
-  mutate(current_antibiotics = fct_relevel(
-    current_antibiotics, "pre", "current", "post")) %>%
-  mutate(host_frac = host / (host + non_host)) %>%
-  mutate(new_sample_id = paste(
-    "farmm", SubjectID, sprintf("%02d", study_day), sep = "."))
+  "data-raw/farmm_samples.tsv",
+  show_col_types = FALSE
+) |>
+  filter(SampleType %in% "Feces") |>
+  filter(Keep) |>
+  filter(!(study_day %in% "PS")) |>
+  mutate(study_day = as.integer(study_day)) |>
+  mutate(study_group = fct_relevel(study_group, "Omnivore", "Vegan", "EEN")) |>
+  mutate(
+    current_antibiotics = fct_recode(
+      current_antibiotics,
+      pre = "Pre Antibiotics",
+      current = "Antibiotics Treatment",
+      post = "Post Antibiotics"
+    )
+  ) |>
+  mutate(
+    current_antibiotics = fct_relevel(
+      current_antibiotics,
+      "pre",
+      "current",
+      "post"
+    )
+  ) |>
+  mutate(host_frac = host / (host + non_host)) |>
+  mutate(
+    new_sample_id = paste(
+      "farmm",
+      SubjectID,
+      sprintf("%02d", study_day),
+      sep = "."
+    )
+  )
 
 farmm_new_sample_ids <- farmm_samples$new_sample_id
 names(farmm_new_sample_ids) <- farmm_samples$SampleID
 
 farmm_bc <- read_tsv(
-  "data-raw/farmm_bc.tsv", show_col_types = FALSE) %>%
-  rename(sample_id = `...1`) %>%
-  column_to_rownames("sample_id") %>%
-  as.matrix() %>%
-  `[`(farmm_samples$SampleID, farmm_samples$SampleID)
+  "data-raw/farmm_bc.tsv",
+  show_col_types = FALSE
+) |>
+  rename(sample_id = `...1`) |>
+  column_to_rownames("sample_id") |>
+  as.matrix() |>
+  (\(x) x[farmm_samples$SampleID, farmm_samples$SampleID])()
 
 colnames(farmm_bc) <- farmm_new_sample_ids[colnames(farmm_bc)]
 rownames(farmm_bc) <- farmm_new_sample_ids[rownames(farmm_bc)]
 
-farmm_samples <- farmm_samples %>%
+farmm_samples <- farmm_samples |>
   select(
     sample_id = new_sample_id,
     subject_id = SubjectID,
@@ -42,11 +63,12 @@ farmm_samples <- farmm_samples %>%
     age = Age,
     bacterial_16S_copies = copy_num_per_gram_feces,
     num_reads = non_host,
-    host_frac) %>%
+    host_frac
+  ) |>
   arrange(diet, subject_id, study_day)
 
-farmm_bc <- farmm_bc %>%
-  as.dist() %>%
+farmm_bc <- farmm_bc |>
+  as.dist() |>
   dist_subset(farmm_samples$sample_id)
 
 usethis::use_data(farmm_samples, overwrite = TRUE)
diff --git a/tools/readme/farmm_final_pcoa-1.png b/man/figures/README-farmm_final_pcoa-1.png
similarity index 100%
rename from tools/readme/farmm_final_pcoa-1.png
rename to man/figures/README-farmm_final_pcoa-1.png
diff --git a/man/figures/README-farmm_pcoa-1.png b/man/figures/README-farmm_pcoa-1.png
new file mode 100644
index 0000000..fa1476a
Binary files /dev/null and b/man/figures/README-farmm_pcoa-1.png differ
diff --git a/man/pipe.Rd b/man/pipe.Rd
deleted file mode 100644
index 1f8f237..0000000
--- a/man/pipe.Rd
+++ /dev/null
@@ -1,20 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils-pipe.R
-\name{\%>\%}
-\alias{\%>\%}
-\title{Pipe operator}
-\usage{
-lhs \%>\% rhs
-}
-\arguments{
-\item{lhs}{A value or the magrittr placeholder.}
-
-\item{rhs}{A function call using the magrittr semantics.}
-}
-\value{
-The result of calling `rhs(lhs)`.
-}
-\description{
-See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details.
-}
-\keyword{internal}
diff --git a/tests/testthat/test-adonis-custom.R b/tests/testthat/test-adonis-custom.R
index 7f27a37..082f17e 100644
--- a/tests/testthat/test-adonis-custom.R
+++ b/tests/testthat/test-adonis-custom.R
@@ -1,7 +1,5 @@
 example_data <- tibble::tibble(
-  SampleID = c(
-    "C1", "C2", "C3", "C4", "C5", "P1", "P2", "P3", "P4", "P5", "P6"
-  ),
+  SampleID = paste0(c(rep("C", 5), rep("P", 6)), c(1:5, 1:6)),
   time_point = c(rep("C", 5), rep("P", 6)),
   SubjectID = paste0("S", c(1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 6)),
   study_group = rep(c("G1", "G2", "G1", "G2"), c(3, 2, 3, 3)),
@@ -10,22 +8,82 @@ example_data <- tibble::tibble(
 
 example_dist <- structure(
   c(
-    0.367, 0.495, 0.425, 0.432, 0.698, 0.812, 0.698, 0.605, 0.41, 0.522, 0.461,
-    0.417, 0.411, 0.687, 0.809, 0.655, 0.618, 0.385, 0.47, 0.348, 0.233, 0.416,
-    0.84, 0.472, 0.497, 0.34, 0.412, 0.393, 0.559, 0.781, 0.499, 0.484, 0.409,
-    0.467, 0.474, 0.804, 0.598, 0.495, 0.257, 0.403, 0.82, 0.59, 0.586, 0.525,
-    0.594, 0.877, 0.832, 0.794, 0.83, 0.498, 0.6, 0.556, 0.447, 0.483, 0.348
+    0.367,
+    0.495,
+    0.425,
+    0.432,
+    0.698,
+    0.812,
+    0.698,
+    0.605,
+    0.41,
+    0.522,
+    0.461,
+    0.417,
+    0.411,
+    0.687,
+    0.809,
+    0.655,
+    0.618,
+    0.385,
+    0.47,
+    0.348,
+    0.233,
+    0.416,
+    0.84,
+    0.472,
+    0.497,
+    0.34,
+    0.412,
+    0.393,
+    0.559,
+    0.781,
+    0.499,
+    0.484,
+    0.409,
+    0.467,
+    0.474,
+    0.804,
+    0.598,
+    0.495,
+    0.257,
+    0.403,
+    0.82,
+    0.59,
+    0.586,
+    0.525,
+    0.594,
+    0.877,
+    0.832,
+    0.794,
+    0.83,
+    0.498,
+    0.6,
+    0.556,
+    0.447,
+    0.483,
+    0.348
   ),
   Labels = c("C1", "C2", "C3", "C4", "C5", "P1", "P2", "P3", "P4", "P5", "P6"),
-  Size = 11L, call = quote(as.dist.default(m = cp_dist)), class = "dist",
-  Diag = FALSE, Upper = FALSE
+  Size = 11L,
+  call = quote(as.dist.default(m = cp_dist)),
+  class = "dist",
+  Diag = FALSE,
+  Upper = FALSE
 )
 
 test_that("adonispost works", {
-  example_data_3groups <- example_data %>%
-    dplyr::mutate(study_group = rep(c("G1", "G2", "G3", "G1", "G2", "G3"), c(2, 2, 1, 2, 2, 2)))
+  example_data_3groups <- example_data |>
+    dplyr::mutate(
+      study_group = rep(
+        c("G1", "G2", "G3", "G1", "G2", "G3"),
+        c(2, 2, 1, 2, 2, 2)
+      )
+    )
   observed <- adonispost(
-    example_data_3groups, example_dist, distmat ~ study_group,
+    example_data_3groups,
+    example_dist,
+    distmat ~ study_group,
     which = study_group
   )
   expected <- tibble::tibble(
@@ -33,13 +91,22 @@ test_that("adonispost works", {
     term = c("study_group", "study_group", "study_group", "study_group"),
     df = c(2, 1, 1, 1),
     sumsq = c(
-      0.396595356060606, 0.2534925, 0.172570726190476, 0.160944833333333
+      0.396595356060606,
+      0.2534925,
+      0.172570726190476,
+      0.160944833333333
     ),
     r.squared = c(
-      0.2457711729374, 0.18721007197831, 0.162938895956797, 0.26470972004068
+      0.2457711729374,
+      0.18721007197831,
+      0.162938895956797,
+      0.26470972004068
     ),
     statistic = c(
-      1.30343028067264, 1.38198123911777, 0.973279580007741, 1.80003549112144
+      1.30343028067264,
+      1.38198123911777,
+      0.973279580007741,
+      1.80003549112144
     ),
     p.value = c(0.154, 0.173, 0.478, 0.092)
   )
@@ -68,7 +135,9 @@ test_that("adonisplus fails if LHS is not 'distmat'", {
 
 test_that("adonisplus works for one group with restricted permutations", {
   observed <- adonisplus(
-    example_data, example_dist, distmat ~ study_group,
+    example_data,
+    example_dist,
+    distmat ~ study_group,
     rep_meas_var = SubjectID,
     shuffle = c(study_group = "between"),
     permutations = 9
@@ -89,27 +158,42 @@ test_that("adonisplus works for one group with restricted permutations", {
 
 test_that("adonisplus works for two groups", {
   observed <- adonisplus(
-    example_data, example_dist, distmat ~ study_group * time_point,
+    example_data,
+    example_dist,
+    distmat ~ study_group * time_point,
     rep_meas_var = SubjectID,
     shuffle = c(study_group = "between", time_point = "within"),
     permutations = 9
   )
   expected <- tibble::tibble(
     term = c(
-      "study_group", "time_point", "study_group:time_point",
-      "Residual", "Total"
+      "study_group",
+      "time_point",
+      "study_group:time_point",
+      "Residual",
+      "Total"
     ),
     df = c(1, 1, 1, 7, 10),
     sumsq = c(
-      0.141236772727273, 0.2546555, 0.161871833333333,
-      1.05591316666667, 1.61367727272727
+      0.141236772727273,
+      0.2546555,
+      0.161871833333333,
+      1.05591316666667,
+      1.61367727272727
     ),
     r.squared = c(
-      0.0875247951460386, 0.157810675222318, 0.100312395835974,
-      0.654352133795669, 1
+      0.0875247951460386,
+      0.157810675222318,
+      0.100312395835974,
+      0.654352133795669,
+      1
     ),
     statistic = c(
-      0.936305598131642, 1.68819610956015, 1.07310228634646, NA, NA
+      0.936305598131642,
+      1.68819610956015,
+      1.07310228634646,
+      NA,
+      NA
     ),
     p.value = c(0.4, 0.1, 0.3, NA, NA)
   )
@@ -118,28 +202,46 @@ test_that("adonisplus works for two groups", {
 
 test_that("adonisplus works for two groups with covariate", {
   observed <- adonisplus(
-    example_data, example_dist, distmat ~ age + study_group * time_point,
+    example_data,
+    example_dist,
+    distmat ~ age + study_group * time_point,
     rep_meas_var = SubjectID,
     shuffle = c(study_group = "between", time_point = "within"),
     permutations = 9
   )
   expected <- tibble::tibble(
     term = c(
-      "age", "study_group", "time_point", "study_group:time_point",
-      "Residual", "Total"
+      "age",
+      "study_group",
+      "time_point",
+      "study_group:time_point",
+      "Residual",
+      "Total"
     ),
     df = c(1, 1, 1, 1, 6, 10),
     sumsq = c(
-      0.147625953372434, 0.12677620170778, 0.270223650980392,
-      0.148155646153846, 0.920895820512821, 1.61367727272727
+      0.147625953372434,
+      0.12677620170778,
+      0.270223650980392,
+      0.148155646153846,
+      0.920895820512821,
+      1.61367727272727
     ),
     r.squared = c(
-      0.0914841869978943, 0.0785635417009472, 0.167458298847878,
-      0.0918124389912542, 0.570681533462026, 1
+      0.0914841869978943,
+      0.0785635417009472,
+      0.167458298847878,
+      0.0918124389912542,
+      0.570681533462026,
+      1
     ),
     statistic = c(
-      0.961841394547052, 0.825997027354399, 1.76061381729312,
-      0.965292552232516, NA, NA
+      0.961841394547052,
+      0.825997027354399,
+      1.76061381729312,
+      0.965292552232516,
+      NA,
+      NA
     ),
     p.value = c(0.7, 0.5, 0.1, 0.5, NA, NA)
   )
diff --git a/tests/testthat/test-farmm.R b/tests/testthat/test-farmm.R
index 374d046..74f6858 100644
--- a/tests/testthat/test-farmm.R
+++ b/tests/testthat/test-farmm.R
@@ -1,9 +1,21 @@
 test_that("farmm_samples data is available", {
   expect_equal(nrow(farmm_samples), 414)
-  expect_equal(colnames(farmm_samples), c(
-    "sample_id", "subject_id", "study_day", "diet", "antibiotics", "height",
-    "weight", "age", "bacterial_16S_copies", "num_reads", "host_frac"
-  ))
+  expect_equal(
+    colnames(farmm_samples),
+    c(
+      "sample_id",
+      "subject_id",
+      "study_day",
+      "diet",
+      "antibiotics",
+      "height",
+      "weight",
+      "age",
+      "bacterial_16S_copies",
+      "num_reads",
+      "host_frac"
+    )
+  )
   expect_equal(levels(farmm_samples$diet), c("Omnivore", "Vegan", "EEN"))
   expect_equal(levels(farmm_samples$antibiotics), c("pre", "current", "post"))
 })
diff --git a/tests/testthat/test-ordination.R b/tests/testthat/test-ordination.R
index 9f9cbfa..c7b9670 100644
--- a/tests/testthat/test-ordination.R
+++ b/tests/testthat/test-ordination.R
@@ -5,7 +5,17 @@ ax <- function(x) {
 
 example_data <- tibble::tibble(
   SampleID = c(
-    "C1", "C2", "C3", "C4", "C5", "P1", "P2", "P3", "P4", "P5", "P6"
+    "C1",
+    "C2",
+    "C3",
+    "C4",
+    "C5",
+    "P1",
+    "P2",
+    "P3",
+    "P4",
+    "P5",
+    "P6"
   ),
   time_point = c(rep("C", 5), rep("P", 6)),
   SubjectID = paste0("S", c(1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 6)),
@@ -15,50 +25,138 @@ example_data <- tibble::tibble(
 
 example_dist <- structure(
   c(
-    0.367, 0.495, 0.425, 0.432, 0.698, 0.812, 0.698, 0.605, 0.41, 0.522, 0.461,
-    0.417, 0.411, 0.687, 0.809, 0.655, 0.618, 0.385, 0.47, 0.348, 0.233, 0.416,
-    0.84, 0.472, 0.497, 0.34, 0.412, 0.393, 0.559, 0.781, 0.499, 0.484, 0.409,
-    0.467, 0.474, 0.804, 0.598, 0.495, 0.257, 0.403, 0.82, 0.59, 0.586, 0.525,
-    0.594, 0.877, 0.832, 0.794, 0.83, 0.498, 0.6, 0.556, 0.447, 0.483, 0.348
+    0.367,
+    0.495,
+    0.425,
+    0.432,
+    0.698,
+    0.812,
+    0.698,
+    0.605,
+    0.41,
+    0.522,
+    0.461,
+    0.417,
+    0.411,
+    0.687,
+    0.809,
+    0.655,
+    0.618,
+    0.385,
+    0.47,
+    0.348,
+    0.233,
+    0.416,
+    0.84,
+    0.472,
+    0.497,
+    0.34,
+    0.412,
+    0.393,
+    0.559,
+    0.781,
+    0.499,
+    0.484,
+    0.409,
+    0.467,
+    0.474,
+    0.804,
+    0.598,
+    0.495,
+    0.257,
+    0.403,
+    0.82,
+    0.59,
+    0.586,
+    0.525,
+    0.594,
+    0.877,
+    0.832,
+    0.794,
+    0.83,
+    0.498,
+    0.6,
+    0.556,
+    0.447,
+    0.483,
+    0.348
   ),
   Labels = c("C1", "C2", "C3", "C4", "C5", "P1", "P2", "P3", "P4", "P5", "P6"),
-  Size = 11L, call = quote(as.dist.default(m = cp_dist)), class = "dist",
-  Diag = FALSE, Upper = FALSE
+  Size = 11L,
+  call = quote(as.dist.default(m = cp_dist)),
+  class = "dist",
+  Diag = FALSE,
+  Upper = FALSE
 )
 
 expected_axes <- tibble::tibble(
   Axis.1 = c(
-    0.0218046657139025, 0.0373054791505106, 0.13833146386329,
-    0.0482957523422068, 0.0906680281431401, 0.0120057411305378,
-    -0.680446513436419, 0.0915777190805797, 0.0658719664029247,
-    0.0777538535042316, 0.0968318441050968
+    0.0218046657139025,
+    0.0373054791505106,
+    0.13833146386329,
+    0.0482957523422068,
+    0.0906680281431401,
+    0.0120057411305378,
+    -0.680446513436419,
+    0.0915777190805797,
+    0.0658719664029247,
+    0.0777538535042316,
+    0.0968318441050968
   ),
   Axis.2 = c(
-    -0.303086502259418, -0.277649583508124, 0.0538409488164232,
-    -0.0371164252609219, -0.0685894632025572, 0.275689960626523,
-    0.020089905908166, 0.306076031692592, 0.162663666052037,
-    -0.100614869652573, -0.0313036692121472
+    -0.303086502259418,
+    -0.277649583508124,
+    0.0538409488164232,
+    -0.0371164252609219,
+    -0.0685894632025572,
+    0.275689960626523,
+    0.020089905908166,
+    0.306076031692592,
+    0.162663666052037,
+    -0.100614869652573,
+    -0.0313036692121472
   ),
   Axis.3 = c(
-    0.0373672324641792, 0.0446548183955713, -0.106676969942398,
-    0.0726191880827978, -0.148599519606932, -0.30153337102521,
-    0.0100963735032305, 0.236709673732038, 0.169041261896952,
-    -0.0652297761452144, 0.0515510886449863
+    0.0373672324641792,
+    0.0446548183955713,
+    -0.106676969942398,
+    0.0726191880827978,
+    -0.148599519606932,
+    -0.30153337102521,
+    0.0100963735032305,
+    0.236709673732038,
+    0.169041261896952,
+    -0.0652297761452144,
+    0.0515510886449863
   )
 )
 
 expected_nmds_axes <- tibble::tibble(
   MDS1 = c(
-    0.836706003503475, -0.478138656124125, -0.873168624865947,
-    0.831612773313092, -0.873168624865947, 1.12730524434323,
-    0.347992717293982, 0.213938532476495, -0.476699598859632,
-    -0.873168624865126, 0.216788858650501
+    0.836706003503475,
+    -0.478138656124125,
+    -0.873168624865947,
+    0.831612773313092,
+    -0.873168624865947,
+    1.12730524434323,
+    0.347992717293982,
+    0.213938532476495,
+    -0.476699598859632,
+    -0.873168624865126,
+    0.216788858650501
   ),
   MDS2 = c(
-    0.724394134805183, -1.0267096839492, 0.000665297749562208,
-    -0.72769854123764, 0.000665297749561982, -0.00467674105569537,
-    0.00496749526862173, -1.0525703201631, 1.02721958016895,
-    0.00066529774956238, 1.05307818291419
+    0.724394134805183,
+    -1.0267096839492,
+    0.000665297749562208,
+    -0.72769854123764,
+    0.000665297749561982,
+    -0.00467674105569537,
+    0.00496749526862173,
+    -1.0525703201631,
+    1.02721958016895,
+    0.00066529774956238,
+    1.05307818291419
   )
 )
 
diff --git a/tools/readme/farmm_pcoa-1.png b/tools/readme/farmm_pcoa-1.png
deleted file mode 100644
index 34dd6d5..0000000
Binary files a/tools/readme/farmm_pcoa-1.png and /dev/null differ