From 03cc34fe19d69cecd447efcda3b76c957325ed12 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 13:37:07 +0100
Subject: [PATCH 01/36] run checks monthly
---
.github/workflows/R-CMD-check.yaml | 3 +--
1 file changed, 1 insertion(+), 2 deletions(-)
diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index 6117b1c..ad0c65b 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -4,7 +4,7 @@ on:
branches:
- master
schedule:
- - cron: '0 8 * * 5'
+ - cron: '0 8 1 * *'
name: R-CMD-check
@@ -21,7 +21,6 @@ jobs:
config:
- { os: windows-latest, r: 'release'}
- { os: macOS-latest, r: 'release'}
- #- { os: ubuntu-16.04, r: 'release', bioc: 'release', cran: "https://demo.rstudiopm.com/all/__linux__/xenial/latest"}
- { os: ubuntu-latest, r: 'release'}
env:
From 40a360c25786cc73a4f8dd2f1097917fc81ec673 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 13:37:28 +0100
Subject: [PATCH 02/36] rename branch
---
.github/workflows/R-CMD-check.yaml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index ad0c65b..ad35fcf 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -2,7 +2,7 @@ on:
push:
pull_request:
branches:
- - master
+ - main
schedule:
- cron: '0 8 1 * *'
From fb162483f6c6c77ff43f3740d09dc3fd6b27ab4c Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 13:59:29 +0100
Subject: [PATCH 03/36] add basic unit tests
---
DESCRIPTION | 4 +++-
tests/testthat.R | 12 ++++++++++
tests/testthat/test-equality.R | 42 +++++++++++++++++++++++++++++++++
tests/testthat/test-existence.R | 35 +++++++++++++++++++++++++++
4 files changed, 92 insertions(+), 1 deletion(-)
create mode 100644 tests/testthat.R
create mode 100644 tests/testthat/test-equality.R
create mode 100644 tests/testthat/test-existence.R
diff --git a/DESCRIPTION b/DESCRIPTION
index 4c63dfc..2680af9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -37,6 +37,8 @@ Imports:
magrittr
Suggests:
knitr,
- rmarkdown
+ rmarkdown,
+ testthat (>= 3.0.0)
VignetteBuilder: knitr
URL: https://stopsack.github.io/batchtma, https://github.com/stopsack/batchtma
+Config/testthat/edition: 3
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..ac0e094
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,12 @@
+# This file is part of the standard setup for testthat.
+# It is recommended that you do not modify it.
+#
+# Where should you do additional test configuration?
+# Learn more about the roles of various files in:
+# * https://r-pkgs.org/testing-design.html#sec-tests-files-overview
+# * https://testthat.r-lib.org/articles/special-files.html
+
+library(testthat)
+library(batchtma)
+
+test_check("batchtma")
diff --git a/tests/testthat/test-equality.R b/tests/testthat/test-equality.R
new file mode 100644
index 0000000..6ca325a
--- /dev/null
+++ b/tests/testthat/test-equality.R
@@ -0,0 +1,42 @@
+test_that("equivalent approaches give same result", {
+ df <- data.frame(
+ tma = rep(1:2, times = 10),
+ biomarker = rep(1:2, times = 10) +
+ runif(max = 5, n = 20),
+ confounder = rep(0:1, times = 10) +
+ runif(max = 10, n = 20),
+ unity = 1
+ )
+
+ df_adj2 <- adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = simple,
+ suffix = "adj"
+ ) |> # drop all attributes
+ data.frame()
+
+ df_adj3 <- adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = standardize,
+ confounders = unity,
+ suffix = "adj"
+ ) |>
+ data.frame()
+
+ df_adj4 <- adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = ipw,
+ confounders = unity,
+ suffix = "adj"
+ ) |>
+ data.frame()
+
+ expect_equal(df_adj2, df_adj3)
+ expect_equal(df_adj2, df_adj4)
+})
diff --git a/tests/testthat/test-existence.R b/tests/testthat/test-existence.R
new file mode 100644
index 0000000..826d553
--- /dev/null
+++ b/tests/testthat/test-existence.R
@@ -0,0 +1,35 @@
+test_that("new variables are generated", {
+ df <- data.frame(
+ tma = rep(1:2, times = 10),
+ biomarker = rep(1:2, times = 10) +
+ runif(max = 5, n = 20),
+ confounder = rep(0:1, times = 10) +
+ runif(max = 10, n = 20)
+ )
+
+ df_adj5 <- adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = quantreg,
+ suffix = "adj"
+ )
+
+ expect_equal(
+ object = sum(!is.na(df_adj5$biomarkeradj)),
+ expected = sum(!is.na(df_adj5$biomarker))
+ )
+
+ df_adj6 <- adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = quantnorm,
+ suffix = "adj"
+ )
+
+ expect_equal(
+ object = sum(!is.na(df_adj6$biomarkeradj)),
+ expected = sum(!is.na(df_adj6$biomarker))
+ )
+})
From d6d47361ec133017598843eaa5efa21a237fc46c Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 13:59:49 +0100
Subject: [PATCH 04/36] increase dev version no
---
DESCRIPTION | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 2680af9..fbffdf3 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: batchtma
Title: Batch Effect Adjustments
-Version: 0.1.7
+Version: 0.1.7-9000
Authors@R:
c(person(given = "Konrad",
family = "Stopsack",
From 705ad41891693bb6b4348e52a67c29bf37723143 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 14:04:09 +0100
Subject: [PATCH 05/36] move functions to separate .R files
---
R/adjust_batch.R | 414 --------------------------------------
R/batch_quantnorm.R | 31 +++
R/batch_rq.R | 103 ++++++++++
R/batchmean_ipw.R | 154 ++++++++++++++
R/batchmean_simple.R | 29 +++
R/batchmean_standardize.R | 74 +++++++
R/factor_drop.R | 14 ++
7 files changed, 405 insertions(+), 414 deletions(-)
create mode 100644 R/batch_quantnorm.R
create mode 100644 R/batch_rq.R
create mode 100644 R/batchmean_ipw.R
create mode 100644 R/batchmean_simple.R
create mode 100644 R/batchmean_standardize.R
create mode 100644 R/factor_drop.R
diff --git a/R/adjust_batch.R b/R/adjust_batch.R
index 035cd6f..1b07e0e 100644
--- a/R/adjust_batch.R
+++ b/R/adjust_batch.R
@@ -2,420 +2,6 @@
## quiets concerns of R CMD check re: the .'s that appear in pipelines
if (getRversion() >= "2.15.1") utils::globalVariables(c("."))
-#' Drop empty factor levels
-#'
-#' @description
-#' Avoids predict() issues. This is \code{forcats::fct_drop()}
-#' without bells and whistles.
-#'
-#' @param f factor
-#'
-#' @return Factor
-#' @noRd
-factor_drop <- function(f) {
- factor_levels <- levels(f)
- factor(f, levels = setdiff(factor_levels, factor_levels[table(f) == 0]))
-}
-
-
-#' Batch means for approach 2: Unadjusted means
-#'
-#' @param data Data set
-#' @param markers Biomarkers to adjust
-#'
-#' @importFrom magrittr %>%
-#' @importFrom rlang .data
-#' @return Tibble of means per marker and batch
-#' @noRd
-batchmean_simple <- function(data, markers) {
- values <- data %>%
- dplyr::select(.data$.id, .data$.batchvar, {{ markers }}) %>%
- dplyr::group_by(.data$.batchvar) %>%
- dplyr::summarize_at(
- .vars = dplyr::vars(-.data$.id),
- .funs = mean,
- na.rm = TRUE
- ) %>%
- dplyr::mutate_at(
- .vars = dplyr::vars(-.data$.batchvar),
- .funs = ~ . - mean(., na.rm = TRUE)
- ) %>%
- tidyr::pivot_longer(
- col = c(-.data$.batchvar),
- names_to = "marker",
- values_to = "batchmean"
- )
- return(list(list(values = values, models = NULL)))
-}
-
-#' Batch means for approach 3: Marginal standardization
-#'
-#' @param data Data set
-#' @param markers Vector of variables to batch-adjust
-#' @param confounders Confounders: features that differ
-#' between batches that should be retained
-#'
-#' @return Tibble with conditional means per marker and batch
-#' @noRd
-batchmean_standardize <- function(data, markers, confounders) {
- res <- data %>%
- tidyr::pivot_longer(
- cols = {{ markers }},
- names_to = "marker",
- values_to = "value"
- ) %>%
- dplyr::filter(!is.na(.data$value)) %>%
- dplyr::group_by(.data$marker) %>%
- tidyr::nest(data = c(-.data$marker)) %>%
- dplyr::mutate(
- data = purrr::map(
- .x = .data$data,
- .f = ~ .x %>%
- dplyr::mutate(.batchvar = factor_drop(.data$.batchvar))
- ),
- model = purrr::map(
- .x = .data$data,
- .f = ~ stats::lm(
- formula = stats::as.formula(paste0(
- "value ~ .batchvar +",
- paste(
- confounders,
- collapse = " + ",
- sep = " + "
- )
- )),
- data = .x
- )
- ),
- .batchvar = purrr::map(
- .x = .data$data,
- .f = ~ .x %>%
- dplyr::pull(.data$.batchvar) %>%
- levels()
- )
- )
-
- values <- res %>%
- tidyr::unnest(cols = .data$.batchvar) %>%
- dplyr::mutate(
- data = purrr::map2(
- .x = .data$data,
- .y = .data$.batchvar,
- .f = ~ .x %>% dplyr::mutate(.batchvar = .y)
- ),
- pred = purrr::map2(.x = .data$model, .y = .data$data, .f = stats::predict)
- ) %>%
- dplyr::select(.data$marker, .data$.batchvar, .data$pred) %>%
- tidyr::unnest(cols = .data$pred) %>%
- dplyr::group_by(.data$marker, .data$.batchvar) %>%
- dplyr::summarize(batchmean = mean(.data$pred, na.rm = TRUE)) %>%
- dplyr::group_by(.data$marker) %>%
- dplyr::mutate(markermean = mean(.data$batchmean)) %>%
- dplyr::ungroup() %>%
- dplyr::transmute(
- marker = .data$marker,
- .batchvar = .data$.batchvar,
- batchmean = .data$batchmean - .data$markermean
- )
- return(list(list(
- models = res %>% dplyr::ungroup() %>% dplyr::pull("model"),
- values = values
- )))
-}
-
-#' Batch means for approach 4: IPW
-#'
-#' @param data Data set
-#' @param markers Variables to batch-adjust
-#' @param confounders Confounders: features that differ
-#' @param truncate Lower and upper extreme quantiles to
-#' truncate stabilized weights at. Defaults to c(0.025, 0.975).
-#'
-#' @return Tibble of batch means per batch and marker
-#' @noRd
-batchmean_ipw <- function(
- data,
- markers,
- confounders,
- truncate = c(0.025, 0.975)
-) {
- ipwbatch <- function(data, variable, confounders, truncate) {
- data <- data %>%
- dplyr::rename(variable = dplyr::one_of(variable)) %>%
- dplyr::filter(!is.na(.data$variable)) %>%
- dplyr::mutate(.batchvar = factor_drop(.data$.batchvar))
-
- res <- data %>%
- tidyr::nest(data = dplyr::everything()) %>%
- dplyr::mutate(
- num = purrr::map(
- .x = .data$data,
- .f = ~ nnet::multinom(
- formula = .batchvar ~ 1,
- data = .x,
- trace = FALSE
- )
- ),
- den = purrr::map(
- .x = .data$data,
- .f = ~ nnet::multinom(
- formula = stats::as.formula(
- paste(".batchvar ~", confounders)
- ),
- data = .x,
- trace = FALSE
- )
- )
- )
-
- values <- res %>%
- dplyr::mutate_at(
- .vars = dplyr::vars(.data$num, .data$den),
- .funs = ~ purrr::map(.x = ., .f = stats::predict, type = "probs") %>%
- purrr::map(.x = ., .f = tibble::as_tibble) %>%
- purrr::map2(
- .x = .,
- .y = .data$data,
- .f = ~ .x %>%
- dplyr::mutate(.batchvar = .y %>%
- purrr::pluck(".batchvar"))
- )
- )
-
- # multinom()$fitted.values is just a vector of probabilities for
- # the 2nd outcome level if there are only two levels
- if (length(levels(factor(data$.batchvar))) == 2) {
- values <- values %>%
- dplyr::mutate_at(
- .vars = dplyr::vars(.data$num, .data$den),
- .funs = ~ purrr::map(.x = ., .f = ~ .x %>%
- dplyr::mutate(
- probs = dplyr::if_else(
- .data$.batchvar ==
- levels(factor(.data$.batchvar))[1],
- true = 1 - .data$value,
- false = .data$value
- )
- ) %>%
- dplyr::pull(.data$probs))
- )
- # otherwise probabilities are a data frame
- } else {
- values <- values %>%
- dplyr::mutate_at(
- .vars = dplyr::vars(.data$num, .data$den),
- .funs =
- ~ purrr::map(
- .x = .,
- .f = ~ .x %>%
- tidyr::pivot_longer(
- -.data$.batchvar,
- names_to = "batch",
- values_to = "prob"
- ) %>%
- dplyr::filter(.data$batch == .data$.batchvar) %>%
- dplyr::pull(.data$prob)
- )
- )
- }
-
- values <- values %>%
- tidyr::unnest(cols = c(.data$data, .data$num, .data$den)) %>%
- dplyr::mutate(
- sw = .data$num / .data$den,
- trunc = dplyr::case_when(
- .data$sw < stats::quantile(.data$sw, truncate[1]) ~
- stats::quantile(.data$sw, truncate[1]),
- .data$sw > stats::quantile(.data$sw, truncate[2]) ~
- stats::quantile(.data$sw, truncate[2]),
- TRUE ~ .data$sw
- )
- )
-
- xlev <- unique(data %>% dplyr::pull(.data$.batchvar))
-
- values <- geepack::geeglm(
- formula = variable ~ .batchvar,
- data = values,
- weights = values$trunc,
- id = values$.id,
- corstr = "independence"
- ) %>%
- broom::tidy() %>%
- dplyr::filter(!stringr::str_detect(
- string = .data$term,
- pattern = "(Intercept)"
- )) %>%
- dplyr::mutate(term = as.character(
- stringr::str_remove_all(
- string = .data$term,
- pattern = ".batchvar"
- )
- )) %>%
- dplyr::full_join(tibble::tibble(term = as.character(xlev)), by = "term") %>%
- dplyr::mutate(
- estimate = dplyr::if_else(
- is.na(.data$estimate),
- true = 0,
- false = .data$estimate
- ),
- estimate = .data$estimate - mean(.data$estimate),
- marker = variable,
- term = .data$term
- ) %>%
- dplyr::arrange(.data$term) %>%
- dplyr::select(.data$marker, .batchvar = .data$term, batchmean = .data$estimate)
- list(values = values, models = res %>% dplyr::pull(.data$den))
- }
-
- purrr::map(
- .x = data %>% dplyr::select({{ markers }}) %>% names(),
- .f = ipwbatch,
- data = data %>%
- dplyr::filter(dplyr::across(dplyr::all_of(confounders), ~ !is.na(.x))),
- truncate = truncate,
- confounders = paste(confounders, sep = " + ", collapse = " + ")
- )
-}
-
-
-#' Quantiles for approach 5: Quantile regression
-#'
-#' @param data Data set
-#' @param variable Single variable to batch-adjust
-#' @param confounders Confounders: features that differ
-#' @param tau Quantiles to use for scaling
-#' @param rq_method Algorithmic method to fit quantile regression.
-#'
-#' @return Tibble of quantiles per batch
-#' @noRd
-batchrq <- function(data, variable, confounders, tau, rq_method) {
- res <- data %>%
- dplyr::rename(variable = {{ variable }}) %>%
- dplyr::filter(!is.na(.data$variable)) %>%
- dplyr::mutate(.batchvar = factor_drop(.data$.batchvar)) %>%
- tidyr::nest(data = dplyr::everything()) %>%
- dplyr::mutate(
- un = purrr::map(
- .x = .data$data,
- .f = ~ quantreg::rq(
- formula = variable ~ .batchvar,
- data = .x,
- tau = tau,
- method = rq_method
- )
- ),
- ad = purrr::map(
- .x = .data$data,
- .f = ~ quantreg::rq(
- formula = stats::reformulate(
- response = "variable",
- termlabels = c(".batchvar", confounders)
- ),
- data = .x,
- tau = tau,
- method = rq_method
- )
- ),
- .batchvar = purrr::map(
- .x = .data$data,
- .f = ~ .x %>%
- dplyr::pull(.data$.batchvar) %>%
- levels()
- )
- )
-
- values <- res %>%
- tidyr::unnest(cols = .data$.batchvar) %>%
- dplyr::mutate(
- data = purrr::map2(
- .x = .data$data,
- .y = .data$.batchvar,
- .f = ~ .x %>% dplyr::mutate(.batchvar = .y)
- ),
- un = purrr::map2(.x = .data$un, .y = .data$data, .f = stats::predict),
- ad = purrr::map2(.x = .data$ad, .y = .data$data, .f = stats::predict),
- un = purrr::map(
- .x = .data$un,
- .f = tibble::as_tibble,
- .name_repair = ~ c("un_lo", "un_hi")
- ),
- ad = purrr::map(
- .x = .data$ad,
- .f = tibble::as_tibble,
- .name_repair = ~ c("ad_lo", "ad_hi")
- ),
- all_lo = purrr::map_dbl(
- .x = .data$data,
- .f = ~ stats::quantile(.x$variable, probs = 0.25)
- ),
- all_hi = purrr::map_dbl(
- .x = .data$data,
- .f = ~ stats::quantile(.x$variable, probs = 0.75)
- ),
- all_iq = .data$all_hi - .data$all_lo
- ) %>%
- dplyr::select(
- .data$.batchvar,
- .data$un,
- .data$ad,
- .data$all_lo,
- .data$all_hi,
- .data$all_iq
- ) %>%
- tidyr::unnest(cols = c(.data$un, .data$ad)) %>%
- dplyr::group_by(.data$.batchvar) %>%
- dplyr::summarize(
- un_lo = stats::quantile(.data$un_lo, probs = 0.25),
- ad_lo = stats::quantile(.data$ad_lo, probs = 0.25),
- un_hi = stats::quantile(.data$un_hi, probs = 0.75),
- ad_hi = stats::quantile(.data$ad_hi, probs = 0.75),
- all_lo = stats::median(.data$all_lo),
- all_iq = stats::median(.data$all_iq)
- ) %>%
- dplyr::mutate(
- un_iq = .data$un_hi - .data$un_lo,
- ad_iq = .data$ad_hi - .data$ad_lo,
- marker = {{ variable }}
- )
-
- models <- res %>% dplyr::pull(.data$ad)
- return(tibble::lst(values, models))
-}
-
-#
-#' Helper function for approach 6: Quantile normalize
-#'
-#' @param var Single variable to quantile-normalize
-#' @param batch Variable indicating batch
-#'
-#' @return Tibble of means per batch for one variable
-#' @noRd
-batch_quantnorm <- function(var, batch) {
- tibble::tibble(var, batch) %>%
- tibble::rowid_to_column() %>%
- tidyr::pivot_wider(names_from = batch, values_from = var) %>%
- dplyr::select(-.data$rowid) %>%
- dplyr::select_if(~ !all(is.na(.))) %>%
- as.matrix() %>%
- limma::normalizeQuantiles() %>%
- tibble::as_tibble() %>%
- dplyr::transmute(
- result = purrr::pmap_dbl(
- .l = .,
- .f = function(...) {
- mean(c(...), na.rm = TRUE)
- }
- ),
- result = dplyr::if_else(
- is.nan(.data$result),
- true = NA_real_,
- false = .data$result
- )
- ) %>%
- dplyr::pull(.data$result)
-}
-
#' Adjust for batch effects
#'
#' @description
diff --git a/R/batch_quantnorm.R b/R/batch_quantnorm.R
new file mode 100644
index 0000000..aa8c1bf
--- /dev/null
+++ b/R/batch_quantnorm.R
@@ -0,0 +1,31 @@
+#' Helper function for approach 6: Quantile normalize
+#'
+#' @param var Single variable to quantile-normalize
+#' @param batch Variable indicating batch
+#'
+#' @return Tibble of means per batch for one variable
+#' @noRd
+batch_quantnorm <- function(var, batch) {
+ tibble::tibble(var, batch) %>%
+ tibble::rowid_to_column() %>%
+ tidyr::pivot_wider(names_from = batch, values_from = var) %>%
+ dplyr::select(-.data$rowid) %>%
+ dplyr::select_if(~ !all(is.na(.))) %>%
+ as.matrix() %>%
+ limma::normalizeQuantiles() %>%
+ tibble::as_tibble() %>%
+ dplyr::transmute(
+ result = purrr::pmap_dbl(
+ .l = .,
+ .f = function(...) {
+ mean(c(...), na.rm = TRUE)
+ }
+ ),
+ result = dplyr::if_else(
+ is.nan(.data$result),
+ true = NA_real_,
+ false = .data$result
+ )
+ ) %>%
+ dplyr::pull(.data$result)
+}
diff --git a/R/batch_rq.R b/R/batch_rq.R
new file mode 100644
index 0000000..f847946
--- /dev/null
+++ b/R/batch_rq.R
@@ -0,0 +1,103 @@
+#' Quantiles for approach 5: Quantile regression
+#'
+#' @param data Data set
+#' @param variable Single variable to batch-adjust
+#' @param confounders Confounders: features that differ
+#' @param tau Quantiles to use for scaling
+#' @param rq_method Algorithmic method to fit quantile regression.
+#'
+#' @return Tibble of quantiles per batch
+#' @noRd
+batchrq <- function(data, variable, confounders, tau, rq_method) {
+ res <- data %>%
+ dplyr::rename(variable = {{ variable }}) %>%
+ dplyr::filter(!is.na(.data$variable)) %>%
+ dplyr::mutate(.batchvar = factor_drop(.data$.batchvar)) %>%
+ tidyr::nest(data = dplyr::everything()) %>%
+ dplyr::mutate(
+ un = purrr::map(
+ .x = .data$data,
+ .f = ~ quantreg::rq(
+ formula = variable ~ .batchvar,
+ data = .x,
+ tau = tau,
+ method = rq_method
+ )
+ ),
+ ad = purrr::map(
+ .x = .data$data,
+ .f = ~ quantreg::rq(
+ formula = stats::reformulate(
+ response = "variable",
+ termlabels = c(".batchvar", confounders)
+ ),
+ data = .x,
+ tau = tau,
+ method = rq_method
+ )
+ ),
+ .batchvar = purrr::map(
+ .x = .data$data,
+ .f = ~ .x %>%
+ dplyr::pull(.data$.batchvar) %>%
+ levels()
+ )
+ )
+
+ values <- res %>%
+ tidyr::unnest(cols = .data$.batchvar) %>%
+ dplyr::mutate(
+ data = purrr::map2(
+ .x = .data$data,
+ .y = .data$.batchvar,
+ .f = ~ .x %>% dplyr::mutate(.batchvar = .y)
+ ),
+ un = purrr::map2(.x = .data$un, .y = .data$data, .f = stats::predict),
+ ad = purrr::map2(.x = .data$ad, .y = .data$data, .f = stats::predict),
+ un = purrr::map(
+ .x = .data$un,
+ .f = tibble::as_tibble,
+ .name_repair = ~ c("un_lo", "un_hi")
+ ),
+ ad = purrr::map(
+ .x = .data$ad,
+ .f = tibble::as_tibble,
+ .name_repair = ~ c("ad_lo", "ad_hi")
+ ),
+ all_lo = purrr::map_dbl(
+ .x = .data$data,
+ .f = ~ stats::quantile(.x$variable, probs = 0.25)
+ ),
+ all_hi = purrr::map_dbl(
+ .x = .data$data,
+ .f = ~ stats::quantile(.x$variable, probs = 0.75)
+ ),
+ all_iq = .data$all_hi - .data$all_lo
+ ) %>%
+ dplyr::select(
+ .data$.batchvar,
+ .data$un,
+ .data$ad,
+ .data$all_lo,
+ .data$all_hi,
+ .data$all_iq
+ ) %>%
+ tidyr::unnest(cols = c(.data$un, .data$ad)) %>%
+ dplyr::group_by(.data$.batchvar) %>%
+ dplyr::summarize(
+ un_lo = stats::quantile(.data$un_lo, probs = 0.25),
+ ad_lo = stats::quantile(.data$ad_lo, probs = 0.25),
+ un_hi = stats::quantile(.data$un_hi, probs = 0.75),
+ ad_hi = stats::quantile(.data$ad_hi, probs = 0.75),
+ all_lo = stats::median(.data$all_lo),
+ all_iq = stats::median(.data$all_iq)
+ ) %>%
+ dplyr::mutate(
+ un_iq = .data$un_hi - .data$un_lo,
+ ad_iq = .data$ad_hi - .data$ad_lo,
+ marker = {{ variable }}
+ )
+
+ models <- res %>% dplyr::pull(.data$ad)
+ return(tibble::lst(values, models))
+}
diff --git a/R/batchmean_ipw.R b/R/batchmean_ipw.R
new file mode 100644
index 0000000..e0d8216
--- /dev/null
+++ b/R/batchmean_ipw.R
@@ -0,0 +1,154 @@
+#' Batch means for approach 4: IPW
+#'
+#' @param data Data set
+#' @param markers Variables to batch-adjust
+#' @param confounders Confounders: features that differ
+#' @param truncate Lower and upper extreme quantiles to
+#' truncate stabilized weights at. Defaults to c(0.025, 0.975).
+#'
+#' @return Tibble of batch means per batch and marker
+#' @noRd
+batchmean_ipw <- function(
+ data,
+ markers,
+ confounders,
+ truncate = c(0.025, 0.975)
+) {
+ ipwbatch <- function(data, variable, confounders, truncate) {
+ data <- data %>%
+ dplyr::rename(variable = dplyr::one_of(variable)) %>%
+ dplyr::filter(!is.na(.data$variable)) %>%
+ dplyr::mutate(.batchvar = factor_drop(.data$.batchvar))
+
+ res <- data %>%
+ tidyr::nest(data = dplyr::everything()) %>%
+ dplyr::mutate(
+ num = purrr::map(
+ .x = .data$data,
+ .f = ~ nnet::multinom(
+ formula = .batchvar ~ 1,
+ data = .x,
+ trace = FALSE
+ )
+ ),
+ den = purrr::map(
+ .x = .data$data,
+ .f = ~ nnet::multinom(
+ formula = stats::as.formula(
+ paste(".batchvar ~", confounders)
+ ),
+ data = .x,
+ trace = FALSE
+ )
+ )
+ )
+
+ values <- res %>%
+ dplyr::mutate_at(
+ .vars = dplyr::vars(.data$num, .data$den),
+ .funs = ~ purrr::map(.x = ., .f = stats::predict, type = "probs") %>%
+ purrr::map(.x = ., .f = tibble::as_tibble) %>%
+ purrr::map2(
+ .x = .,
+ .y = .data$data,
+ .f = ~ .x %>%
+ dplyr::mutate(.batchvar = .y %>%
+ purrr::pluck(".batchvar"))
+ )
+ )
+
+ # multinom()$fitted.values is just a vector of probabilities for
+ # the 2nd outcome level if there are only two levels
+ if (length(levels(factor(data$.batchvar))) == 2) {
+ values <- values %>%
+ dplyr::mutate_at(
+ .vars = dplyr::vars(.data$num, .data$den),
+ .funs = ~ purrr::map(.x = ., .f = ~ .x %>%
+ dplyr::mutate(
+ probs = dplyr::if_else(
+ .data$.batchvar ==
+ levels(factor(.data$.batchvar))[1],
+ true = 1 - .data$value,
+ false = .data$value
+ )
+ ) %>%
+ dplyr::pull(.data$probs))
+ )
+ # otherwise probabilities are a data frame
+ } else {
+ values <- values %>%
+ dplyr::mutate_at(
+ .vars = dplyr::vars(.data$num, .data$den),
+ .funs =
+ ~ purrr::map(
+ .x = .,
+ .f = ~ .x %>%
+ tidyr::pivot_longer(
+ -.data$.batchvar,
+ names_to = "batch",
+ values_to = "prob"
+ ) %>%
+ dplyr::filter(.data$batch == .data$.batchvar) %>%
+ dplyr::pull(.data$prob)
+ )
+ )
+ }
+
+ values <- values %>%
+ tidyr::unnest(cols = c(.data$data, .data$num, .data$den)) %>%
+ dplyr::mutate(
+ sw = .data$num / .data$den,
+ trunc = dplyr::case_when(
+ .data$sw < stats::quantile(.data$sw, truncate[1]) ~
+ stats::quantile(.data$sw, truncate[1]),
+ .data$sw > stats::quantile(.data$sw, truncate[2]) ~
+ stats::quantile(.data$sw, truncate[2]),
+ TRUE ~ .data$sw
+ )
+ )
+
+ xlev <- unique(data %>% dplyr::pull(.data$.batchvar))
+
+ values <- geepack::geeglm(
+ formula = variable ~ .batchvar,
+ data = values,
+ weights = values$trunc,
+ id = values$.id,
+ corstr = "independence"
+ ) %>%
+ broom::tidy() %>%
+ dplyr::filter(!stringr::str_detect(
+ string = .data$term,
+ pattern = "(Intercept)"
+ )) %>%
+ dplyr::mutate(term = as.character(
+ stringr::str_remove_all(
+ string = .data$term,
+ pattern = ".batchvar"
+ )
+ )) %>%
+ dplyr::full_join(tibble::tibble(term = as.character(xlev)), by = "term") %>%
+ dplyr::mutate(
+ estimate = dplyr::if_else(
+ is.na(.data$estimate),
+ true = 0,
+ false = .data$estimate
+ ),
+ estimate = .data$estimate - mean(.data$estimate),
+ marker = variable,
+ term = .data$term
+ ) %>%
+ dplyr::arrange(.data$term) %>%
+ dplyr::select(.data$marker, .batchvar = .data$term, batchmean = .data$estimate)
+ list(values = values, models = res %>% dplyr::pull(.data$den))
+ }
+
+ purrr::map(
+ .x = data %>% dplyr::select({{ markers }}) %>% names(),
+ .f = ipwbatch,
+ data = data %>%
+ dplyr::filter(dplyr::across(dplyr::all_of(confounders), ~ !is.na(.x))),
+ truncate = truncate,
+ confounders = paste(confounders, sep = " + ", collapse = " + ")
+ )
+}
diff --git a/R/batchmean_simple.R b/R/batchmean_simple.R
new file mode 100644
index 0000000..ff8af67
--- /dev/null
+++ b/R/batchmean_simple.R
@@ -0,0 +1,29 @@
+#' Batch means for approach 2: Unadjusted means
+#'
+#' @param data Data set
+#' @param markers Biomarkers to adjust
+#'
+#' @importFrom magrittr %>%
+#' @importFrom rlang .data
+#' @return Tibble of means per marker and batch
+#' @noRd
+batchmean_simple <- function(data, markers) {
+ values <- data %>%
+ dplyr::select(.data$.id, .data$.batchvar, {{ markers }}) %>%
+ dplyr::group_by(.data$.batchvar) %>%
+ dplyr::summarize_at(
+ .vars = dplyr::vars(-.data$.id),
+ .funs = mean,
+ na.rm = TRUE
+ ) %>%
+ dplyr::mutate_at(
+ .vars = dplyr::vars(-.data$.batchvar),
+ .funs = ~ . - mean(., na.rm = TRUE)
+ ) %>%
+ tidyr::pivot_longer(
+ col = c(-.data$.batchvar),
+ names_to = "marker",
+ values_to = "batchmean"
+ )
+ return(list(list(values = values, models = NULL)))
+}
diff --git a/R/batchmean_standardize.R b/R/batchmean_standardize.R
new file mode 100644
index 0000000..208b186
--- /dev/null
+++ b/R/batchmean_standardize.R
@@ -0,0 +1,74 @@
+#' Batch means for approach 3: Marginal standardization
+#'
+#' @param data Data set
+#' @param markers Vector of variables to batch-adjust
+#' @param confounders Confounders: features that differ
+#' between batches that should be retained
+#'
+#' @return Tibble with conditional means per marker and batch
+#' @noRd
+batchmean_standardize <- function(data, markers, confounders) {
+ res <- data %>%
+ tidyr::pivot_longer(
+ cols = {{ markers }},
+ names_to = "marker",
+ values_to = "value"
+ ) %>%
+ dplyr::filter(!is.na(.data$value)) %>%
+ dplyr::group_by(.data$marker) %>%
+ tidyr::nest(data = c(-.data$marker)) %>%
+ dplyr::mutate(
+ data = purrr::map(
+ .x = .data$data,
+ .f = ~ .x %>%
+ dplyr::mutate(.batchvar = factor_drop(.data$.batchvar))
+ ),
+ model = purrr::map(
+ .x = .data$data,
+ .f = ~ stats::lm(
+ formula = stats::as.formula(paste0(
+ "value ~ .batchvar +",
+ paste(
+ confounders,
+ collapse = " + ",
+ sep = " + "
+ )
+ )),
+ data = .x
+ )
+ ),
+ .batchvar = purrr::map(
+ .x = .data$data,
+ .f = ~ .x %>%
+ dplyr::pull(.data$.batchvar) %>%
+ levels()
+ )
+ )
+
+ values <- res %>%
+ tidyr::unnest(cols = .data$.batchvar) %>%
+ dplyr::mutate(
+ data = purrr::map2(
+ .x = .data$data,
+ .y = .data$.batchvar,
+ .f = ~ .x %>% dplyr::mutate(.batchvar = .y)
+ ),
+ pred = purrr::map2(.x = .data$model, .y = .data$data, .f = stats::predict)
+ ) %>%
+ dplyr::select(.data$marker, .data$.batchvar, .data$pred) %>%
+ tidyr::unnest(cols = .data$pred) %>%
+ dplyr::group_by(.data$marker, .data$.batchvar) %>%
+ dplyr::summarize(batchmean = mean(.data$pred, na.rm = TRUE)) %>%
+ dplyr::group_by(.data$marker) %>%
+ dplyr::mutate(markermean = mean(.data$batchmean)) %>%
+ dplyr::ungroup() %>%
+ dplyr::transmute(
+ marker = .data$marker,
+ .batchvar = .data$.batchvar,
+ batchmean = .data$batchmean - .data$markermean
+ )
+ return(list(list(
+ models = res %>% dplyr::ungroup() %>% dplyr::pull("model"),
+ values = values
+ )))
+}
diff --git a/R/factor_drop.R b/R/factor_drop.R
new file mode 100644
index 0000000..d663fb3
--- /dev/null
+++ b/R/factor_drop.R
@@ -0,0 +1,14 @@
+#' Drop empty factor levels
+#'
+#' @description
+#' Avoids predict() issues. This is \code{forcats::fct_drop()}
+#' without bells and whistles.
+#'
+#' @param f factor
+#'
+#' @return Factor
+#' @noRd
+factor_drop <- function(f) {
+ factor_levels <- levels(f)
+ factor(f, levels = setdiff(factor_levels, factor_levels[table(f) == 0]))
+}
From 801698146e6dd46df04ac9a0becad15b071217fa Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 14:09:36 +0100
Subject: [PATCH 06/36] let air reformat code
---
R/adjust_batch.R | 80 +++++++++++++++++++++++++++++------------
R/batchmean_ipw.R | 90 +++++++++++++++++++++++++++--------------------
R/plot_batch.R | 24 +++++++++----
3 files changed, 127 insertions(+), 67 deletions(-)
diff --git a/R/adjust_batch.R b/R/adjust_batch.R
index 1b07e0e..ff32bd8 100644
--- a/R/adjust_batch.R
+++ b/R/adjust_batch.R
@@ -1,6 +1,8 @@
# as per https://github.com/jennybc/googlesheets/blob/master/R/googlesheets.R:
## quiets concerns of R CMD check re: the .'s that appear in pipelines
-if (getRversion() >= "2.15.1") utils::globalVariables(c("."))
+if (getRversion() >= "2.15.1") {
+ utils::globalVariables(c("."))
+}
#' Adjust for batch effects
#'
@@ -208,10 +210,13 @@ adjust_batch <- function(
))
}
- if (method %in% c("simple", "quantnorm") &
- data %>%
- dplyr::select({{ confounders }}) %>%
- ncol() > 0) {
+ if (
+ method %in% c("simple", "quantnorm") &
+ data %>%
+ dplyr::select({{ confounders }}) %>%
+ ncol() >
+ 0
+ ) {
message(paste0(
"Batch effect correction via 'method = ",
method,
@@ -226,10 +231,14 @@ adjust_batch <- function(
# Mean-based methods
if (method %in% c("simple", "standardize", "ipw")) {
- if (method %in% c("standardize", "ipw") &
- data %>%
- dplyr::select({{ confounders }}) %>%
- ncol() == 0) {
+ if (
+ method %in%
+ c("standardize", "ipw") &
+ data %>%
+ dplyr::select({{ confounders }}) %>%
+ ncol() ==
+ 0
+ ) {
message(paste0(
"Batch effect correction via 'method = ",
method,
@@ -254,7 +263,10 @@ adjust_batch <- function(
truncate = ipw_truncate
)
)
- adjust_parameters <- purrr::map_dfr(.x = res, .f = ~ purrr::pluck(.x, "values"))
+ adjust_parameters <- purrr::map_dfr(
+ .x = res,
+ .f = ~ purrr::pluck(.x, "values")
+ )
method_indices <- c("simple" = 2, "standardize" = 3, "ipw" = 4)
if (suffix == "_adjX") {
suffix <- paste0("_adj", method_indices[method[1]])
@@ -281,7 +293,10 @@ adjust_batch <- function(
.x = data %>% dplyr::select({{ markers }}) %>% names(),
.f = batchrq,
data = data %>%
- dplyr::filter(dplyr::across(dplyr::all_of({{ confounders }}), ~ !is.na(.x))),
+ dplyr::filter(dplyr::across(
+ dplyr::all_of({{ confounders }}),
+ ~ !is.na(.x)
+ )),
confounders = dplyr::if_else(
dplyr::enexpr(confounders) != "",
true = paste0(
@@ -297,7 +312,10 @@ adjust_batch <- function(
tau = quantreg_tau,
rq_method = quantreg_method
)
- adjust_parameters <- purrr::map_dfr(.x = res, .f = ~ purrr::pluck(.x, "values"))
+ adjust_parameters <- purrr::map_dfr(
+ .x = res,
+ .f = ~ purrr::pluck(.x, "values")
+ )
if (suffix == "_adjX") {
suffix <- "_adj5"
}
@@ -316,8 +334,12 @@ adjust_batch <- function(
dplyr::group_by(.data$marker) %>%
dplyr::mutate(
value_adjusted = (.data$value - .data$un_lo) /
- .data$un_iq * .data$all_iq * (.data$un_iq / .data$ad_iq) +
- .data$all_lo - .data$ad_lo + .data$un_lo,
+ .data$un_iq *
+ .data$all_iq *
+ (.data$un_iq / .data$ad_iq) +
+ .data$all_lo -
+ .data$ad_lo +
+ .data$un_lo,
marker = paste0(.data$marker, suffix)
) %>%
dplyr::select(
@@ -348,15 +370,20 @@ adjust_batch <- function(
) %>%
dplyr::mutate(marker = paste0(.data$marker, suffix)) %>%
dplyr::group_by(.data$marker) %>%
- dplyr::mutate(value_adjusted = batch_quantnorm(
- var = .data$value,
- batch = .data$.batchvar
- )) %>%
+ dplyr::mutate(
+ value_adjusted = batch_quantnorm(
+ var = .data$value,
+ batch = .data$.batchvar
+ )
+ ) %>%
dplyr::ungroup() %>%
dplyr::select(-.data$value)
res <- list(list(res = NULL, models = NULL))
- adjust_parameters <- tibble::tibble(marker = data %>%
- dplyr::select({{ markers }}) %>% names())
+ adjust_parameters <- tibble::tibble(
+ marker = data %>%
+ dplyr::select({{ markers }}) %>%
+ names()
+ )
}
# Dataset to return
@@ -372,12 +399,19 @@ adjust_batch <- function(
# Meta-data to return as attribute
attr_list <- list(
adjust_method = method,
- markers = data_orig %>% dplyr::select({{ markers }}) %>% names(),
+ markers = data_orig %>%
+ dplyr::select({{ markers }}) %>%
+ names(),
suffix = suffix,
- batchvar = data_orig %>% dplyr::select({{ batch }}) %>% names(),
+ batchvar = data_orig %>%
+ dplyr::select({{ batch }}) %>%
+ names(),
confounders = dplyr::enexpr(confounders),
adjust_parameters = adjust_parameters,
- model_fits = purrr::map(.x = res, .f = ~ purrr::pluck(.x, "models"))
+ model_fits = purrr::map(
+ .x = res,
+ .f = ~ purrr::pluck(.x, "models")
+ )
)
class(attr_list) <- c("batchtma", class(res))
attr(values, which = ".batchtma") <- attr_list
diff --git a/R/batchmean_ipw.R b/R/batchmean_ipw.R
index e0d8216..141325f 100644
--- a/R/batchmean_ipw.R
+++ b/R/batchmean_ipw.R
@@ -9,10 +9,10 @@
#' @return Tibble of batch means per batch and marker
#' @noRd
batchmean_ipw <- function(
- data,
- markers,
- confounders,
- truncate = c(0.025, 0.975)
+ data,
+ markers,
+ confounders,
+ truncate = c(0.025, 0.975)
) {
ipwbatch <- function(data, variable, confounders, truncate) {
data <- data %>%
@@ -52,8 +52,10 @@ batchmean_ipw <- function(
.x = .,
.y = .data$data,
.f = ~ .x %>%
- dplyr::mutate(.batchvar = .y %>%
- purrr::pluck(".batchvar"))
+ dplyr::mutate(
+ .batchvar = .y %>%
+ purrr::pluck(".batchvar")
+ )
)
)
@@ -63,34 +65,35 @@ batchmean_ipw <- function(
values <- values %>%
dplyr::mutate_at(
.vars = dplyr::vars(.data$num, .data$den),
- .funs = ~ purrr::map(.x = ., .f = ~ .x %>%
- dplyr::mutate(
- probs = dplyr::if_else(
- .data$.batchvar ==
- levels(factor(.data$.batchvar))[1],
- true = 1 - .data$value,
- false = .data$value
- )
- ) %>%
- dplyr::pull(.data$probs))
+ .funs = ~ purrr::map(
+ .x = .,
+ .f = ~ .x %>%
+ dplyr::mutate(
+ probs = dplyr::if_else(
+ .data$.batchvar == levels(factor(.data$.batchvar))[1],
+ true = 1 - .data$value,
+ false = .data$value
+ )
+ ) %>%
+ dplyr::pull(.data$probs)
+ )
)
# otherwise probabilities are a data frame
} else {
values <- values %>%
dplyr::mutate_at(
.vars = dplyr::vars(.data$num, .data$den),
- .funs =
- ~ purrr::map(
- .x = .,
- .f = ~ .x %>%
- tidyr::pivot_longer(
- -.data$.batchvar,
- names_to = "batch",
- values_to = "prob"
- ) %>%
- dplyr::filter(.data$batch == .data$.batchvar) %>%
- dplyr::pull(.data$prob)
- )
+ .funs = ~ purrr::map(
+ .x = .,
+ .f = ~ .x %>%
+ tidyr::pivot_longer(
+ -.data$.batchvar,
+ names_to = "batch",
+ values_to = "prob"
+ ) %>%
+ dplyr::filter(.data$batch == .data$.batchvar) %>%
+ dplyr::pull(.data$prob)
+ )
)
}
@@ -117,17 +120,24 @@ batchmean_ipw <- function(
corstr = "independence"
) %>%
broom::tidy() %>%
- dplyr::filter(!stringr::str_detect(
- string = .data$term,
- pattern = "(Intercept)"
- )) %>%
- dplyr::mutate(term = as.character(
- stringr::str_remove_all(
+ dplyr::filter(
+ !stringr::str_detect(
string = .data$term,
- pattern = ".batchvar"
+ pattern = "(Intercept)"
)
- )) %>%
- dplyr::full_join(tibble::tibble(term = as.character(xlev)), by = "term") %>%
+ ) %>%
+ dplyr::mutate(
+ term = as.character(
+ stringr::str_remove_all(
+ string = .data$term,
+ pattern = ".batchvar"
+ )
+ )
+ ) %>%
+ dplyr::full_join(
+ tibble::tibble(term = as.character(xlev)),
+ by = "term"
+ ) %>%
dplyr::mutate(
estimate = dplyr::if_else(
is.na(.data$estimate),
@@ -139,7 +149,11 @@ batchmean_ipw <- function(
term = .data$term
) %>%
dplyr::arrange(.data$term) %>%
- dplyr::select(.data$marker, .batchvar = .data$term, batchmean = .data$estimate)
+ dplyr::select(
+ .data$marker,
+ .batchvar = .data$term,
+ batchmean = .data$estimate
+ )
list(values = values, models = res %>% dplyr::pull(.data$den))
}
diff --git a/R/plot_batch.R b/R/plot_batch.R
index 2bb6b97..0aa1128 100644
--- a/R/plot_batch.R
+++ b/R/plot_batch.R
@@ -73,8 +73,15 @@
#' color = confounder
#' ) +
#' ggplot2::labs(y = "Biomarker (variable 'noisy')")
-plot_batch <- function(data, marker, batch, color = NULL, maxlevels = 15, title = NULL, ...) {
-
+plot_batch <- function(
+ data,
+ marker,
+ batch,
+ color = NULL,
+ maxlevels = 15,
+ title = NULL,
+ ...
+) {
# Set levels to number of discrete entries for `color` or 101 if NULL.
# If distinct count of `color` is only 1, that indicates NULL was passed.
nlevels <- dplyr::n_distinct(dplyr::select(data, {{ color }}))
@@ -115,15 +122,20 @@ plot_batch <- function(data, marker, batch, color = NULL, maxlevels = 15, title
mapping = ggplot2::aes(color = {{ color }}, shape = {{ color }})
) +
ggplot2::scale_shape_manual(name = dplyr::enexpr(color), values = 15:30) +
- ggplot2::scale_color_viridis_d(name = dplyr::enexpr(color), option = "cividis")
- }
- else {
+ ggplot2::scale_color_viridis_d(
+ name = dplyr::enexpr(color),
+ option = "cividis"
+ )
+ } else {
myplot +
ggplot2::geom_jitter(
width = 0.2,
height = 0,
mapping = ggplot2::aes(color = {{ color }})
) +
- ggplot2::scale_color_viridis_c(name = dplyr::enexpr(color), option = "cividis")
+ ggplot2::scale_color_viridis_c(
+ name = dplyr::enexpr(color),
+ option = "cividis"
+ )
}
}
From 814cc01cabca7d81c5507070bc2e6683ec6a9302 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 19:55:03 +0100
Subject: [PATCH 07/36] switch to native pipe
---
DESCRIPTION | 2 +
NAMESPACE | 1 -
R/adjust_batch.R | 80 +++++++++++++++++++--------------------
R/batch_quantnorm.R | 18 ++++-----
R/batch_rq.R | 32 ++++++++--------
R/batchmean_ipw.R | 60 ++++++++++++++---------------
R/batchmean_simple.R | 11 +++---
R/batchmean_standardize.R | 40 ++++++++++----------
R/batchtma.R | 2 +-
9 files changed, 123 insertions(+), 123 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index fbffdf3..a9d369c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -21,6 +21,8 @@ Encoding: UTF-8
Roxygen: list(markdown = TRUE)
RoxygenNote: 7.2.3
biocViews:
+Depends:
+ R (>= 4.1.0)
Imports:
broom (>= 0.7.0),
dplyr (>= 1.0.0),
diff --git a/NAMESPACE b/NAMESPACE
index 38231f7..078d0a7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,6 +5,5 @@ export(adjust_batch)
export(diagnose_models)
export(plot_batch)
importFrom(ggplot2,labs)
-importFrom(magrittr,"%>%")
importFrom(rlang,":=")
importFrom(rlang,.data)
diff --git a/R/adjust_batch.R b/R/adjust_batch.R
index ff32bd8..9e58007 100644
--- a/R/adjust_batch.R
+++ b/R/adjust_batch.R
@@ -185,17 +185,17 @@ adjust_batch <- function(
) {
method <- as.character(dplyr::enexpr(method))
allmethods <- c("simple", "standardize", "ipw", "quantreg", "quantnorm")
- data_orig <- data %>%
+ data_orig <- data |>
dplyr::mutate(.id = dplyr::row_number())
- data <- data_orig %>%
- dplyr::rename(.batchvar = {{ batch }}) %>%
+ data <- data_orig |>
+ dplyr::rename(.batchvar = {{ batch }}) |>
dplyr::mutate(
.batchvar = factor(.data$.batchvar),
.batchvar = factor_drop(.data$.batchvar)
- ) %>%
+ ) |>
dplyr::select(.data$.id, .data$.batchvar, {{ markers }}, {{ confounders }})
- confounders <- data %>%
- dplyr::select({{ confounders }}) %>%
+ confounders <- data |>
+ dplyr::select({{ confounders }}) |>
names()
# Check inputs: method and confounders
@@ -212,8 +212,8 @@ adjust_batch <- function(
if (
method %in% c("simple", "quantnorm") &
- data %>%
- dplyr::select({{ confounders }}) %>%
+ data |>
+ dplyr::select({{ confounders }}) |>
ncol() >
0
) {
@@ -225,7 +225,7 @@ adjust_batch <- function(
paste(dplyr::enexpr(confounders), sep = ", ", collapse = ", "),
"). They will be ignored."
))
- data <- data %>% dplyr::select(-dplyr::any_of({{ confounders }}))
+ data <- data |> dplyr::select(-dplyr::any_of({{ confounders }}))
confounders <- NULL
}
@@ -234,8 +234,8 @@ adjust_batch <- function(
if (
method %in%
c("standardize", "ipw") &
- data %>%
- dplyr::select({{ confounders }}) %>%
+ data |>
+ dplyr::select({{ confounders }}) |>
ncol() ==
0
) {
@@ -272,27 +272,27 @@ adjust_batch <- function(
suffix <- paste0("_adj", method_indices[method[1]])
}
- values <- data %>%
- dplyr::select(-dplyr::any_of({{ confounders }})) %>%
+ values <- data |>
+ dplyr::select(-dplyr::any_of({{ confounders }})) |>
tidyr::pivot_longer(
cols = c(-.data$.id, -.data$.batchvar),
names_to = "marker",
values_to = "value"
- ) %>%
- dplyr::left_join(adjust_parameters, by = c("marker", ".batchvar")) %>%
+ ) |>
+ dplyr::left_join(adjust_parameters, by = c("marker", ".batchvar")) |>
dplyr::mutate(
value_adjusted = .data$value - .data$batchmean,
marker = paste0(.data$marker, suffix)
- ) %>%
+ ) |>
dplyr::select(-.data$batchmean, -.data$value)
}
# Quantile regression
if (method == "quantreg") {
res <- purrr::map(
- .x = data %>% dplyr::select({{ markers }}) %>% names(),
+ .x = data |> dplyr::select({{ markers }}) |> names(),
.f = batchrq,
- data = data %>%
+ data = data |>
dplyr::filter(dplyr::across(
dplyr::all_of({{ confounders }}),
~ !is.na(.x)
@@ -320,7 +320,7 @@ adjust_batch <- function(
suffix <- "_adj5"
}
- values <- data %>%
+ values <- data |>
tidyr::pivot_longer(
cols = c(
-.data$.id,
@@ -329,9 +329,9 @@ adjust_batch <- function(
),
names_to = "marker",
values_to = "value"
- ) %>%
- dplyr::left_join(adjust_parameters, by = c("marker", ".batchvar")) %>%
- dplyr::group_by(.data$marker) %>%
+ ) |>
+ dplyr::left_join(adjust_parameters, by = c("marker", ".batchvar")) |>
+ dplyr::group_by(.data$marker) |>
dplyr::mutate(
value_adjusted = (.data$value - .data$un_lo) /
.data$un_iq *
@@ -341,7 +341,7 @@ adjust_batch <- function(
.data$ad_lo +
.data$un_lo,
marker = paste0(.data$marker, suffix)
- ) %>%
+ ) |>
dplyr::select(
-dplyr::any_of({{ confounders }}),
-.data$value,
@@ -361,50 +361,50 @@ adjust_batch <- function(
if (suffix == "_adjX") {
suffix <- "_adj6"
}
- values <- data %>%
- dplyr::select(-dplyr::any_of({{ confounders }})) %>%
+ values <- data |>
+ dplyr::select(-dplyr::any_of({{ confounders }})) |>
tidyr::pivot_longer(
cols = c(-.data$.id, -.data$.batchvar),
names_to = "marker",
values_to = "value"
- ) %>%
- dplyr::mutate(marker = paste0(.data$marker, suffix)) %>%
- dplyr::group_by(.data$marker) %>%
+ ) |>
+ dplyr::mutate(marker = paste0(.data$marker, suffix)) |>
+ dplyr::group_by(.data$marker) |>
dplyr::mutate(
value_adjusted = batch_quantnorm(
var = .data$value,
batch = .data$.batchvar
)
- ) %>%
- dplyr::ungroup() %>%
+ ) |>
+ dplyr::ungroup() |>
dplyr::select(-.data$value)
res <- list(list(res = NULL, models = NULL))
adjust_parameters <- tibble::tibble(
- marker = data %>%
- dplyr::select({{ markers }}) %>%
+ marker = data |>
+ dplyr::select({{ markers }}) |>
names()
)
}
# Dataset to return
- values <- values %>%
+ values <- values |>
tidyr::pivot_wider(
names_from = .data$marker,
values_from = .data$value_adjusted
- ) %>%
- dplyr::select(-.data$.batchvar) %>%
- dplyr::left_join(x = data_orig, by = ".id") %>%
+ ) |>
+ dplyr::select(-.data$.batchvar) |>
+ dplyr::left_join(x = data_orig, by = ".id") |>
dplyr::select(-.data$.id)
# Meta-data to return as attribute
attr_list <- list(
adjust_method = method,
- markers = data_orig %>%
- dplyr::select({{ markers }}) %>%
+ markers = data_orig |>
+ dplyr::select({{ markers }}) |>
names(),
suffix = suffix,
- batchvar = data_orig %>%
- dplyr::select({{ batch }}) %>%
+ batchvar = data_orig |>
+ dplyr::select({{ batch }}) |>
names(),
confounders = dplyr::enexpr(confounders),
adjust_parameters = adjust_parameters,
diff --git a/R/batch_quantnorm.R b/R/batch_quantnorm.R
index aa8c1bf..481c729 100644
--- a/R/batch_quantnorm.R
+++ b/R/batch_quantnorm.R
@@ -6,14 +6,14 @@
#' @return Tibble of means per batch for one variable
#' @noRd
batch_quantnorm <- function(var, batch) {
- tibble::tibble(var, batch) %>%
- tibble::rowid_to_column() %>%
- tidyr::pivot_wider(names_from = batch, values_from = var) %>%
- dplyr::select(-.data$rowid) %>%
- dplyr::select_if(~ !all(is.na(.))) %>%
- as.matrix() %>%
- limma::normalizeQuantiles() %>%
- tibble::as_tibble() %>%
+ tibble::tibble(var, batch) |>
+ tibble::rowid_to_column() |>
+ tidyr::pivot_wider(names_from = batch, values_from = var) |>
+ dplyr::select(-.data$rowid) |>
+ dplyr::select_if(~ !all(is.na(.))) |>
+ as.matrix() |>
+ limma::normalizeQuantiles() |>
+ tibble::as_tibble() |>
dplyr::transmute(
result = purrr::pmap_dbl(
.l = .,
@@ -26,6 +26,6 @@ batch_quantnorm <- function(var, batch) {
true = NA_real_,
false = .data$result
)
- ) %>%
+ ) |>
dplyr::pull(.data$result)
}
diff --git a/R/batch_rq.R b/R/batch_rq.R
index f847946..a6855b4 100644
--- a/R/batch_rq.R
+++ b/R/batch_rq.R
@@ -9,11 +9,11 @@
#' @return Tibble of quantiles per batch
#' @noRd
batchrq <- function(data, variable, confounders, tau, rq_method) {
- res <- data %>%
- dplyr::rename(variable = {{ variable }}) %>%
- dplyr::filter(!is.na(.data$variable)) %>%
- dplyr::mutate(.batchvar = factor_drop(.data$.batchvar)) %>%
- tidyr::nest(data = dplyr::everything()) %>%
+ res <- data |>
+ dplyr::rename(variable = {{ variable }}) |>
+ dplyr::filter(!is.na(.data$variable)) |>
+ dplyr::mutate(.batchvar = factor_drop(.data$.batchvar)) |>
+ tidyr::nest(data = dplyr::everything()) |>
dplyr::mutate(
un = purrr::map(
.x = .data$data,
@@ -38,19 +38,19 @@ batchrq <- function(data, variable, confounders, tau, rq_method) {
),
.batchvar = purrr::map(
.x = .data$data,
- .f = ~ .x %>%
- dplyr::pull(.data$.batchvar) %>%
+ .f = ~ .x |>
+ dplyr::pull(.data$.batchvar) |>
levels()
)
)
- values <- res %>%
- tidyr::unnest(cols = .data$.batchvar) %>%
+ values <- res |>
+ tidyr::unnest(cols = .data$.batchvar) |>
dplyr::mutate(
data = purrr::map2(
.x = .data$data,
.y = .data$.batchvar,
- .f = ~ .x %>% dplyr::mutate(.batchvar = .y)
+ .f = ~ .x |> dplyr::mutate(.batchvar = .y)
),
un = purrr::map2(.x = .data$un, .y = .data$data, .f = stats::predict),
ad = purrr::map2(.x = .data$ad, .y = .data$data, .f = stats::predict),
@@ -73,7 +73,7 @@ batchrq <- function(data, variable, confounders, tau, rq_method) {
.f = ~ stats::quantile(.x$variable, probs = 0.75)
),
all_iq = .data$all_hi - .data$all_lo
- ) %>%
+ ) |>
dplyr::select(
.data$.batchvar,
.data$un,
@@ -81,9 +81,9 @@ batchrq <- function(data, variable, confounders, tau, rq_method) {
.data$all_lo,
.data$all_hi,
.data$all_iq
- ) %>%
- tidyr::unnest(cols = c(.data$un, .data$ad)) %>%
- dplyr::group_by(.data$.batchvar) %>%
+ ) |>
+ tidyr::unnest(cols = c(.data$un, .data$ad)) |>
+ dplyr::group_by(.data$.batchvar) |>
dplyr::summarize(
un_lo = stats::quantile(.data$un_lo, probs = 0.25),
ad_lo = stats::quantile(.data$ad_lo, probs = 0.25),
@@ -91,13 +91,13 @@ batchrq <- function(data, variable, confounders, tau, rq_method) {
ad_hi = stats::quantile(.data$ad_hi, probs = 0.75),
all_lo = stats::median(.data$all_lo),
all_iq = stats::median(.data$all_iq)
- ) %>%
+ ) |>
dplyr::mutate(
un_iq = .data$un_hi - .data$un_lo,
ad_iq = .data$ad_hi - .data$ad_lo,
marker = {{ variable }}
)
- models <- res %>% dplyr::pull(.data$ad)
+ models <- res |> dplyr::pull(.data$ad)
return(tibble::lst(values, models))
}
diff --git a/R/batchmean_ipw.R b/R/batchmean_ipw.R
index 141325f..3f1c866 100644
--- a/R/batchmean_ipw.R
+++ b/R/batchmean_ipw.R
@@ -15,13 +15,13 @@ batchmean_ipw <- function(
truncate = c(0.025, 0.975)
) {
ipwbatch <- function(data, variable, confounders, truncate) {
- data <- data %>%
- dplyr::rename(variable = dplyr::one_of(variable)) %>%
- dplyr::filter(!is.na(.data$variable)) %>%
+ data <- data |>
+ dplyr::rename(variable = dplyr::one_of(variable)) |>
+ dplyr::filter(!is.na(.data$variable)) |>
dplyr::mutate(.batchvar = factor_drop(.data$.batchvar))
- res <- data %>%
- tidyr::nest(data = dplyr::everything()) %>%
+ res <- data |>
+ tidyr::nest(data = dplyr::everything()) |>
dplyr::mutate(
num = purrr::map(
.x = .data$data,
@@ -43,17 +43,17 @@ batchmean_ipw <- function(
)
)
- values <- res %>%
+ values <- res |>
dplyr::mutate_at(
.vars = dplyr::vars(.data$num, .data$den),
- .funs = ~ purrr::map(.x = ., .f = stats::predict, type = "probs") %>%
- purrr::map(.x = ., .f = tibble::as_tibble) %>%
+ .funs = ~ purrr::map(.x = ., .f = stats::predict, type = "probs") |>
+ purrr::map(.x = ., .f = tibble::as_tibble) |>
purrr::map2(
.x = .,
.y = .data$data,
- .f = ~ .x %>%
+ .f = ~ .x |>
dplyr::mutate(
- .batchvar = .y %>%
+ .batchvar = .y |>
purrr::pluck(".batchvar")
)
)
@@ -62,43 +62,43 @@ batchmean_ipw <- function(
# multinom()$fitted.values is just a vector of probabilities for
# the 2nd outcome level if there are only two levels
if (length(levels(factor(data$.batchvar))) == 2) {
- values <- values %>%
+ values <- values |>
dplyr::mutate_at(
.vars = dplyr::vars(.data$num, .data$den),
.funs = ~ purrr::map(
.x = .,
- .f = ~ .x %>%
+ .f = ~ .x |>
dplyr::mutate(
probs = dplyr::if_else(
.data$.batchvar == levels(factor(.data$.batchvar))[1],
true = 1 - .data$value,
false = .data$value
)
- ) %>%
+ ) |>
dplyr::pull(.data$probs)
)
)
# otherwise probabilities are a data frame
} else {
- values <- values %>%
+ values <- values |>
dplyr::mutate_at(
.vars = dplyr::vars(.data$num, .data$den),
.funs = ~ purrr::map(
.x = .,
- .f = ~ .x %>%
+ .f = ~ .x |>
tidyr::pivot_longer(
-.data$.batchvar,
names_to = "batch",
values_to = "prob"
- ) %>%
- dplyr::filter(.data$batch == .data$.batchvar) %>%
+ ) |>
+ dplyr::filter(.data$batch == .data$.batchvar) |>
dplyr::pull(.data$prob)
)
)
}
- values <- values %>%
- tidyr::unnest(cols = c(.data$data, .data$num, .data$den)) %>%
+ values <- values |>
+ tidyr::unnest(cols = c(.data$data, .data$num, .data$den)) |>
dplyr::mutate(
sw = .data$num / .data$den,
trunc = dplyr::case_when(
@@ -110,7 +110,7 @@ batchmean_ipw <- function(
)
)
- xlev <- unique(data %>% dplyr::pull(.data$.batchvar))
+ xlev <- unique(data |> dplyr::pull(.data$.batchvar))
values <- geepack::geeglm(
formula = variable ~ .batchvar,
@@ -118,14 +118,14 @@ batchmean_ipw <- function(
weights = values$trunc,
id = values$.id,
corstr = "independence"
- ) %>%
- broom::tidy() %>%
+ ) |>
+ broom::tidy() |>
dplyr::filter(
!stringr::str_detect(
string = .data$term,
pattern = "(Intercept)"
)
- ) %>%
+ ) |>
dplyr::mutate(
term = as.character(
stringr::str_remove_all(
@@ -133,11 +133,11 @@ batchmean_ipw <- function(
pattern = ".batchvar"
)
)
- ) %>%
+ ) |>
dplyr::full_join(
tibble::tibble(term = as.character(xlev)),
by = "term"
- ) %>%
+ ) |>
dplyr::mutate(
estimate = dplyr::if_else(
is.na(.data$estimate),
@@ -147,20 +147,20 @@ batchmean_ipw <- function(
estimate = .data$estimate - mean(.data$estimate),
marker = variable,
term = .data$term
- ) %>%
- dplyr::arrange(.data$term) %>%
+ ) |>
+ dplyr::arrange(.data$term) |>
dplyr::select(
.data$marker,
.batchvar = .data$term,
batchmean = .data$estimate
)
- list(values = values, models = res %>% dplyr::pull(.data$den))
+ list(values = values, models = res |> dplyr::pull(.data$den))
}
purrr::map(
- .x = data %>% dplyr::select({{ markers }}) %>% names(),
+ .x = data |> dplyr::select({{ markers }}) |> names(),
.f = ipwbatch,
- data = data %>%
+ data = data |>
dplyr::filter(dplyr::across(dplyr::all_of(confounders), ~ !is.na(.x))),
truncate = truncate,
confounders = paste(confounders, sep = " + ", collapse = " + ")
diff --git a/R/batchmean_simple.R b/R/batchmean_simple.R
index ff8af67..fe1f6cc 100644
--- a/R/batchmean_simple.R
+++ b/R/batchmean_simple.R
@@ -3,23 +3,22 @@
#' @param data Data set
#' @param markers Biomarkers to adjust
#'
-#' @importFrom magrittr %>%
#' @importFrom rlang .data
#' @return Tibble of means per marker and batch
#' @noRd
batchmean_simple <- function(data, markers) {
- values <- data %>%
- dplyr::select(.data$.id, .data$.batchvar, {{ markers }}) %>%
- dplyr::group_by(.data$.batchvar) %>%
+ values <- data |>
+ dplyr::select(.data$.id, .data$.batchvar, {{ markers }}) |>
+ dplyr::group_by(.data$.batchvar) |>
dplyr::summarize_at(
.vars = dplyr::vars(-.data$.id),
.funs = mean,
na.rm = TRUE
- ) %>%
+ ) |>
dplyr::mutate_at(
.vars = dplyr::vars(-.data$.batchvar),
.funs = ~ . - mean(., na.rm = TRUE)
- ) %>%
+ ) |>
tidyr::pivot_longer(
col = c(-.data$.batchvar),
names_to = "marker",
diff --git a/R/batchmean_standardize.R b/R/batchmean_standardize.R
index 208b186..f12953b 100644
--- a/R/batchmean_standardize.R
+++ b/R/batchmean_standardize.R
@@ -8,19 +8,19 @@
#' @return Tibble with conditional means per marker and batch
#' @noRd
batchmean_standardize <- function(data, markers, confounders) {
- res <- data %>%
+ res <- data |>
tidyr::pivot_longer(
cols = {{ markers }},
names_to = "marker",
values_to = "value"
- ) %>%
- dplyr::filter(!is.na(.data$value)) %>%
- dplyr::group_by(.data$marker) %>%
- tidyr::nest(data = c(-.data$marker)) %>%
+ ) |>
+ dplyr::filter(!is.na(.data$value)) |>
+ dplyr::group_by(.data$marker) |>
+ tidyr::nest(data = c(-.data$marker)) |>
dplyr::mutate(
data = purrr::map(
.x = .data$data,
- .f = ~ .x %>%
+ .f = ~ .x |>
dplyr::mutate(.batchvar = factor_drop(.data$.batchvar))
),
model = purrr::map(
@@ -39,36 +39,36 @@ batchmean_standardize <- function(data, markers, confounders) {
),
.batchvar = purrr::map(
.x = .data$data,
- .f = ~ .x %>%
- dplyr::pull(.data$.batchvar) %>%
+ .f = ~ .x |>
+ dplyr::pull(.data$.batchvar) |>
levels()
)
)
- values <- res %>%
- tidyr::unnest(cols = .data$.batchvar) %>%
+ values <- res |>
+ tidyr::unnest(cols = .data$.batchvar) |>
dplyr::mutate(
data = purrr::map2(
.x = .data$data,
.y = .data$.batchvar,
- .f = ~ .x %>% dplyr::mutate(.batchvar = .y)
+ .f = ~ .x |> dplyr::mutate(.batchvar = .y)
),
pred = purrr::map2(.x = .data$model, .y = .data$data, .f = stats::predict)
- ) %>%
- dplyr::select(.data$marker, .data$.batchvar, .data$pred) %>%
- tidyr::unnest(cols = .data$pred) %>%
- dplyr::group_by(.data$marker, .data$.batchvar) %>%
- dplyr::summarize(batchmean = mean(.data$pred, na.rm = TRUE)) %>%
- dplyr::group_by(.data$marker) %>%
- dplyr::mutate(markermean = mean(.data$batchmean)) %>%
- dplyr::ungroup() %>%
+ ) |>
+ dplyr::select(.data$marker, .data$.batchvar, .data$pred) |>
+ tidyr::unnest(cols = .data$pred) |>
+ dplyr::group_by(.data$marker, .data$.batchvar) |>
+ dplyr::summarize(batchmean = mean(.data$pred, na.rm = TRUE)) |>
+ dplyr::group_by(.data$marker) |>
+ dplyr::mutate(markermean = mean(.data$batchmean)) |>
+ dplyr::ungroup() |>
dplyr::transmute(
marker = .data$marker,
.batchvar = .data$.batchvar,
batchmean = .data$batchmean - .data$markermean
)
return(list(list(
- models = res %>% dplyr::ungroup() %>% dplyr::pull("model"),
+ models = res |> dplyr::ungroup() |> dplyr::pull("model"),
values = values
)))
}
diff --git a/R/batchtma.R b/R/batchtma.R
index 296f59d..6c057da 100644
--- a/R/batchtma.R
+++ b/R/batchtma.R
@@ -14,7 +14,7 @@
#'
#' \code{\link[batchtma]{plot_batch}}: Plot biomarkers by batch
#'
-#' @docType package
+#' @docType _PACKAGE
#' @name batchtma
#' @seealso \url{https://stopsack.github.io/batchtma/}
#' @references
From 39acc07c0fb58776a2466d3d8e1aaf8530def1e3 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 19:55:22 +0100
Subject: [PATCH 08/36] r(e)oxygenize
---
DESCRIPTION | 2 +-
man/batchtma.Rd | 2 +-
2 files changed, 2 insertions(+), 2 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index a9d369c..b64d9fa 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -19,7 +19,7 @@ Description: Different adjustment methods for batch effects in biomarker data,
License: GPL-3
Encoding: UTF-8
Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.3
biocViews:
Depends:
R (>= 4.1.0)
diff --git a/man/batchtma.Rd b/man/batchtma.Rd
index fba57e3..df225bc 100644
--- a/man/batchtma.Rd
+++ b/man/batchtma.Rd
@@ -1,6 +1,6 @@
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/batchtma.R
-\docType{package}
+\docType{_PACKAGE}
\name{batchtma}
\alias{batchtma}
\title{batchtma: Methods to address batch effects}
From b7c4d37a6e9e74153f7ea773770e5f3c089a9cf6 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 20:41:49 +0100
Subject: [PATCH 09/36] replace . obj; ~ fn definition; filter(across())
---
DESCRIPTION | 2 +-
R/adjust_batch.R | 16 ++---
R/batch_quantnorm.R | 20 +++---
R/batch_rq.R | 79 +++++++++++++++---------
R/batchmean_ipw.R | 124 ++++++++++++++++++++++++--------------
R/batchmean_simple.R | 2 +-
R/batchmean_standardize.R | 49 +++++++++------
7 files changed, 185 insertions(+), 107 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index b64d9fa..bd78d79 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -25,7 +25,7 @@ Depends:
R (>= 4.1.0)
Imports:
broom (>= 0.7.0),
- dplyr (>= 1.0.0),
+ dplyr (>= 1.0.4),
geepack,
ggplot2,
limma,
diff --git a/R/adjust_batch.R b/R/adjust_batch.R
index 9e58007..d766f96 100644
--- a/R/adjust_batch.R
+++ b/R/adjust_batch.R
@@ -265,7 +265,7 @@ adjust_batch <- function(
)
adjust_parameters <- purrr::map_dfr(
.x = res,
- .f = ~ purrr::pluck(.x, "values")
+ .f = \(x) purrr::pluck(x, "values")
)
method_indices <- c("simple" = 2, "standardize" = 3, "ipw" = 4)
if (suffix == "_adjX") {
@@ -293,10 +293,12 @@ adjust_batch <- function(
.x = data |> dplyr::select({{ markers }}) |> names(),
.f = batchrq,
data = data |>
- dplyr::filter(dplyr::across(
- dplyr::all_of({{ confounders }}),
- ~ !is.na(.x)
- )),
+ dplyr::filter(
+ dplyr::if_all(
+ dplyr::all_of({{ confounders }}),
+ \(x) !is.na(x)
+ )
+ ),
confounders = dplyr::if_else(
dplyr::enexpr(confounders) != "",
true = paste0(
@@ -314,7 +316,7 @@ adjust_batch <- function(
)
adjust_parameters <- purrr::map_dfr(
.x = res,
- .f = ~ purrr::pluck(.x, "values")
+ .f = \(x) purrr::pluck(x, "values")
)
if (suffix == "_adjX") {
suffix <- "_adj5"
@@ -410,7 +412,7 @@ adjust_batch <- function(
adjust_parameters = adjust_parameters,
model_fits = purrr::map(
.x = res,
- .f = ~ purrr::pluck(.x, "models")
+ .f = \(x) purrr::pluck(x, "models")
)
)
class(attr_list) <- c("batchtma", class(res))
diff --git a/R/batch_quantnorm.R b/R/batch_quantnorm.R
index 481c729..856caf5 100644
--- a/R/batch_quantnorm.R
+++ b/R/batch_quantnorm.R
@@ -6,17 +6,22 @@
#' @return Tibble of means per batch for one variable
#' @noRd
batch_quantnorm <- function(var, batch) {
- tibble::tibble(var, batch) |>
+ df <- tibble::tibble(var, batch) |>
tibble::rowid_to_column() |>
- tidyr::pivot_wider(names_from = batch, values_from = var) |>
+ tidyr::pivot_wider(
+ names_from = batch,
+ values_from = var
+ ) |>
dplyr::select(-.data$rowid) |>
- dplyr::select_if(~ !all(is.na(.))) |>
+ dplyr::select(dplyr::where(\(x) !all(is.na(x)))) |>
as.matrix() |>
limma::normalizeQuantiles() |>
- tibble::as_tibble() |>
- dplyr::transmute(
+ tibble::as_tibble()
+
+ df |>
+ dplyr::mutate(
result = purrr::pmap_dbl(
- .l = .,
+ .l = df,
.f = function(...) {
mean(c(...), na.rm = TRUE)
}
@@ -25,7 +30,8 @@ batch_quantnorm <- function(var, batch) {
is.nan(.data$result),
true = NA_real_,
false = .data$result
- )
+ ),
+ .keep = "none"
) |>
dplyr::pull(.data$result)
}
diff --git a/R/batch_rq.R b/R/batch_rq.R
index a6855b4..268e8bf 100644
--- a/R/batch_rq.R
+++ b/R/batch_rq.R
@@ -17,30 +17,36 @@ batchrq <- function(data, variable, confounders, tau, rq_method) {
dplyr::mutate(
un = purrr::map(
.x = .data$data,
- .f = ~ quantreg::rq(
- formula = variable ~ .batchvar,
- data = .x,
- tau = tau,
- method = rq_method
- )
+ .f = \(x) {
+ quantreg::rq(
+ formula = variable ~ .batchvar,
+ data = x,
+ tau = tau,
+ method = rq_method
+ )
+ }
),
ad = purrr::map(
.x = .data$data,
- .f = ~ quantreg::rq(
- formula = stats::reformulate(
- response = "variable",
- termlabels = c(".batchvar", confounders)
- ),
- data = .x,
- tau = tau,
- method = rq_method
- )
+ .f = \(x) {
+ quantreg::rq(
+ formula = stats::reformulate(
+ response = "variable",
+ termlabels = c(".batchvar", confounders)
+ ),
+ data = x,
+ tau = tau,
+ method = rq_method
+ )
+ }
),
.batchvar = purrr::map(
.x = .data$data,
- .f = ~ .x |>
- dplyr::pull(.data$.batchvar) |>
- levels()
+ .f = \(x) {
+ x |>
+ dplyr::pull(.data$.batchvar) |>
+ levels()
+ }
)
)
@@ -50,27 +56,46 @@ batchrq <- function(data, variable, confounders, tau, rq_method) {
data = purrr::map2(
.x = .data$data,
.y = .data$.batchvar,
- .f = ~ .x |> dplyr::mutate(.batchvar = .y)
+ .f = \(x, y) {
+ x |>
+ dplyr::mutate(.batchvar = y)
+ }
+ ),
+ un = purrr::map2(
+ .x = .data$un,
+ .y = .data$data,
+ .f = stats::predict
+ ),
+ ad = purrr::map2(
+ .x = .data$ad,
+ .y = .data$data,
+ .f = stats::predict
),
- un = purrr::map2(.x = .data$un, .y = .data$data, .f = stats::predict),
- ad = purrr::map2(.x = .data$ad, .y = .data$data, .f = stats::predict),
un = purrr::map(
.x = .data$un,
- .f = tibble::as_tibble,
- .name_repair = ~ c("un_lo", "un_hi")
+ .f = \(x) {
+ tibble::as_tibble(
+ x,
+ .name_repair = ~ c("un_lo", "un_hi")
+ )
+ }
),
ad = purrr::map(
.x = .data$ad,
- .f = tibble::as_tibble,
- .name_repair = ~ c("ad_lo", "ad_hi")
+ .f = \(x) {
+ tibble::as_tibble(
+ x,
+ .name_repair = ~ c("ad_lo", "ad_hi")
+ )
+ }
),
all_lo = purrr::map_dbl(
.x = .data$data,
- .f = ~ stats::quantile(.x$variable, probs = 0.25)
+ .f = \(x) stats::quantile(x$variable, probs = 0.25)
),
all_hi = purrr::map_dbl(
.x = .data$data,
- .f = ~ stats::quantile(.x$variable, probs = 0.75)
+ .f = \(x) stats::quantile(x$variable, probs = 0.75)
),
all_iq = .data$all_hi - .data$all_lo
) |>
diff --git a/R/batchmean_ipw.R b/R/batchmean_ipw.R
index 3f1c866..9fcc38b 100644
--- a/R/batchmean_ipw.R
+++ b/R/batchmean_ipw.R
@@ -25,38 +25,57 @@ batchmean_ipw <- function(
dplyr::mutate(
num = purrr::map(
.x = .data$data,
- .f = ~ nnet::multinom(
- formula = .batchvar ~ 1,
- data = .x,
- trace = FALSE
- )
+ .f = \(x) {
+ nnet::multinom(
+ formula = .batchvar ~ 1,
+ data = x,
+ trace = FALSE
+ )
+ }
),
den = purrr::map(
.x = .data$data,
- .f = ~ nnet::multinom(
- formula = stats::as.formula(
- paste(".batchvar ~", confounders)
- ),
- data = .x,
- trace = FALSE
- )
+ .f = \(x) {
+ nnet::multinom(
+ formula = stats::as.formula(
+ paste(".batchvar ~", confounders)
+ ),
+ data = x,
+ trace = FALSE
+ )
+ }
)
)
values <- res |>
dplyr::mutate_at(
.vars = dplyr::vars(.data$num, .data$den),
- .funs = ~ purrr::map(.x = ., .f = stats::predict, type = "probs") |>
- purrr::map(.x = ., .f = tibble::as_tibble) |>
- purrr::map2(
- .x = .,
- .y = .data$data,
- .f = ~ .x |>
- dplyr::mutate(
- .batchvar = .y |>
- purrr::pluck(".batchvar")
+ .funs = \(num_den) {
+ purrr::map(
+ .x = num_den,
+ .f = \(x) {
+ stats::predict(
+ x,
+ type = "probs"
)
- )
+ }
+ ) |>
+ purrr::map(
+ .x = _,
+ .f = tibble::as_tibble
+ ) |>
+ purrr::map2(
+ .x = _,
+ .y = .data$data,
+ .f = \(x, y) {
+ x |>
+ dplyr::mutate(
+ .batchvar = y |>
+ purrr::pluck(".batchvar")
+ )
+ }
+ )
+ }
)
# multinom()$fitted.values is just a vector of probabilities for
@@ -65,35 +84,43 @@ batchmean_ipw <- function(
values <- values |>
dplyr::mutate_at(
.vars = dplyr::vars(.data$num, .data$den),
- .funs = ~ purrr::map(
- .x = .,
- .f = ~ .x |>
- dplyr::mutate(
- probs = dplyr::if_else(
- .data$.batchvar == levels(factor(.data$.batchvar))[1],
- true = 1 - .data$value,
- false = .data$value
- )
- ) |>
- dplyr::pull(.data$probs)
- )
+ .funs = \(num_den) {
+ purrr::map(
+ .x = num_den,
+ .f = \(x) {
+ x |>
+ dplyr::mutate(
+ probs = dplyr::if_else(
+ .data$.batchvar == levels(factor(.data$.batchvar))[1],
+ true = 1 - .data$value,
+ false = .data$value
+ )
+ ) |>
+ dplyr::pull(.data$probs)
+ }
+ )
+ }
)
# otherwise probabilities are a data frame
} else {
values <- values |>
dplyr::mutate_at(
.vars = dplyr::vars(.data$num, .data$den),
- .funs = ~ purrr::map(
- .x = .,
- .f = ~ .x |>
- tidyr::pivot_longer(
- -.data$.batchvar,
- names_to = "batch",
- values_to = "prob"
- ) |>
- dplyr::filter(.data$batch == .data$.batchvar) |>
- dplyr::pull(.data$prob)
- )
+ .funs = \(num_den) {
+ purrr::map(
+ .x = num_den,
+ .f = \(x) {
+ x |>
+ tidyr::pivot_longer(
+ -.data$.batchvar,
+ names_to = "batch",
+ values_to = "prob"
+ ) |>
+ dplyr::filter(.data$batch == .data$.batchvar) |>
+ dplyr::pull(.data$prob)
+ }
+ )
+ }
)
}
@@ -161,7 +188,12 @@ batchmean_ipw <- function(
.x = data |> dplyr::select({{ markers }}) |> names(),
.f = ipwbatch,
data = data |>
- dplyr::filter(dplyr::across(dplyr::all_of(confounders), ~ !is.na(.x))),
+ dplyr::filter(
+ dplyr::if_all(
+ .cols = dplyr::all_of(confounders),
+ .fns = \(x) !is.na(x)
+ )
+ ),
truncate = truncate,
confounders = paste(confounders, sep = " + ", collapse = " + ")
)
diff --git a/R/batchmean_simple.R b/R/batchmean_simple.R
index fe1f6cc..74ee06b 100644
--- a/R/batchmean_simple.R
+++ b/R/batchmean_simple.R
@@ -17,7 +17,7 @@ batchmean_simple <- function(data, markers) {
) |>
dplyr::mutate_at(
.vars = dplyr::vars(-.data$.batchvar),
- .funs = ~ . - mean(., na.rm = TRUE)
+ .funs = \(x) x - mean(x, na.rm = TRUE)
) |>
tidyr::pivot_longer(
col = c(-.data$.batchvar),
diff --git a/R/batchmean_standardize.R b/R/batchmean_standardize.R
index f12953b..bd49e55 100644
--- a/R/batchmean_standardize.R
+++ b/R/batchmean_standardize.R
@@ -20,28 +20,34 @@ batchmean_standardize <- function(data, markers, confounders) {
dplyr::mutate(
data = purrr::map(
.x = .data$data,
- .f = ~ .x |>
- dplyr::mutate(.batchvar = factor_drop(.data$.batchvar))
+ .f = \(x) {
+ x |>
+ dplyr::mutate(.batchvar = factor_drop(.data$.batchvar))
+ }
),
model = purrr::map(
.x = .data$data,
- .f = ~ stats::lm(
- formula = stats::as.formula(paste0(
- "value ~ .batchvar +",
- paste(
- confounders,
- collapse = " + ",
- sep = " + "
- )
- )),
- data = .x
- )
+ .f = \(x) {
+ stats::lm(
+ formula = stats::as.formula(paste0(
+ "value ~ .batchvar +",
+ paste(
+ confounders,
+ collapse = " + ",
+ sep = " + "
+ )
+ )),
+ data = x
+ )
+ }
),
.batchvar = purrr::map(
.x = .data$data,
- .f = ~ .x |>
- dplyr::pull(.data$.batchvar) |>
- levels()
+ .f = \(x) {
+ x |>
+ dplyr::pull(.data$.batchvar) |>
+ levels()
+ }
)
)
@@ -51,9 +57,16 @@ batchmean_standardize <- function(data, markers, confounders) {
data = purrr::map2(
.x = .data$data,
.y = .data$.batchvar,
- .f = ~ .x |> dplyr::mutate(.batchvar = .y)
+ .f = \(x, y) {
+ x |>
+ dplyr::mutate(.batchvar = y)
+ }
),
- pred = purrr::map2(.x = .data$model, .y = .data$data, .f = stats::predict)
+ pred = purrr::map2(
+ .x = .data$model,
+ .y = .data$data,
+ .f = stats::predict
+ )
) |>
dplyr::select(.data$marker, .data$.batchvar, .data$pred) |>
tidyr::unnest(cols = .data$pred) |>
From 4c3cda0192e4550d91eef1fad3d0c6a1f10dec3d Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 20:50:27 +0100
Subject: [PATCH 10/36] dots should no longer appear in pipes
---
R/adjust_batch.R | 6 ------
1 file changed, 6 deletions(-)
diff --git a/R/adjust_batch.R b/R/adjust_batch.R
index d766f96..866a214 100644
--- a/R/adjust_batch.R
+++ b/R/adjust_batch.R
@@ -1,9 +1,3 @@
-# as per https://github.com/jennybc/googlesheets/blob/master/R/googlesheets.R:
-## quiets concerns of R CMD check re: the .'s that appear in pipelines
-if (getRversion() >= "2.15.1") {
- utils::globalVariables(c("."))
-}
-
#' Adjust for batch effects
#'
#' @description
From ce42136574d3ac832c4ca726db4949bffd9295b3 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Mon, 2 Feb 2026 21:50:50 +0100
Subject: [PATCH 11/36] modernize tidyselect
---
R/adjust_batch.R | 73 +++++++++++++++++++++++----------------
R/batch_quantnorm.R | 2 +-
R/batch_rq.R | 16 ++++-----
R/batchmean_ipw.R | 31 +++++++++--------
R/batchmean_simple.R | 8 ++---
R/batchmean_standardize.R | 12 ++++---
6 files changed, 80 insertions(+), 62 deletions(-)
diff --git a/R/adjust_batch.R b/R/adjust_batch.R
index 866a214..b608389 100644
--- a/R/adjust_batch.R
+++ b/R/adjust_batch.R
@@ -187,7 +187,7 @@ adjust_batch <- function(
.batchvar = factor(.data$.batchvar),
.batchvar = factor_drop(.data$.batchvar)
) |>
- dplyr::select(.data$.id, .data$.batchvar, {{ markers }}, {{ confounders }})
+ dplyr::select(".id", ".batchvar", {{ markers }}, {{ confounders }})
confounders <- data |>
dplyr::select({{ confounders }}) |>
names()
@@ -219,7 +219,7 @@ adjust_batch <- function(
paste(dplyr::enexpr(confounders), sep = ", ", collapse = ", "),
"). They will be ignored."
))
- data <- data |> dplyr::select(-dplyr::any_of({{ confounders }}))
+ data <- data |> dplyr::select(!dplyr::any_of({{ confounders }}))
confounders <- NULL
}
@@ -267,24 +267,29 @@ adjust_batch <- function(
}
values <- data |>
- dplyr::select(-dplyr::any_of({{ confounders }})) |>
+ dplyr::select(!dplyr::any_of({{ confounders }})) |>
tidyr::pivot_longer(
- cols = c(-.data$.id, -.data$.batchvar),
+ cols = !c(".id", ".batchvar"),
names_to = "marker",
values_to = "value"
) |>
- dplyr::left_join(adjust_parameters, by = c("marker", ".batchvar")) |>
+ dplyr::left_join(
+ adjust_parameters,
+ by = c("marker", ".batchvar")
+ ) |>
dplyr::mutate(
value_adjusted = .data$value - .data$batchmean,
marker = paste0(.data$marker, suffix)
) |>
- dplyr::select(-.data$batchmean, -.data$value)
+ dplyr::select(!c("batchmean", "value"))
}
# Quantile regression
if (method == "quantreg") {
res <- purrr::map(
- .x = data |> dplyr::select({{ markers }}) |> names(),
+ .x = data |>
+ dplyr::select({{ markers }}) |>
+ names(),
.f = batchrq,
data = data |>
dplyr::filter(
@@ -318,15 +323,18 @@ adjust_batch <- function(
values <- data |>
tidyr::pivot_longer(
- cols = c(
- -.data$.id,
- -.data$.batchvar,
- -dplyr::any_of({{ confounders }})
+ cols = !c(
+ ".id",
+ ".batchvar",
+ dplyr::any_of({{ confounders }})
),
names_to = "marker",
values_to = "value"
) |>
- dplyr::left_join(adjust_parameters, by = c("marker", ".batchvar")) |>
+ dplyr::left_join(
+ adjust_parameters,
+ by = c("marker", ".batchvar")
+ ) |>
dplyr::group_by(.data$marker) |>
dplyr::mutate(
value_adjusted = (.data$value - .data$un_lo) /
@@ -339,16 +347,18 @@ adjust_batch <- function(
marker = paste0(.data$marker, suffix)
) |>
dplyr::select(
- -dplyr::any_of({{ confounders }}),
- -.data$value,
- -.data$un_lo,
- -.data$un_hi,
- -.data$ad_lo,
- -.data$ad_hi,
- -.data$un_iq,
- -.data$ad_iq,
- -.data$all_iq,
- -.data$all_lo
+ !c(
+ dplyr::any_of({{ confounders }}),
+ "value",
+ "un_lo",
+ "un_hi",
+ "ad_lo",
+ "ad_hi",
+ "un_iq",
+ "ad_iq",
+ "all_iq",
+ "all_lo"
+ )
)
}
@@ -358,9 +368,9 @@ adjust_batch <- function(
suffix <- "_adj6"
}
values <- data |>
- dplyr::select(-dplyr::any_of({{ confounders }})) |>
+ dplyr::select(!dplyr::any_of({{ confounders }})) |>
tidyr::pivot_longer(
- cols = c(-.data$.id, -.data$.batchvar),
+ cols = !c(".id", ".batchvar"),
names_to = "marker",
values_to = "value"
) |>
@@ -373,7 +383,7 @@ adjust_batch <- function(
)
) |>
dplyr::ungroup() |>
- dplyr::select(-.data$value)
+ dplyr::select(!"value")
res <- list(list(res = NULL, models = NULL))
adjust_parameters <- tibble::tibble(
marker = data |>
@@ -385,12 +395,15 @@ adjust_batch <- function(
# Dataset to return
values <- values |>
tidyr::pivot_wider(
- names_from = .data$marker,
- values_from = .data$value_adjusted
+ names_from = "marker",
+ values_from = "value_adjusted"
+ ) |>
+ dplyr::select(!".batchvar") |>
+ dplyr::left_join(
+ x = data_orig,
+ by = ".id"
) |>
- dplyr::select(-.data$.batchvar) |>
- dplyr::left_join(x = data_orig, by = ".id") |>
- dplyr::select(-.data$.id)
+ dplyr::select(!".id")
# Meta-data to return as attribute
attr_list <- list(
diff --git a/R/batch_quantnorm.R b/R/batch_quantnorm.R
index 856caf5..c904d4e 100644
--- a/R/batch_quantnorm.R
+++ b/R/batch_quantnorm.R
@@ -12,7 +12,7 @@ batch_quantnorm <- function(var, batch) {
names_from = batch,
values_from = var
) |>
- dplyr::select(-.data$rowid) |>
+ dplyr::select(!"rowid") |>
dplyr::select(dplyr::where(\(x) !all(is.na(x)))) |>
as.matrix() |>
limma::normalizeQuantiles() |>
diff --git a/R/batch_rq.R b/R/batch_rq.R
index 268e8bf..008072a 100644
--- a/R/batch_rq.R
+++ b/R/batch_rq.R
@@ -51,7 +51,7 @@ batchrq <- function(data, variable, confounders, tau, rq_method) {
)
values <- res |>
- tidyr::unnest(cols = .data$.batchvar) |>
+ tidyr::unnest(cols = ".batchvar") |>
dplyr::mutate(
data = purrr::map2(
.x = .data$data,
@@ -100,14 +100,14 @@ batchrq <- function(data, variable, confounders, tau, rq_method) {
all_iq = .data$all_hi - .data$all_lo
) |>
dplyr::select(
- .data$.batchvar,
- .data$un,
- .data$ad,
- .data$all_lo,
- .data$all_hi,
- .data$all_iq
+ ".batchvar",
+ "un",
+ "ad",
+ "all_lo",
+ "all_hi",
+ "all_iq"
) |>
- tidyr::unnest(cols = c(.data$un, .data$ad)) |>
+ tidyr::unnest(cols = c("un", "ad")) |>
dplyr::group_by(.data$.batchvar) |>
dplyr::summarize(
un_lo = stats::quantile(.data$un_lo, probs = 0.25),
diff --git a/R/batchmean_ipw.R b/R/batchmean_ipw.R
index 9fcc38b..6bf35ca 100644
--- a/R/batchmean_ipw.R
+++ b/R/batchmean_ipw.R
@@ -49,7 +49,7 @@ batchmean_ipw <- function(
values <- res |>
dplyr::mutate_at(
- .vars = dplyr::vars(.data$num, .data$den),
+ .vars = c("num", "den"),
.funs = \(num_den) {
purrr::map(
.x = num_den,
@@ -57,16 +57,13 @@ batchmean_ipw <- function(
stats::predict(
x,
type = "probs"
- )
+ ) |>
+ tibble::as_tibble()
}
) |>
- purrr::map(
- .x = _,
- .f = tibble::as_tibble
- ) |>
purrr::map2(
.x = _,
- .y = .data$data,
+ .y = res$data,
.f = \(x, y) {
x |>
dplyr::mutate(
@@ -83,7 +80,7 @@ batchmean_ipw <- function(
if (length(levels(factor(data$.batchvar))) == 2) {
values <- values |>
dplyr::mutate_at(
- .vars = dplyr::vars(.data$num, .data$den),
+ .vars = dplyr::vars("num", "den"),
.funs = \(num_den) {
purrr::map(
.x = num_den,
@@ -125,7 +122,7 @@ batchmean_ipw <- function(
}
values <- values |>
- tidyr::unnest(cols = c(.data$data, .data$num, .data$den)) |>
+ tidyr::unnest(cols = c("data", "num", "den")) |>
dplyr::mutate(
sw = .data$num / .data$den,
trunc = dplyr::case_when(
@@ -177,15 +174,21 @@ batchmean_ipw <- function(
) |>
dplyr::arrange(.data$term) |>
dplyr::select(
- .data$marker,
- .batchvar = .data$term,
- batchmean = .data$estimate
+ "marker",
+ .batchvar = "term",
+ batchmean = "estimate"
)
- list(values = values, models = res |> dplyr::pull(.data$den))
+ list(
+ values = values,
+ models = res |>
+ dplyr::pull(.data$den)
+ )
}
purrr::map(
- .x = data |> dplyr::select({{ markers }}) |> names(),
+ .x = data |>
+ dplyr::select({{ markers }}) |>
+ names(),
.f = ipwbatch,
data = data |>
dplyr::filter(
diff --git a/R/batchmean_simple.R b/R/batchmean_simple.R
index 74ee06b..5d9fcbe 100644
--- a/R/batchmean_simple.R
+++ b/R/batchmean_simple.R
@@ -8,19 +8,19 @@
#' @noRd
batchmean_simple <- function(data, markers) {
values <- data |>
- dplyr::select(.data$.id, .data$.batchvar, {{ markers }}) |>
+ dplyr::select(".id", ".batchvar", {{ markers }}) |>
dplyr::group_by(.data$.batchvar) |>
dplyr::summarize_at(
- .vars = dplyr::vars(-.data$.id),
+ .vars = dplyr::vars(!".id"),
.funs = mean,
na.rm = TRUE
) |>
dplyr::mutate_at(
- .vars = dplyr::vars(-.data$.batchvar),
+ .vars = dplyr::vars(!".batchvar"),
.funs = \(x) x - mean(x, na.rm = TRUE)
) |>
tidyr::pivot_longer(
- col = c(-.data$.batchvar),
+ col = !".batchvar",
names_to = "marker",
values_to = "batchmean"
)
diff --git a/R/batchmean_standardize.R b/R/batchmean_standardize.R
index bd49e55..02d3067 100644
--- a/R/batchmean_standardize.R
+++ b/R/batchmean_standardize.R
@@ -16,7 +16,7 @@ batchmean_standardize <- function(data, markers, confounders) {
) |>
dplyr::filter(!is.na(.data$value)) |>
dplyr::group_by(.data$marker) |>
- tidyr::nest(data = c(-.data$marker)) |>
+ tidyr::nest(data = !"marker") |>
dplyr::mutate(
data = purrr::map(
.x = .data$data,
@@ -52,7 +52,7 @@ batchmean_standardize <- function(data, markers, confounders) {
)
values <- res |>
- tidyr::unnest(cols = .data$.batchvar) |>
+ tidyr::unnest(cols = ".batchvar") |>
dplyr::mutate(
data = purrr::map2(
.x = .data$data,
@@ -68,8 +68,8 @@ batchmean_standardize <- function(data, markers, confounders) {
.f = stats::predict
)
) |>
- dplyr::select(.data$marker, .data$.batchvar, .data$pred) |>
- tidyr::unnest(cols = .data$pred) |>
+ dplyr::select("marker", ".batchvar", "pred") |>
+ tidyr::unnest(cols = "pred") |>
dplyr::group_by(.data$marker, .data$.batchvar) |>
dplyr::summarize(batchmean = mean(.data$pred, na.rm = TRUE)) |>
dplyr::group_by(.data$marker) |>
@@ -81,7 +81,9 @@ batchmean_standardize <- function(data, markers, confounders) {
batchmean = .data$batchmean - .data$markermean
)
return(list(list(
- models = res |> dplyr::ungroup() |> dplyr::pull("model"),
+ models = res |>
+ dplyr::ungroup() |>
+ dplyr::pull("model"),
values = values
)))
}
From be3cb030e008bea6092f76739783093ea678ce16 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Tue, 3 Feb 2026 08:45:28 +0100
Subject: [PATCH 12/36] _at() -> across()
---
R/batchmean_ipw.R | 123 ++++++++++++++++++++++---------------------
R/batchmean_simple.R | 20 ++++---
2 files changed, 77 insertions(+), 66 deletions(-)
diff --git a/R/batchmean_ipw.R b/R/batchmean_ipw.R
index 6bf35ca..05e4dd2 100644
--- a/R/batchmean_ipw.R
+++ b/R/batchmean_ipw.R
@@ -48,76 +48,81 @@ batchmean_ipw <- function(
)
values <- res |>
- dplyr::mutate_at(
- .vars = c("num", "den"),
- .funs = \(num_den) {
- purrr::map(
- .x = num_den,
- .f = \(x) {
- stats::predict(
- x,
- type = "probs"
- ) |>
- tibble::as_tibble()
- }
- ) |>
- purrr::map2(
- .x = _,
- .y = res$data,
- .f = \(x, y) {
- x |>
- dplyr::mutate(
- .batchvar = y |>
- purrr::pluck(".batchvar")
- )
+ dplyr::mutate(
+ dplyr::across(
+ .cols = c("num", "den"),
+ .fns = \(num_den) {
+ purrr::map(
+ .x = num_den,
+ .f = \(x) {
+ stats::predict(
+ x,
+ type = "probs"
+ ) |>
+ tibble::as_tibble()
}
- )
- }
+ ) |>
+ purrr::map2( # .x = _ # would require R 4.2.0
+ .y = res$data,
+ .f = \(x, y) {
+ x |>
+ dplyr::mutate(
+ .batchvar = y |>
+ purrr::pluck(".batchvar")
+ )
+ }
+ )
+ }
+ )
)
# multinom()$fitted.values is just a vector of probabilities for
# the 2nd outcome level if there are only two levels
if (length(levels(factor(data$.batchvar))) == 2) {
values <- values |>
- dplyr::mutate_at(
- .vars = dplyr::vars("num", "den"),
- .funs = \(num_den) {
- purrr::map(
- .x = num_den,
- .f = \(x) {
- x |>
- dplyr::mutate(
- probs = dplyr::if_else(
- .data$.batchvar == levels(factor(.data$.batchvar))[1],
- true = 1 - .data$value,
- false = .data$value
- )
- ) |>
- dplyr::pull(.data$probs)
- }
- )
- }
+ dplyr::mutate(
+ dplyr::across(
+ .cols = c("num", "den"),
+ .fns = \(num_den) {
+ purrr::map(
+ .x = num_den,
+ .f = \(x) {
+ x |>
+ dplyr::mutate(
+ probs = dplyr::if_else(
+ .data$.batchvar == levels(factor(.data$.batchvar))[1],
+ true = 1 - .data$value,
+ false = .data$value
+ )
+ ) |>
+ dplyr::pull(.data$probs)
+ }
+ )
+ }
+ )
)
# otherwise probabilities are a data frame
} else {
values <- values |>
- dplyr::mutate_at(
- .vars = dplyr::vars(.data$num, .data$den),
- .funs = \(num_den) {
- purrr::map(
- .x = num_den,
- .f = \(x) {
- x |>
- tidyr::pivot_longer(
- -.data$.batchvar,
- names_to = "batch",
- values_to = "prob"
- ) |>
- dplyr::filter(.data$batch == .data$.batchvar) |>
- dplyr::pull(.data$prob)
- }
- )
- }
+ dplyr::mutate(
+ dplyr::across(
+ .cols = c("num", "den"),
+ .fns = \(num_den) {
+ purrr::map(
+ .x = num_den,
+ .f = \(x) {
+ x |>
+ tidyr::pivot_longer(
+ cols = c(!".batchvar"),
+ names_to = "batch",
+ values_to = "prob"
+ ) |>
+ dplyr::filter(.data$batch == .data$.batchvar) |>
+ dplyr::pull(.data$prob)
+ }
+ )
+ }
+ )
)
}
diff --git a/R/batchmean_simple.R b/R/batchmean_simple.R
index 5d9fcbe..aca2139 100644
--- a/R/batchmean_simple.R
+++ b/R/batchmean_simple.R
@@ -10,14 +10,20 @@ batchmean_simple <- function(data, markers) {
values <- data |>
dplyr::select(".id", ".batchvar", {{ markers }}) |>
dplyr::group_by(.data$.batchvar) |>
- dplyr::summarize_at(
- .vars = dplyr::vars(!".id"),
- .funs = mean,
- na.rm = TRUE
+ dplyr::summarize(
+ dplyr::across(
+ .cols = !".id",
+ .fns = \(x) mean(
+ x,
+ na.rm = TRUE
+ )
+ )
) |>
- dplyr::mutate_at(
- .vars = dplyr::vars(!".batchvar"),
- .funs = \(x) x - mean(x, na.rm = TRUE)
+ dplyr::mutate(
+ dplyr::across(
+ .cols = !".batchvar",
+ .fns = \(x) x - mean(x, na.rm = TRUE)
+ )
) |>
tidyr::pivot_longer(
col = !".batchvar",
From 70dfc767e04fb544ff3989f10324a5f7001da91c Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Tue, 3 Feb 2026 08:45:55 +0100
Subject: [PATCH 13/36] transmute() has been superseded
---
R/batchmean_standardize.R | 5 +++--
1 file changed, 3 insertions(+), 2 deletions(-)
diff --git a/R/batchmean_standardize.R b/R/batchmean_standardize.R
index 02d3067..3d0de84 100644
--- a/R/batchmean_standardize.R
+++ b/R/batchmean_standardize.R
@@ -75,10 +75,11 @@ batchmean_standardize <- function(data, markers, confounders) {
dplyr::group_by(.data$marker) |>
dplyr::mutate(markermean = mean(.data$batchmean)) |>
dplyr::ungroup() |>
- dplyr::transmute(
+ dplyr::mutate(
marker = .data$marker,
.batchvar = .data$.batchvar,
- batchmean = .data$batchmean - .data$markermean
+ batchmean = .data$batchmean - .data$markermean,
+ .keep = "none"
)
return(list(list(
models = res |>
From f7438fe3aeb09b41ae6bd1a3eb44adade3dff8d6 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Tue, 3 Feb 2026 08:47:34 +0100
Subject: [PATCH 14/36] change with roxygen 7.0.0 in 2019
---
R/batchtma.R | 4 +---
man/batchtma.Rd | 12 +++++++++++-
2 files changed, 12 insertions(+), 4 deletions(-)
diff --git a/R/batchtma.R b/R/batchtma.R
index 6c057da..283fb33 100644
--- a/R/batchtma.R
+++ b/R/batchtma.R
@@ -14,7 +14,6 @@
#'
#' \code{\link[batchtma]{plot_batch}}: Plot biomarkers by batch
#'
-#' @docType _PACKAGE
#' @name batchtma
#' @seealso \url{https://stopsack.github.io/batchtma/}
#' @references
@@ -23,5 +22,4 @@
#' Mucci LA+ (+ equal contribution). Extent, impact, and mitigation of
#' batch effects in tumor biomarker studies using tissue microarrays.
#' eLife 2021;10:e71265. doi: https://doi.org/10.7554/elife.71265
-NULL
-# > NULL
+"_PACKAGE"
diff --git a/man/batchtma.Rd b/man/batchtma.Rd
index df225bc..a889c21 100644
--- a/man/batchtma.Rd
+++ b/man/batchtma.Rd
@@ -1,7 +1,8 @@
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/batchtma.R
-\docType{_PACKAGE}
+\docType{package}
\name{batchtma}
+\alias{batchtma-package}
\alias{batchtma}
\title{batchtma: Methods to address batch effects}
\description{
@@ -31,3 +32,12 @@ eLife 2021;10:e71265. doi: https://doi.org/10.7554/elife.71265
\seealso{
\url{https://stopsack.github.io/batchtma/}
}
+\author{
+\strong{Maintainer}: Konrad Stopsack \email{stopsack@post.harvard.edu} (\href{https://orcid.org/0000-0002-0722-1311}{ORCID})
+
+Authors:
+\itemize{
+ \item Travis Gerke \email{gerket@mskcc.org} (\href{https://orcid.org/0000-0002-9500-8907}{ORCID})
+}
+
+}
From 206d0ed02d08360329a13a6a78750624cd7ef5bb Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Tue, 3 Feb 2026 08:48:09 +0100
Subject: [PATCH 15/36] %>% pipe and . variables are gone
---
DESCRIPTION | 3 +--
1 file changed, 1 insertion(+), 2 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index bd78d79..9d99606 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -35,8 +35,7 @@ Imports:
rlang (>= 0.4.0),
stringr,
tibble,
- tidyr (>= 1.1.0),
- magrittr
+ tidyr (>= 1.1.0)
Suggests:
knitr,
rmarkdown,
From d30d5edfa2cfca9a8dd0dc2db4325864ea6693d2 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Tue, 3 Feb 2026 08:54:41 +0100
Subject: [PATCH 16/36] unify syntax
---
R/batchmean_ipw.R | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/R/batchmean_ipw.R b/R/batchmean_ipw.R
index 05e4dd2..5a03b12 100644
--- a/R/batchmean_ipw.R
+++ b/R/batchmean_ipw.R
@@ -113,7 +113,7 @@ batchmean_ipw <- function(
.f = \(x) {
x |>
tidyr::pivot_longer(
- cols = c(!".batchvar"),
+ cols = !".batchvar",
names_to = "batch",
values_to = "prob"
) |>
From 0c531e1a6607eb4d9b6cc7d68fa53ca278889d21 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Tue, 3 Feb 2026 08:55:11 +0100
Subject: [PATCH 17/36] misnamed argument, apparently inconsequential
---
R/batchmean_simple.R | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/R/batchmean_simple.R b/R/batchmean_simple.R
index aca2139..e19777d 100644
--- a/R/batchmean_simple.R
+++ b/R/batchmean_simple.R
@@ -26,7 +26,7 @@ batchmean_simple <- function(data, markers) {
)
) |>
tidyr::pivot_longer(
- col = !".batchvar",
+ cols = !".batchvar",
names_to = "marker",
values_to = "batchmean"
)
From 4aec1a765b95320e273f5e0616d5234a78480501 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Fri, 6 Feb 2026 23:56:26 +0100
Subject: [PATCH 18/36] add covr action from rifttable
---
.github/workflows/R-CMD-check.yaml | 48 ++++++++++++++++++++++--------
1 file changed, 35 insertions(+), 13 deletions(-)
diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index ad35fcf..f6d2dcc 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -100,11 +100,6 @@ jobs:
BiocCheck::BiocCheck(".")
shell: Rscript {0}
- - name: Show testthat output
- if: always()
- run: find check -name 'testthat.Rout*' -exec cat '{}' \; || true
- shell: bash
-
- name: Upload check results
if: failure()
uses: actions/upload-artifact@HEAD
@@ -112,19 +107,46 @@ jobs:
name: ${{ runner.os }}-r${{ matrix.config.r }}-bioc-${{ matrix.config.bioc }}-results
path: check
+ - name: pkgdown dependencies
+ if: github.event_name == 'push' && matrix.config.os == 'macOS-latest'
+ uses: r-lib/actions/setup-r-dependencies@v2
+ with:
+ extra-packages: any::pkgdown, local::., any::covr, any::xml2
+ needs: website, coverage
+
- name: Test coverage
- if: matrix.config.os == 'macOS-latest'
run: |
- install.packages("covr")
- covr::codecov(token = "${{secrets.CODECOV_TOKEN}}")
+ 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}
- - name: pkgdown dependencies
- if: github.event_name == 'push' && matrix.config.os == 'macOS-latest'
- uses: r-lib/actions/setup-r-dependencies@v2
+ - 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:
- extra-packages: any::pkgdown, local::.
- needs: website
+ name: coverage-test-failures
+ path: ${{ runner.temp }}/package
- name: Build site
if: github.event_name == 'push' && matrix.config.os == 'macOS-latest'
From c04b3c1ffb9e28c439165c371bb1980c27c630e9 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Fri, 6 Feb 2026 23:57:42 +0100
Subject: [PATCH 19/36] trigger new action
---
.github/workflows/R-CMD-check.yaml | 2 ++
1 file changed, 2 insertions(+)
diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index f6d2dcc..d83fe20 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -1,5 +1,7 @@
on:
push:
+ branches:
+ - dev
pull_request:
branches:
- main
From e9e9711cb9c737a98e54ef626fb65b20256723af Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Fri, 20 Feb 2026 20:12:57 +0100
Subject: [PATCH 20/36] give intelligible error for empty method argument
---
R/adjust_batch.R | 8 ++++++++
1 file changed, 8 insertions(+)
diff --git a/R/adjust_batch.R b/R/adjust_batch.R
index b608389..5058f44 100644
--- a/R/adjust_batch.R
+++ b/R/adjust_batch.R
@@ -193,6 +193,14 @@ adjust_batch <- function(
names()
# Check inputs: method and confounders
+ if(method[1] == "c") {
+ stop(paste0(
+ "No valid argument 'method =' provided. Available methods: ",
+ paste(allmethods, collapse = ", "),
+ ".\n",
+ "See: help(\"adjust_batch\")."
+ ))
+ }
if (!(method[1] %in% allmethods)) {
stop(paste0(
"Method '",
From 94fcf6df87b28ff23d2de64fe9d132aea5b9142a Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Fri, 20 Feb 2026 20:13:15 +0100
Subject: [PATCH 21/36] add some missing tests
---
tests/testthat/test-diagnostics.R | 21 +++++++++++++++
tests/testthat/test-equality.R | 33 +++++++++++++++++++++++
tests/testthat/test-input.R | 45 +++++++++++++++++++++++++++++++
tests/testthat/test-plot.R | 17 ++++++++++++
4 files changed, 116 insertions(+)
create mode 100644 tests/testthat/test-diagnostics.R
create mode 100644 tests/testthat/test-input.R
create mode 100644 tests/testthat/test-plot.R
diff --git a/tests/testthat/test-diagnostics.R b/tests/testthat/test-diagnostics.R
new file mode 100644
index 0000000..bad385f
--- /dev/null
+++ b/tests/testthat/test-diagnostics.R
@@ -0,0 +1,21 @@
+test_that("diagnostics run", {
+ df <- data.frame(
+ tma = rep(1:2, times = 10),
+ biomarker = rep(1:2, times = 10) +
+ runif(max = 5, n = 20),
+ confounder = rep(0:1, times = 10) +
+ runif(max = 10, n = 20)
+ )
+ expect_no_error(
+ object = print(
+ diagnose_models(
+ adjust_batch(
+ data = df,
+ marker = biomarker,
+ batch = tma,
+ method = simple
+ )
+ )
+ )
+ )
+})
diff --git a/tests/testthat/test-equality.R b/tests/testthat/test-equality.R
index 6ca325a..bf9d93c 100644
--- a/tests/testthat/test-equality.R
+++ b/tests/testthat/test-equality.R
@@ -40,3 +40,36 @@ test_that("equivalent approaches give same result", {
expect_equal(df_adj2, df_adj3)
expect_equal(df_adj2, df_adj4)
})
+
+
+test_that("multinomial model works", {
+ df <- data.frame(
+ tma = rep(1:3, times = 10),
+ biomarker = rep(1:2, times = 15) +
+ runif(max = 5, n = 30),
+ confounder = rep(0:1, times = 15) +
+ runif(max = 10, n = 30),
+ unity = 1
+ )
+
+ df_adj2 <- adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = simple,
+ suffix = "adj"
+ ) |> # drop all attributes
+ data.frame()
+
+ df_adj4 <- adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = ipw,
+ confounders = unity,
+ suffix = "adj"
+ ) |>
+ data.frame()
+
+ expect_equal(df_adj2, df_adj4)
+})
diff --git a/tests/testthat/test-input.R b/tests/testthat/test-input.R
new file mode 100644
index 0000000..0207d87
--- /dev/null
+++ b/tests/testthat/test-input.R
@@ -0,0 +1,45 @@
+test_that("input errors are detected", {
+ df <- data.frame(
+ tma = rep(1:2, times = 10),
+ biomarker = rep(1:2, times = 10) +
+ runif(max = 5, n = 20),
+ confounder = rep(0:1, times = 10) +
+ runif(max = 10, n = 20)
+ )
+ expect_error(
+ object = adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma
+ ),
+ regexp = "No valid argument \'method"
+ )
+ expect_error(
+ object = adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = inexistant
+ ),
+ regexp = "is not implemented"
+ )
+ expect_message(
+ object = adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = simple,
+ confounders = confounder
+ ),
+ regexp = "does not support adjustment"
+ )
+ expect_message(
+ object = adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = standardize
+ ),
+ regexp = "no valid confounders"
+ )
+})
diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R
new file mode 100644
index 0000000..1179701
--- /dev/null
+++ b/tests/testthat/test-plot.R
@@ -0,0 +1,17 @@
+test_that("ggplot runs", {
+ df <- data.frame(
+ tma = rep(1:2, times = 10),
+ biomarker = rep(1:2, times = 10) +
+ runif(max = 5, n = 20),
+ confounder = rep(0:1, times = 10) +
+ runif(max = 10, n = 20)
+ )
+ expect_no_error(
+ object = plot_batch(
+ data = df,
+ marker = biomarker,
+ batch = tma,
+ color = confounder
+ )
+ )
+})
From 0a2855c2f7bb6400b5cd0a596ae691b5333e4c40 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Fri, 20 Feb 2026 20:16:07 +0100
Subject: [PATCH 22/36] catch pipe in vignette
---
vignettes/batchtma.Rmd | 58 +++++++++++++++++++++---------------------
1 file changed, 29 insertions(+), 29 deletions(-)
diff --git a/vignettes/batchtma.Rmd b/vignettes/batchtma.Rmd
index 0d79ea0..300d903 100644
--- a/vignettes/batchtma.Rmd
+++ b/vignettes/batchtma.Rmd
@@ -63,7 +63,7 @@ df
We plot the biomarker values (*y*-axis) by batch (*x*-axis), using the `plot_batch()` function. Color/shape symbolizes which participant/tumor the measurements came from. Boxes span from the 25th to the 75th percentile (interquartile range); thick lines indicate medians; asterisks indicate means.
```{r}
-df %>% plot_batch(marker = noisy, batch = batch, color = person)
+df |> plot_batch(marker = noisy, batch = batch, color = person)
```
@@ -72,14 +72,14 @@ df %>% plot_batch(marker = noisy, batch = batch, color = person)
We add systematic differences between batches such that there is differential measurement error between batches in terms of mean and variance. As shown above, true values were the same beyond random error.
```{r}
-df <- df %>%
+df <- df |>
# Multiply by batch number to differentially change variance by batch,
# divide by mean batch number to keep overall variance the same:
mutate(noisy_batch = noisy * batchnum / mean(batchnum) +
# Similarly, change mean value per batch, keeping overall mean the same:
batchnum * 3 - mean(batchnum) * 3)
-df %>% plot_batch(marker = noisy_batch, batch = batch, color = person)
+df |> plot_batch(marker = noisy_batch, batch = batch, color = person)
```
@@ -90,9 +90,9 @@ df %>% plot_batch(marker = noisy_batch, batch = batch, color = person)
`method = simple` calculates the mean for each batch and subtracts the difference between this mean and the grand mean, such that all batches end up having a mean equivalent to the grand mean. Differences in variance between batches will remain, if they exist (as in this example).
```{r}
-df %>%
+df |>
adjust_batch(markers = noisy_batch, batch = batch,
- method = simple) %>%
+ method = simple) |>
plot_batch(marker = noisy_batch_adj2, batch = batch, color = person)
```
@@ -102,9 +102,9 @@ df %>%
`method = standardize` performs marginal standardization by fitting a linear regression model for biomarker values with `batch` and `confounders` as predictors, and obtains the marginal means per batch if they had the same distribution of `confounders`. Differences between these marginal means and the grand mean are subtracted as in `method = simple`. In this example, the confounder is a random variable, and the results are essentially the same as for `method = simple`.
```{r}
-df %>%
+df |>
adjust_batch(markers = noisy_batch, batch = batch,
- method = standardize, confounders = random) %>%
+ method = standardize, confounders = random) |>
plot_batch(marker = noisy_batch_adj3, batch = batch, color = person)
```
@@ -114,9 +114,9 @@ df %>%
`method = ipw` predicts the probability of a measurement being from a specific batch, given the `confounders`. Mean differences between batches are obtained from a marginal structural model with stabilized inverse-probability weights and then used as in the two preceding methods. Again, the confounder is merely a random variable in this example.
```{r}
-df %>%
+df |>
adjust_batch(markers = noisy_batch, batch = batch,
- method = ipw, confounders = random) %>%
+ method = ipw, confounders = random) |>
plot_batch(marker = noisy_batch_adj4, batch = batch, color = person)
```
@@ -126,9 +126,9 @@ df %>%
`method = quantreg`, unlike the three preceding mean-only methods, addresses two distinct properties of batches: the offset values (a lower quantile), potentially reflective of background signal, and an inter-quantile range, potentially reflective of the dynamic range of the measurement. By default, the first and third qua**r**tile are used.
```{r}
-df %>%
+df |>
adjust_batch(markers = noisy_batch, batch = batch,
- method = quantreg, confounders = random) %>%
+ method = quantreg, confounders = random) |>
plot_batch(marker = noisy_batch_adj5, batch = batch, color = person)
```
@@ -137,9 +137,9 @@ df %>%
`method = quantnorm` performs quantile normalization: values are ranked within each batch, and then each rank is assigned the mean per rank across batches. Quantile normalization ensures that all batches have near-identical biomarker distributions. However, quantile normalization does not allow for accounting for `confounders`.
```{r}
-df %>%
+df |>
adjust_batch(markers = noisy_batch, batch = batch,
- method = quantnorm) %>%
+ method = quantnorm) |>
plot_batch(marker = noisy_batch_adj6, batch = batch, color = person)
```
@@ -157,25 +157,25 @@ All following plots show color/symbol shape by `confounder`.
```{r}
set.seed(123) # for reproducibility
-df <- df %>%
+df <- df |>
# Make confounder associated with batch:
mutate(confounder = round(batchnum + runif(n = 200, max = 2)),
# Make biomarker values associated with confounder:
noisy_conf = noisy + confounder * 3 - mean(confounder) * 3)
-df %>% plot_batch(marker = noisy_conf, batch = batch, color = confounder)
+df |> plot_batch(marker = noisy_conf, batch = batch, color = confounder)
```
Second, batch effects for mean and variance are added.
```{r}
-df <- df %>%
+df <- df |>
# Add batch effects to confounded biomarker values:
mutate(noisy_conf_batch = noisy_conf * batchnum / mean(batchnum) +
batchnum * 3 - mean(batchnum) * 3)
-df %>% plot_batch(marker = noisy_conf_batch, batch = batch, color = confounder)
+df |> plot_batch(marker = noisy_conf_batch, batch = batch, color = confounder)
```
The following examples show to what extent batch effects can be removed from `noisy_conf_batch` to recover the "ground truth," `noisy_conf`.
@@ -186,9 +186,9 @@ The following examples show to what extent batch effects can be removed from `no
`method = standardize` reduces batch effects while allowing batch means to differ because of confounders. Batches with higher true biomarker values associated with the confounders, like batch B, retain higher means after batch effect adjustment.
```{r}
-df %>%
+df |>
adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = standardize, confounders = confounder) %>%
+ method = standardize, confounders = confounder) |>
plot_batch(marker = noisy_conf_batch_adj3, batch = batch, color = confounder)
```
@@ -198,9 +198,9 @@ df %>%
`method = ipw` also reduces batch effects while allowing batch means to differ because of confounders. As with marginal standardization, differences in variance between batches are not addressed.
```{r}
-df %>%
+df |>
adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = ipw, confounders = confounder) %>%
+ method = ipw, confounders = confounder) |>
plot_batch(marker = noisy_conf_batch_adj4, batch = batch, color = confounder)
```
@@ -210,9 +210,9 @@ df %>%
`method = quantreg` reduces for batch effects of offset and dynamic range separately, accounting for differences in confounders. In the example, batch B with its high levels of the `confounder` retains somewhat higher values and a higher variance, unlike with the two preceding mean-only methods.
```{r}
-df %>%
+df |>
adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = quantreg, confounders = confounder) %>%
+ method = quantreg, confounders = confounder) |>
plot_batch(marker = noisy_conf_batch_adj5, batch = batch, color = confounder)
```
@@ -224,9 +224,9 @@ df %>%
`method = simple` leads to all batches having the same mean after adjustment. In the example, this method ignores that batch B should have higher values because of higher values of the `confounder`.
```{r}
-df %>%
+df |>
adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = simple, confounders = confounder) %>%
+ method = simple, confounders = confounder) |>
plot_batch(marker = noisy_conf_batch_adj2, batch = batch, color = confounder)
```
@@ -236,9 +236,9 @@ df %>%
`method = quantreg` also ignores the confounder.
```{r}
-df %>%
+df |>
adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = quantnorm, confounders = confounder) %>%
+ method = quantnorm, confounders = confounder) |>
plot_batch(marker = noisy_conf_batch_adj6, batch = batch, color = confounder)
```
@@ -264,7 +264,7 @@ To obtain an overview of model diagnostics, `diagnose_models()` is called on the
```{r}
# df2 is the new dataset that also contains "noisy_conf_batch_adj2":
-df2 <- df %>% adjust_batch(markers = noisy_conf_batch, batch = batch,
+df2 <- df |> adjust_batch(markers = noisy_conf_batch, batch = batch,
method = standardize, confounders = confounder)
# Show overview of model diagnostics:
@@ -287,7 +287,7 @@ Diagnostics for regression models can be performed on fitted models just as for
```{r}
tibble(fitted = fitted.values(fit),
- residuals = residuals(fit)) %>%
+ residuals = residuals(fit)) |>
ggplot(mapping = aes(x = fitted, y = residuals)) +
geom_point() +
theme_minimal()
From 6c78e36ea61e333fd4e3272bf1a1deaf11999445 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Fri, 20 Feb 2026 20:49:31 +0100
Subject: [PATCH 23/36] match install instructions with index.Rmd
---
vignettes/batchtma.Rmd | 14 +++++++-------
1 file changed, 7 insertions(+), 7 deletions(-)
diff --git a/vignettes/batchtma.Rmd b/vignettes/batchtma.Rmd
index 300d903..8015609 100644
--- a/vignettes/batchtma.Rmd
+++ b/vignettes/batchtma.Rmd
@@ -16,17 +16,17 @@ knitr::opts_chunk$set(
# Installation
-Batchtma can be installed from [GitHub](https://github.com/) using:
+batchtma can be installed from [CRAN](https://cran.r-project.org/) using:
-```{r, eval = FALSE}
-# install.packages("remotes") # The "remotes" package needs to be installed
-remotes::install_github("stopsack/batchtma")
+```{r}
+install.packages("batchtma")
```
-To have vignettes like this current site be available offline via, *e.g.*, `vignette("batchtma")`, modify the last command:
+To install a potentially newer version from [GitHub](https://github.com/), use:
-```{r, eval = FALSE}
-remotes::install_github("stopsack/batchtma", build_vignettes = TRUE)
+```{r}
+# install.packages("remotes") # The "remotes" package needs to be installed
+remotes::install_github("stopsack/batchtma")
```
From 1f5b45c75ce680d6cb1de70e4ada89b5b0c557ad Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Fri, 20 Feb 2026 20:49:39 +0100
Subject: [PATCH 24/36] reformat code
---
vignettes/batchtma.Rmd | 242 ++++++++++++++++++++++++++++++++---------
1 file changed, 189 insertions(+), 53 deletions(-)
diff --git a/vignettes/batchtma.Rmd b/vignettes/batchtma.Rmd
index 8015609..3a0c3b9 100644
--- a/vignettes/batchtma.Rmd
+++ b/vignettes/batchtma.Rmd
@@ -46,16 +46,34 @@ We construct a toy dataset of 10 individuals (*e.g.,* tumors), each with 40 meas
set.seed(123) # for reproducibility
df <- tibble(
# Batches:
- batch = rep(paste0("batch", LETTERS[1:4]), times = 100),
- batchnum = rep(c(1, 5, 2, 3), times = 100),
+ batch = rep(
+ paste0("batch", LETTERS[1:4]),
+ times = 100
+ ),
+ batchnum = rep(
+ c(1, 5, 2, 3),
+ times = 100
+ ),
# Participants:
- person = rep(letters[1:10], each = 40),
+ person = rep(
+ letters[1:10],
+ each = 40
+ ),
# Instead of a confounder, we will use a random variable for now:
- random = runif(n = 400, min = -2, max = 2),
+ random = runif(
+ n = 400,
+ min = -2,
+ max = 2
+ ),
# The true (usually unobservable biomarker value):
- true = rep(c(2, 2.5, 3, 5, 6, 8, 10, 12, 15, 12), each = 40),
+ true = rep(
+ c(2, 2.5, 3, 5, 6, 8, 10, 12, 15, 12),
+ each = 40
+ ),
# The observed biomarker value with random error ("noise"):
- noisy = true + runif(max = true / 3, n = 400) * 4)
+ noisy = true +
+ runif(max = true / 3, n = 400) * 4
+)
df
```
@@ -63,7 +81,12 @@ df
We plot the biomarker values (*y*-axis) by batch (*x*-axis), using the `plot_batch()` function. Color/shape symbolizes which participant/tumor the measurements came from. Boxes span from the 25th to the 75th percentile (interquartile range); thick lines indicate medians; asterisks indicate means.
```{r}
-df |> plot_batch(marker = noisy, batch = batch, color = person)
+df |>
+ plot_batch(
+ marker = noisy,
+ batch = batch,
+ color = person
+ )
```
@@ -75,11 +98,19 @@ We add systematic differences between batches such that there is differential me
df <- df |>
# Multiply by batch number to differentially change variance by batch,
# divide by mean batch number to keep overall variance the same:
- mutate(noisy_batch = noisy * batchnum / mean(batchnum) +
- # Similarly, change mean value per batch, keeping overall mean the same:
- batchnum * 3 - mean(batchnum) * 3)
+ mutate(
+ noisy_batch = noisy * batchnum / mean(batchnum) +
+ # Similarly, change mean value per batch, keeping overall mean the same:
+ batchnum * 3 -
+ mean(batchnum) * 3
+ )
-df |> plot_batch(marker = noisy_batch, batch = batch, color = person)
+df |>
+ plot_batch(
+ marker = noisy_batch,
+ batch = batch,
+ color = person
+ )
```
@@ -91,9 +122,16 @@ df |> plot_batch(marker = noisy_batch, batch = batch, color = person)
```{r}
df |>
- adjust_batch(markers = noisy_batch, batch = batch,
- method = simple) |>
- plot_batch(marker = noisy_batch_adj2, batch = batch, color = person)
+ adjust_batch(
+ markers = noisy_batch,
+ batch = batch,
+ method = simple
+ ) |>
+ plot_batch(
+ marker = noisy_batch_adj2,
+ batch = batch,
+ color = person
+ )
```
@@ -103,9 +141,17 @@ df |>
```{r}
df |>
- adjust_batch(markers = noisy_batch, batch = batch,
- method = standardize, confounders = random) |>
- plot_batch(marker = noisy_batch_adj3, batch = batch, color = person)
+ adjust_batch(
+ markers = noisy_batch,
+ batch = batch,
+ method = standardize,
+ confounders = random
+ ) |>
+ plot_batch(
+ marker = noisy_batch_adj3,
+ batch = batch,
+ color = person
+ )
```
@@ -115,9 +161,17 @@ df |>
```{r}
df |>
- adjust_batch(markers = noisy_batch, batch = batch,
- method = ipw, confounders = random) |>
- plot_batch(marker = noisy_batch_adj4, batch = batch, color = person)
+ adjust_batch(
+ markers = noisy_batch,
+ batch = batch,
+ method = ipw,
+ confounders = random
+ ) |>
+ plot_batch(
+ marker = noisy_batch_adj4,
+ batch = batch,
+ color = person
+ )
```
@@ -127,9 +181,17 @@ df |>
```{r}
df |>
- adjust_batch(markers = noisy_batch, batch = batch,
- method = quantreg, confounders = random) |>
- plot_batch(marker = noisy_batch_adj5, batch = batch, color = person)
+ adjust_batch(
+ markers = noisy_batch,
+ batch = batch,
+ method = quantreg,
+ confounders = random
+ ) |>
+ plot_batch(
+ marker = noisy_batch_adj5,
+ batch = batch,
+ color = person
+ )
```
## Quantile normalization
@@ -138,9 +200,16 @@ df |>
```{r}
df |>
- adjust_batch(markers = noisy_batch, batch = batch,
- method = quantnorm) |>
- plot_batch(marker = noisy_batch_adj6, batch = batch, color = person)
+ adjust_batch(
+ markers = noisy_batch,
+ batch = batch,
+ method = quantnorm
+ ) |>
+ plot_batch(
+ marker = noisy_batch_adj6,
+ batch = batch,
+ color = person
+ )
```
@@ -159,11 +228,18 @@ All following plots show color/symbol shape by `confounder`.
set.seed(123) # for reproducibility
df <- df |>
# Make confounder associated with batch:
- mutate(confounder = round(batchnum + runif(n = 200, max = 2)),
- # Make biomarker values associated with confounder:
- noisy_conf = noisy + confounder * 3 - mean(confounder) * 3)
+ mutate(
+ confounder = round(batchnum + runif(n = 200, max = 2)),
+ # Make biomarker values associated with confounder:
+ noisy_conf = noisy + confounder * 3 - mean(confounder) * 3
+ )
-df |> plot_batch(marker = noisy_conf, batch = batch, color = confounder)
+df |>
+ plot_batch(
+ marker = noisy_conf,
+ batch = batch,
+ color = confounder
+ )
```
@@ -172,10 +248,18 @@ Second, batch effects for mean and variance are added.
```{r}
df <- df |>
# Add batch effects to confounded biomarker values:
- mutate(noisy_conf_batch = noisy_conf * batchnum / mean(batchnum) +
- batchnum * 3 - mean(batchnum) * 3)
+ mutate(
+ noisy_conf_batch = noisy_conf * batchnum / mean(batchnum) +
+ batchnum * 3 -
+ mean(batchnum) * 3
+ )
-df |> plot_batch(marker = noisy_conf_batch, batch = batch, color = confounder)
+df |>
+ plot_batch(
+ marker = noisy_conf_batch,
+ batch = batch,
+ color = confounder
+ )
```
The following examples show to what extent batch effects can be removed from `noisy_conf_batch` to recover the "ground truth," `noisy_conf`.
@@ -187,9 +271,17 @@ The following examples show to what extent batch effects can be removed from `no
```{r}
df |>
- adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = standardize, confounders = confounder) |>
- plot_batch(marker = noisy_conf_batch_adj3, batch = batch, color = confounder)
+ adjust_batch(
+ markers = noisy_conf_batch,
+ batch = batch,
+ method = standardize,
+ confounders = confounder
+ ) |>
+ plot_batch(
+ marker = noisy_conf_batch_adj3,
+ batch = batch,
+ color = confounder
+ )
```
@@ -199,9 +291,17 @@ df |>
```{r}
df |>
- adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = ipw, confounders = confounder) |>
- plot_batch(marker = noisy_conf_batch_adj4, batch = batch, color = confounder)
+ adjust_batch(
+ markers = noisy_conf_batch,
+ batch = batch,
+ method = ipw,
+ confounders = confounder
+ ) |>
+ plot_batch(
+ marker = noisy_conf_batch_adj4,
+ batch = batch,
+ color = confounder
+ )
```
@@ -211,9 +311,17 @@ df |>
```{r}
df |>
- adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = quantreg, confounders = confounder) |>
- plot_batch(marker = noisy_conf_batch_adj5, batch = batch, color = confounder)
+ adjust_batch(
+ markers = noisy_conf_batch,
+ batch = batch,
+ method = quantreg,
+ confounders = confounder
+ ) |>
+ plot_batch(
+ marker = noisy_conf_batch_adj5,
+ batch = batch,
+ color = confounder
+ )
```
@@ -225,9 +333,17 @@ df |>
```{r}
df |>
- adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = simple, confounders = confounder) |>
- plot_batch(marker = noisy_conf_batch_adj2, batch = batch, color = confounder)
+ adjust_batch(
+ markers = noisy_conf_batch,
+ batch = batch,
+ method = simple,
+ confounders = confounder
+ ) |>
+ plot_batch(
+ marker = noisy_conf_batch_adj2,
+ batch = batch,
+ color = confounder
+ )
```
@@ -237,9 +353,17 @@ df |>
```{r}
df |>
- adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = quantnorm, confounders = confounder) |>
- plot_batch(marker = noisy_conf_batch_adj6, batch = batch, color = confounder)
+ adjust_batch(
+ markers = noisy_conf_batch,
+ batch = batch,
+ method = quantnorm,
+ confounders = confounder
+ ) |>
+ plot_batch(
+ marker = noisy_conf_batch_adj6,
+ batch = batch,
+ color = confounder
+ )
```
@@ -264,8 +388,13 @@ To obtain an overview of model diagnostics, `diagnose_models()` is called on the
```{r}
# df2 is the new dataset that also contains "noisy_conf_batch_adj2":
-df2 <- df |> adjust_batch(markers = noisy_conf_batch, batch = batch,
- method = standardize, confounders = confounder)
+df2 <- df |>
+ adjust_batch(
+ markers = noisy_conf_batch,
+ batch = batch,
+ method = standardize,
+ confounders = confounder
+ )
# Show overview of model diagnostics:
diagnose_models(df2)
@@ -286,9 +415,16 @@ diagnose_models(df2)$adjust_parameters
Diagnostics for regression models can be performed on fitted models just as for any other regression model in R. For example, residuals *vs.* fitted can be plotted:
```{r}
-tibble(fitted = fitted.values(fit),
- residuals = residuals(fit)) |>
- ggplot(mapping = aes(x = fitted, y = residuals)) +
+tibble(
+ fitted = fitted.values(fit),
+ residuals = residuals(fit)
+) |>
+ ggplot(
+ mapping = aes(
+ x = fitted,
+ y = residuals
+ )
+ ) +
geom_point() +
theme_minimal()
```
From 92d0b18070400463060d6b9fee8c185e4999313b Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Fri, 20 Feb 2026 20:53:05 +0100
Subject: [PATCH 25/36] do not actually install pkg in vignette
---
vignettes/batchtma.Rmd | 4 ++--
1 file changed, 2 insertions(+), 2 deletions(-)
diff --git a/vignettes/batchtma.Rmd b/vignettes/batchtma.Rmd
index 3a0c3b9..3b5406c 100644
--- a/vignettes/batchtma.Rmd
+++ b/vignettes/batchtma.Rmd
@@ -18,13 +18,13 @@ knitr::opts_chunk$set(
batchtma can be installed from [CRAN](https://cran.r-project.org/) using:
-```{r}
+```{r install_cran, eval = FALSE}
install.packages("batchtma")
```
To install a potentially newer version from [GitHub](https://github.com/), use:
-```{r}
+```{r install_github, eval = FALSE}
# install.packages("remotes") # The "remotes" package needs to be installed
remotes::install_github("stopsack/batchtma")
```
From 27390195abbb2c97a87271d690af16e7a1e43983 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Fri, 20 Feb 2026 20:53:37 +0100
Subject: [PATCH 26/36] release v0.2.0 (also fixes #3)
---
DESCRIPTION | 2 +-
NEWS.md | 14 +++++++++++++-
2 files changed, 14 insertions(+), 2 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 9d99606..b8efd68 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: batchtma
Title: Batch Effect Adjustments
-Version: 0.1.7-9000
+Version: 0.2.0
Authors@R:
c(person(given = "Konrad",
family = "Stopsack",
diff --git a/NEWS.md b/NEWS.md
index f0722f3..5102098 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,15 @@
+# batchtma 0.2.0
+
+* Braking change: Require R >= 4.1
+* No other user-visible changes.
+* Internal changes:
+ + Replace long-deprecated `filter(across())` with `if_all()` to be {dplyr}
+ 1.2.0-compatible (#3, thanks to @DavisVaughan).
+ + Split code into multiple files.
+ + Use modern tidyselect, pipe, and anonymous functions.
+ + Add some unit tests.
+
+
# batchtma 0.1.7
* Separately call tidyverse packages
@@ -16,7 +28,7 @@
# batchtma 0.1.4
-* Bioconductor limma dependency added for `remotes::install_github()``
+* Bioconductor limma dependency added for `remotes::install_github()`
* Code review, updates to curly-curly, and tidyverse style changes,
thanks to @tgerke
From 2693a83658f71f99745856e7e918468f1a180493 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Fri, 20 Feb 2026 20:54:55 +0100
Subject: [PATCH 27/36] update action version
---
.github/workflows/R-CMD-check.yaml | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index d83fe20..55ff19c 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -64,7 +64,7 @@ jobs:
- name: Cache R packages
if: runner.os != 'Windows' && matrix.config.image == null
- uses: actions/cache@v1
+ uses: actions/cache@v4
with:
path: ${{ env.R_LIBS_USER }}
key: ${{ runner.os }}-r-${{ matrix.config.r }}-bioc-${{ matrix.config.bioc }}-${{ hashFiles('depends.Rds') }}
From 3d84e30af2b5dc721f7a88330692ff9a1542c947 Mon Sep 17 00:00:00 2001
From: "Konrad H. Stopsack" <35557324+stopsack@users.noreply.github.com>
Date: Sun, 22 Feb 2026 22:41:54 +0100
Subject: [PATCH 28/36] code formatting only
---
index.Rmd | 30 ++++++++++++----
index.md | 59 +++++++++++++++++++++-----------
man/figures/index-example-1.png | Bin 25505 -> 25165 bytes
man/figures/index-example-2.png | Bin 26483 -> 26740 bytes
4 files changed, 63 insertions(+), 26 deletions(-)
diff --git a/index.Rmd b/index.Rmd
index 11984f3..82826e2 100644
--- a/index.Rmd
+++ b/index.Rmd
@@ -48,16 +48,34 @@ library(batchtma)
Define example data with batch effects:
```{r eval=TRUE}
-df <- data.frame(tma = rep(1:2, times = 30),
- biomarker = rep(1:2, times = 30) + runif(max = 3, n = 60))
+df <- data.frame(
+ tma = rep(1:2, times = 30),
+ biomarker = rep(1:2, times = 30) +
+ runif(max = 3, n = 60)
+)
```
Run the `adjust_batch()` function to adjust for batch effects:
```{r example, eval=TRUE, fig.show="hold", out.width='40%'}
-df_adjust <- adjust_batch(data = df, markers = biomarker, batch = tma, method = simple)
-
-plot_batch(data = df, marker = biomarker, batch = tma, title = "Raw data")
-plot_batch(data = df_adjust, marker = biomarker_adj2, batch = tma, title = "Adjusted data")
+df_adjust <- adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = simple
+)
+
+plot_batch(
+ data = df,
+ marker = biomarker,
+ batch = tma,
+ title = "Raw data"
+)
+plot_batch(
+ data = df_adjust,
+ marker = biomarker_adj2,
+ batch = tma,
+ title = "Adjusted data"
+)
```
* The ["Get Started"](articles/batchtma.html) vignette has details on how to use the package.
diff --git a/index.md b/index.md
index 39c9b5e..ca02d24 100644
--- a/index.md
+++ b/index.md
@@ -4,6 +4,7 @@
# batchtma R package: Methods to address batch effects
+
The goal of the batchtma package is to provide functions for batch
@@ -40,23 +41,41 @@ library(batchtma)
Define example data with batch effects:
``` r
-df <- data.frame(tma = rep(1:2, times = 30),
- biomarker = rep(1:2, times = 30) + runif(max = 3, n = 60))
+df <- data.frame(
+ tma = rep(1:2, times = 30),
+ biomarker = rep(1:2, times = 30) +
+ runif(max = 3, n = 60)
+)
```
Run the `adjust_batch()` function to adjust for batch effects:
``` r
-df_adjust <- adjust_batch(data = df, markers = biomarker, batch = tma, method = simple)
-
-plot_batch(data = df, marker = biomarker, batch = tma, title = "Raw data")
-plot_batch(data = df_adjust, marker = biomarker_adj2, batch = tma, title = "Adjusted data")
+df_adjust <- adjust_batch(
+ data = df,
+ markers = biomarker,
+ batch = tma,
+ method = simple
+)
+
+plot_batch(
+ data = df,
+ marker = biomarker,
+ batch = tma,
+ title = "Raw data"
+)
+plot_batch(
+ data = df_adjust,
+ marker = biomarker_adj2,
+ batch = tma,
+ title = "Adjusted data"
+)
```
-
+
-- The [“Get Started”](articles/batchtma.html) vignette has details on
- how to use the package.
+- The [“Get Started”](articles/batchtma.html) vignette has details on
+ how to use the package.
## Methodology
@@ -69,17 +88,17 @@ characteristics that are expected to lead to differences in biomarker
levels. They should ideally be retained when performing batch effect
adjustments.
-| \# | `method =` | Approach | Addressed | Retains “true” between-batch differences |
-|-----|---------------|-------------------------------|------------------------|------------------------------------------|
-| 1 | — | Unadjusted | — | Yes |
-| 2 | `simple` | Simple means | Means | No |
-| 3 | `standardize` | Standardized batch means | Means | Yes |
-| 4 | `ipw` | Inverse-probability weighting | Means | Yes |
-| 5 | `quantreg` | Quantile regression | Low and high quantiles | Yes |
-| 6 | `quantnorm` | Quantile normalization | All ranks | No |
-
-- [“Get Started”](articles/batchtma.html) shows general examples on
- the different methods in absence and presence of confounding.
+| \# | `method =` | Approach | Addressed | Retains “true” between-batch differences |
+|----|----|----|----|----|
+| 1 | — | Unadjusted | — | Yes |
+| 2 | `simple` | Simple means | Means | No |
+| 3 | `standardize` | Standardized batch means | Means | Yes |
+| 4 | `ipw` | Inverse-probability weighting | Means | Yes |
+| 5 | `quantreg` | Quantile regression | Low and high quantiles | Yes |
+| 6 | `quantnorm` | Quantile normalization | All ranks | No |
+
+- [“Get Started”](articles/batchtma.html) shows general examples on the
+ different methods in absence and presence of confounding.
## Reference
diff --git a/man/figures/index-example-1.png b/man/figures/index-example-1.png
index 5e8d0f7378c80450fad5837b7b9220b522ba2959..20477cad4dbe17dd44a1e739765b6fe9ec8197e2 100644
GIT binary patch
literal 25165
zcmeFZXHZm8(=9xVgCvOxC`ynhSqVyJKr#qO&Y+SJP?BVZq#zO`gCxmGKnY3^Fq4!l
ziVCRYEIIS_0eznL=dJpxZq>bYe~8vNoH=Lj-rcL$T0P;K>WZYqbi^nWid0!iP78%X
zXQ5DdBZNoc74HYbLMRl5%<^_*&;t(5}hpvowtoK*thw-;>#oS)`b&>jY@p4Bca3!$go&Np0X>e|V>HW~a
z+Vk~&o0&3AMZB;dp?+?;>|52^5&dC7nbtRix)A@7dV53H
zr{T<`TT$_)?q=?Sz7!jZ;b2eqR}yXS=jC>)FPNthyZdXU9UCWSK#lz3;;ryrXJM$g
zV#gKM9;t78@xxNs#*1+pYqZ}yuZ^A9K$mEq^?ezh`H0_3lzIhHZ}p2Zl4N(yd7qzn
z&n9`-EbW+E%X#;m)h!SG^c2Ehv7)pjbPw40^{W5Kwgl)`bL3K(*SG~l|N7o?3US|Dg*C{tv8#;Slx0b=abeGQk*X9!g73LF{jE?OQ_#X+F
z`KI?=ct)s0kJzuy+4!yEvZT7DR35MH&1;o8*@nV;TWH31lwYl9J)zA}
z7hZ{TmZ6fb-@LX|)tIEd8$xuA-`nF-$?-^9)B2+yDx57#{7P}OMXs|>rInw5O6w>j
z9_8@O|Ggl1MBV9k2#d1xDdteedhd|A)qD?sRms8f!1u0yZc5-^_ciwwk~$&%`d09%
z;fThovw!+rrk=Hv{d!=RUR9kseg94C=1)t+mN3q&8hbji4WYjHzBJq4-u3odKGHXgdXcq%WS`BG!c_@kWsUq%!~Z9fiW8l;vb}Jn`m!
z5G0?F-`e;3B@&$w*F(9SUJT|
zS&5Jdhe99zc&o%Dsxa_V;D5Yy5)Z8@Mfl$rE>JU}A~`S7UPj?zD6q&6Ce)bVf3Cnp
zp^b;`N_KI*@b`b=EpUM&|MLPZr+Z66Z#@krAa5Jq-@B737#kPICM>L5tY6M|=ZAaJ
z2~MrI6UmJ(``PJ8|I-{o*~nAlP=^2vp!eZC|U
zm+*LAV?8e2;?cX3?$hz5#?RG1WlHm}PWQ+0H*Uys=@gFmOAqlMC1>DaQFss|X!A*F
zsNvOz^~`|1T8YKxBhT(E40-AY&anJ@oU7D?Pr~_%<^N=EzFZ&TZcpOB6Y#{aQekcO
zE6e43WC^+@dMOs2`zdRK^|yY?in>k4yd=J2M|rR;SX7`1A!b-EDWW6Cy-l9hlr{J$
zsX%n?#HaKb*OVXiciOjBXU4wQy1Klh*7vN=avE+lC8p+L%~MVC$h@<$cun-y2+mIV
zV-kO+PZL8)^tNjJxt~KxyPO`4d=cZluclW86>#Vc3lH)cA{%uypL)@4#_K0DI$Zc7PBkxV?J(RcC=EiPn3bxPeeG*?#s&yI5^_UxUY3F_3qN{IewUi<0QzN?e_#Q)S
z!2b5I%i64=cGN=_#njzbeIF*S&;DI9JPZ$?#U@H;TgyIcYOCk-Ygi4H0`+uO7M9EN
z^$Sv${tWM}UE185DIaU+Z;EG8jA|)w+T)V+UWw6=-qp!hPq&ZiMJ-Nb2GEeso7fFC
zc#Y>}%&5Rl$4huEDSe2dhnmyqkma@hk?3bw_{f>(`<-q$#OOC8Z^
z&y#sllP0y)!ZKcQXV_2Qdwt$<>%lZO_w41n$;kb>(>MMnL~1CzWR&#O40&Tx48uhI
zsrgJqmNzFdJ3H)!h2C{~f6qwik4@Pef8^|~%i(Z7u9d~%R?OF?+rJemCB2p`%jd>F
zKC<~-c=e5+Q@zKWfOEqig>{cL+9L#H9Ub2BnJLs1J}=!r3J)ajRHAJv_e;tf+h&Je
z`dEJ+VJ`pu?C2RU`hli^{f<%@)%hDiPolRmX8$?O7bz=z88!VL_Lm+Ay<#)tK>>Rv+LE4IDQURbr$Y>p+FTfVtVnqRFGE$X$DE4s4O&Yx6@)gJFi
zxnaa!y{gR^J#%BRX>YYW`$>v?1dT9n5uB)OZlBDB68-Xv(cdMvW2a><4f|G}IzJ}%
z7~YX>KDboG`+V?{``lnPdw`@u`#ZuU!)y1U;Quy-R!
zckV;cuDr04nu>}*;A1k0p612R#?J1Ual4^M-wZkn`SjT3NW`s4&h-{VPu}pD`(3nf
z$+oL`<;ssb_Y%(Fi9Z?T!B3(mL>mgwH3oHP)w$Ef
zcsWr_<_(@t?QZwkrDphV*~Nvgo$~H;ad#9E*8NQ=8Ot%Udf=GA2T`7(16eBmZS%
z1Jh;?yuzUAkU6nIQ^q$As`>d=?TM{_et*bw9&YTZveWKG7jHLhj-TGO&bJicf3|1$
z{l!xvkyXliGKLZMKrhws5|3I%=Sb6lOWL(Bf3Js5eg^@VHv6D>H
z+BsJ!QT8HN`G}@p)!qlj0F(2me#PDX
zV%@lJ&P}F`+w?%gzf1Cs&G%9C4hpe0_-^CLr6YoZ4pCPs`AMV%
zXD)cllV*pD$#0{P1rzm%I?-`d&l&!T}|cI-BeXe(li09
zK3}K$xr>$j@mm-!j)k-TYKR_H#jMfdPq}$RnM}huC4IH~@;z)q?ipI?DCq*Xtrfdo
zs^TL_=qHqLj{OCRhh2kq4F)(0#3~s7jOt6KeIB#gYGang1{`VY-iovR1b;_ibQ?xNhO+
zJ;t@c4@c32avC_4kax}K10+RQTnNljj8?I?FV-l!e}a=4>Ic=Xx-2Sil4OxkYU{lT
z9@7r47_pzGoJ%(A?A!&zNGk@z7l+Pj6Ep
zP!LCt-xklw!ln|x$hr0N&AmK;T~gams(u4pS$ajlsO0TkmvCR+kfJ34g=bs#MTF?!
zqoUB8OJi~TO*=g&Uf%lodR}q<^u?PLe$gQ)^i#@~(~sJ3Q6I);=~>N}-tJ$U0-Dqw
zUHo`rGSS$H((x2qDEoORNhxCP#U^gu@&a}}teobJzP6`sRpsClV(I<1izMLgM$`gh
z*;Lzt$)zSc+7z=sUwL#fVXQqVVokMQ=dmol;nVA(2qvLyN#{vpYG(1|_--{Fc92uC
zAs8AylXI$xJhurl%We`)no{|aE5MSTWQ#gj_=?EpXz9Ikke>v1#odkaK0dniEq~VI
z+$2EH9kS8Emp;$w?`xdmWroW}d?&mgh+H;C9MA8Vgy*7#;EO_or`MO}o)3Fl>C}5T
zetmH}N2ge))j|4n?eRkW@{BYozrh(x`euNQkl$`HQ7GYPy0Bsm_YStx3vVt!#Jc08
zVWoAX#$r+E07_S1j9Fd*Yn=2Un$cN++n^#VJ04x?RGS@sVtk2HGYdO;
z>h|3m$h!K>6p0}56o8km&Ka<`gLkB>qeHDSC}!8e;R#(|k+$MWLUE&S-p_I!HDxFR_p_
ztu_{avIobj@;PnDo#`
zjvk9S?cC7aaiwUc;m%Z#JTwgQ1n*-Yx>F(NeQw~7Cq8&9Ch#aJU6}tv_zwoOM9}jO
ze%bN=S_SNDLJAyv$CF&K(1lcJsSHTv)le`8GWY=6EE
zeX~bj>%D5f39z|CD^KipKvmCFZ^0GAN`j8<^#$^JvTiQH&G!n;v3RyOi$Cce)`(jf
z%^QDyk3Ids5a+wO)QThto6d}U`D&Sz>i)>^H8{K1B!0xK}ki?A0J6}gPIP~MTEdjJi|z)Sim+9yWUv5}`Qr3NbE
zlhDNb%-3>Q3^yy8XYQ18Vo~-Vw`GeCA0zAqE~;EfMP;e=(%^4CzNQ_QZhk5>A0ET2
zde)7ogLT3N5lymqIR13YgNv(1(g}a(jiGzM4jHl_(|=klWP$Ll9ApOe4igL#TByJ;
z#1gB0+AiJlwBc(n8)l!is|9HAzUs+SMZj2WVj?>v3)(qaqekD&GbHQq&@K5_9?2DG
zFrTI149qR=O-(sqSwYN3Z(?qT^_%6P}ZLUrAD(0)DFmtGgVDceb+NYT@~EsJ4sJp#F$Q-z4UG+VDz<6Xsyt~=Eq6_G&G+n|A#?PJT$4b5kZjbC8duH?FoG^0MurcVmGdBjxOoJs-
z-<0xY9l7h4)x9?N?Fvbn`ks`_sFEOoZH3kT#_#AcXwF@uvj-@b2l?HgWy#klDX
zcrp2E1H-03hwRGss%)F!m*Ebk&{$@6^AL>%N_8{Hyr|8|$IE`GbJe8g-nc%EvB6sG~X1`W)6y6Xe#1kh3+Nr(k
z@#Qi210?^UT-fmt-tPbHNu1!*hYZNu{Z3p
z8E{TN>U6C;GkAZ0i0*t!@pTyL1{JReQJ@^fm!~)G2|0XLTS;n=A_Zs{%41rmvqCON
zh1R3!%B*DtjMIQsQMrHPcgoJ@G65#i8@1K1Z{7+tQtpipQTa@e%3~n$YVQ4s4p0Jq
zWtWCeR8-WWZ+4y4gW2M0;5VoGD^GY~r^Ps9hw#_Nek-(b1R=#YH#24L?2xOP?Kmy-MRjjeiuqxE}Df^R7JNhX?#t
z5bjRpgX8vU^ps&wvEAA7^2VyezROoceK#tq?ECNYox9LTU_yojf3pz+-s%|R*=6uo
zC>uTntIGauOCKgn2NBT!Js|MLukglCjr)BJ|6U5KMMaO?Q2&Ha*De
z#AB=?&Gq-+@BvIfQYgGYlIN1u44c4I?J$tIJzuX;Ry$7DO*E^Fl$a`vn|lp3II4h6
z^W0}XF9CC?FrzgZfkPmP6zeROX${S-Ulz}2TIXr%`ToHvYoc#bep{Y%Hb}<|ojI^b
zR7zi;ZvhanN$eo!fP6r}AQ3vI;m#}H$GJxdu@9zL0
zW0{!Il{TG;@VO}?+JYS&9YUMa)1j=G>t~@$C$TCEY{thB)AG$8$Js~7{rUBd{}#J5
zl5QW9%>Jm?cZ>Qp8_pI~8Bsyci-+#AZH?tepCv?i?AX}YU)7ccHA<7Cs}Q7ErXKWQ
z=g!Jhui00cbDG$)#x47U_Tzv2nk|Vf^mJ1`y2U!2MTV96r3OzCzMMjq@63G^TAj8_
zcR(SMj_sxQD)9!vO@K?59myB&swNSA{@kUKz|{&}Lb|%r?hV#vrreV
zl-;kC4sPHQdaF+^Gk3H-Ndo}S*}=#OI423_0lWQ{&uuyx5S071yz-{RO4nuU_mQC_
zv?@=pzsq$PC^OTz=AwP{lMIsqJfcOO$U)$s@JL7+7GIO(yV3dD_uy96hXay0weyvF
z^RJBE=Pskzy{x|u0H*htY~6HmiCmKs;aek|;^I@93Gc&NEGL64&=m7oE7ku5lyDqG
zM(T+6um~KTNvSK2-)a)H3)E8vc7m~`P>pZ!ep)?@?XVy$AkbWkS3zc?92^{~(Cvb5
z(eKdYWNB#`7Il~7ECsX+U7m%mR}M26GtYV3E+&LgojOu%I5WabZ;%`P&AmTo78tPr
z?F&}YQI42wWVcCmmyjk7_Ub)aEHakQ;(}!h#njNtdW3BUs7%wl;JPx|EtPjl+=Bmo5x1i0H6}L=L89y>h0`%<@S>(27{?z_`UoEwW|IUk)8Nuz%
zx&YX&uXidVae2_X@w$tjJ9c=jS45*9AL*}f`-p)C12Bze+o^`KA|8cTnCC?Bk5m!C
zGo7-9u#E*sTYfwvV9)C)vsZY;Y>GVnTkY*^o=to8nDDJR1lpe_0{G~|5TLk)+(nMl
zr3CqZFa6?&4d3^{65Tv^3mP0~vZ6^<)gK;0rm!a1r=0~nwu+Qy?(zm^V0;tm+-Fi_
z8AWb2TMK8sIE>)8`p}#-?wS0sM%N83Rdz~oBee#KPD&5?84*nZ0fC!($L}AOTI>!i
z)^D561jlpdufHWc?(0osIx+G{&?ahv?;n{8)f(M}5@`sXg!kCOJK
zMA3hMBSW?Y4`Yi*ihE=Hcdb-Lps5w`-M9se%#AWWuknskr-bd}OHFD8^e0@Y&`(c*
zrkQtf;C@{*NP(P}?;d-|k?E_I$Y)BoI&&LY0wygqbnk-8$A)jcAzhj&Z@L@*K3HS4XeTe
zgOICVZk8RObjf*0Vv
z*Fdk4t<91yHqb-?_R_Sssr#sI#1
zmL^n%cz*Mt>WBqcOq>~ziJ=p=--u5c>Ae?1f03%f5>+e{a+H=8
zvqPkoeEwDG8-LLGmTUY6#?G%2+*zvvB-m{iE2J6e<}y~qTUIvku<_{TSqhTC*?<@C
zex;hQ*aW5pX?6+{N=XWrZV$FLF;?N7>$-U%krdmRDee42;{;noNi(3hx75CM#J{{T
z`%p|DVf$7|bj3&Yi*=N<|9HWb2V7oyB-}%87%PI_3+kd|QAp-lorql99DV8|^7^T$
z{z8Iwfbf^-aZJ)YWf9t^-NA67DrErPem0p*X
z$*Z;NLq2&>8py;k#ZZ?J(`ro$#}3quIm$@9%;3U?&}
z_IK}W`EJ=%_6-H$lPqnYy(fOX6(w}<*|XDl(Ah;!-q-+&=+WHkP!iRL%<`itHW?Sb
zL8!Qh{4`$Du
zOuHQ!&MKs(9HYcC9YrI8IA8?)n|S?F3u4Uv5XAFTdHMNm&=w#9v&X{ndNz=LdmfMc
zK>NQaQ|dsG7gaOBp|-!UzuR=C^Cu~WnGFsQr}9#Gcr)m^tk6xnlFJ%BvkbS*>;f8X
zS>H!(9z`pI3REJV{JqK`={=g9=TKFDG;^y;aYDH9*M05|mvpgR6LoE&vE;`ePM67C
zDKn)DvmaoeqsbazI!kc{Z}4kQk?Ig7?PO101TZBX?kn|WxD^19iaSP@2v7W0yBVTO
zp}wA>7CK*%+}3H=;$X6v){GT9s$M2*{cGZMP{+D8**Xh}Y@1Ugr*y~_qt3K8K5~Ra
z(OiT68mvLTX$Ck-n{i-r<2%o{DI%BUkC(P3MB7?g7R9GtPIxsaFi}CljSgYP0`cQ#tgc@&VpRcD5zH8Ep)y8(gG
z^lQ4*l?z^=-hYgJvy~BbzZ@9mDb-@ZAj~5Z?bfNR98kQImESxPvb(M@K=6p*!zQhW
z$87Tcr7tryN^SO%ZEntg(0CCfljP%7RYHbM{x9RtU6rj$z$!XUPkuDuaCGN|tgV(J
zAP@D2W~SX$62r`%wVwpa-ifC>ErfWOJq-&05IuG74k9~*UFd0oP}9%UT}^!fn`5+w
zRkjap+id_E=iU7J>@)9qwttQhMqpGKsXV^OFfs6Ec-(sYfN)d+JK!+??RmMx-p21*
z(R5d4)(X!**QW-nOH90{@({%#EGlMQt!&-hOe4LiS8(c$n(2efoy+Tu6HoWdW+bwz
z?fV6yf30>mNX0^;Xqps`0Qq>~#2P-u7f^YtZhjSLrfiT$wb>vw_|^*?nSJ+%#RP7B
zj~8;jCFQ%K)y5Cig(ntTKLN+3hUZSE8L&Rxw7*;KHrXZP7T?gcexzSLX|n^o7`2TI
z1u}*Hz&QD2$Bo6w9$!3DKrxu|Ni-GJC%T@R%8ZaL1_oH;?D{(k88cMVj8tN}rcSYB
zw|Yhf1gA<_$^tbqw!I>;teE7xoK*Km4>SZE_5nty>w(+!m!~;)cpo#$=K)KS;sL5J
zjO942X_PcC0L*s$-V~pI!SY%(!Zu_5yK8jobLdO6qc^u#rf#e$?p(JF!zMy&w0ZIn
zomd<>&`LT|1atIyl;hZy7lvNSZ>Q=NX}#;
z;k;YVZ2-^P4kh^Vr{Lph+;EZaJ{lhn;Z93JGvvNJsDswZmIEe_q3~Du0*zJM79E$y
z0dM^ew@g8(I`ZSr@T$)qSOO;RNQ+ap@L^G@rXk3&L++mNTuf~2_0JLO5=#0mTp*Sp
z=q4ry5?7Q^!{I8U(0iDz^#y*qTi)!Z`@UUSG8y%rk=Q6W&!Rnxf58Lc#__K3>48$i
z8+vM7TDf6}0OAi3&8b@)-}L6|JzO+yJ(FL`x$-}n=_2(D$ULKqKi@PuR}cg;A?O1g
ze)+#xni-yiPE^(9Uw;cXD+~CIkd^FE??Y}sEDhU*Hwo2%<&N(Ec;l^7c-hv-DAjjw
zb0R7<`xl|XRcJu)wVqy=?-t!$arjo_sK>)@b$ID!Ri?tLiwG11JpGuuuRV$s^NXtj
z2pYTY>`?tvpnqHKl|zt!22tt|9T_My6{4y%AIk~=yW_Rg$}t8KSPZypOx!ZLry9(K
z?fY&yEuvMJ4-Ef(~tEuRR#zqg^lG>OtNmI5|tmqudS#l)Zd
zx=j~)=JOL^M`SXa18-9z*{gOQ7~qfhnVxdJMSKuIVE(FWoNyOLJ%JEXYya)#SJG5w
zR9*bE(p!|15h68a5-mfJwMC+_9i*NqWqQ%+`#Ex{jjY-4wc69*?`Un6iMCR4>kR8#
zsGBK^gL{`}I*-w2NE#E0PDbCJiz3-O}w-ChJ5MmaSc!1
zEb}IR7jWiy&bc~101`!NbM$y+N_@O!hBZ>tvsI@=&pQ67OW4;?f6?7XiHEt#Lg9ug
z#=D9DNx)yRC1rJOI+_IH-~G4e4m8>Qyg>nxt;aT!m^K!($X^K!oY{(k3`%
zC#?|>%;461%ny6Rz5a!I1AbM&JV+^OjHaYg!gd4kv#(4W+NZpBA(q3*)R)w@P7We
zx$v1ce)qI1ncQAz;s*vW37Q&_uQ7GvO4x1LOi|bOCGY7vo)|oR?U0#p5uAmrh%V0=
zyZ02BN*O^alT3S1A(__Ergd%)!Iu@mM+^+RdDx!0>uyv>Mu~@4_p5s8y^Wkw(|SDy
zpV;^5Y1)q!a7a;KBS2H0AS@;?G}9n3n&htl_E7J0xjFwJSMf6G^Z>A5UTL(i9?9x`
zruGD!f&kokJ~PcY%#zRA?4QBU
zkqBARf5e^Za%XOnVl{X4S~vUI5uv~=W2}4C6Dj+KpLb7>BE~;hV)LG!9?LHk*Ln&x
z7z=8zn%^U~4ba$SnFTZ&k5~)=D)R1b{CNC4+@6KPJ?LQNFg8>sJo(kBA@VHAviz89z8=x;WyqogI-2%sQf(
zzFq-JeRlL2_pKCKGfx?#nwy-NerqWvF0E{|RlOgmq<2^1>c=L7IBr4BM+zl(AxUbwcp19{|1
zS5bPw{_a9RJ3x!mIv3pBMnFTn$*vawRz-ONxk=VRgh_)+G0QVTM?$Gaucf6mTJz0`
zWLz^!hUq6h`e+u?qL#e$UM)|fM$TQu`2Y3+(>qAf2OyN
z_yHQ$1|2Y2dlg@3<8dB)0}pcsmp+*!$~?SHariiiRD>sS(@T%Ux*M8u0}Hu-cuk})
zbt^K~UVPa&)68ai&nWIw<_kb-$HTno;~g&-RCiqsR=-k~N0n~0k6neLne>rSaBx^f
zG1I=JL7hlJmgywkZf&d|tIPQbwyg_09m*@>YP-MKMn*G52d%WX6({21F@SG%D?j~2
z=9)S;!&ueua<-Um+^5sU(mVWB=n%6RVEHo&#HWRuj_nLP&Alzp`6N6gHk)cxWqbDY
zC8rhf*L`Q2r3oQaC*HS88%vyPFk2G*5D^T`1^dD0Fp3S}6;6bMJ{S*Xm^R+$-(yQo
z=R~VW@9(%)UY9F3duh3~GG%Qh6p%6|ln~l<TeVab_R
zC@67D5JJfx((5`C+^FwCtNsEu&VjqQj3MWeOwyn471x3{tpH#rdDDZth?c6$iWSH%
z&Mn1A?kN**&=2l7viGfn$u?|Hx_irBn%R_0KWVuaCq$=QegO!ij94(b)ER`-X%Sbm_t;cHUF{A{6}^h1umX1
z9m^$wlVAiIh4Ug@pU(|BO}SONIun*Rf7Dff=-xA@^4+>~jURw!u}%>)2-~aE
zK2KO)vQ$m^+?1yfnVTzYCOEitw;Yx9LI@N
z?PfXPnm$U6vzCqqklSK7BE7#EjW~*`yr&uc2DcCwa|GFBgj?Nss)O44$vN$rl;j)z
z+Q9>aOcZ4}KE>kR3CE#^uR_*SCnC4*xdW&E`eW*{BYVu5n_ljF@1T>oJl>h<{J9Rw
zS5b<^9%2KhHTS*nY>{pGf~>gi(2kVme37AwJh8ukG1
z`uUf-LKU$01SFFi(O|@Chngi9tdSo07?^Z2n%hqO!S;ClA7<-zegZ<458#|8P$KZ(
znR5;h^n0)H;N;ozhNb(i)J3A`t5jW1++_)Gv)zd9)F8Pq4-^3F`$D^Ko`F?ue%%v6
z@gdw#o6B$-O0*{lE(;6G(!0>*$zvK~{FkD-u`(dk?
zW+0HeYqYk#`H_>fh+8eJ3P)K-TXLsmFsWIQ4xakn5zZd&2OCcvrLEDx^=kmiC$SuKf0LB
z`*-9exe2zL02(K=a)7%TR4dB;I92zSohZfaU(C+t8mWY-cD3zDEtt~HTR+;8eQ-nS>fo2$rp`OX!*^ww-oKpz_-
z)6IH-aBDjF&`UYfc#8(?HEhy9`6gD*@5hzXran$daR4_EjYl|V*!jiJqrJp=yjRQk
z>FG6OL5$kc_gk~NxEh_+8polyr8SE9-@d3@fN4np{BOPIay2;I{}7sDoNEc~jmfVFz&7Rg9a|m0i&((PERN
z-hSyJWkLf}ra~Ic-CHv*kDd>Cm!xEgiHaIZG!asj`aha8lybZ_AZp`2NT$!FnT4~9
zyyLixsE)UaaeKafnc*4A1!6s$b%$30kN1AVM263(`TkbYsh8HPDpwodQglhF4NiQa
z4w=?+6gwSUHh3B7!_xR}4g0O9aKx4=Fp)~RThq~xUY2d1M(;Zq)nS#ZHIzX`QbT1O
ziK@cyOK>yQ-wLb~+F($|OV7UnN9$onYy%Yfz$^^7fv#^E*f*(|ZZ14jWFVSuFrd#3
zoB;lFbgCdT5um6MKvK)Ll4ZwOii$ddx9K8gZ~crUR!Zzo!>7%n=(7O5OtA%G-7;P~0k^g|
z0=GDdabyj=I)yPBt5xQRo7C^17nG$XEXLQVq!)85JZNE{D@v{&mfUxpo+cTTjt{|R
zlT5xQ0R}J*3bXb~YqR)50}ITnTv@T#md2c6%R0h<9B#vE`vm=$FqJYItJ<*k6f#8Z
zaR_L?m2Sn4_Z7h^wOZSgerJf2(QoBpr@seNj1s!;MB$#+k{a-#@nD5XVY8E5U#0_G
z#hn%G14aPH_mZ~fi*wqmb=^SRtu?v^jgI^?u+lA8$_ALnQed&HNT=%SMokx>yzd%I
zq2v|_t10KE)3}E3%9g%jxE=ZfN?MtzuVQl|8Ucw28J=hrt@ro?obmTklM7P2zn7uK
zKDMdF=)dup$E@MY&YlG5k+v*RQdd{fl-WnsL|-kT=>#95Nubf3wtQZ9aUJlp^XGs8
zy-aBRR8lWFXOuJ%9-)o6Ls^l5`&2%W>fBlBR+kxq`+;wbRB`$HWw^55p#qJ}Nwbje
zlof#p*o+DqNU#$s9`M3~2;(r$=l}^35%sd;`s2WBq?(E3pCuSqQDDF95*~(CP`G7
zsgOUl7w+Icqw5LDV-3-Tm)jm23@ywWoDprDvI(+)6pz2H495u
zF{eJgEZA#~>XcZ1pn)phO>Np;Grm_WQZgv!M}}8eIbW-5Un1SPy3PtEtpPvC8&h5|
zi91P;>qA;fe_l4#M4{s@4(G)zno$cTUuB!Gqs1aFujK~FfkhS}oe+D2Of}B_DdLgb
zgmi!ina})py?YDtz%}Rs`}4lB&xA6N0#mH_gbTQ!7-UXFGU-Xp^mAZULrg>tJ)6Mo
zup!umwofk?*o8;sO4K-Cim=dTP$lx*rh2knJ*Ek@O!Yiq+n>G=I?$+)BzF@HZe!3m
zQ$F?7ecFZrqPX7{EPkw+i7<}l02546xH5~}Bja#Mo5X(ZtJxv$#EQy?_kJrG+lOOXFhaH&EpWmO})PmWGYS35ruF
z-NG}?@Bp{txH2|KKVkYT4tEdtsqnFO{Zzvj9{D-EH>wy}3P|Wzu3QmzJvW_|JF7e}
zPzFx^wI4TGViwVL4{Yv{(IeeBk-0t>C?ESA72wD2zMBkYmlpav-_MHa{ep&k5lVd+
z8gCm{2B)bbz0eEh`g9q*MkzO8E+H2J710g5{U6_W|KOorZ+@t!A(>+H9(s^CyLI))
zo^bY2z#vN4(#@SDw6XOaRD$^*;!5+IZnBwI7l2PZGcuCxnA
z0wY#dKQGp43R1~g%*uuOy`WY{;`N>HW=YEswM%HRquG#`U0_b~&iGbc9FV=X+_r0O(|w%Rq-JNifO0cnvLLYqy7frr$9h863S5obP
zlLgBD)shxCeb_at?R8Timm?D)^LiZ?L7Jpe-<~@=`^&rp=f~1WF!}q^4igYtvboQL
z8nt}Ve`&OJoT&A8v^1aJ)=K!id@Fq*<^A+ciYk`>$1Y$6pofIx(`;#-U;)rq6}LhbM94v
z`@-2n$)+L^$Hw@|A)p>eT?gfN6+&Xgv6j4e`V`ZfB15j@1f3(UfGVndjb_$rRdaHZ
zQYyPYz9M-F=idpN%0SDiHE6lQ}e&D6O3+}tnYuBF;Kl7HSF&%ZhaW{mM@
zS7X0l*~SK{N3)E^ZjGBg;4@rtz2ovDr~LE~wM1}Nk#<392h5JCR-R#eW#d`%RK$4S
z`OiK#O^*o+9lETu%v9h^NO@+a}^BlT`{+;ca|wx9UL4=JzZYFYjp5;C|)xclWbPR1$a$
zdkWaMDiZ>Oqo_x{M0*m(i*-w^^%8qVbQn9!jjG#Lik_>`pHp0cr_Qla
z{Y1nXBe~I6ln8)suI%x|q^BvHy!QEak2=8-d_)$&p~bMvd)o`^#)78v8^9@z+Y}wm
z60xnCYMF?cRcsULh7uJG)yygb$~O!ILJV
zgm}_wC=Eb-t%BS;N*yIcm@a;&j4yO_5QU;IM!pNc$D#*}sK51W<}|>&f4-MoHOaHz
zPQ5C#|1!DciRS{huN~GILhR--rkG%yN5(+jY4CK@FD<1)g%P5kSc%Uqfs=(;d!MNBOh3!y7?PMQESuTs#D&-Dnji&YghAcV|{~XvfkIOe#m*&sZ
z>&~n~=h6mR;KvgjVuM0Jjk(QhK?f|SVFA`IWXcn!!U9(KWIqvpF#s#o#YV3al8tC&
z#tY%C0fo{I2VrinB3wH%=fkvxuG0+tL{M6(KKt#t24bE8TZHNge+M;1#Y6CLUaT_=
z*V0M%+!*p%5auSTcKoIn#UN?_t=JU4YCsDtir*45!NZQV#Fqm8@05^j%o$Oax0uXS
zADG9HBftRO2XkB%Zm%;){{n@G-&R%1F;fm%Y1i=wFqX=%1s)}szCRrudSj%M4eqB7
zeVjvmr53kV|GDfu$kIz*(e#8Ox;Prmf}u^3{_VUZvsj2B7(S24(mpX&g%vXtL%>Ky
zJr>H0UvU;*;of(av$*Dc2;_a6S0yJ25wm9%b)TS}>Ll2crfMUA-ofE`@sauC%seM|
zyfA=jU!O^xpaWt7ym+dKzPrSe@z^~h^x$*ng>_F4M3nZRLiA&~q{S=h>NJj|{WAB(
zmP7h&+AHENSY0`Ww}oSW?^kt&w@T3exzmqP7-G9va?I|P@VU9}oO_PbSQ&(_!lJLF
z3Wd)Q+XlA1_aQWp2PYiYo!2sBUFuuF7ACWS?
zh<^^dPsg*!3zz3bg%5_Vl`Qbmcl^QIKxA^cIFS(4*zjAlMMn@48uuID?JlM1?IMr@
zT!5);19gQ*3=LflHH6j^PyKVp54q~I#C=bU^#KmCz+4skV&BvpOyky2xcuH;^phFU`4+ExU>wKT3G#gPtLG45yLuRpk0V2E-SUaP+C{2dL<0N|m98?0!QP$FP
z;CYVa(8!pr|Ej<=4|jLF3ID`j?6`;l8XhrwYh_{a?kLd3+naX*DOtA;)`}7^;uwvo
zYmyya@HIuUfYp6lX-!ifDe#{=fr%UO?VJ>K;X-#Gw_xa+{2D8G_%IiacQw00eaMdf
zQN3SKb050`mm1&3x=X0Gl^iE4-Gik?$MGg+w;~
zuqN`p`c=w`D){X@XY$Yk4@Ov0GicglFnuj}#0PLQ8z8SZNKc-mE8+TIU+|$lXe}sW
ze=IPK9Zc?>N7J2Qkn$B_rZcX&`4Hw+jD;5aBv!Xm1Z_q^42V?PvkZdHh7g;+N0haP
z?Ab?rH>$NqPIawnqB;Ph0Tf2U7?D(|kstWAvc0_6NwGE`Q^yVhIYt!^=1yPVraRbe
zVT?)=A20ZVtItS!_nlBm$~pEtPixzEOZ5TIdckWcj4ZeeNq;;LJa||h<{0R*T3cJC
zcrzA2Xf|M*?Rp%Xo0sSHgfDjhu(-Furi%cYa3>QL9J9fCkhs1WWr2L29sN_X7uYAs
zORlS1@JMlxPud&2mOJca5fKoKmje1;^1oP3$0iY<*GQC{|7ajq0v_^rQD^l5wtr1k
z2s$$nF;dpq%?CP!m>hv8Fxae9_Xictl%MyHsqr;sLe?Nu=BV6@h`zJg=)u>EK#7H~
zv9u2rQ+1O(JOz1A&SL(1E&F62&1>*Mv_n#26SVFKQyQ{_Xziv12f0|;sJ(U|+?8sM
zLZl_|L)qo47tG-AF>pgz8Hzs|(!QU*fu?ht-XD1Bv!3eS1z!~Ma3dlAF5wpwNIQ>)
zvG*Ab1A!=NfiE2SXxUpPWeG-U!6#fG--XY?5iwq
zDh2|3f{XzA-Ou3AV3-vfF&GdC2opK*A
zov4%^6bG^1F+TutuUP>HWc_STz^W~|(W$%f;^yJR5e08xGA2vxwsVvJF49y9xJ_yT
zeo?EMzT8L{!j@H1(D3K|X{|!dvv;J#n$N8w$QJCoubQPOk6ZlkX
z4X6WW4xfwSc3>g61=_U~DacP>RP$b&S7iIHX
z0LUTwGJSCy3=hW}qdW>Nz{nxr^ee@wZi>A-yh3<#MGZWpE8j=mV~`2%X94?rOwr@#
zkz5k^op>^txU(Sq>vFg^(sL^&ry^f;Q{K2$TxHj*0ET?4wXqcq38P0YTqN6+VM(~7
zsYrq*I?kT#gfb!{Re|n=vYh$rR^66sXeMt>W9SJmaZt%IvEvmQ0rT4?ad+VJOOT!j
zMtk*oA)y_WV1GC~N^8kH{G)4T4aQqr!FkF5=%QmRe3&v!l3eZ#5qyBmd(iJzhDPJT
zhFef08&oIzPFq8g4J6hW&;{-v1&W1&&@2n+=WO%8f(?}neC%Y**~_Dd3=EPHd(bKO
z|F+MBmzk#VCV}ka0ueyi93k@R3&C?8C+0r=;fcK!JVgYsO?zk?{(FgO4}$uC)x7=p%Ks-^#_P8^MrGtgAoK=p
zh|l33Ib-7Z$49C_L#JuEPalr^pwJBP^b98(9%M8_YDb(~@O?;T-PjdCp91dFI<+8e
zXWzN%^6z_6REZY}hhm?z{ULTozU2of8a5V|H<(PLH_#NcgJ1i2`LjM~Tn`osO-zA1
zI{NdY5>cs9wF(^WxdhWdG`=MiIkW9Ae>XeuEoApAJ3cqV=@L1$-xr*L_B9@czBd2s
zGG|4MsMm+h+xwh;y$z3zk00gJ0A(T%hPZl)b-CB2{(b)qP7v@@`t8i$X6-PM-UboD
ze(C!yq&T#Wyb87gKkaRPTKs<>{DtH^vDDgl)XQINNhZuNuM~M>M&sRI4rgQ^3|(V!
zhlvl#e*X*KarOWEOS1kixB&8<78hPJ)7EEUVc2vO5+1EY)DghQLWm-z1SF`eczDji
zA~W;FT>2#%d+iGvid-xhNgWF3J~)}UR6=*)LlXY
z877SkiTsl=qBOp2A9`B$)IoUPWNAL<`RzLwvAu2}47ay}W46>tqAe7;LmX}lLXr25
z>bV$Gj7d(>aKOHwvZ`tv@a1p#ChBuh#sTa;`wkv@`jh?J!nX|qIxk}alYVp1ymmJt~Y##SL@Bx5aE
z^Igw5)c4Q#-}m45@B8C@&-=Xhb3gZe-Pd)!Z$6MY-MPx5IXzN)sLR=RY%RGJ?l%ns
zwhzxM80=IYVxIALFcIpbc3@3an7^h7bxN^TyRN8A2u{su;j9O0gAa|DvMmrdhZ=wS
z3c&H(Zlom7+kCr2Sth^Gb@&bQ9(BHiwsn@
zEqF*Y9;l-QJeq93w#4XwxjT3_K#!hYu{V>!+jwCRXs4Wcrtx=@64Ea{SQZ?-+wlkD
z21W*~cOfrSjw*R&`qS-SuJEl#@Nbmd^aT@TCl&nVb=9mxnmkI+7|Zp2tmAS+fk3#-
zMmfH1uP|#=PZk~zqc(<{UNxn0_-sf6op&@%p8u_j-T~UU9_>;ktBgZ92GhIPM#lxH
zkN5^bzpcBW-fY{2kFy&t6HYBZn~t0*T0V?i!E$!FGDnxK;B7DT%K@%-_$G|hFu9e6
zMaNokXU0;`6LK7HWw84i?WkNz`7x+q8g`<+k)PFPmG5j6SJG^}J8P0?!-;{u5wMJb
zyfvqjY5vBH0*hl3%}}mqfO|)}3gTiW@t}Tt4g$v(&T}z+b4zsA;8zH&e|)QDaw6%#
zRqJl)*!zfr7+OVa^30&YA{KHD_g@d60+j3owyPJ)Gj+7E^YfDpuO0o$z1MEmkN`#y
z3j(+-Mk*jmo1u^hDxdt8)!$Hl;gN7!XsUl0uoY~|a?Tr?eheZ(s7X@GvH00!SIKlQ
z&-&OdP!P@1^e#B$a-O>UD8QvBK?^rEPChf$F$}E)hzoUE|r3Ar^;mRT%5cWgG`o7-bc?McUP1!FqBvr@Purqbgr@^
zA1Fq1;=4_OE@y!Yos8%hIyy@Vn;ev@Z%zrjvx0=z9UuQxHzixdl;W0J1b(h@GOE#N
zSmq?Qxph7qe-W^$wo21)Sy`)1GJ-wWMLX*H$b$lD?7JxoB>sGX%n}
zm6C=>b+=8b-71cJQ%JPDGqignk2Up}pW9O$p&ijC+G1ZTJJ{;|!(js2kd?MvpTfbE
z!~Li@_DnqTdTXVCqV*p52=;?4EVx8%wV
ztkR8snd3O+%Jy&LFI^#;D814fR|7r_ISue^)((7<@aV03&|DIfLtV0eWhFBo`3>?z
zUdQV6C!GBf9+*(oLXYenkhH+cAcqV98^z^s;Ai=EAj&+#(EX~8sOGU0mWXTxK
zBfEC7G90$IGxn>-7vBRM04k`4{)e}D%Alg{iFYrV`JI0T?YFO-h7&7U08)Str|ih>
z6gJNGJ>(?#&@2!Y58rJKitfxxBL0gYei;}JS#SgH}U
z(1dN@g`!>hoxd1xQZAkxZkdkWr&qne&hV9pR)7U^&ZTg>^L0z7(|Nk_3X=s6p||!e
zYaCyp6UvKdFLA
z3Y#I0fGNI4!%>VQKLiMi+-dqa_-CcG{2AdOBH#I3gZi!yH)1P*B$p6>^ylsah3u(d
zdK6)9c6|j=W`RR>pSX{{-?=nrEaZ9WxCJ&$Z5QkPh|gQXPS?20LObPBlD3>i{~ugE
zCPig1wJHh2C5V2a;Qmq0+>oqL2s>|M_P;_3vff^W>C>2xjabbnn!)T77g~aDkG*@JHf-iBc?7>qg(_
z$Y)oOh*ER!;Z+?g^t>CsW@PM!^L9p6PcG{llMpiH4X+W+kcmxFSSpA~XDf#t)cm#K
z{m+*Qxt}`n7ukBss4WE^&eeXhvx{Wa3I52sA3N+~zbG@w-J?|bb?A0pC3r&
zBMw)#v6=>Qo!9>jVfY+bN)E9{F54zJRT!a4c?0a9R7+(i@wzFP
zT{&hP^{|J4hsA&`u={NAH}LZ0I+z*3r2sKdnOgi^R_uq-2`mJ^h6$i#HP@`{5;k5-
z6NrEp%@6PY$(mpWqov?=ntv_?ZdG3eDuDkrKlQ>r^BL$MFsQnCsCYVdEv1xdT-hK0
z)6J!Xx7z{rAelMmgE>unbVhO@mT4E_$l&Mf*Ax{&_%5@rZi;vK$V?&i^#gW7b`wkFkLH*ZOFz3`vrJ4GPBpexqD2c*|Kcl4%KU%9!tB
znL5w|z#4*>iT3OEpyeg?UDR=
zQCO~jch>A{cHZ!FWymXzs2&Gu2@eGnp94_+i4*!1GIR|_rm}>k1!C88&~+UF!IQCO
z4f+;2j8{G5F@-Zv`=#q8$m2u@c`19b*zE<_Ow6i!9cBnSOh+U>;1!=RMf^jp_gu5>
z0gTI9W#8c10aGfs^?QiT8TiYNq0gbuc;?-pB~2zmP)s~WjWeY5c;!vPqax)i(d-iZ
z_~G$Cby3QWYKIA8J?LVbzwtYZz0@DV(s&9IQ6Ng>Dq2`6k&HTIAb@L#7qvgs8~bMp
zeBO$Q1dfHwV?Ic9Rx>DN_WySGUJBT)@gT!K!!dXpQrCfLprc7~&4AYS0m$gHGvp=#
zs<&6vJwruIg!%F8i6z|=KeSi}p?xqOyQ8@e=v-GboB0uj&g`H#;jbAd(>mC5vGY_$
zVW(OV_&51ikQrgPC*yZD03DOB-1A4(;0|SF+9GP;}Xn`tIl?MC(7H_OX-jtEFBJ
zeLV3h^#p)I%@D#dpXzS86Ymw{quxLa&P>Pi)o*!|ZKH|Otn;;!sR466$;yO~m0>-6
z>0I%}v65G5n?`u4%u94C_KQ^G9~pAXSAsW_X2~O7Mtu2vilSo+EWI`B-yBv!o2vOz
zvXmw=^qk%R
literal 25505
zcmeFZRalh&7d<*N3{ny*p$LK?T`DCZ37SFAjlQ3$*Cg{s00K8
zIY@j2K5>6BAcR1mDJ-sBQM0y&m`Zbmp9b9K;~85Okb=8yWsFnvIUnEwG7p#
z^uf@2vr-+lTUN#=d?PclzYjLhjYV#}oLP`!zfA^zi4Z?%8?f@BLkC#q0ZK
z(>ZDiNdG>+Us;#xS%ypydoDgQ?4f1s&nNcdN?H_%&eTwA!tiJa6ss`ns
z_o5ykC|Qf(b_oCTlq&w#NUxKBF26!P;QVWtrR)l_eg`WfIqSV*>w>&BsCY&MWjHRLO%{?2y!i9yU!
zr`q$*JF8m{wPGWQheE{Y$mn0N3EVDRl&$sAD&$C}GAwfP2^#9HRqhF`N5)sh(~qi!
zGKJ3Gk@`yW=4Bn`(_^<8n)5~{Kt;fWi_G^uGvhad1h1XQpN^AG=gjCW`J8IIJ7OK0_dgqscjXz5i|QWTBl0HjneMq=EIcjL
zc$?I#%uesC!t!MmW2qFrTh@2;-zMn@-`+woH6Xl7UCW8hj_B=-tk+40SkZX>yErh<`AEt<@jocQAA?e8A)wuVC5XHZNo~m;3fT>6k^t&E9xf-Cul?oW_2a
zKY#gbtfV?Y^W-tfT>h%5cNlTT!
zIg@F6RnOIytsn8r^P|E3_O83vGC{`xq9kFTV4uyE?A!gH4G|$+6I4P1!N0$egZTeEZrIa{KwuDxax$8($k{$3PwLLCg~g}z
z$C+1HF?jMOA_&_wM{dPk&cxoLd=?D)Ok5XNl1`rDh%S_7w!@Xbx^x=
z$x#$BGZu09>nauz#)SUgAtj_#*zd2S1J3?CfS3xqB>3Myn5ndosFo53#|!^`1&LDg
zKJxDpNc0f5+Ql3#p3mPN&iC_Oc}hVp@z!;5NWIe2P3+I(6C=-ga$Y%%*84Gr)-r0O
zsl2vmPpx}Qazd^_H>|*MTqQ%d(DJ6rAMJlt2UegO(Uh6A!
z-pzYT!5AXFzoQYzAjZO{R}>~txgl%Wl{4!j-70X#b3x{muw@v7l&ALf@E4Mm`!o5?
zBi~Qkg)vGzJXc0~hw0y#e;%&`)-qylOL?Q_}fEt-)VB~dtbiIw666|
zv#v^D#gG47q5U6Op8)HL?MfS)QvFhglq|!l{yPc{uR}ua_5Z}`HY8ZUAgCwG`XyM=
zYBD(<-~VquDURoPk(tYc6Hzd{s>VKO2{!u
zMtOg$tE$lhx7qL5e7P+_hG27Nxp}u-cq2BfU+=?*4>!c(T!%>V
zr9-wJ>%U<9H)z-&gw{-0OT+8rLSBa5x|dw$?o9s9PfJcoQH+=J+WA`C_b}RXz0YN-
zQ_HZf|D%V=?++$z#de09EA6VIyXym@H5Q1C?fI&JBQBLIv&92Z5$-txUYny4>Pd1W
zr4va3p5>1!trFZiZ)rsiYUptm@{{l-P`{91Q4X&|HzHvR%`u7L1AgL{L
z72C{2F7LVFwV2-U#+V2Lcdm5tQs2lMXG(p28~12cJvKsHP
zw>duUG}Td8#kh5^_SIRj(x4?r2ZEu5;_f0FtEFc0m>+wQf@WjA_Fu@wN__(!1vW6+
zQ#{xk`!xXE7
zVH;P__drG~UD~cJ9ag0X{@0@o5z#-@FOLm=@*QhoSB|=7o4{r0vs-8p+3{1=OEM?u
zQ*Py!W$Vwzjtjb#{OI}C!yRC8~jc#o}u=Ke}emsDN
z#}svi(KQ)&H`n0f&iLcX79;KLmrvMZMYzvf47B42h@MNCc1gu;c}I)b+!*L6wCakw
zh(5!#A>Ouyj{9Qu->}W$hmD@pwVo@jHosZQN>{Y-TN|CW<$xObNoVw?s`Re0E`xJa
zD$yC&YU4i=DxyM&&P@?#caGS-YXmJ*ote=xvfY%Lt&3Var+WFhdVKY=>&w5fm-v@{^#f%3qp8QU
z;@$d8U#yf;Q{9%w6t7-=a+F@^nu_$UQ#Dlf=L1dcLy{=DJ
z=&_bdQwvb%9+b8-ZLZJv6bo<^&zq~vy>;2Id}J`9!W$pj!_@i0
zplf=xi7+Zu!Y`z*OO({RJylukqsMwv9r^hmd&$YkF+FxwmpO`k(tbD1p0?xfxWBU;bjpwE(kJhV=}<$)UqSzd{;O>DYos}^Q7ys0po6KVN7k$%3!nJj
zIrcYwBhgv3%!oqWH#`6P0%KbmUovF{tY1&@o{DZ`H<27ab`$#EsFzYyK6lra1DpWZ*&b`D5dWVE>=84px)E@
zDr7(M7g>+|JWKU$ZI%`p?hb
zfQ)SX=@H>vCasJzFCPwN)@tHItEzX_WFgz=Z_xJnz3YGKA
zM-G>CD{+iVzgtLbEThxAy_fF8FGWXs9h=F!o8kZE1k?b^acSk;sXX5)Ao)N9MuZiuMJNUrIw
z+=%+uJf~yYQk0Zl(hJKm`BGD3Qt3JRYNyudShSspOqxhpuh)9t2ak2zcI{yS?OL+g?`waaGS!QX$DHW;@Pyh|8GVn%
zXj3wQ9x1)tD0mTYipHQxLwFs;QC5JQFBA!D(FfZVca1R4XV2|~XE>L})s=>I(RUEV
zgH0cF)MOic)ZJR2x7$jzRLjWDWm5_dtMuAF}cCf>H&%TO@u|`SH(ifREpP7WrSYXp5lJi=xz&uzJz+!>x7JpVym(Sj&
z*+A-W6_cc^+SL~)&vySvfZ%`8vV)WBX42!~X?enX98?iE7b=HR0b^NOS|+j%detod
z#W%$yFaZkuC4(qMRaMG&y1QGeCXi)hwo#%5g_?Rc94ay9TaS)Znx%D#f59G!7T-+@
z*U@tku5PNWrP@VcWHE%iK4_-eEP?=1`Vb~p-rU3JWdF80(kzavo31lrWj^VMwV-L^
z!mPtr7e>q~MEVApJKy^4uP0RBp5Hw%7{~gU(^Q0Jm!g(75%fs-9DXHd&
z_hf0lydLAZ@QG;b!OpMu_fu}Zz50i!+bKvzC*{6`6XNwl6W6WDG?RtFPpJi`MXZVG
zM^YpT<}%SE!;zkN`Hs3?2q!(n1!mPRr%B^s(Ii3hEV4tSo>obj0d*wuK>zS|;g
zlVgKYoZ02d0PESo6U>NJ#s6?1Da{5T&@X|l;cBz83HL7pexhEftaRB~tN{emd%Pf?
z+>Y#PA!T}TM#kmCjnqV!4t(a9D+Qx@vhNG$B+uA=iSv`y8INvDAC(O6>OTe-X!@9%
zdKd~A3bAYJ9k-PS67G>p6$Xu0;|RWli**J}qJ*}jNN^Ykz!-7^)&nl~DySM3l5u?)
zR}z8$l5U$|eEXf^@(07;P5rLK4lLuTQUiZ^qojKmW+lohha`GlgNVy{uJ^Yyh33Qr
zdV;{+$vfr%hYTWdvpvPpG&@d@R%2oq1ef%*ge-sMB@$UO8GP$5v}y=pQ~KP0^0Iq@
zC!p7*E<+!}kD?xH)}_E99<9CqV}sPieDvBsd}Pd76ute=^qW^88d=^JXrsNXtY7|+
ze0pB{+U*=4y>UVPVq1Of%Y}9mO|cSR*7CoUqtFTcAJ;n=@EeI4g=GAlO<*BJMSi3F
z;cN)(-y?HF&N&zZCEKatO^AGde6fZwld7NoLg3GYPE;@Ls+9-J>EDD47)sgG
zq}}jmJEx}9jRbAm`0zm$zPnHJ{-xgIeD)J9XWQwV@rdg4x?IH5X;yn7l;k66`pZr7
z|0QhB4+x}raNQ;3L2$wU8Z^hxn2>u+06y&uTHv}2Jd|%|>$24_`Bbb-w%4@NYM#L*
z0ad4!54eA@{}|ZLn>LQHBz1i)L>`WHO1W#%AU|G0j|En9>-pDg(14YU)|SIlv&R*@r?
zI&`NuoP;mIw)|}e{w_I_uu>wn1C<|R9Gk@CCEF{X3g{Kx==a(85~MUS;u5~ZAY!dO
zn3HbM93!StxjiRQS6{yrpS$e)?k*6J-EG%D%1HfnZ{QnC>jOR~uw%=t7&3FI^E3Du
zQNfasxjFH5O8h~IUfbA@=oFACLgM1)a@8}`QE&b>-AW8*isLhd;rFUnMp!YXmuBu{n6Cwncz
zCW(1To%i?rJePlnCdo_yQ|9Xy5V)+=z10i82!i1+mfG6
z(192d{Ae{Xg7aVgN_37QZ3IA9)ap)&MdVcI_qz9R?=22~=EDPs
znT-ShME*aX#)zK%f5%`m2`LyZG2;Wdbzj*napx&x_9Yi*zl=L^u=?CS%d@QtW30i*z*E0&)MCsQ=KRNHs`{)3Ys=Z);8WPPt)VG%fFk^Cm7+ztTkVY45v5Ap9pK{c5`1#>rL2i;?KHa10w
zR>=xMXW||%-W&Cme&Tu)D}ER@u$Vh&h)1SdJ`|@e$`x7vW&zuC3#uLf{`
zYh_GfM-Crci|iZ+N5jVAu)e+mw^j)vXpl^Fi1>dH{3KqRqYBOX9|k9-%Nu3g0EtOZ
zRr}Ak@q=S=Z)sXA?5`0sQ+eQ4R-AqpT4~J)KEwoecgj{RDHR>XD1~vJea)lQl-oFS
zz!LluF^fZzM+44s=$sv-yQZwn7Pn}Ht{N;>$s^MHw*n9-S~4oE?$H}-|NG<<-2SIC
z)`zozUpv~o7=-18^S}Q-JA8tMFCFp7^E((?{4b;T2xR)VYUBT?4>J{fiGJo2?Z1_U
z#eReH(h`LJO`>o;OdbMh{{P?1|35PS*2w?kn(?Z$`;LOpN4|MjtL*
zZQfU*p{((qOV29Sx~)EVZaM_ynpVy1&+qT4(t~UfE8!9jM5GQxLyIa|^`sKguou&v
z4gw;O03i#!3SKJHi)?!8VUx=lrlqAdNX%ETY8iSt_cnz1lT8im
z>e+qlyyHIK|8Z>FkIsM-!t>J@r#8iBoX5KgaF{!K7@v(lL~bL&3dUVI`XcG#@`qf$
z+8qY~FG``=!oq@mql2d(B*72B6sHAjeR<5M&CIz0+vD+YV)^g2SD0Q@3*2~=0wu_f
zD^vov4$#CxSv4ruqoDT590lUJwA?IBH0wD1@VeF*7ltac3Ypcu0ni*u5q;9S*fZ%Ue3-4{cfw;<$mutJ3Yc@=>6x
zVU3M9s&JKO&KqHUda?zC!o@r{^7Jybj6Oc66X4&HVXdN$K&-qCk_stvo{p-F-=sNb
z@AL3sk$C{uL3yIu9YV0xOO6|Fn1PbMV-R!rdM6zebx^noCp$AE3h9az=dyJRLrw{r
zUaJXYWMgBK54^mljhBQ|V6sHG5}m)b*wPuiG;_*QEp-@O`c^23=v4IrNt$k#@}xzZ
z|8FwsH^J`q4EPX~f_Z=MTUKNLvRkojpYH9PJ6w3JYJ!=85i9s`Mq9PiX`Rd>bmBoVPK1u<@|PIEvm5sD}N(e%n^CE~2t
z6)cPmSVuE!R^?HOrQc}^}h2;gb>br2QD
zd%n5+7wnOgBwuZ)E1urS;QqJhvA_Vri?L4(^b|e^|1Sqr@LtDJ
zmpf0hgAC?i&c`W&|1LIyaHt_6>#ovKC(9H1db$xCtbY?J
zHDLKzhcr$OR4J-^DMLdKtzumPr4YC2?BWF(cqVJs1g5}dv2ljUj*L*xwcln%Ty%u*
zgCL}P%@MLacvB7F#dxq9chRoDqU5zFs3A#Y!f1BBns7cnlnd~9G}W~?*qlXgd`+Tw
zz?(;UXTf)6=C=m$*5JG)`@0^H;<}eB)wqd>C)6DTKCXoD8dgf$rq}?`dz_P*HEQ0P
zB;Q+K?FY|~+d(`1*alVQqWN#-%MVv3T0hwI6z!-NAS9!`cK!PE2#a(Ne@KPI#L)-gr-}B`?n?U?
zZu|Y+^@fxP18!lfpGQ4{7~j*FwF%&O3a`
zhsz3w>TXL7>L&=5qW{U&5OW4=IXzd$5zjkX|C%f-Zq^_zHTC-b-mZ-nKdCVjjLq+H
z1GVG$Rf;=F+yN>LFP*+jmP)Mn-R;z`x(!*$!l2S_!j0Q4Xu
zc(jW6b7-ns!@*EhsS9>f*2l{1eN#a)x(-+UVU@P`HmFc|vr?k)OxS58&EPRfX(HlZ
z+D3~*-%MrAW))H_U=5fysWRTx{cvvI%f@q>Qx`=TD81(jG8xGUc7Jc9@}cKEjrXU~
zw9k}Su^4{@U&2{tYZee{V~5-6JS|I|G;>wPb|%x}yEEQ|5rVpeFW;SUpE&CVY@YFi
z=iX|Up=n;zn+x}omhP0;C%Wlgd^qzfy=nX$dWZ{i;<{vEDJ@E#nwargVj7c}gWTF&
zpWfG)OTc1h_!Crx+uGYrCR!55e9dldxmzyaAy_Na58jC!M?IuVAncfP@n3S&9ekIW
zL+a6ylg!?oO@q0ARWA536B!w({32F+;KUS+X7^1`<43(qn6lWyO12>y}
z_U>OD-b(eaTmLHA-!l06F%XbrEzC~c1r`pVAw+-D)LjM{brVS}wM$)LpAOSQ@$mt5
zPd`dr*9#<#y+IA7Xb}1Y&6>iX?9fxb5Jmy9=WI*j6`~W5r`{=abKnbqBBNX5@jg2Q
zfd$PRL47t4$~onZ!Nw0D{{DwpjwK;zz*(%iB`h{pf23w1(Jc?uj{alBM`fYl<*e~Q
zKSoxe{Di0b_;~e8;n8Cir{}qyua{tRs@t;lMDcad9|elK%(KmZ25&Ize>D=GbcC^R0oSZmsWb@d&1K8H+{}pSuI9`k!Jf
zs^r!nUTt+>or*MB4hapF*D~7#4XjUBe+q!{P2!oJV%3%4vZi@GwYD^FoZ@!46|=5(
zTZJO=ubWtfg@+Y9VRbYp$lxEGr*G>NSw}X;OB;G`j?xKM#q)40uz_;Ij@mk@A3@gX
zJ-q3)yEUC_N-8oy==|B8T)S{osLIiEd-V>iJ!;5E!}^Iga*n4PR#k!`P^g^;go+H;GAQWuE~UD4QC_D=kk*8i5`FJ10DUs~r2~Lt
z?}z0-dBJHyO>lk;@Y8Yiw9k)J4HJoFy^edjB*7XzF;8!DAlLbqR^uaNwq9}g`};NZ
zuqTsEe?)x%Q^6aEhbFP{H-6_(>rX3xW`ucf&x#6O
z4#{y%RV4UwobCGKGg0$JO(bQY-{xz7OOg*50b!a#{I-^K+X}>(rytX{gCyoI>3lGL
zX7LfP0yxug$)L9GKGwujo=;L0_-83Gz;K8e7Td;dI(V4`Xv>hZoBfEo^)-LbV*~Kr
zTX9uA*;SsZgm$zs#yf}S1_oMEOX;I4bMnY0RM3>c!IsX?fDk{n|1=)c^LPV
ziL+i?E5V?@%EwZ2=n5q0ix!X_;qrlXvSZ_lS7sI-*N=&bG@-iF-oO9=|6Ti|
ze+Y?4XL;G{+BMR^qQEv;rH2MiC;v2yb4!hJpUNoWX7U46WJe79abRb!)jL0GYi*V%
zTGTH5exTU*(}**K)A4Ga(G1}?#A}rfT!=v>=deqH=#jp@U
z2ASSpn2xTzZgzhH><~M!iJ$S(cLj!G`3Tf(5*3w{0!|3GAW{1ja&`z_C$5D&a+TbD4%(){r;gaT!F@Rc)
zWWe3cLP%4-Q|6n5Hcq&a5#!d}oe;Q`zJ!*~d5AnZ7>;*8lY@d{mX_b0&>=b!uj`%t
zbnn6h06*IhNzMe;5=G*#8VL{s?^Xq@Y(+*!F7iGZ?Yge4?A$m&)NnzTzuGRAib^Dw
zFTru%f;lDj_xjB5;hB1EHTPM)l%P_S2JFKLCW7NITT}u#2#$4FB*7!y+UpVgd>W}Y
zo@pH`QF%mY{Tmf!bb^0g4JTsLT{~a3Z!0BZ70^V)y?c(s^1Reiip{SkB99h`jm%w|I%(|!ON0!O#>66K)EU3yrnNMGKntoRH
zq##E}Sr)$lG0Znn1z7g$ans*56QE#}l+S^NBGx&rjUa7=4ecnsr0lGDJ(3@Ql~W&B@-<{2hoWAX&1t
z?`uT3{6+a#>__y*v~}*bj76Y2hCpo>P%XrSUmDCjJW98-w9lHCS%@R_}lt
z3OB@&{a4H^MTr|}YUxW1zYt7@5YU#RU!J$Ojd8mI65{FgNZmsT@s%1fmOU4iJHA}7
zNk`i1JTPB2o%1+6Ab{g4^AIAe@X-lJ-;%$~1~|#k3GvmaCepSfd9wXGb**DALgSNn
zegKy}S}#lY)pdKe7%$Y3bh{(LCO?O9e-9ILGEf3rDj3H_z*u_YVeyj3z54c_gp`00`9=lr5L~S4>>j_RN)oWrx6R
zkiIxAssaOPGS=JiNt%m`%QTTb;<|L^E-S04t&qc<4peDTpIP#ji
zds*al8J}0HC&)=~73jTxyT;+jt$A%6=Fn2k@6&JLA1^cUgN6&jlPkXm*!k{+Qw*6+~J{
zduq7uj!W-fi2t6Z8gF1!=q#5^;n}_?BkoLmpc^o|Mo^(DPS
z4WmnKfkpeywkr;+%A6Y<-u?25^w*KsPVWTtYxMPg7X~~{4P=yv+Ea34Y?ca!U7@cc
zF)IXe#tTUp$Fco!*l-DWhE0X3{)xXzW!?St`2nWekuZjW&Cg7}P#8ZQM+>6wkn_f$I22TlhK&bW1Mo#NeIm}#g>|kGK1-z&1L2AZtg88E!2V`GZ~)T@6MCr
zkDTWp49T@3)VCd`B66^ad=JHm?ieU`kx;~1$P^1xs+VAkPEIpa;1oZuTscdE6d;2<
z^ZTfK3>-?LGHB3KNl(yoP6n;)dtaXvioThIk
zwdeC7GmlsGLenZ*_7r7T?C)(Sg>qPKmmcTLV)>f4)wj^dP=|AW!NMTW%|*(pP@MWnL5+6Pw4R5qo3$I0qsKon(q
z&$_=U4}!t=Jd?RoOLS4ebcbQ`4q-ujj{=j|xy$Y=ufR>&04cQTL&>J&VjvM(x4`_E
zh-l%%eF&1iDTPepPUhNszkmN;g6DraTl|G$HY!0i#q+*IK~}f2%=jkzSt^KpdYGCZ
zCOdJ5ZC$-^UY&c7TA-Lv-yl%M5cI);JcxMWt;a5lJiGm5O2l(-(5W!VTDb-Xysui7
zKvkbt-Y0&dqVqQU0PB9WmAP7);lmqSqEHvyX%M%#3==^t$AF#UCKXRoV_KZV5XCNy
z;=LSP8mYDKYtEd;U*kVSn1+h?T8rd0sL(Fh);8E}Ap4Q*Pt*%i0|f7wg!}^M_h3hn
z*`E23*cY?%-opr+pAqtOK6QFRi#UmE@)wN~0=w*Yp~uZ+&!k79{Kj&j1JL6@mouP75hTLSvoX=59Dr
z*bYCXR$8vyRWLXFE!ZsCeoLC@?&}xr6tV0l?#dd^`}z5??BagQs?Rj)EVGrE=dp;sx}{>?
z^N7w#yQOH2^9$|mqY;MuOyGOo)B3}BaE>0*iJc~ihckILo1DLV2%WrJA;A82!nGr#
zTiLzbC!^bYBA!YGR)GR=Qgepj*{5~D9{E{LM?W~2&=MX=Afw}7Q9YWoe8KYRT!p&#
z_Of`m*8!DSJKkdPXh!_jtB=XY_2>sMX;H#f^$%(o>74I6-90r+JJ_Hq*jAZPwr{@l
zvJ4Y5p@_RzO7Pr)C1dj$M;xcBR<7>=KNbaoObz{_O|}GX}gS}_&7JDj6ee;a9Jhsj97OUJVUHZ
zh}C`igWAH0@$vmYf6yzorHOqNR~b)Hl|jMlEr(3cd7GIXkI(NDwK-t-cypV{&}%hR
z6TE_6(Vjs#Y|*15Yxyr_!i{lQGr|Ia@+0{9+la>`XLjt)A2m2jC4uutseO#$Z%v7R
zbA+6p`l0iYvwp$<83F(UG6AMhv)OFSEl*cbQT>NUU1&lmm_^*V^JZOW*M#cT^O5jakE-TPBZ}UDU?nQ8rCjX}W&*-JN
zUxGpTPmZ2q!E^n_Hy5rrneLZ5j@RQWtZ(iI2DiRd7=Li2{2+}V{vxLHPbbu5%AeLG
zFUzq3T%8P#7WMb{t|}zGR^egALr!{tD&<)4$>5V&04<;Kue`gaqja(7cA=&~t*?9+1S7De`^x5~
zoRIMW;`jknBiW!)#9;uy-{3l_+zPQl=DyGVx{pJvJe?xc_T`BfC+kV@TXVX`8#HP{
zQZIlo>_)Q?}oR^lDbTt+c>xcsV2P3h{#9rhg;E3XVX)Z+0LnNYkn
zu4}Y}0(F6=CEifEo9b)xgBuogjVPpk^0jjNu?Q9l3_c;cZTxwv%+%7-$fEOH@bYE-
zoYFD+<#H?`A;;!7t!g^#)PRP-{>vA+;#vW>
z6Y{j#lW5g!=DJ2jJXWVn*b@b=-(Bh5-x-$%Dk0JO8O-sD2A+%Zn^3Dt0mF}*X-W;v
z?FC6yr#cmZu0FQnSC+pt6F~b(GCYTl;T+ulJv#pDnM=^cr?3#Bwa526YQ~$TwzJIJ
zeQt(`R8_9DDJ=nHh-wFS5Pw;mrwl>F^APEp|4sZQC1$vx>^9aQNfvy6Lu8;L5yPg`
zvq=vCA{R2xO-{EYXS?A&dL*)jZlg4c*RNk>9dZ_C!}C@=Zvq0cj)p1mSedx{(@T22
z*HM1^CSfX|YeL&6AAdfZE~L>ULQM1~9)|XVr8`e5M~$@{EHSP-u2bfeA!O4nWXboa
zIeqdHo~<70V=pmsUKnBf4youyiLt0{75RVb4ELBst^@ThrTz5_{*No_rQ>lPjotG&
z4bEG+27N1(%G!J-=l^$a_wn~OUouRQKhz3a=`%N7o!rV@j`#mONxXFUUW$Xhip*mM
z&Yia!W$%|>qVt!7{ovG7WMj(JRB^HM*?(KsNB{>_iAb?xXFnNzO;b?0O%-QA8S673
zOc9DZ;eQEoQmB(#BaB#0*1kc$sj@}qyVo($0zt@MI3jkF`^Djyk)tV?qWs
z@>g!$!%G-o1A4RW^yv*mq$tEHI5r3E+~1s($fClZ_2vR_jX-my^d^EbsK-E4NPyFP
zzdq>N8m)Cx|00SJB~rgkgXq}w4H29oK|83&cuTHwS6RD?#1CGo(ATd`=%ZHVD*|Vk
z&q5tg?(JY9LbT5*6i*sNYIxV5nh$exgQ;R=Ltrfec6eRLlP)O{{%@a2xFQK<2&ZFmerQ8-ZNUbL^9Y
z28k%o@D6yb$I5?W)dOYXc{Ly>UMLtBPL8Jpb%JeUP#@`aaTDa9r^pH@KAFm#B9|-l
zh+QKmzkzqY!lodc8xlrcp!1(NR>5B$85$ZoF1`R&B08HQ4k%FFBZRQY3Z~FB28coh
z?k`RWewS3qR7;e-NA|t*7>4;S>axdL7Jy5LQRQL>0QZJRWVNxnHXsFQpGXSsdVI2{4JQB}1a0b1t+h
zFx5RGB)hSYk!H!13l7id{e1bx=U}mzg=|!*bsBDw5$s2u-
zo$3LiJ6Z8W`{#@IzaJG`xTYtsS7gm4ZUX62l-N9T573O+RX?iSrxLCH{H^Qyv0F
zCt*(i>7L?Ti&X$G-Z!2nAd{NT-K_N7d?Ph9)%y17=J^d+YS9js`+Q!qi@>wLCifa2
z%C;C@Leis{@gLD@z(t7YowZV<8wvbWSPN^^6r_GCM
z0oofmsr+!R*Fsmvu4;Gfjmq6M@P|dt5df3S?|64h37`M7)W>q-Z(8C3M2UPtHDU$2
z&?@jZCO)L)Jde{tzvaGDtGp;>jAmB?;C>*wv?zPWMl6!xkdoQCSWvW`I&mLA-2IcWP)60=|Dj~TgdAEDzJsV&>@Wc8AL*el8
z$3)*??YEa-_#NK$nPZ%Jge$-(Lc#iZnl!e16sYfOsD#}rbD6jO8S~RS@xWibRL61V
zPg2xF+LI#IZzQ*8@!t@8d7$DpG
zQVE<<*H@>znpX%`cdJebnZMekqEj}AIYhVE(`eWPg@i)6HBS?1x;Bq15}#AYNTW^6
z%(ep(e0`hS>4FpB=y~`vQBWm9elP1*W
zW=Rn0?RszS2jzR48fR;FFizR7Togdjmds9|uQU)Zppr(fh0k4d17dh{a6+kNX`rCn3Le
z_w+%A$9*NC!m>BT@#rxp@Er96tbFb_6n`gdiw@oM7-ks|en}pUfA~o}aeo7fr(k=*
zw`f}NU(5LmJk*Sls1}yc2ajS2!A*J@T(!r%&drm_0XuKD)PGoM5k`KXCpn@BRlbaf
zk#o)m@LYbIG|GkpdInr1Gp@&*59Djj}zZ_IOkFy0n$f|~tatc6d{2HS7
zpR!bho>w=Ou908Y)$I=tW?-*oKKHcm{PWYLnA^9*0ZU$iNX{(bGAHyYt;l<7WaZ!J
zUyDH?&4*@6KYT;+Xyt@;baW(SPHZhWjx}D$PL#3B``3^2-~^N*PYcsMFUEjh!r090
zss0)Wz^m;Z#KEqODW|k(OB3H=;Mg#=-;^a>nKXiKyQRxba6u{dY%nZT`aRq}
zEp9jpu+F#|i_FP<^2YfOlQ&+enZdF+oKGZ)SL=aK-YfS`HXSxn0BG%
z>y@d_h8tC*qXT?vrL1A+=#hp9-6CrpRRHSW1082Lj2aLf-ZvKcf<<(uS7ZPft9uGC
zJk8A3elP2kk_7kxt#1P~7{4O%y+AE{a<=9z4&P|BYVmFbXxk*`d9E^_oB^%gZv8yj
z)6QM`W>VMg`0TEZY406wO$4s@!YO8EDpe%86nLk<@Mi>tT@T|tRwNxP;na4}OSrW%
z$qL@}tHIBcK8@>`#j*jD3Wcr!+Q~Zb
zsyJx7xUteeQb+d~kET=QH5#&lMz5UsWR_0cb9nFw9~jO1`Df%gnPL?Xh%vV)
zX$j0`Ka~8uh(}ndPCk}`$UT-06|FXOylAhzFO}yQL5=O-pvQhvLW2m;h!0$(
z`Z@1`SAX!{btec94!+CoXOL%vD*%%i<||;ZTL#dyiwc5UI7nW-{1714ByQ|JKbNrV
z&3~r^dR!+ymAlM~82ZJ!Li^2Uia%b5hu#r@L$~y95C*WEhgRZJ5K_oMoZ*`;0ULwQ
z&ocK;yoXv`y)r=)I@urfk6niERTrV~whurqKP4ASr*q2jZNWu+JcRjG5Px^)kFmYT
z$H-XuO`ESP@t4X8ZOr`4j@~rs1&hrll&MDnaEmHE8xxk-h{>!Tg~D!cS^H?>INEY+pX5Q!+{bvY|EZ_-DJ;{e$HZ#
zuTW9B;K1--R~yXA1_d$)2Zs>|DC!@K74ba}!v?hc6KBA8u->a=c`p%TzL4>@&sjyy
zkd%QjzDqK(a(E*UFxDkyqsg|E-;_WJ}cYi2cS2=z(k53jYLupn<-(@
zZT$=s(8`IKkUFTuk|AR=SfTcfj)!Rd;M7sr2~2fgpPd*5j`tc&iVgn9Sj~|yzx`OT
zKU08p-x$+aeMAVg91cLx!2nxYv+SuA7bY88=^k4ueS#MtO!fT0*TR@Bk&HWU6~|Bf
zJT^h^IYn&`dbGcIb5hK&Dj%QrDRGOobNHUYn<}H)=+9_#BsuF{BGf(&-rr^dsSi9#
zuhXvy5%(3@=xLAdre7R}%&<^V4R`eAsFUUvR4>>B1w-A$zBv85YXlyWv`F(j7-3Z$JWM3+8m_A&K{JAw{Vepz
zAKEU{G6_waZ
zG=8WS7Y_0EC%++n(HnHPwoEMnyuif(HlbsK$13#a(4z}k)sR~rP30){Pp1FGi;K3U
zqu~kb2kpjA$G5&brHp>|2g+0ekaa?7`E=2rCek;SMjN6BHOI+26;|%o1PU`b2-~xr
zePs>{SZ)BBGw*iL19}QCXG&AwwL8#4YiU}r4qzMI)RhE4Xllt$4*7SU+J?^y?06;s
z3pV3wP7{As0yETN$e1_;xfU>T=0_MqK?kuh?baMWXY#5z}tMv4ek`7TLW9yS;{YdYA>#Q4g
zUK)2zg{7tsPB3HkBu|l%zufu^GWJxY})9N_r
z&0{$K>;07|Q9FK+nU-@IbKg!DN+A9U))KJQx;N+4MD$3*IB&}0Yt{9)L3M`qtLGCj
z9%WBFBtx5M-pSaRM_qbwZ-W0XYzF6+?<;F(3CEoa0+1au@uVGUaLfg&3_^7e
z&_~f5%1_|FypZ6
z-bd@HY3_ba-M{OEeax;eZwke(B%Fo!BxKN}1~PZyZyv&ox#Ajj^CZU5;Z_sgpD8B%
zFtj@F?ib^wkvdb-OS^jS@(m}1?7Qwbl1mTRmbFskDgmO~{N0#KA5N6ue%Ejz12+l^
zu>g!o1aw}EX|dfN>%_7xJf!eWF!CzbjZtgYxV>|BcSNIu33cO9I7uFHKhcFX+PK
zkXrOm8|vNeExb_~NUoZ3$`I!vP{ISNE>i6?SbpeW9BYcTs+8IoB#<6~POua(2JPAW
z?_~nSu2*1wAY`R3fzjOkO6Xt54PJ?+=BA2}O7IJRM1Vix;i)BA0fR;ADa-#|G6i8(q)xm1PLT_ux0(?oyo+v2
z#&xx|ui`Fy)I-7Wm?1#s%Zyd|Z#lRb^jA7ecbX_p|U2GSR
z=SHcN>0U}V>JIL{mmOxVNaDH|l|Lo{5p|#YJ!v;zCfEn89mk&?0fySE~@H-1{e{g71sK6i%#TQO>W73
zqX4Tlp-uv!T5p#SBCt~x3!PjNhpvd8h<7QJ_r$1T2gdWV
z7Z=cK;_{kwv_jQIN1~Q7f&|pm)Ck$uThbrnm(ED;w2W2}Pybe(%n^Q!4`9eV%Vslu
z_o|Q9(RNK2zjh_?v5y=Lv$Y(Tz(}+c%=^{LJ{G}>L^8c#A|riP8j(%YZ-WYKTA#7=
zsXh}m?2?ZYYav~FDhlrji&s^cyJj<=ccsz4yqAhaAE0^b{q>^HtB0^unjLL-ss@6u
zU_pIX6UJAMlPPbyZ8uzDxud3B4{GYo^*Y*naW^S%dyOAZ4
z7{sJa0P2-O(hnz+KEzg~(~!Yml&mSD;r2CX(?d!hKMbc}BsI+QP*s(ITzAoK+XeTl
z>a(BM*1oyh?7GfwLEeXlgc|}{-H*3&@2c*;v7P3=Uv_8`QT+yZvsicY(!vvTGZiyu
z?!jE1N+`a*wm*!H4pTprcOG3##XLqBYRWv$120vg;eG%$tnO(7~R&rBE
z@Y*DiRRZc~8#Q?9JLZ{GPicjSAUj6c&Sh^W$c~ONx7}SYlx*u6Us@mVd}MFfGVs(`
zJ9peur`GAaV^|<`!1q~k%=(+VFQ1&X?m2BM47Y?7EIH%T)=g59eA;nM@U{xAetM++
zqgT?=3yqZby6#=ZlgX3=BZ2J1-QBUBhoT2NQfOV+ONhBA^cEjrc|v#hd)K@xX7=NO
zf-!PUU$h7;hIqKiU$jm%WM|-;5dENlVMe{*rE8tIKcX!j*3E>Uv6MeGQPpD^wc9D<
zcs2nnjG|TKrcR2KLw~2m6iKrVRVy6BX#I4sR_xy%+)j-nt!uGki^f_`=X!(GI_^!&oDV6xqIJb{F}^eAj|NeDDxn@v1z)VT*vHc{
zZo3KhebNh%@qm8Pq=a}657m9!h(HHjGkbmGDJAZ^$X^P<@l+C&^_h8v;eaaSAbOXZ
z7f9({OYEMO#C@Ep-TZ<0&?()1`-I)oHi!j@2&fh&Kt#Ctq05-a>v;JSWKe=K#|;lIphmDUVRw`Wt3*L!
zR*u@dueHd#w$S5{7Pgz={cR51t2xRerED|r-^_OAL-%$8Q{LFNXq-_Pf)R_x5>}W5
zH|R@k@Ecm&ZT}qLBgjFs>jQJ_NR-HWFjW$zv2xv4cm=&JNW3!r{w;XFWTh>P;Z$xz
zl6VlqMiJTWqnWGRA8RK|!eGB1hVoRhq^0`%iu;37IFsx_@I%p{Uq0k*p9G>$@~sF)
z*kko2c12!yWHeHY!eR*63x!!K5X*@8cm*KK;n+(4bWB%`Oa25OhjTsLxS9h|NYGvf
zh;^al5p(y$M6wU#FR05QNfUP(so+GcwJWst`_g9wL!DTHpb2oeh9@xiH<~xJ$7EBy|Ar-^g;T&;YN!a6PH}!-Av7s;Lpj;-CuHgPV&T7k2MdfS4?;RY!O_m
zuuXZ#?1nj5FmJce>JTTIB$lg~oJG_tzj1TY5Pa2}?pZko1l;gypyqxMu+Zr%C1ovh
z;GB!THnw>U4T)USiTj9#x|WgHCM|0Rz{=h+Uof%2nyE+Um6IH|MX-Qbjz{H{
z0dc-M^b}1Z{^+|g$&K9G)
ztTQIOcC%J0mGs=K-S&$XL6)n|20X`2E9n#u3-JGa7`1r~kb{ShjThsMr#%W6b56Zft*DyOA%%pcLBD(UHDyB7y`J
zw!vUE>s(ek4RfzygSM!dQRwnJ-A-P1=?9=PB&MTm=>ekTB`HD(`7Uq1RBLD=8);SP
zamehi#yu}{jX*&VV{I01GcDBcfkhV++FL_S#2onyI{ne?;>VIr_j*AJ%NbSJ*K5s?
zu=k=7EMQsV&J1f4QiMJ;;H;sQZwd`4_gek)$rub!gj8miz4xm%U;nNXT%w#UOdx{A
z%C=OXH)>(+s6{SYANjPL`OK{V91Y^oMyfjTa~zCLVI?&Cop$^=5i7hhJ^C^3vYCgs
znF2$zdWMXs-e?kv%OkIL!
z6`vxV`qaz+PU$?nEOa^F+s=W1Fr6`H@;;jY@@X7_nakp_n;!af9xLz@$J7&T2%I~6@Dn)CTSpW{T)EyCd)
z)$k%bcvV1vF8`qXHZa%sBi8$CwuoL5E#5RS7)c~I!&IqMeJLSEgjK>Ld5V`k2eWB-
z=kv6xx^z0&Q%YAycyJCf#cEzu1Ed*KK{Tx;-f$Ro#>r7k!)ev-4z;qkp%e|hu(`f`
zSwIANM2Z9jULG|;_9p1aBf&3TcVN>G6ipMPPXc34O>tg(InlfA!y%T8FY9yX?H>#a
ze6487(=}>iEHL$;i=;l4uMQ(Df3uY}kX3nItE_uQB2W2`0&DKhcm*HG0uPAe7tz6J
zf!NSuCGXO95Vynfm7l_tdy9i71CE3h0R4Xv)Wcb=K(B8Wjx(qL?q#k|pZrQbZy{9?
zI%p)^BgU5Gv(y`7GU1Z3^w)${BFol)^7jRc@!Rck)GCk<9?%q+30NM+4Xgq#1Yd~!
zJ3x~cVP@MEI`G{q3-QW-^B%^01`I5;Fo>7iYfPo6%}B3X3>!&xO`hLa>~nAZkycx*
zrsqF~J0DW2*EZa6n~}@V?gbc5Q8Zw!YZJBcC3yfn_&pvC6HERV!R3sWCs@a(MUBv6W=~UyvOdK5);KLVobL7r<~HC+cArt&pL4&
zTZJk@RwBqyLvIe*o0Pp4mK!Y%1gRbFc?G
z25-Ad98yyI$+L&_T^Ua*0{!q<)h*gPGO!STXkz@>4_3MfYz?sdUp?vMMMmoTw18C2
z7B`9XSx|x7Iw0z6Z~wAoV9DhD*P-i9g`*Q9dUrK!`RIQTkpW?Pk(Li0Pvsb;%x-7g
zN+M3u*sWKhvJ*S4x1%@^mt&qo*K3LjvC!T~h)lMnrb+uU)ug9tCc}-0{zg5)GByTg
z;2db#@)Imb2CfRmH68m5X;e9yu6+OuWG_x&QhYLwWti%#F`_}moeAT#9Ft4O&cAU-
z5}G<%hXS-D{tx!tahz<>_KnoL3ivt2du}2ST+*t3^_}WiTo*#0Z|)gc2u#0(|K!7c
z78VxP`{Uqq`vg8>1ZpL^BJcNl=7YTIbBTGbGWmqzhR6(jqMiSw(iH-q=jG)|OO$SK
zGZ{bQQv0lzloao$Z-&n@{Luy&lA~65$r^#aQ0`)Uyo=DP7!1
diff --git a/man/figures/index-example-2.png b/man/figures/index-example-2.png
index 3e841052f1c40cdddf0b6115f4eedbcfbd87316c..093fbf939f11e7a0d5719c42691ec3472dee93f6 100644
GIT binary patch
literal 26740
zcmd?RbySq^7d<*N(u#tVBBdZ5N=gq1A|N3l(x`MxNr%`dEr_IagQPSjDIL-%-7O_>
z&kOqb{{Fas-nH(!>;9ojVdkCpeV+52v(Mi9%yU&`IbtGeA`}WmEPq`_4TZuapit;h
z0zCMM`^yo26bk30#nr2-R#)Y&+F9EQ;UCBYHH}oeSz>IYqj=A{nPs!@a*tHxJC_
za#iKf&xRiV$iC8e;g<<&pW