diff --git a/README.Rmd b/README.Rmd
index 4b4aa44f..f31d029b 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -207,12 +207,12 @@ plot(contrasts, estimate_means(model)) +
 ```{r}
 model <- lm(Sepal.Width ~ Species * Petal.Length, data = iris)
 
-estimate_contrasts(model, at = "Petal.Length", length = 3)
+estimate_contrasts(model, by = "Petal.Length", length = 3)
 ```
 
 ```{r}
 # Recompute contrasts with a higher precision (for a smoother plot)
-contrasts <- estimate_contrasts(model, at = "Petal.Length", length = 20)
+contrasts <- estimate_contrasts(model, by = "Petal.Length", length = 20)
 
 # Add Contrast column by concatenating
 contrasts$Contrast <- paste(contrasts$Level1, "-", contrasts$Level2)
@@ -302,7 +302,7 @@ model <- mgcv::gam(Sepal.Width ~ s(Petal.Length), data = iris)
 # 1. Compute derivatives
 deriv <- estimate_slopes(model,
   trend = "Petal.Length",
-  at = "Petal.Length",
+  by = "Petal.Length",
   length = 100
 )
 
@@ -363,7 +363,7 @@ For instance, the plot below shows that the effect of `hp` (the y-axis) is signi
 ```{r}
 model <- lm(mpg ~ hp * wt, data = mtcars)
 
-slopes <- estimate_slopes(model, trend = "hp", at = "wt")
+slopes <- estimate_slopes(model, trend = "hp", by = "wt")
 
 plot(slopes)
 ```
diff --git a/README.md b/README.md
index 4fa0d68f..1ce9cc8b 100644
--- a/README.md
+++ b/README.md
@@ -163,7 +163,10 @@ means
 ## versicolor | 2.77 | 0.05 | [2.68, 2.86]
 ## virginica  | 2.97 | 0.05 | [2.88, 3.07]
 ## 
-## Marginal means estimated at Species
+## Marginal means estimated at
+```
+
+``` r
 
 # 3. Plot
 ggplot(iris, aes(x = Species, y = Sepal.Width)) +
@@ -233,7 +236,7 @@ contrasts
 ``` r
 model <- lm(Sepal.Width ~ Species * Petal.Length, data = iris)
 
-estimate_contrasts(model, at = "Petal.Length", length = 3)
+estimate_contrasts(model, by = "Petal.Length", length = 3)
 ## Marginal Contrasts Analysis
 ## 
 ## Level1     |     Level2 | Petal.Length | Difference |        95% CI |   SE | t(144) |      p
@@ -254,7 +257,7 @@ estimate_contrasts(model, at = "Petal.Length", length = 3)
 
 ``` r
 # Recompute contrasts with a higher precision (for a smoother plot)
-contrasts <- estimate_contrasts(model, at = "Petal.Length", length = 20)
+contrasts <- estimate_contrasts(model, by = "Petal.Length", length = 20)
 
 # Add Contrast column by concatenating
 contrasts$Contrast <- paste(contrasts$Level1, "-", contrasts$Level2)
@@ -306,6 +309,9 @@ head(pred1, n = 5)
 ## 5.00         |      2.19 | 0.10 | [1.99, 2.39] |     -0.79 |         1.40
 ## 
 ## Variable predicted: Petal.Length
+```
+
+``` r
 
 # Same for model 2
 model2 <- lm(Petal.Length ~ Sepal.Length * Species, data = iris)
@@ -354,6 +360,9 @@ random
 ## cyl   |     6 |        drat |       -0.09 | 0.54 | [-1.15,  0.98]
 ## cyl   |     8 | (Intercept) |        3.32 | 0.73 | [ 1.89,  4.74]
 ## cyl   |     8 |        drat |       -2.15 | 0.47 | [-3.07, -1.23]
+```
+
+``` r
 
 plot(random)
 ```
@@ -391,7 +400,7 @@ model <- mgcv::gam(Sepal.Width ~ s(Petal.Length), data = iris)
 # 1. Compute derivatives
 deriv <- estimate_slopes(model,
   trend = "Petal.Length",
-  at = "Petal.Length",
+  by = "Petal.Length",
   length = 100
 )
 
@@ -465,7 +474,7 @@ is significantly negative only when `wt` is low (`< ~4`).
 ``` r
 model <- lm(mpg ~ hp * wt, data = mtcars)
 
-slopes <- estimate_slopes(model, trend = "hp", at = "wt")
+slopes <- estimate_slopes(model, trend = "hp", by = "wt")
 
 plot(slopes)
 ```
diff --git a/man/figures/unnamed-chunk-10-1.png b/man/figures/unnamed-chunk-10-1.png
index f9edfe30..02771cd0 100644
Binary files a/man/figures/unnamed-chunk-10-1.png and b/man/figures/unnamed-chunk-10-1.png differ
diff --git a/man/figures/unnamed-chunk-12-1.png b/man/figures/unnamed-chunk-12-1.png
index e2e1f938..15df600a 100644
Binary files a/man/figures/unnamed-chunk-12-1.png and b/man/figures/unnamed-chunk-12-1.png differ
diff --git a/man/figures/unnamed-chunk-14-1.png b/man/figures/unnamed-chunk-14-1.png
index be9dad67..8bf298bc 100644
Binary files a/man/figures/unnamed-chunk-14-1.png and b/man/figures/unnamed-chunk-14-1.png differ
diff --git a/man/figures/unnamed-chunk-15-1.png b/man/figures/unnamed-chunk-15-1.png
index bc080a1b..e53708b2 100644
Binary files a/man/figures/unnamed-chunk-15-1.png and b/man/figures/unnamed-chunk-15-1.png differ
diff --git a/man/figures/unnamed-chunk-16-1.png b/man/figures/unnamed-chunk-16-1.png
index 46caa64a..0431b604 100644
Binary files a/man/figures/unnamed-chunk-16-1.png and b/man/figures/unnamed-chunk-16-1.png differ
diff --git a/man/figures/unnamed-chunk-17-1.png b/man/figures/unnamed-chunk-17-1.png
index f7be5b17..87bf15e7 100644
Binary files a/man/figures/unnamed-chunk-17-1.png and b/man/figures/unnamed-chunk-17-1.png differ
diff --git a/man/figures/unnamed-chunk-3-1.png b/man/figures/unnamed-chunk-3-1.png
index c33bd5fa..6fd5b8b0 100644
Binary files a/man/figures/unnamed-chunk-3-1.png and b/man/figures/unnamed-chunk-3-1.png differ
diff --git a/man/figures/unnamed-chunk-4-1.png b/man/figures/unnamed-chunk-4-1.png
index e0555599..056f0037 100644
Binary files a/man/figures/unnamed-chunk-4-1.png and b/man/figures/unnamed-chunk-4-1.png differ
diff --git a/man/figures/unnamed-chunk-6-1.png b/man/figures/unnamed-chunk-6-1.png
index 0b063d12..4e5e9527 100644
Binary files a/man/figures/unnamed-chunk-6-1.png and b/man/figures/unnamed-chunk-6-1.png differ
diff --git a/man/figures/unnamed-chunk-8-1.png b/man/figures/unnamed-chunk-8-1.png
index cdf07ef7..c2510ab5 100644
Binary files a/man/figures/unnamed-chunk-8-1.png and b/man/figures/unnamed-chunk-8-1.png differ
diff --git a/man/figures/unnamed-chunk-9-1.png b/man/figures/unnamed-chunk-9-1.png
index 60dc7e93..3caeb4e0 100644
Binary files a/man/figures/unnamed-chunk-9-1.png and b/man/figures/unnamed-chunk-9-1.png differ
diff --git a/tests/testthat/test-attributes_visualisation.R b/tests/testthat/test-attributes_visualisation.R
index 7cc27833..6751932e 100644
--- a/tests/testthat/test-attributes_visualisation.R
+++ b/tests/testthat/test-attributes_visualisation.R
@@ -3,15 +3,15 @@ test_that("attributes_means", {
   model <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
 
   estim <- suppressMessages(estimate_means(model))
-  expect_identical(attributes(estim)$at, "Species")
-  expect_identical(attributes(estim)$fixed, NULL)
+  expect_identical(attributes(estim)$by, "Species")
+  expect_null(attributes(estim)$fixed)
 
   estim <- suppressMessages(estimate_means(model, fixed = "Sepal.Width"))
-  expect_identical(attributes(estim)$at, "Species")
+  expect_identical(attributes(estim)$by, "Species")
   expect_identical(attributes(estim)$fixed, "Sepal.Width")
 
-  estim <- suppressMessages(estimate_means(model, at = "all"))
-  expect_identical(attributes(estim)$at, c("Species", "Sepal.Width"))
+  estim <- suppressMessages(estimate_means(model, by = "all"))
+  expect_identical(attributes(estim)$by, c("Species", "Sepal.Width"))
 })
 
 
@@ -22,13 +22,13 @@ test_that("attributes_contrasts", {
 
   estim <- suppressMessages(estimate_contrasts(model))
   expect_identical(attributes(estim)$contrast, "Species")
-  expect_identical(attributes(estim)$at, NULL)
-  expect_identical(attributes(estim)$fixed, NULL)
+  expect_null(attributes(estim)$by)
+  expect_null(attributes(estim)$fixed)
 
   estim <- suppressMessages(estimate_contrasts(model, fixed = "Sepal.Width"))
   expect_identical(attributes(estim)$contrast, "Species")
   expect_identical(attributes(estim)$fixed, "Sepal.Width")
-  expect_identical(attributes(estim)$modulate, NULL)
+  expect_null(attributes(estim)$modulate)
 })
 
 
diff --git a/tests/testthat/test-brms.R b/tests/testthat/test-brms.R
index 2350ae10..a71177b7 100644
--- a/tests/testthat/test-brms.R
+++ b/tests/testthat/test-brms.R
@@ -5,7 +5,7 @@ test_that("estimate_means - brms", {
   skip_if_not_installed("emmeans")
   model <- brms::brm(Sepal.Length ~ Species * Sepal.Width, data = iris, refresh = 0, iter = 1000)
   estim <- estimate_means(model)
-  expect_equal(dim(estim), c(3, 5))
+  expect_identical(dim(estim), c(3L, 5L))
 })
 
 test_that("estimate_relation - brms", {
@@ -13,7 +13,7 @@ test_that("estimate_relation - brms", {
   skip_if_not_installed("emmeans")
   model <- brms::brm(Sepal.Length ~ Species * Sepal.Width, data = iris, refresh = 0, iter = 1000)
   estim <- estimate_relation(model, preserve_range = FALSE)
-  expect_equal(dim(estim), c(30, 6))
+  expect_identical(dim(estim), c(30L, 6L))
 
   # estim <- estimate_relation(model, preserve_range=FALSE, iterations = 10)
   # expect_equal(dim(estim), c(30, 6))
@@ -23,6 +23,6 @@ test_that("estimate_means - brms", {
   skip_if_not_installed("brms")
   skip_if_not_installed("emmeans")
   model <- brms::brm(Sepal.Length ~ Species * Sepal.Width, data = iris, refresh = 0, iter = 1000)
-  estim <- estimate_slopes(model, at = "Species")
-  expect_equal(dim(estim), c(3, 5))
+  estim <- estimate_slopes(model, by = "Species")
+  expect_identical(dim(estim), c(3L, 5L))
 })
diff --git a/tests/testthat/test-estimate_means.R b/tests/testthat/test-estimate_means.R
index 5eee7896..e686d873 100644
--- a/tests/testthat/test-estimate_means.R
+++ b/tests/testthat/test-estimate_means.R
@@ -13,26 +13,26 @@ test_that("estimate_means() - core", {
   # Simple
   model <- lm(vs ~ cyl, data = dat)
   estim1 <- suppressMessages(estimate_means(model))
-  expect_equal(dim(estim1), c(3, 5))
+  expect_identical(dim(estim1), c(3L, 5L))
   estim2 <- suppressMessages(estimate_means(model, backend = "marginaleffects"))
-  expect_equal(dim(estim2), c(3, 5))
-  expect_true(max(estim1$Mean - estim2$Mean) < 1e-10)
+  expect_identical(dim(estim2), c(3L, 5L))
+  expect_lt(max(estim1$Mean - estim2$Mean), 1e-10)
 
   # Interaction (factor * continuous)
   model <- lm(mpg ~ wt * gear, data = dat)
   estim1 <- suppressMessages(estimate_means(model))
-  expect_equal(dim(estim1), c(3, 5))
+  expect_identical(dim(estim1), c(3L, 5L))
   estim2 <- suppressMessages(estimate_means(model, backend = "marginaleffects"))
-  expect_equal(dim(estim2), c(3, 6))
-  expect_true(max(estim1$Mean - estim2$Mean) < 1e-10)
+  expect_identical(dim(estim2), c(3L, 6L))
+  expect_lt(max(estim1$Mean - estim2$Mean), 1e-10)
 
   # At specific levels
   model <- lm(Sepal.Width ~ Species, data = iris)
-  estim1 <- suppressMessages(estimate_means(model, at = "Species=c('versicolor', 'virginica')"))
-  expect_equal(dim(estim1), c(2, 5))
-  estim2 <- suppressMessages(estimate_means(model, at = "Species=c('versicolor', 'virginica')", backend = "marginaleffects"))
-  expect_equal(dim(estim2), c(2, 5))
-  expect_true(max(estim1$Mean - estim2$Mean) < 1e-10)
+  estim1 <- suppressMessages(estimate_means(model, by = "Species=c('versicolor', 'virginica')"))
+  expect_identical(dim(estim1), c(2L, 5L))
+  estim2 <- suppressMessages(estimate_means(model, by = "Species=c('versicolor', 'virginica')", backend = "marginaleffects"))
+  expect_identical(dim(estim2), c(2L, 5L))
+  expect_lt(max(estim1$Mean - estim2$Mean), 1e-10)
 
   # Interactions between factors
   dat <- iris
@@ -40,24 +40,24 @@ test_that("estimate_means() - core", {
   dat <<- dat
 
   model <- lm(Sepal.Width ~ Species * Petal.Length_factor, data = dat)
-  estim1 <- suppressMessages(estimate_means(model, at = "all"))
-  expect_equal(dim(estim1), c(6, 6))
-  estim2 <- suppressWarnings(suppressMessages(estimate_means(model, at = "all", backend = "marginaleffects")))
-  expect_equal(dim(estim2), c(6, 6))
+  estim1 <- suppressMessages(estimate_means(model, by = "all"))
+  expect_identical(dim(estim1), c(6L, 6L))
+  estim2 <- suppressWarnings(suppressMessages(estimate_means(model, by = "all", backend = "marginaleffects")))
+  expect_identical(dim(estim2), c(6L, 6L))
 
   # No interaction (two factors)
   model <- lm(Petal.Length ~ Sepal.Width + Species, data = iris)
   estim1 <- suppressMessages(estimate_means(model))
-  expect_equal(dim(estim1), c(3, 5))
+  expect_identical(dim(estim1), c(3L, 5L))
   estim2 <- suppressMessages(estimate_means(model, backend = "marginaleffects"))
-  expect_equal(dim(estim2), c(3, 6))
-  expect_true(max(estim1$Mean - estim2$Mean) < 1e-10)
+  expect_identical(dim(estim2), c(3L, 6L))
+  expect_lt(max(estim1$Mean - estim2$Mean), 1e-10)
 
 
   # At specific levels of continuous
-  estim1 <- suppressMessages(estimate_means(model, at = "Sepal.Width"))
+  estim1 <- suppressMessages(estimate_means(model, by = "Sepal.Width"))
   expect_equal(dim(estim1), c(10, 5))
-  estim2 <- suppressMessages(estimate_means(model, at = "Sepal.Width", backend = "marginaleffects"))
+  estim2 <- suppressMessages(estimate_means(model, by = "Sepal.Width", backend = "marginaleffects"))
   expect_equal(dim(estim2), c(10, 6))
   # Note that the absolute values are different here... for unclear reasons
   expect_true(max(diff(estim1$Mean) - diff(estim2$Mean)) < 1e-10)
@@ -69,44 +69,44 @@ test_that("estimate_means() - core", {
   model <- glm(y ~ Species, family = "binomial", data = dat)
 
   estim <- suppressMessages(estimate_means(model))
-  expect_equal(dim(estim), c(3, 5))
+  expect_identical(dim(estim), c(3L, 5L))
   estim <- suppressMessages(estimate_means(model, transform = "response"))
-  expect_equal(dim(estim), c(3, 5))
+  expect_identical(dim(estim), c(3L, 5L))
   expect_true(all(estim$Probability >= 0) & all(estim$Probability <= 1))
 
 
   model <- lm(Petal.Length ~ Sepal.Width + Species, data = iris)
   estim <- suppressMessages(estimate_means(model))
-  expect_equal(dim(estim), c(3, 5))
+  expect_identical(dim(estim), c(3L, 5L))
 
-  estim <- suppressMessages(estimate_means(model, at = "all"))
+  estim <- suppressMessages(estimate_means(model, by = "all"))
   expect_equal(dim(estim), c(30, 6))
 
   # In formula modification
   # FIXME: this got broken but it seems just to tedious to fix. Don't use in formula transforms.
   # model <- lm(mpg ~ wt * as.factor(gear), data = mtcars)
   # estim <- suppressMessages(estimate_means(model))
-  # expect_equal(dim(estim), c(3, 5))
+  # expect_equal(dim(estim), c(3L, 5L))
 
   # One continuous and one factor
   model <- lm(Petal.Length ~ Species * Sepal.Width, data = iris)
 
   estim <- suppressMessages(estimate_means(model))
-  expect_equal(dim(estim), c(3, 5))
+  expect_identical(dim(estim), c(3L, 5L))
   estim <- suppressMessages(estimate_means(model, fixed = "Sepal.Width"))
-  expect_equal(dim(estim), c(3, 6))
-  estim <- suppressMessages(estimate_means(model, at = c("Species", "Sepal.Width"), length = 2))
-  expect_equal(dim(estim), c(6, 6))
-  estim <- suppressMessages(estimate_means(model, at = "Species=c('versicolor', 'setosa')"))
-  expect_equal(dim(estim), c(2, 5))
-  estim <- suppressMessages(estimate_means(model, at = "Sepal.Width=c(2, 4)"))
-  expect_equal(dim(estim), c(2, 5))
-  estim <- suppressMessages(estimate_means(model, at = c("Species", "Sepal.Width=0")))
-  expect_equal(dim(estim), c(3, 6))
-  estim <- suppressMessages(estimate_means(model, at = "Sepal.Width", length = 5))
+  expect_identical(dim(estim), c(3L, 6L))
+  estim <- suppressMessages(estimate_means(model, by = c("Species", "Sepal.Width"), length = 2))
+  expect_identical(dim(estim), c(6L, 6L))
+  estim <- suppressMessages(estimate_means(model, by = "Species=c('versicolor', 'setosa')"))
+  expect_identical(dim(estim), c(2L, 5L))
+  estim <- suppressMessages(estimate_means(model, by = "Sepal.Width=c(2, 4)"))
+  expect_identical(dim(estim), c(2L, 5L))
+  estim <- suppressMessages(estimate_means(model, by = c("Species", "Sepal.Width=0")))
+  expect_identical(dim(estim), c(3L, 6L))
+  estim <- suppressMessages(estimate_means(model, by = "Sepal.Width", length = 5))
   expect_equal(dim(estim), c(5, 5))
-  estim <- suppressMessages(estimate_means(model, at = "Sepal.Width=c(2, 4)"))
-  expect_equal(dim(estim), c(2, 5))
+  estim <- suppressMessages(estimate_means(model, by = "Sepal.Width=c(2, 4)"))
+  expect_identical(dim(estim), c(2L, 5L))
 
   # Two factors
   dat <- iris
@@ -114,12 +114,12 @@ test_that("estimate_means() - core", {
   dat <<- dat
   model <- lm(Petal.Length ~ Species * Petal.Length_factor, data = dat)
 
-  estim <- suppressMessages(estimate_means(model, at = "all"))
-  expect_equal(dim(estim), c(6, 6))
-  estim <- suppressMessages(estimate_means(model, at = "Petal.Length_factor"))
-  expect_equal(dim(estim), c(2, 5))
-  estim <- suppressMessages(estimate_means(model, at = "Petal.Length_factor='B'"))
-  expect_true(as.character(unique(estim$Petal.Length_factor)) == "B")
+  estim <- suppressMessages(estimate_means(model, by = "all"))
+  expect_identical(dim(estim), c(6L, 6L))
+  estim <- suppressMessages(estimate_means(model, by = "Petal.Length_factor"))
+  expect_identical(dim(estim), c(2L, 5L))
+  estim <- suppressMessages(estimate_means(model, by = "Petal.Length_factor='B'"))
+  expect_identical(as.character(unique(estim$Petal.Length_factor)), "B")
 
 
   # Three factors
@@ -132,7 +132,7 @@ test_that("estimate_means() - core", {
   expect_equal(dim(estim), c(12, 7))
   estim <- suppressMessages(estimate_means(model, fixed = "am"))
   expect_equal(dim(estim), c(6, 7))
-  estim <- suppressMessages(estimate_means(model, at = c("gear='5'", "vs")))
+  estim <- suppressMessages(estimate_means(model, by = c("gear='5'", "vs")))
   expect_equal(dim(estim), c(2, 7))
 
   dat <- iris
@@ -156,13 +156,13 @@ test_that("estimate_means() - core", {
   model <- glm(Petal.Length_factor ~ Species, data = dat, family = "binomial")
 
   estim <- suppressMessages(estimate_means(model))
-  expect_equal(dim(estim), c(3, 5))
+  expect_identical(dim(estim), c(3L, 5L))
   estim <- suppressMessages(estimate_means(model, transform = "none"))
-  expect_equal(dim(estim), c(3, 5))
+  expect_identical(dim(estim), c(3L, 5L))
 
   model <- glm(Petal.Length ~ Species, data = iris, family = "Gamma")
   estim <- suppressMessages(estimate_means(model))
-  expect_equal(dim(estim), c(3, 5))
+  expect_identical(dim(estim), c(3L, 5L))
 })
 
 test_that("estimate_means() - mixed models", {
@@ -176,15 +176,15 @@ test_that("estimate_means() - mixed models", {
   model <- lme4::lmer(Sepal.Width ~ Species + (1 | Petal.Length_factor), data = dat)
 
   estim1 <- suppressMessages(estimate_means(model))
-  expect_equal(dim(estim1), c(3, 5))
+  expect_identical(dim(estim1), c(3L, 5L))
   estim2 <- suppressMessages(estimate_means(model, backend = "marginaleffects"))
-  expect_equal(dim(estim2), c(3, 5))
-  expect_true(max(estim1$Mean - estim2$Mean) < 1e-10)
+  expect_identical(dim(estim2), c(3L, 5L))
+  expect_lt(max(estim1$Mean - estim2$Mean), 1e-10)
 
   model <- lme4::glmer(Sepal.Width ~ Species + (1 | Petal.Length_factor), data = dat, family = "Gamma")
   estim1 <- suppressMessages(estimate_means(model))
-  expect_equal(dim(estim1), c(3, 5))
+  expect_identical(dim(estim1), c(3L, 5L))
   estim2 <- suppressMessages(estimate_means(model, backend = "marginaleffects"))
-  expect_equal(dim(estim2), c(3, 5))
-  expect_true(max(estim1$Mean - estim2$Mean) < 1e-10)
+  expect_identical(dim(estim2), c(3L, 5L))
+  expect_lt(max(estim1$Mean - estim2$Mean), 1e-10)
 })
diff --git a/tests/testthat/test-estimate_predicted.R b/tests/testthat/test-estimate_predicted.R
index 530887a2..17ee21e5 100644
--- a/tests/testthat/test-estimate_predicted.R
+++ b/tests/testthat/test-estimate_predicted.R
@@ -137,7 +137,7 @@ test_that("estimate_response - Frequentist", {
   expect_equal(dim(estim), c(32, 4))
 
   model <- glm(vs ~ wt + cyl, data = mtcars, family = "binomial")
-  estim <- estimate_link(model, at = "wt")
+  estim <- estimate_link(model, by = "wt")
   expect_equal(dim(estim), c(10, 6))
 
 
diff --git a/tests/testthat/test-estimate_slopes.R b/tests/testthat/test-estimate_slopes.R
index 56dd3a79..301783fb 100644
--- a/tests/testthat/test-estimate_slopes.R
+++ b/tests/testthat/test-estimate_slopes.R
@@ -6,21 +6,21 @@ test_that("estimate_slopes", {
   estim <- suppressMessages(estimate_slopes(model))
   expect_equal(dim(estim), c(1, 8))
 
-  estim <- suppressMessages(estimate_slopes(model, at = "Species"))
+  estim <- suppressMessages(estimate_slopes(model, by = "Species"))
   expect_equal(dim(estim), c(3, 9))
 
-  estim <- suppressMessages(estimate_slopes(model, at = "Petal.Length"))
+  estim <- suppressMessages(estimate_slopes(model, by = "Petal.Length"))
   expect_equal(dim(estim), c(10, 9))
 
-  estim <- suppressMessages(estimate_slopes(model, at = c("Species", "Petal.Length")))
+  estim <- suppressMessages(estimate_slopes(model, by = c("Species", "Petal.Length")))
   expect_equal(dim(estim), c(30, 10))
 
   model <- lm(Petal.Length ~ poly(Sepal.Width, 4), data = iris)
 
-  estim <- suppressMessages(estimate_slopes(model, at = "Sepal.Width"))
+  estim <- suppressMessages(estimate_slopes(model, by = "Sepal.Width"))
   expect_equal(dim(estim), c(10, 9))
-  estim <- suppressMessages(estimate_slopes(model, at = "Sepal.Width", length = 5))
+  estim <- suppressMessages(estimate_slopes(model, by = "Sepal.Width", length = 5))
   expect_equal(dim(estim), c(5, 9))
-  estim <- suppressMessages(estimate_slopes(model, at = "Sepal.Width = c(1, 2, 3)"))
+  estim <- suppressMessages(estimate_slopes(model, by = "Sepal.Width = c(1, 2, 3)"))
   expect_equal(dim(estim), c(3, 9))
 })
diff --git a/tests/testthat/test-get_marginaleffects.R b/tests/testthat/test-get_marginaleffects.R
index d0e8826c..15dbbc9a 100644
--- a/tests/testthat/test-get_marginaleffects.R
+++ b/tests/testthat/test-get_marginaleffects.R
@@ -4,7 +4,7 @@ test_that("get_marginaleffects", {
 
   model <- lm(Sepal.Width ~ Species * Petal.Length, data = iris)
 
-  expect_equal(nrow(get_marginaleffects(model, trend = "Petal.Length", at = "Species")), 3L)
+  expect_equal(nrow(get_marginaleffects(model, trend = "Petal.Length", by = "Species")), 3L)
 
-  # get_marginaleffects(model, trend = "Petal.Length", at = "Species", length = 10)
+  # get_marginaleffects(model, trend = "Petal.Length", by = "Species", length = 10)
 })
diff --git a/tests/testthat/test-glmmTMB.R b/tests/testthat/test-glmmTMB.R
index 61d72c8b..da0988a8 100644
--- a/tests/testthat/test-glmmTMB.R
+++ b/tests/testthat/test-glmmTMB.R
@@ -41,7 +41,7 @@ test_that("estimate_contrasts - glmmTMB", {
 test_that("estimate_slope - glmmTMB", {
   estim <- suppressMessages(estimate_slopes(model2,
     trend = "cover",
-    at = "mined", regrid = "response"
+    by = "mined", regrid = "response"
   ))
   estim2 <- as.data.frame(emmeans::emtrends(model2, "mined", var = "cover", regrid = "response"))
   expect_equal(estim$Coefficient, estim2$cover.trend, tolerance = 1e-2)
diff --git a/vignettes/derivatives.Rmd b/vignettes/derivatives.Rmd
index c7b72897..af66bd34 100644
--- a/vignettes/derivatives.Rmd
+++ b/vignettes/derivatives.Rmd
@@ -143,7 +143,7 @@ We know that the slope is `12.75` (from our parameters analysis). Does it change
 To compute the derivative, we can use the `estimate_slopes()` function, and specify that we want to know: the trend of *x* **over the course of** ("at") itself.
 
 ```{r}
-deriv <- modelbased::estimate_slopes(model_lm, trend = "x", at = "x")
+deriv <- modelbased::estimate_slopes(model_lm, trend = "x", by = "x")
 
 plot(deriv) + # add a dashed line at 0 to show absence of effect
   geom_hline(yintercept = 0, linetype = "dashed")
@@ -162,7 +162,7 @@ Such as GAMs. Lets' to the same for our GAM model:
 <!-- TODO: currently fails with emmeans 1.8.0 //-->
 
 ```{r, eval=FALSE}
-deriv <- modelbased::estimate_slopes(model_gam, trend = "x", at = "x")
+deriv <- modelbased::estimate_slopes(model_gam, trend = "x", by = "x")
 
 plot(deriv, line = list(color = "blue")) +
   geom_hline(yintercept = 0, linetype = "dashed")
@@ -205,7 +205,7 @@ Instead, we can rely on our *good ol'* derivatives to obtain the "linear slope"
 
 
 ```{r}
-deriv <- modelbased::estimate_slopes(model_poly, trend = "x", at = "x")
+deriv <- modelbased::estimate_slopes(model_poly, trend = "x", by = "x")
 
 plot(deriv) +
   geom_hline(yintercept = 0, linetype = "dashed")
@@ -228,7 +228,7 @@ model_gam2 <- mgcv::gam(y2 ~ s(x), data = data)
 plot(modelbased::estimate_relation(model_gam2, length = 100), line = list(color = "blue"))
 
 # Increase precision
-deriv <- modelbased::estimate_slopes(model_gam2, trend = "x", at = "x", length = 100)
+deriv <- modelbased::estimate_slopes(model_gam2, trend = "x", by = "x", length = 100)
 
 plot(deriv, line = list(color = "blue")) +
   geom_hline(yintercept = 0, linetype = "dashed")
@@ -250,7 +250,7 @@ model_gam3 <- mgcv::gam(y3 ~ s(x), data = data)
 
 plot(modelbased::estimate_relation(model_gam3, length = 100), line = list(color = "blue"))
 
-deriv <- modelbased::estimate_slopes(model_gam3, trend = "x", at = "x", length = 100)
+deriv <- modelbased::estimate_slopes(model_gam3, trend = "x", by = "x", length = 100)
 
 plot(deriv, line = list(color = "blue")) +
   geom_hline(yintercept = 0, linetype = "dashed")
diff --git a/vignettes/describe_nonlinear.Rmd b/vignettes/describe_nonlinear.Rmd
index 06c3c760..9b7ee988 100644
--- a/vignettes/describe_nonlinear.Rmd
+++ b/vignettes/describe_nonlinear.Rmd
@@ -159,7 +159,7 @@ We can visualize the link between the Delay and the Memory score by using the
 ```{r}
 library(modelbased)
 
-estim <- estimate_relation(model, at = "Delay", ci = c(0.50, 0.69, 0.89, 0.97))
+estim <- estimate_relation(model, by = "Delay", ci = c(0.50, 0.69, 0.89, 0.97))
 
 ggplot(estim, aes(x = Delay, y = Predicted)) +
   geom_jitter(data = df, aes(y = Memory), width = 0.2, height = 0.2) +
diff --git a/vignettes/estimate_contrasts.Rmd b/vignettes/estimate_contrasts.Rmd
index 2d36d06e..785e66ce 100644
--- a/vignettes/estimate_contrasts.Rmd
+++ b/vignettes/estimate_contrasts.Rmd
@@ -111,7 +111,7 @@ continuous variable. Based on the model above (including the interaction with
 `Petal.Width`, that we will then visualise.
 
 ```{r, fig}
-contrasts <- estimate_contrasts(model, at = "Petal.Width", length = 100)
+contrasts <- estimate_contrasts(model, by = "Petal.Width", length = 100)
 
 # Create a variable with the two levels concatenated
 contrasts$Contrast <- paste(contrasts$Level1, "-", contrasts$Level2)
diff --git a/vignettes/estimate_slopes.Rmd b/vignettes/estimate_slopes.Rmd
index a8589493..a365d937 100644
--- a/vignettes/estimate_slopes.Rmd
+++ b/vignettes/estimate_slopes.Rmd
@@ -85,7 +85,7 @@ We can see that the effect of `Petal.Length`, **marginalized over Species**, is
 
 
 ```{r}
-slopes <- estimate_slopes(model, trend = "Petal.Length", at = "Species")
+slopes <- estimate_slopes(model, trend = "Petal.Length", by = "Species")
 
 slopes
 plot(slopes)
@@ -105,7 +105,7 @@ library(modelbased)
 
 model <- lm(mpg ~ hp * wt, data = mtcars)
 
-slopes <- estimate_slopes(model, trend = "hp", at = "wt")
+slopes <- estimate_slopes(model, trend = "hp", by = "wt")
 
 plot(slopes)
 ```
@@ -144,7 +144,7 @@ Because the meaning of these parameters is somewhat disconnected with our need f
 # Compute derivative
 deriv <- estimate_slopes(model,
   trend = "Petal.Length",
-  at = "Petal.Length",
+  by = "Petal.Length",
   length = 100
 )
 
diff --git a/vignettes/modelisation_approach.Rmd b/vignettes/modelisation_approach.Rmd
index 39a83fad..869390a3 100644
--- a/vignettes/modelisation_approach.Rmd
+++ b/vignettes/modelisation_approach.Rmd
@@ -160,7 +160,7 @@ What can we conclude from that? Absolutely **nothing**! We need to investigate i
 ## Post-hoc comparison tests
 
 ```{r}
-posthoc <- model_emmeans(results, at = c("Condition", "Category")) %>%
+posthoc <- model_emmeans(results, by = c("Condition", "Category")) %>%
   pairs()
 parameters(posthoc)
 ```