diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 0000000..c2caa95
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,23 @@
+name: R Tests
+
+on:
+  push:
+  pull_request:
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+      - name: Install dependencies
+        run: |
+          install.packages(c(
+            "testthat", "tidyverse", "rstatix", "ggstatsplot", "ggsignif",
+            "see", "ggpmisc", "writexl", "readxl", "emoa", "effectsize"
+          ))
+      - name: Run tests
+        run: |
+          Rscript -e 'testthat::test_dir("tests/testthat", reporter = "summary")'
diff --git a/r_functionality.R b/r_functionality.R
index c1a6e41..ecc44a8 100644
--- a/r_functionality.R
+++ b/r_functionality.R
@@ -23,43 +23,38 @@ BibTeX:
 ")
 
 
-if (!require("httr")) install.packages("httr")
-if (!require("easystats")) install.packages("easystats")
-if (!require("tidyverse")) install.packages("tidyverse")
-if (!require("Cairo")) install.packages("Cairo")
-if (!require("rstatix")) install.packages("rstatix")
-if (!require("afex")) install.packages("afex")
-if (!require("Hmisc")) install.packages("Hmisc")
-if (!require("FSA")) install.packages("FSA")
-if (!require("PMCMRplus")) install.packages("PMCMRplus")
-if (!require("psych")) install.packages("psych")
-if (!require("pals")) install.packages("pals")
-if (!require("wesanderson")) install.packages("wesanderson")
-if (!require("ggstatsplot")) install.packages("ggstatsplot")
-if (!require("ARTool")) install.packages("ARTool")
-if (!require("pastecs")) install.packages("pastecs")
-if (!require("rstantools")) install.packages("rstantools")
-if (!require("styler")) install.packages("styler")
-if (!require("assertthat")) install.packages("assertthat")
-if (!require("reporttools")) install.packages("reporttools")
-if (!require("stargazer")) install.packages("stargazer")
-if (!require("writexl")) install.packages("writexl")
-if (!require("cli")) install.packages("cli")
-if (!require("DT")) install.packages("DT")
-if (!require("flexdashboard")) install.packages("flexdashboard")
-if (!require("sjPlot")) install.packages("sjPlot")
-if (!require("emmeans")) install.packages("emmeans")
-if (!require("stringr")) install.packages("stringr")
-if (!require("ggpmisc")) install.packages("ggpmisc")
-if (!require("ggtext")) install.packages("ggtext")
-if (!require("marginaleffects")) install.packages("marginaleffects")
-if (!require("scales")) install.packages("scales")
-if (!require("conflicted")) install.packages("conflicted")
-if (!require("curl")) install.packages("curl")
-
-
-
-library(easystats)
+auto_install_packages <- isTRUE(getOption("rcode.auto_install", TRUE))
+
+load_or_install_package <- function(pkg) {
+  is_available <- requireNamespace(pkg, quietly = TRUE)
+  if (!is_available && auto_install_packages) {
+    install.packages(pkg)
+    is_available <- requireNamespace(pkg, quietly = TRUE)
+  }
+
+  if (is_available) {
+    suppressPackageStartupMessages(library(pkg, character.only = TRUE))
+    return(TRUE)
+  }
+
+  warning(sprintf("Package '%s' is not installed. Some functionality may be unavailable.", pkg))
+  FALSE
+}
+
+rcode_packages <- c(
+  "httr", "easystats", "tidyverse", "Cairo", "rstatix", "afex", "Hmisc", "FSA",
+  "PMCMRplus", "psych", "pals", "wesanderson", "ggstatsplot", "ARTool",
+  "pastecs", "rstantools", "styler", "assertthat", "reporttools", "stargazer",
+  "writexl", "cli", "DT", "flexdashboard", "sjPlot", "emmeans", "stringr",
+  "ggpmisc", "ggtext", "marginaleffects", "scales", "conflicted", "curl",
+  "clipr", "car", "dunn.test", "xtable", "readxl", "BayesFactor", "bayestestR",
+  "foreign", "see", "marginaleffects", "effectsize", "ggsignif", "emoa",
+  "RColorBrewer"
+)
+
+invisible(lapply(unique(rcode_packages), load_or_install_package))
+
+library(easystats)
 
 ## Keep easystats current
 try(easystats::easystats_update(ask = FALSE), silent = TRUE)
@@ -71,43 +66,7 @@ try(easystats::easystats_update(ask = FALSE), silent = TRUE)
 
 # afex: necessary for ggstatsplot
 # Hmisc: necessary for mean_cl_normal --> 95% confidence intervals
-library(clipr)
-library(tidyverse)
-library(Cairo)
-library(rstatix)
-#library(nparLD)
-library(afex)
-library(Hmisc)
-library(FSA)
-library(PMCMRplus)
-library(psych)
-library(RColorBrewer)
-library(pals)
-library(wesanderson)
-library(ggstatsplot)
-library(styler)
-library(pastecs)
-library(car)
-library(dunn.test)
-library(xtable)
-library(rstantools)
-library(ARTool)
-#library(esquisse)
-library(assertthat)
-library(stargazer)
-library(reporttools)
-library(readxl)
-library(BayesFactor)
-library(bayestestR)
-library(writexl)
-library(foreign)
-library(see)
-library(sjPlot)
-library(emmeans)
-library(stringr)
-library(ggpmisc)
-library(ggtext)
-library(scales)
+library(tidyverse)
 
 
 # JANUARY 2025: no longer available
@@ -119,63 +78,63 @@ options(digits = 10)
 options(digits.secs = 3) # Set to 3 for millisecond precision
 
 ## Prefer dplyr over plyr or Hmisc in name clashes
-conflicted::conflict_prefer("mutate",   "dplyr")
-conflicted::conflict_prefer("filter",   "dplyr")
-conflicted::conflict_prefer("select",   "dplyr")
-conflicted::conflict_prefer("summarise","dplyr")
-conflicted::conflict_prefer("summarize","dplyr")
-conflicted::conflict_prefer("rename",   "dplyr")
-conflicted::conflict_prefer("arrange",  "dplyr")
-
-conflicts_prefer(scales::alpha)
-conflicts_prefer(ggplot2::annotate)
-conflicts_prefer(brms::ar)
-conflicts_prefer(car::bootCase)
-conflicts_prefer(httr::cache_info)
-conflicts_prefer(effectsize::cohens_d)
-conflicts_prefer(scales::col_factor)
-conflicts_prefer(correlation::cor_test)
-conflicts_prefer(brms::cs)
-conflicts_prefer(ggplot2::`%+%`)
-conflicts_prefer(scales::alpha)
-conflicts_prefer(ggplot2::annotate)
-conflicts_prefer(brms::ar)
-conflicts_prefer(car::bootCase)
-conflicts_prefer(httr::cache_info)
-conflicts_prefer(effectsize::cohens_d)
-conflicts_prefer(scales::col_factor)
-conflicts_prefer(correlation::cor_test)
-conflicts_prefer(brms::cs)
-conflicts_prefer(psych::describe)
-conflicts_prefer(purrr::discard)
-conflicts_prefer(report::display)
-conflicts_prefer(effectsize::eta_squared)
-conflicts_prefer(tidyr::expand)
-conflicts_prefer(tidyr::extract)
-conflicts_prefer(dplyr::filter)
-conflicts_prefer(dplyr::first)
-conflicts_prefer(insight::format_error)
-conflicts_prefer(insight::format_message)
-conflicts_prefer(insight::format_warning)
-conflicts_prefer(modelbased::get_emmeans)
-conflicts_prefer(tibble::has_name)
-conflicts_prefer(psych::headtail)
-conflicts_prefer(xtable::label)
-conflicts_prefer(Hmisc::`label<-`)
-conflicts_prefer(dplyr::lag)
-conflicts_prefer(dplyr::last)
-conflicts_prefer(lme4::lmer)
-conflicts_prefer(psych::logit)
-conflicts_prefer(lme4::ngrps)
-conflicts_prefer(tidyr::pack)
-conflicts_prefer(effectsize::phi)
-conflicts_prefer(dplyr::recode)
-conflicts_prefer(datawizard::rescale)
-conflicts_prefer(purrr::some)
-conflicts_prefer(dplyr::src)
-conflicts_prefer(dplyr::summarize)
-conflicts_prefer(devtools::test)
-conflicts_prefer(tidyr::unpack)
+try(conflicted::conflict_prefer("mutate",   "dplyr"), silent = TRUE)
+try(conflicted::conflict_prefer("filter",   "dplyr"), silent = TRUE)
+try(conflicted::conflict_prefer("select",   "dplyr"), silent = TRUE)
+try(conflicted::conflict_prefer("summarise","dplyr"), silent = TRUE)
+try(conflicted::conflict_prefer("summarize","dplyr"), silent = TRUE)
+try(conflicted::conflict_prefer("rename",   "dplyr"), silent = TRUE)
+try(conflicted::conflict_prefer("arrange",  "dplyr"), silent = TRUE)
+
+try(conflicts_prefer(scales::alpha), silent = TRUE)
+try(conflicts_prefer(ggplot2::annotate), silent = TRUE)
+try(conflicts_prefer(brms::ar), silent = TRUE)
+try(conflicts_prefer(car::bootCase), silent = TRUE)
+try(conflicts_prefer(httr::cache_info), silent = TRUE)
+try(conflicts_prefer(effectsize::cohens_d), silent = TRUE)
+try(conflicts_prefer(scales::col_factor), silent = TRUE)
+try(conflicts_prefer(correlation::cor_test), silent = TRUE)
+try(conflicts_prefer(brms::cs), silent = TRUE)
+try(conflicts_prefer(ggplot2::`%+%`), silent = TRUE)
+try(conflicts_prefer(scales::alpha), silent = TRUE)
+try(conflicts_prefer(ggplot2::annotate), silent = TRUE)
+try(conflicts_prefer(brms::ar), silent = TRUE)
+try(conflicts_prefer(car::bootCase), silent = TRUE)
+try(conflicts_prefer(httr::cache_info), silent = TRUE)
+try(conflicts_prefer(effectsize::cohens_d), silent = TRUE)
+try(conflicts_prefer(scales::col_factor), silent = TRUE)
+try(conflicts_prefer(correlation::cor_test), silent = TRUE)
+try(conflicts_prefer(brms::cs), silent = TRUE)
+try(conflicts_prefer(psych::describe), silent = TRUE)
+try(conflicts_prefer(purrr::discard), silent = TRUE)
+try(conflicts_prefer(report::display), silent = TRUE)
+try(conflicts_prefer(effectsize::eta_squared), silent = TRUE)
+try(conflicts_prefer(tidyr::expand), silent = TRUE)
+try(conflicts_prefer(tidyr::extract), silent = TRUE)
+try(conflicts_prefer(dplyr::filter), silent = TRUE)
+try(conflicts_prefer(dplyr::first), silent = TRUE)
+try(conflicts_prefer(insight::format_error), silent = TRUE)
+try(conflicts_prefer(insight::format_message), silent = TRUE)
+try(conflicts_prefer(insight::format_warning), silent = TRUE)
+try(conflicts_prefer(modelbased::get_emmeans), silent = TRUE)
+try(conflicts_prefer(tibble::has_name), silent = TRUE)
+try(conflicts_prefer(psych::headtail), silent = TRUE)
+try(conflicts_prefer(xtable::label), silent = TRUE)
+try(conflicts_prefer(Hmisc::`label<-`), silent = TRUE)
+try(conflicts_prefer(dplyr::lag), silent = TRUE)
+try(conflicts_prefer(dplyr::last), silent = TRUE)
+try(conflicts_prefer(lme4::lmer), silent = TRUE)
+try(conflicts_prefer(psych::logit), silent = TRUE)
+try(conflicts_prefer(lme4::ngrps), silent = TRUE)
+try(conflicts_prefer(tidyr::pack), silent = TRUE)
+try(conflicts_prefer(effectsize::phi), silent = TRUE)
+try(conflicts_prefer(dplyr::recode), silent = TRUE)
+try(conflicts_prefer(datawizard::rescale), silent = TRUE)
+try(conflicts_prefer(purrr::some), silent = TRUE)
+try(conflicts_prefer(dplyr::src), silent = TRUE)
+try(conflicts_prefer(dplyr::summarize), silent = TRUE)
+try(conflicts_prefer(devtools::test), silent = TRUE)
+try(conflicts_prefer(tidyr::unpack), silent = TRUE)
 
 
 # not in 
@@ -192,17 +151,46 @@ na.zero <- function (x) {
 #  Converting a Windows path to the format that works in R
 # No need for an argument. The path is printed to your console correctly and written to your clipboard for easy pasting to a script
 # From: https://stackoverflow.com/questions/8425409/file-path-issues-in-r-using-windows-hex-digits-in-character-string-error
-pathPrep <- function(path = "clipboard") {
-  y <- if (path == "clipboard") {
-    readClipboard()
-  } else {
-    cat("Please enter the path:\n\n")
-    readline()
-  }
-  x <- chartr("\\", "/", y)
-  writeClipboard(x)
-  return(x)
-}
+pathPrep <- function(path = "clipboard", read_fn = NULL, write_fn = NULL) {
+  get_clip_reader <- function() {
+    if (!is.null(read_fn)) {
+      return(read_fn)
+    }
+    if (requireNamespace("clipr", quietly = TRUE) && clipr::clipr_available()) {
+      return(clipr::read_clip)
+    }
+    if (exists("readClipboard", mode = "function")) {
+      return(get("readClipboard", mode = "function"))
+    }
+    stop("Clipboard is not available. Provide a custom `read_fn` or a direct path.")
+  }
+
+  get_clip_writer <- function() {
+    if (!is.null(write_fn)) {
+      return(write_fn)
+    }
+    if (requireNamespace("clipr", quietly = TRUE) && clipr::clipr_available()) {
+      return(clipr::write_clip)
+    }
+    if (exists("writeClipboard", mode = "function")) {
+      return(get("writeClipboard", mode = "function"))
+    }
+    return(function(...) invisible(NULL))
+  }
+
+  writer <- get_clip_writer()
+
+  y <- if (identical(path, "clipboard")) {
+    reader <- get_clip_reader()
+    reader()
+  } else {
+    path
+  }
+
+  x <- chartr("\\", "/", y)
+  writer(x)
+  return(x)
+}
 
 # for label of number of data points
 n_fun <- function(x){
@@ -292,51 +280,75 @@ normalize <- function(x_vector, old_min, old_max, new_min, new_max) {
 }
 
 
-#' Checking the version of R and effectsize as well as ggstatsplot. If not appropriate, a message for the user is generated.
-#'
-#' @return
-#' @export
-#'
-#' @examples
-checkPackageVersions <- function() {
-  required_r_version <- package_version("4.5.2")
-
-  if (getRversion() >= required_r_version) {
-    message("R Version OK!")
-  } else {
-    message("updateR()!")
-    message("Attention: novel version of RTools is required!")
-  }
-
-  check_pkg_version <- function(pkg, required_version, ok_msg, update_msg) {
-    if (!requireNamespace(pkg, quietly = TRUE)) {
-      message(sprintf("Package '%s' is not installed - %s", pkg, update_msg))
-      return(invisible(FALSE))
-    }
-
-    if (packageVersion(pkg) >= package_version(required_version)) {
-      message(ok_msg)
-    } else {
-      message(update_msg)
-    }
-
-    invisible(TRUE)
-  }
-
-  check_pkg_version(
-    pkg = "effectsize",
-    required_version = "1.0.1",
-    ok_msg = "effectsize OK!",
-    update_msg = "update effectsize - highly important!"
-  )
-
-  check_pkg_version(
-    pkg = "ggstatsplot",
-    required_version = "0.13.2",
-    ok_msg = "ggstatsplot OK!",
-    update_msg = "update ggstatsplot - highly important!"
-  )
-}
+get_latest_package_version <- function(pkg) {
+  repos <- getOption("repos")
+  if (is.null(repos) || identical(repos["CRAN"], "@CRAN@")) {
+    repos["CRAN"] <- "https://cloud.r-project.org"
+  }
+
+  tryCatch(
+    {
+      pkg_info <- available.packages(repos = repos)
+      if (pkg %in% rownames(pkg_info)) {
+        return(package_version(pkg_info[pkg, "Version"]))
+      }
+      NA
+    },
+    error = function(...) NA
+  )
+}
+
+#' Checking the version of R and effectsize as well as ggstatsplot. If not appropriate, a message for the user is generated.
+#'
+#' @return
+#' @export
+#'
+#' @examples
+checkPackageVersions <- function() {
+  required_r_version <- package_version("4.5.2")
+
+  if (getRversion() >= required_r_version) {
+    message("R Version OK!")
+  } else {
+    message("updateR()!")
+    message("Attention: novel version of RTools is required!")
+  }
+
+  check_pkg_version <- function(pkg, required_version, ok_msg, update_msg) {
+    if (!requireNamespace(pkg, quietly = TRUE)) {
+      message(sprintf("Package '%s' is not installed - %s", pkg, update_msg))
+      return(invisible(FALSE))
+    }
+
+    installed_version <- packageVersion(pkg)
+    if (installed_version >= package_version(required_version)) {
+      message(ok_msg)
+    } else {
+      message(update_msg)
+    }
+
+    latest_version <- get_latest_package_version(pkg)
+    if (!is.na(latest_version) && latest_version > installed_version) {
+      message(sprintf("A newer CRAN version of %s is available: %s", pkg, latest_version))
+    }
+
+    invisible(TRUE)
+  }
+
+  check_pkg_version(
+    pkg = "effectsize",
+    required_version = "1.0.1",
+    ok_msg = "effectsize OK!",
+    update_msg = "update effectsize - highly important!"
+  )
+
+  check_pkg_version(
+    pkg = "ggstatsplot",
+    required_version = "0.13.2",
+    ok_msg = "ggstatsplot OK!",
+    update_msg = "update ggstatsplot - highly important!"
+  )
+}
 
 
 
@@ -887,10 +899,9 @@ checkAssumptionsForAnova <- function(data, y, factors) {
 }
 
 
-#' Generate the Latex-text based on the NPAV by Lüpsen (see \url{http://www.uni-koeln.de/~luepsen/R/}).
-#' Only significant main and interaction effects are reported.
-#' P-values are rounded for the third digit.
-#' Attention: Effect sizes are not calculated!
+#' Generate the Latex-text based on the NPAV by Lüpsen (see \url{http://www.uni-koeln.de/~luepsen/R/}).
+#' Only significant main and interaction effects are reported.
+#' P-values are rounded for the third digit and partial eta squared values are provided when possible.
 #' Attention: the independent variables of the formula and the term specifying the participant must be factors (i.e., use as.factor()).
 #'
 #' To easily copy and paste the results to your manuscript, the following commands must be defined in Latex:
@@ -954,36 +965,74 @@ reportNPAV <- function(model, dv = "Testdependentvariable", write_to_clipboard =
           }
 
 
-          if (str_detect(model$descriptions[i], "X")) {
-            stringtowrite <- paste0("The NPAV found a significant interaction effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\F{", numeratordf, "}{", denominatordf, "}{", sprintf("%.2f", Fvalue), "}, ", pValue, "). ")
-          } else {
-            stringtowrite <- paste0("The NPAV found a significant main effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\F{", numeratordf, "}{", denominatordf, "}{", sprintf("%.2f", Fvalue), "}, ", pValue, "). ")
-          }
+          if (str_detect(model$descriptions[i], "X")) {
+            stringtowrite <- paste0("The NPAV found a significant interaction effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\F{", numeratordf, "}{", denominatordf, "}{", sprintf("%.2f", Fvalue), "}, ", pValue, ")")
+          } else {
+            stringtowrite <- paste0("The NPAV found a significant main effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\F{", numeratordf, "}{", denominatordf, "}{", sprintf("%.2f", Fvalue), "}, ", pValue, ")")
+          }
+
+          effect_size_text <- ""
+          if (!is.na(denominatordf) && is.finite(denominatordf)) {
+            effect_size <- tryCatch(
+              effectsize::F_to_eta2(
+                f = Fvalue,
+                df = numeratordf,
+                df_error = denominatordf,
+                ci = 0.95
+              ),
+              error = function(e) NULL
+            )
+
+            if (!is.null(effect_size)) {
+              effect_size <- as.data.frame(effect_size)
+              eta_value <- effect_size$Eta2_partial
+              ci_low <- effect_size$CI_low
+              ci_high <- effect_size$CI_high
+              if (!is.null(eta_value) && !is.na(eta_value)) {
+                effect_size_text <- paste0(
+                  ", $\\eta_{p}^{2}=",
+                  sprintf("%.2f", eta_value)
+                )
+                if (!is.null(ci_low) && !is.null(ci_high) && !any(is.na(c(ci_low, ci_high)))) {
+                  effect_size_text <- paste0(
+                    effect_size_text,
+                    " [",
+                    sprintf("%.2f", ci_low),
+                    ", ",
+                    sprintf("%.2f", ci_high),
+                    "]"
+                  )
+                }
+              }
+            }
+          }
+
+          stringtowrite <- paste0(stringtowrite, effect_size_text, ". ")
 
           # gsub backslash needs four \: https://stackoverflow.com/questions/27491986/r-gsub-replacing-backslashes
           # nice format of X in Latex via \times
           # Replace "X" with LaTeX code if preceded by a space
           stringtowrite <- gsub("(?<=\\s)X", "$\\\\times$ \\\\", stringtowrite, perl = TRUE)
 
-          if (write_to_clipboard) {
-            cat(stringtowrite)
-            write_clip(stringtowrite)
-          } else {
-            cat(stringtowrite)
-          }
-        }
-      }
+          if (write_to_clipboard) {
+            cat(stringtowrite)
+            write_clip(stringtowrite)
+          } else {
+            cat(stringtowrite)
+          }
+        }
+      }
     }
   }
 }
 
 
 
-#' Generate the Latex-text based on the NPAV by Lüpsen (see \url{http://www.uni-koeln.de/~luepsen/R/}).
-#' Only significant main and interaction effects are reported.
-#' P-values are rounded for the third digit.
-#' This is the version to report CHI Square values, which is necessary for between-subject studies!
-#' Attention: Effect sizes are not calculated!
+#' Generate the Latex-text based on the NPAV by Lüpsen (see \url{http://www.uni-koeln.de/~luepsen/R/}).
+#' Only significant main and interaction effects are reported.
+#' P-values are rounded for the third digit.
+#' This is the version to report CHI Square values, which is necessary for between-subject studies!
+#' Cohen's $w$ effect sizes are provided when a sample size is supplied.
 #' Attention: the independent variables of the formula and the term specifying the participant must be factors (i.e., use as.factor()).
 #'
 #'
@@ -992,12 +1041,13 @@ reportNPAV <- function(model, dv = "Testdependentvariable", write_to_clipboard =
 #' \newcommand{\p}{\textit{p=}}
 #' \newcommand{\pminor}{\textit{p$<$}}
 #'
-#' @param model the model of the np.anova
-#' @param dv the name of the dependent variable that should be reported
+#' @param model the model of the np.anova
+#' @param dv the name of the dependent variable that should be reported
+#' @param sample_size Optional sample size used to compute Cohen's $w$ from the chi-square statistic.
 #'
 #' @return
 #' @export
-reportNPAVChi <- function(model, dv = "Testdependentvariable", write_to_clipboard = FALSE) {
+reportNPAVChi <- function(model, dv = "Testdependentvariable", write_to_clipboard = FALSE, sample_size = NULL) {
   .Deprecated("ARTool")
   assertthat::not_empty(model)
   assertthat::not_empty(dv)
@@ -1039,13 +1089,50 @@ reportNPAVChi <- function(model, dv = "Testdependentvariable", write_to_clipboar
       }
       
       
-      if (str_detect(model$descriptions[i], "X")) {
-        
-        stringtowrite <- paste0("The NPAV found a significant interaction effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\chisq~(1)=", Chivalue, ", ", pValue, "). ")
-      
-        } else {
-        stringtowrite <- paste0("The NPAV found a significant main effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\chisq~(1)=", Chivalue, ", ",  pValue, "). ")
-      }
+      if (str_detect(model$descriptions[i], "X")) {
+
+        stringtowrite <- paste0("The NPAV found a significant interaction effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\chisq~(1)=", Chivalue, ", ", pValue, ")")
+
+        } else {
+        stringtowrite <- paste0("The NPAV found a significant main effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\chisq~(1)=", Chivalue, ", ",  pValue, ")")
+      }
+
+      effect_size_text <- ""
+      if (!is.null(sample_size) && is.numeric(sample_size) && sample_size > 0) {
+        effect_size <- tryCatch(
+          effectsize::chisq_to_w(
+            chi = Chivalue,
+            n = sample_size,
+            ci = 0.95
+          ),
+          error = function(e) NULL
+        )
+
+        if (!is.null(effect_size)) {
+          effect_size <- as.data.frame(effect_size)
+          w_value <- effect_size$Cohens_w
+          ci_low <- effect_size$CI_low
+          ci_high <- effect_size$CI_high
+          if (!is.null(w_value) && !is.na(w_value)) {
+            effect_size_text <- paste0(
+              ", $w=",
+              sprintf("%.2f", w_value)
+            )
+            if (!is.null(ci_low) && !is.null(ci_high) && !any(is.na(c(ci_low, ci_high)))) {
+              effect_size_text <- paste0(
+                effect_size_text,
+                " [",
+                sprintf("%.2f", ci_low),
+                ", ",
+                sprintf("%.2f", ci_high),
+                "]"
+              )
+            }
+          }
+        }
+      }
+
+      stringtowrite <- paste0(stringtowrite, effect_size_text, ". ")
       
       # gsub backslash needs four \: https://stackoverflow.com/questions/27491986/r-gsub-replacing-backslashes
       # nice format of X in Latex via \times
@@ -1208,9 +1295,8 @@ reportART <- function(model, dv = "Testdependentvariable", write_to_clipboard =
 
 #' Report the model produced by nparLD. The model provided must be the model generated by the command 'nparLD' \code{\link[nparLD]{nparLD}} (see \url{https://cran.r-project.org/web/packages/nparLD/nparLD.pdf}).
 #'
-#' #' Only significant main and interaction effects are reported.
-#' P-values are rounded for the third digit.
-#' Attention: Effect sizes are not calculated!
+#' #' Only significant main and interaction effects are reported.
+#' P-values are rounded for the third digit and relative treatment effects (RTE) are included when available.
 #' Attention: the independent variables of the formula and the term specifying the participant must be factors (i.e., use as.factor()).
 #'
 #' #' To easily copy and paste the results to your manuscript, the following commands must be defined in Latex:
@@ -1256,23 +1342,36 @@ reportNparLD <- function(model, dv = "Testdependentvariable", write_to_clipboard
       }
 
 
-      if (str_detect(model$descriptions[i], "X")) {
-        stringtowrite <- paste0("The NPVA found a significant interaction effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\F{", Fvalue, "}, \\df{", numeratordf, "}, ", pValue, "). ")
-      } else {
-        stringtowrite <- paste0("The NPVA found a significant main effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\F{", Fvalue, "}, \\df{", numeratordf, "}, ", pValue, "). ")
-      }
+      if (str_detect(model$descriptions[i], "X")) {
+        stringtowrite <- paste0("The NPVA found a significant interaction effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\F{", Fvalue, "}, \\df{", numeratordf, "}, ", pValue, ")")
+      } else {
+        stringtowrite <- paste0("The NPVA found a significant main effect of \\", trimws(model$descriptions[i]), " on ", dv, " (\\F{", Fvalue, "}, \\df{", numeratordf, "}, ", pValue, ")")
+      }
+
+      effect_size_text <- ""
+      if ("RTE" %in% names(model)) {
+        rte_value <- model$RTE[i]
+        if (!is.null(rte_value) && !is.na(rte_value)) {
+          effect_size_text <- paste0(
+            ", $RTE=",
+            sprintf("%.2f", rte_value)
+          )
+        }
+      }
+
+      stringtowrite <- paste0(stringtowrite, effect_size_text, ". ")
       
       # gsub backslash needs four \: https://stackoverflow.com/questions/27491986/r-gsub-replacing-backslashes
       # nice format of X in Latex via \times
       # Replace "X" with LaTeX code if preceded by a space
       stringtowrite <- gsub("(?<=\\s)X", "$\\\\times$ \\\\", stringtowrite, perl = TRUE)
       
-      if(write_to_clipboard){
-        cat(stringtowrite)
-        write_clip(stringtowrite)
-      }else{
-        cat(stringtowrite)
-      }
+      if(write_to_clipboard){
+        cat(stringtowrite)
+        write_clip(stringtowrite)
+      }else{
+        cat(stringtowrite)
+      }
     }
   }
 }
@@ -1887,10 +1986,10 @@ reshape_data <- function(input_filepath, sheetName = "Results", marker = "videoi
 #' # Add the Pareto front column
 #' main_df <- add_pareto_emoa_column(data = main_df, objectives)
 #' head(main_df)
-add_pareto_emoa_column <- function(data, objectives) {
-  # Load required library
-  if (!require("emoa")) install.packages("emoa")
-  library(emoa)
+add_pareto_emoa_column <- function(data, objectives) {
+  if (!requireNamespace("emoa", quietly = TRUE)) {
+    stop("Package 'emoa' is required for add_pareto_emoa_column().")
+  }
   
   # Input checks
   assertthat::not_empty(data)
@@ -1906,7 +2005,7 @@ add_pareto_emoa_column <- function(data, objectives) {
   }
   
   # Transpose and convert to matrix as required by the nondominated_points function
-  pareto_points <- emoa::nondominated_points(t(as.matrix(objective_data)))
+  pareto_points <- emoa::nondominated_points(t(as.matrix(objective_data)))
   
   # Convert the Pareto points matrix back to a data frame for comparison
   pareto_df <- as.data.frame(t(pareto_points))
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..e93c4af
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,3 @@
+library(testthat)
+
+testthat::test_dir("tests/testthat", reporter = "summary")
diff --git a/tests/testthat/helper-r_code.R b/tests/testthat/helper-r_code.R
new file mode 100644
index 0000000..96db164
--- /dev/null
+++ b/tests/testthat/helper-r_code.R
@@ -0,0 +1,4 @@
+library(testthat)
+options(rcode.auto_install = FALSE)
+rcode_file <- normalizePath(file.path("..", "..", "r_functionality.R"), mustWork = TRUE)
+suppressWarnings(source(rcode_file, local = FALSE))
diff --git a/tests/testthat/test_r_functionality.R b/tests/testthat/test_r_functionality.R
new file mode 100644
index 0000000..e70a41e
--- /dev/null
+++ b/tests/testthat/test_r_functionality.R
@@ -0,0 +1,309 @@
+library(testthat)
+library(dplyr)
+
+sample_df <- tibble::tibble(
+  ConditionID = rep(c("A", "B"), each = 5),
+  value = c(1:5, 2:6)
+)
+
+sample_model <- data.frame(
+  Df = c(1, 10),
+  `F value` = c(5.2, NA),
+  `Pr(>F)` = c(0.03, NA),
+  check.names = FALSE
+)
+row.names(sample_model) <- c("Factor", "Residuals")
+
+chi_model <- data.frame(
+  Df = 1,
+  ` Chi Sq` = 4.5,
+  ` Pr(>Chi)` = 0.02,
+  check.names = FALSE
+)
+row.names(chi_model) <- "Factor"
+
+nparld_model <- list(
+  ANOVA.test = data.frame(
+    Statistic = 5.1,
+    df = 1,
+    `p-value` = 0.01,
+    RTE = 0.62,
+    check.names = FALSE,
+    row.names = "Factor"
+  )
+)
+
+art_model <- data.frame(
+  Effect = c("Factor"),
+  Df = 1,
+  `F value` = 4.1,
+  `Pr(>F)` = 0.04,
+  Df.res = 12,
+  check.names = FALSE
+)
+
+report_data <- tibble::tibble(
+  group = rep(c("A", "B", "C"), each = 4),
+  score = c(1:4, 2:5, 3:6)
+)
+
+dunn_object <- list(
+  res = data.frame(
+    Comparison = c("A - B", "A - C"),
+    Z = c(2.4, 1.1),
+    P.adj = c(0.01, 0.2)
+  )
+)
+
+posthoc_stats <- list(
+  subtitle_data = data.frame(
+    estimate = 0.45,
+    p.value = 0.04,
+    statistic = 2.1,
+    df = 1,
+    df.error = 10,
+    method = "Paired t-test"
+  ),
+  pairwise_comparisons_data = data.frame(
+    p.value = 0.01,
+    group1 = "A",
+    group2 = "B"
+  )
+)
+
+basic_plot <- ggplot2::ggplot(sample_df, ggplot2::aes(x = ConditionID, y = value)) + ggplot2::geom_point()
+
+with_mock <- testthat::with_mocked_bindings
+
+
+#### basic utilities ---------------------------------------------------------
+
+test_that("basic utility helpers behave", {
+  expect_true(1 %!in% 2:5)
+  expect_false(2 %!in% 2:5)
+  expect_equal(na.zero(c(1, NA, 3)), c(1, 0, 3))
+
+  written <- NULL
+  expect_equal(
+    pathPrep(path = "C:/Temp/sample.txt", read_fn = function() stop("should not be called"), write_fn = function(x) written <<- x),
+    "C:/Temp/sample.txt"
+  )
+  expect_equal(written, "C:/Temp/sample.txt")
+
+  n_result <- n_fun(sample_df$value)
+  expect_s3_class(n_result, "data.frame")
+  expect_equal(n_result$label, paste0("n = ", length(sample_df$value)))
+
+  with_mock(`curl::has_internet` = function(...) TRUE, {
+    expect_true(havingIP())
+  })
+
+  layer <- stat_sum_df(mean)
+  expect_s3_class(layer, "LayerInstance")
+
+  expect_equal(normalize(c(1, 2, 3), 1, 3, 0, 1), c(0, 0.5, 1))
+
+  with_mock(
+    getRversion = function() package_version("4.5.2"),
+    requireNamespace = function(pkg, quietly = TRUE) TRUE,
+    packageVersion = function(pkg) package_version("9.9.9"),
+    get_latest_package_version = function(pkg) package_version("9.9.9"),
+    {
+      expect_invisible(checkPackageVersions())
+    }
+  )
+})
+
+
+#### ggstatsplot wrappers ----------------------------------------------------
+
+test_that("within and between wrappers choose correct type", {
+  skip_if_not_installed("ggstatsplot")
+  skip_if_not_installed("ggsignif")
+  mock_plot <- list()
+  data <- tibble::tibble(group = rep(c("A", "B"), each = 4), value = c(rep(0, 4), rep(1, 4)))
+
+  result <- with_mock(
+    `ggstatsplot::ggwithinstats` = function(..., type) list(type = type),
+    shapiro.test = function(...) list(p.value = 0.2),
+    {
+      ggwithinstatsWithPriorNormalityCheck(data, "group", "value", "Value")
+    }
+  )
+  expect_equal(result$type, "p")
+
+  np_result <- with_mock(
+    `ggstatsplot::ggwithinstats` = function(..., type) list(type = type),
+    shapiro.test = function(...) list(p.value = 0.001),
+    {
+      ggwithinstatsWithPriorNormalityCheck(data, "group", "value", "Value")
+    }
+  )
+  expect_equal(np_result$type, "np")
+
+  between <- with_mock(
+    `ggstatsplot::ggbetweenstats` = function(..., type) list(type = type),
+    shapiro.test = function(...) list(p.value = 0.001),
+    pairwise_comparisons = function(...) data.frame(group1 = "A", group2 = "B", `p.value` = 0.01, stringsAsFactors = FALSE),
+    `ggsignif::geom_signif` = function(...) ggplot2::geom_blank(),
+    {
+      ggbetweenstatsWithPriorNormalityCheck(data, "group", "value", "Value", c("A", "B"))
+    }
+  )
+  expect_equal(between$type, "np")
+
+  expect_s3_class(
+    with_mock(
+      `ggstatsplot::ggbetweenstats` = function(...) ggplot2::ggplot(),
+      pairwise_comparisons = function(...) data.frame(group1 = "A", group2 = "B", `p.value` = 0.01, stringsAsFactors = FALSE),
+      `ggsignif::geom_signif` = function(...) ggplot2::geom_blank(),
+      {
+        ggbetweenstatsWithPriorNormalityCheckAsterisk(data, "group", "value", "Value", c("A", "B"))
+      }
+    ),
+    "ggplot"
+  )
+
+  expect_s3_class(
+    with_mock(
+      `ggstatsplot::ggwithinstats` = function(...) ggplot2::ggplot(),
+      pairwise_comparisons = function(...) data.frame(group1 = "A", group2 = "B", `p.value` = 0.01, stringsAsFactors = FALSE),
+      `ggsignif::geom_signif` = function(...) ggplot2::geom_blank(),
+      shapiro.test = function(...) list(p.value = 0.2),
+      {
+        ggwithinstatsWithPriorNormalityCheckAsterisk(data, "group", "value", "Value", c("A", "B"))
+      }
+    ),
+    "ggplot"
+  )
+})
+
+
+#### effect size helpers -----------------------------------------------------
+
+test_that("effect size helpers print expected summaries", {
+  wilcox_obj <- list(p.value = 0.04, data.name = "Sample")
+  expect_output(rFromWilcox(wilcox_obj, 20), "Effect Size")
+  expect_output(rFromWilcoxAdjusted(wilcox_obj, 20, 2), "Effect Size")
+  expect_output(rFromNPAV(0.02, 30), "\\effectsize")
+})
+
+
+#### debugging and assumptions ----------------------------------------------
+
+test_that("debugging and assumption helpers work", {
+  df <- data.frame(a = c(1, 2, 3), b = c("x", "y", "z"))
+  expect_type(debug_contr_error(df)$nlevels, "integer")
+
+  anova_df <- tibble::tibble(
+    factor1 = rep(c("A", "B"), each = 10),
+    factor2 = rep(c("X", "Y"), times = 10),
+    outcome = rnorm(20)
+  )
+  expect_null(checkAssumptionsForAnova(anova_df, "outcome", c("factor1", "factor2")))
+})
+
+
+#### reporting helpers -------------------------------------------------------
+
+test_that("reporting helpers include effect sizes", {
+  expect_match(capture.output(reportNPAV(sample_model, dv = "score")), "eta")
+  expect_match(capture.output(reportNPAVChi(chi_model, dv = "score", sample_size = 30)), "w=")
+  expect_match(capture.output(reportART(art_model, dv = "score")), "eta")
+  expect_match(capture.output(reportNparLD(nparld_model, dv = "score")), "RTE")
+  expect_output(reportMeanAndSD(sample_df, iv = "ConditionID", dv = "value"), "m")
+
+  expect_output(reportDunnTest(dunn_object, report_data, iv = "group", dv = "score"), "post-hoc")
+  expect_output(reportDunnTestTable(dunn_object, report_data, iv = "group", dv = "score"), "Post-hoc")
+
+  expect_match(
+    capture.output(
+      with_mock(
+        extract_stats = function(...) posthoc_stats,
+        {
+          reportggstatsplot(basic_plot, iv = "group", dv = "score")
+        }
+      )
+    ),
+    "significant"
+  )
+
+  expect_output(
+    with_mock(
+      extract_stats = function(...) posthoc_stats,
+      {
+        reportggstatsplotPostHoc(report_data, basic_plot, iv = "group", dv = "score")
+      }
+    ),
+    "post-hoc"
+  )
+})
+
+
+#### data shaping ------------------------------------------------------------
+
+test_that("data wrangling helpers behave", {
+  skip_if_not_installed("writexl")
+  skip_if_not_installed("readxl")
+  replaced <- replace_values(data.frame(x = c("neg2", "neg1", "0")), c("neg2", "neg1"), c("-2", "-1"))
+  expect_equal(replaced$x[1:2], c("-2", "-1"))
+
+  input_path <- tempfile(fileext = ".xlsx")
+  output_path <- tempfile(fileext = ".xlsx")
+  writexl::write_xlsx(data.frame(
+    ID = 1:2,
+    videoinfo_1 = c("a", "b"),
+    Q1 = 1:2,
+    Q2 = 3:4,
+    videoinfo_2 = c("c", "d"),
+    Q3 = 5:6,
+    Q4 = 7:8
+  ), input_path)
+  reshape_data(input_path, marker = "videoinfo", output_filepath = output_path)
+  expect_true(file.exists(output_path))
+
+  skip_if_not_installed("emoa")
+  pareto_df <- add_pareto_emoa_column(
+    data = data.frame(trust = c(1, 2), predictability = c(2, 1)),
+    objectives = c("trust", "predictability")
+  )
+  expect_true("PARETO_EMOA" %in% names(pareto_df))
+
+  rei <- remove_outliers_REI(data.frame(matrix(sample(1:5, 20, replace = TRUE), ncol = 4)), range = c(1, 5))
+  expect_true("REI" %in% names(rei))
+})
+
+
+#### plotting helpers --------------------------------------------------------
+
+test_that("plotting helpers return ggplot objects", {
+  skip_if_not_installed("see")
+  skip_if_not_installed("ggpmisc")
+  expect_s3_class(generateEffectPlot(sample_df, "ConditionID", "value", "ConditionID"), "ggplot")
+
+  mobo_df <- tibble::tibble(
+    Iteration = rep(1:6, each = 2),
+    trust = rnorm(12),
+    ConditionID = rep(c("A", "B"), times = 6)
+  )
+  expect_s3_class(generateMoboPlot(mobo_df, "Iteration", "trust"), "ggplot")
+
+  mobo_df2 <- tibble::tibble(
+    Iteration = rep(1:6, each = 2),
+    trust = rnorm(12),
+    Phase = rep(c("sampling", "optimization"), each = 6),
+    ConditionID = rep(c("A", "B"), times = 6)
+  )
+  expect_s3_class(generateMoboPlot2(mobo_df2, x = "Iteration", y = "trust", fillColourGroup = "ConditionID"), "ggplot")
+})
+
+
+#### latex and misc ---------------------------------------------------------
+
+test_that("latex helper and np.anova behave", {
+  text <- "- significant effect\n- non-significant effect\nStandardized parameters were obtained by fitting the model"
+  expect_match(latexify_report(text, only_sig = TRUE, remove_std = TRUE), "\\\\item")
+
+  factor_df <- data.frame(y = factor(c("a", "b")), x = factor(c("a", "b")))
+  expect_error(np.anova(y ~ x, data = factor_df), "invalid type")
+})